QGIS API Documentation 3.99.0-Master (357b655ed83)
Loading...
Searching...
No Matches
CloughTocherInterpolator.cpp
Go to the documentation of this file.
1/***************************************************************************
2 CloughTocherInterpolator.cpp
3 ----------------------------
4 copyright : (C) 2004 by Marco Hugentobler
5 email : [email protected]
6 ***************************************************************************/
7
8/***************************************************************************
9 * *
10 * This program is free software; you can redistribute it and/or modify *
11 * it under the terms of the GNU General Public License as published by *
12 * the Free Software Foundation; either version 2 of the License, or *
13 * (at your option) any later version. *
14 * *
15 ***************************************************************************/
16
17#include "CloughTocherInterpolator.h"
18
19#include "MathUtils.h"
20#include "NormVecDecorator.h"
21#include "qgslogger.h"
22
23#include <QString>
24
25using namespace Qt::StringLiterals;
26
31
36
37double CloughTocherInterpolator::calcBernsteinPoly( int n, int i, int j, int k, double u, double v, double w )
38{
39 if ( i < 0 || j < 0 || k < 0 )
40 {
41 QgsDebugError( u"Invalid parameters for Bernstein poly calculation!"_s );
42 return 0;
43 }
44
45 const double result = MathUtils::faculty( n ) * std::pow( u, i ) * std::pow( v, j ) * std::pow( w, k ) / ( MathUtils::faculty( i ) * MathUtils::faculty( j ) * MathUtils::faculty( k ) );
46 return result;
47}
48
49static void normalize( QgsPoint &point )
50{
51 const double length = sqrt( pow( point.x(), 2 ) + pow( point.y(), 2 ) + pow( point.z(), 2 ) );
52 if ( length > 0 )
53 {
54 point.setX( point.x() / length );
55 point.setY( point.y() / length );
56 point.setZ( point.z() / length );
57 }
58}
59
60bool CloughTocherInterpolator::calcNormVec( double x, double y, QgsPoint &result )
61{
62 if ( !result.isEmpty() )
63 {
64 init( x, y );
65 QgsPoint barycoord( 0, 0, 0 ); //barycentric coordinates of (x,y) with respect to the triangle
66 QgsPoint endpointUXY( 0, 0, 0 ); //endpoint of the derivative in u-direction (in xy-coordinates)
67 const QgsPoint endpointV( 0, 0, 0 ); //endpoint of the derivative in v-direction (in barycentric coordinates)
68 QgsPoint endpointVXY( 0, 0, 0 ); //endpoint of the derivative in v-direction (in xy-coordinates)
69
70 //find out, in which subtriangle the point (x,y) is
71 //is the point in the first subtriangle (point1,point2,cp10)?
72 MathUtils::calcBarycentricCoordinates( x, y, &point1, &point2, &cp10, &barycoord );
73 if ( barycoord.x() >= -mEdgeTolerance && barycoord.x() <= ( 1 + mEdgeTolerance ) && barycoord.y() >= -mEdgeTolerance && barycoord.y() <= ( 1 + mEdgeTolerance ) )
74 {
75 const double zu = point1.z() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp1.z() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp2.z() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp3.z() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp4.z() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp7.z() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() );
76 const double zv = cp1.z() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp2.z() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + point2.z() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp4.z() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp5.z() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp8.z() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() );
77 const double zw = cp3.z() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp4.z() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp5.z() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp7.z() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp8.z() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp10.z() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() );
78 MathUtils::BarycentricToXY( barycoord.x() + 1, barycoord.y() - 1, barycoord.z(), &point1, &point2, &cp10, &endpointUXY );
79 endpointUXY.setZ( 3 * ( zu - zv ) );
80 MathUtils::BarycentricToXY( barycoord.x(), barycoord.y() + 1, barycoord.z() - 1, &point1, &point2, &cp10, &endpointVXY );
81 endpointVXY.setZ( 3 * ( zv - zw ) );
82 const Vector3D v1( endpointUXY.x() - x, endpointUXY.y() - y, endpointUXY.z() );
83 const Vector3D v2( endpointVXY.x() - x, endpointVXY.y() - y, endpointVXY.z() );
84 result.setX( v1.getY() * v2.getZ() - v1.getZ() * v2.getY() );
85 result.setY( v1.getZ() * v2.getX() - v1.getX() * v2.getZ() );
86 result.setZ( v1.getX() * v2.getY() - v1.getY() * v2.getX() );
87 normalize( result );
88 return true;
89 }
90 //is the point in the second subtriangle (point2,point3,cp10)?
91 MathUtils::calcBarycentricCoordinates( x, y, &point2, &point3, &cp10, &barycoord );
92 if ( barycoord.x() >= -mEdgeTolerance && barycoord.x() <= ( 1 + mEdgeTolerance ) && barycoord.y() >= -mEdgeTolerance && barycoord.y() <= ( 1 + mEdgeTolerance ) )
93 {
94 const double zu = point2.z() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp9.z() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp16.z() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp5.z() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp13.z() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp8.z() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() );
95 const double zv = cp9.z() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp16.z() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + point3.z() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp13.z() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp15.z() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp12.z() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() );
96 const double zw = cp5.z() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp13.z() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp15.z() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp8.z() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp12.z() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp10.z() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() );
97 MathUtils::BarycentricToXY( barycoord.x() + 1, barycoord.y() - 1, barycoord.z(), &point2, &point3, &cp10, &endpointUXY );
98 endpointUXY.setZ( 3 * ( zu - zv ) );
99 MathUtils::BarycentricToXY( barycoord.x(), barycoord.y() + 1, barycoord.z() - 1, &point2, &point3, &cp10, &endpointVXY );
100 endpointVXY.setZ( 3 * ( zv - zw ) );
101 const Vector3D v1( endpointUXY.x() - x, endpointUXY.y() - y, endpointUXY.z() );
102 const Vector3D v2( endpointVXY.x() - x, endpointVXY.y() - y, endpointVXY.z() );
103 result.setX( v1.getY() * v2.getZ() - v1.getZ() * v2.getY() );
104 result.setY( v1.getZ() * v2.getX() - v1.getX() * v2.getZ() );
105 result.setZ( v1.getX() * v2.getY() - v1.getY() * v2.getX() );
106 normalize( result );
107 return true;
108 }
109 //is the point in the third subtriangle (point3,point1,cp10)?
110 MathUtils::calcBarycentricCoordinates( x, y, &point3, &point1, &cp10, &barycoord );
111 if ( barycoord.x() >= -mEdgeTolerance && barycoord.x() <= ( 1 + mEdgeTolerance ) && barycoord.y() >= -mEdgeTolerance && barycoord.y() <= ( 1 + mEdgeTolerance ) )
112 {
113 const double zu = point3.z() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp14.z() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp6.z() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp15.z() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp11.z() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp12.z() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() );
114 const double zv = cp14.z() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp6.z() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + point1.z() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp11.z() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp3.z() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp7.z() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() );
115 const double zw = cp15.z() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp11.z() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp3.z() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp12.z() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp7.z() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp10.z() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() );
116 MathUtils::BarycentricToXY( barycoord.x() + 1, barycoord.y() - 1, barycoord.z(), &point3, &point1, &cp10, &endpointUXY );
117 endpointUXY.setZ( 3 * ( zu - zv ) );
118 MathUtils::BarycentricToXY( barycoord.x(), barycoord.y() + 1, barycoord.z() - 1, &point3, &point1, &cp10, &endpointVXY );
119 endpointVXY.setZ( 3 * ( zv - zw ) );
120 const Vector3D v1( endpointUXY.x() - x, endpointUXY.y() - y, endpointUXY.z() );
121 const Vector3D v2( endpointVXY.x() - x, endpointVXY.y() - y, endpointVXY.z() );
122 result.setX( v1.getY() * v2.getZ() - v1.getZ() * v2.getY() );
123 result.setY( v1.getZ() * v2.getX() - v1.getX() * v2.getZ() );
124 result.setZ( v1.getX() * v2.getY() - v1.getY() * v2.getX() );
125 normalize( result );
126 return true;
127 }
128
129 //the point is in none of the subtriangles, test if point has the same position as one of the vertices
130 if ( x == point1.x() && y == point1.y() )
131 {
132 result.setX( -der1X );
133 result.setY( -der1Y );
134 result.setZ( 1 );
135 normalize( result );
136 return true;
137 }
138 else if ( x == point2.x() && y == point2.y() )
139 {
140 result.setX( -der2X );
141 result.setY( -der2Y );
142 result.setZ( 1 );
143 normalize( result );
144 return true;
145 }
146 else if ( x == point3.x() && y == point3.y() )
147 {
148 result.setX( -der3X );
149 result.setY( -der3Y );
150 result.setZ( 1 );
151 normalize( result );
152 return true;
153 }
154
155 result.setX( 0 ); //return a vertical normal if failed
156 result.setY( 0 );
157 result.setZ( 1 );
158 return false;
159 }
160 else
161 {
162 QgsDebugError( u"warning, null pointer"_s );
163 return false;
164 }
165}
166
167bool CloughTocherInterpolator::calcPoint( double x, double y, QgsPoint &result )
168{
169 init( x, y );
170 //find out, in which subtriangle the point (x,y) is
171 QgsPoint barycoord( 0, 0, 0 );
172 //is the point in the first subtriangle (point1,point2,cp10)?
173 MathUtils::calcBarycentricCoordinates( x, y, &point1, &point2, &cp10, &barycoord );
174 if ( barycoord.x() >= -mEdgeTolerance && barycoord.x() <= ( 1 + mEdgeTolerance ) && barycoord.y() >= -mEdgeTolerance && barycoord.y() <= ( 1 + mEdgeTolerance ) )
175 {
176 const double z = point1.z() * calcBernsteinPoly( 3, 3, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp1.z() * calcBernsteinPoly( 3, 2, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp2.z() * calcBernsteinPoly( 3, 1, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + point2.z() * calcBernsteinPoly( 3, 0, 3, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp3.z() * calcBernsteinPoly( 3, 2, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp4.z() * calcBernsteinPoly( 3, 1, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp5.z() * calcBernsteinPoly( 3, 0, 2, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp7.z() * calcBernsteinPoly( 3, 1, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() ) + cp8.z() * calcBernsteinPoly( 3, 0, 1, 2, barycoord.x(), barycoord.y(), barycoord.z() ) + cp10.z() * calcBernsteinPoly( 3, 0, 0, 3, barycoord.x(), barycoord.y(), barycoord.z() );
177 result.setX( x );
178 result.setY( y );
179 result.setZ( z );
180 return true;
181 }
182 //is the point in the second subtriangle (point2,point3,cp10)?
183 MathUtils::calcBarycentricCoordinates( x, y, &point2, &point3, &cp10, &barycoord );
184 if ( barycoord.x() >= -mEdgeTolerance && barycoord.x() <= ( 1 + mEdgeTolerance ) && barycoord.y() >= -mEdgeTolerance && barycoord.y() <= ( 1 + mEdgeTolerance ) )
185 {
186 const double z = cp10.z() * calcBernsteinPoly( 3, 0, 0, 3, barycoord.x(), barycoord.y(), barycoord.z() ) + cp8.z() * calcBernsteinPoly( 3, 1, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() ) + cp5.z() * calcBernsteinPoly( 3, 2, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + point2.z() * calcBernsteinPoly( 3, 3, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp12.z() * calcBernsteinPoly( 3, 0, 1, 2, barycoord.x(), barycoord.y(), barycoord.z() ) + cp13.z() * calcBernsteinPoly( 3, 1, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp9.z() * calcBernsteinPoly( 3, 2, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp15.z() * calcBernsteinPoly( 3, 0, 2, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp16.z() * calcBernsteinPoly( 3, 1, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + point3.z() * calcBernsteinPoly( 3, 0, 3, 0, barycoord.x(), barycoord.y(), barycoord.z() );
187 result.setX( x );
188 result.setY( y );
189 result.setZ( z );
190 return true;
191 }
192 //is the point in the third subtriangle (point3,point1,cp10)?
193 MathUtils::calcBarycentricCoordinates( x, y, &point3, &point1, &cp10, &barycoord );
194 if ( barycoord.x() >= -mEdgeTolerance && barycoord.x() <= ( 1 + mEdgeTolerance ) && barycoord.y() >= -mEdgeTolerance && barycoord.y() <= ( 1 + mEdgeTolerance ) )
195 {
196 const double z = point1.z() * calcBernsteinPoly( 3, 0, 3, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp3.z() * calcBernsteinPoly( 3, 0, 2, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp7.z() * calcBernsteinPoly( 3, 0, 1, 2, barycoord.x(), barycoord.y(), barycoord.z() ) + cp10.z() * calcBernsteinPoly( 3, 0, 0, 3, barycoord.x(), barycoord.y(), barycoord.z() ) + cp6.z() * calcBernsteinPoly( 3, 1, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp11.z() * calcBernsteinPoly( 3, 1, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp12.z() * calcBernsteinPoly( 3, 1, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() ) + cp14.z() * calcBernsteinPoly( 3, 2, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp15.z() * calcBernsteinPoly( 3, 2, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + point3.z() * calcBernsteinPoly( 3, 3, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() );
197 result.setX( x );
198 result.setY( y );
199 result.setZ( z );
200 return true;
201 }
202
203 //the point is in none of the subtriangles, test if point has the same position as one of the vertices
204 if ( x == point1.x() && y == point1.y() )
205 {
206 result.setX( x );
207 result.setY( y );
208 result.setZ( point1.z() );
209 return true;
210 }
211 else if ( x == point2.x() && y == point2.y() )
212 {
213 result.setX( x );
214 result.setY( y );
215 result.setZ( point2.z() );
216 return true;
217 }
218 else if ( x == point3.x() && y == point3.y() )
219 {
220 result.setX( x );
221 result.setY( y );
222 result.setZ( point3.z() );
223 return true;
224 }
225 else
226 {
227 result.setX( x );
228 result.setY( y );
229 result.setZ( 0 );
230 }
231
232 return false;
233}
234
235void CloughTocherInterpolator::init( double x, double y ) //version, which has the unintended breaklines within the macrotriangles
236{
237 Vector3D v1, v2, v3; //normals at the three data points
238 int ptn1, ptn2, ptn3; //numbers of the vertex points
239 NormVecDecorator::PointState state1, state2, state3; //states of the vertex points (Normal, BreakLine, EndPoint possible)
240
241 if ( mTIN )
242 {
243 mTIN->getTriangle( x, y, point1, ptn1, &v1, &state1, point2, ptn2, &v2, &state2, point3, ptn3, &v3, &state3 );
244
245 if ( point1 == lpoint1 && point2 == lpoint2 && point3 == lpoint3 ) //if we are in the same triangle as at the last run, we can leave 'init'
246 {
247 return;
248 }
249
250 //calculate the partial derivatives at the data points
251 der1X = -v1.getX() / v1.getZ();
252 der1Y = -v1.getY() / v1.getZ();
253 der2X = -v2.getX() / v2.getZ();
254 der2Y = -v2.getY() / v2.getZ();
255 der3X = -v3.getX() / v3.getZ();
256 der3Y = -v3.getY() / v3.getZ();
257
258 //calculate the control points
259 cp1.setX( point1.x() + ( point2.x() - point1.x() ) / 3 );
260 cp1.setY( point1.y() + ( point2.y() - point1.y() ) / 3 );
261 cp1.setZ( point1.z() + ( cp1.x() - point1.x() ) * der1X + ( cp1.y() - point1.y() ) * der1Y );
262
263 cp2.setX( point2.x() + ( point1.x() - point2.x() ) / 3 );
264 cp2.setY( point2.y() + ( point1.y() - point2.y() ) / 3 );
265 cp2.setZ( point2.z() + ( cp2.x() - point2.x() ) * der2X + ( cp2.y() - point2.y() ) * der2Y );
266
267 cp9.setX( point2.x() + ( point3.x() - point2.x() ) / 3 );
268 cp9.setY( point2.y() + ( point3.y() - point2.y() ) / 3 );
269 cp9.setZ( point2.z() + ( cp9.x() - point2.x() ) * der2X + ( cp9.y() - point2.y() ) * der2Y );
270
271 cp16.setX( point3.x() + ( point2.x() - point3.x() ) / 3 );
272 cp16.setY( point3.y() + ( point2.y() - point3.y() ) / 3 );
273 cp16.setZ( point3.z() + ( cp16.x() - point3.x() ) * der3X + ( cp16.y() - point3.y() ) * der3Y );
274
275 cp14.setX( point3.x() + ( point1.x() - point3.x() ) / 3 );
276 cp14.setY( point3.y() + ( point1.y() - point3.y() ) / 3 );
277 cp14.setZ( point3.z() + ( cp14.x() - point3.x() ) * der3X + ( cp14.y() - point3.y() ) * der3Y );
278
279 cp6.setX( point1.x() + ( point3.x() - point1.x() ) / 3 );
280 cp6.setY( point1.y() + ( point3.y() - point1.y() ) / 3 );
281 cp6.setZ( point1.z() + ( cp6.x() - point1.x() ) * der1X + ( cp6.y() - point1.y() ) * der1Y );
282
283 //set x- and y-coordinates of the central point
284 cp10.setX( ( point1.x() + point2.x() + point3.x() ) / 3 );
285 cp10.setY( ( point1.y() + point2.y() + point3.y() ) / 3 );
286
287 //set the derivatives of the points to new values if they are on a breakline
288 if ( state1 == NormVecDecorator::BreakLine )
289 {
290 Vector3D target1;
291 if ( mTIN->calcNormalForPoint( x, y, ptn1, &target1 ) )
292 {
293 der1X = -target1.getX() / target1.getZ();
294 der1Y = -target1.getY() / target1.getZ();
295
296 if ( state2 == NormVecDecorator::Normal )
297 {
298 //recalculate cp1
299 cp1.setZ( point1.z() + ( cp1.x() - point1.x() ) * der1X + ( cp1.y() - point1.y() ) * der1Y );
300 }
301 if ( state3 == NormVecDecorator::Normal )
302 {
303 //recalculate cp6
304 cp6.setZ( point1.z() + ( cp6.x() - point1.x() ) * der1X + ( cp6.y() - point1.y() ) * der1Y );
305 }
306 }
307 }
308
309 if ( state2 == NormVecDecorator::BreakLine )
310 {
311 Vector3D target2;
312 if ( mTIN->calcNormalForPoint( x, y, ptn2, &target2 ) )
313 {
314 der2X = -target2.getX() / target2.getZ();
315 der2Y = -target2.getY() / target2.getZ();
316
317 if ( state1 == NormVecDecorator::Normal )
318 {
319 //recalculate cp2
320 cp2.setZ( point2.z() + ( cp2.x() - point2.x() ) * der2X + ( cp2.y() - point2.y() ) * der2Y );
321 }
322 if ( state3 == NormVecDecorator::Normal )
323 {
324 //recalculate cp9
325 cp9.setZ( point2.z() + ( cp9.x() - point2.x() ) * der2X + ( cp9.y() - point2.y() ) * der2Y );
326 }
327 }
328 }
329
330 if ( state3 == NormVecDecorator::BreakLine )
331 {
332 Vector3D target3;
333 if ( mTIN->calcNormalForPoint( x, y, ptn3, &target3 ) )
334 {
335 der3X = -target3.getX() / target3.getZ();
336 der3Y = -target3.getY() / target3.getZ();
337
338 if ( state1 == NormVecDecorator::Normal )
339 {
340 //recalculate cp14
341 cp14.setZ( point3.z() + ( cp14.x() - point3.x() ) * der3X + ( cp14.y() - point3.y() ) * der3Y );
342 }
343 if ( state2 == NormVecDecorator::Normal )
344 {
345 //recalculate cp16
346 cp16.setZ( point3.z() + ( cp16.x() - point3.x() ) * der3X + ( cp16.y() - point3.y() ) * der3Y );
347 }
348 }
349 }
350
351
352 cp3.setX( point1.x() + ( cp10.x() - point1.x() ) / 3 );
353 cp3.setY( point1.y() + ( cp10.y() - point1.y() ) / 3 );
354 cp3.setZ( point1.z() + ( cp3.x() - point1.x() ) * der1X + ( cp3.y() - point1.y() ) * der1Y );
355
356 cp5.setX( point2.x() + ( cp10.x() - point2.x() ) / 3 );
357 cp5.setY( point2.y() + ( cp10.y() - point2.y() ) / 3 );
358 cp5.setZ( point2.z() + ( cp5.x() - point2.x() ) * der2X + ( cp5.y() - point2.y() ) * der2Y );
359
360 cp15.setX( point3.x() + ( cp10.x() - point3.x() ) / 3 );
361 cp15.setY( point3.y() + ( cp10.y() - point3.y() ) / 3 );
362 cp15.setZ( point3.z() + ( cp15.x() - point3.x() ) * der3X + ( cp15.y() - point3.y() ) * der3Y );
363
364
365 cp4.setX( ( point1.x() + cp10.x() + point2.x() ) / 3 );
366 cp4.setY( ( point1.y() + cp10.y() + point2.y() ) / 3 );
367
368 const QgsPoint midpoint3( ( cp1.x() + cp2.x() ) / 2, ( cp1.y() + cp2.y() ) / 2, ( cp1.z() + cp2.z() ) / 2 );
369 const Vector3D cp1cp2( cp2.x() - cp1.x(), cp2.y() - cp1.y(), cp2.z() - cp1.z() );
370 Vector3D odir3( 0, 0, 0 ); //direction orthogonal to the line from point1 to point2
371 if ( ( point2.y() - point1.y() ) != 0 ) //avoid division through zero
372 {
373 odir3.setX( 1 );
374 odir3.setY( -( point2.x() - point1.x() ) / ( point2.y() - point1.y() ) );
375 odir3.setZ( ( der1X + der2X ) / 2 + ( der1Y + der2Y ) / 2 * odir3.getY() ); //take the linear interpolated cross-derivative
376 }
377 else
378 {
379 odir3.setY( 1 );
380 odir3.setX( -( point2.y() - point1.y() ) / ( point2.x() - point1.x() ) );
381 odir3.setZ( ( der1X + der2X ) / 2 * odir3.getX() + ( der1Y + der2Y ) / 2 );
382 }
383 Vector3D midpoint3cp4( 0, 0, 0 );
384 MathUtils::derVec( &cp1cp2, &odir3, &midpoint3cp4, cp4.x() - midpoint3.x(), cp4.y() - midpoint3.y() );
385 cp4.setZ( midpoint3.z() + midpoint3cp4.getZ() );
386
387 cp13.setX( ( point2.x() + cp10.x() + point3.x() ) / 3 );
388 cp13.setY( ( point2.y() + cp10.y() + point3.y() ) / 3 );
389 //find the point in the middle of the bezier curve between point2 and point3
390 const QgsPoint midpoint1( ( cp9.x() + cp16.x() ) / 2, ( cp9.y() + cp16.y() ) / 2, ( cp9.z() + cp16.z() ) / 2 );
391 const Vector3D cp9cp16( cp16.x() - cp9.x(), cp16.y() - cp9.y(), cp16.z() - cp9.z() );
392 Vector3D odir1( 0, 0, 0 ); //direction orthogonal to the line from point2 to point3
393 if ( ( point3.y() - point2.y() ) != 0 ) //avoid division through zero
394 {
395 odir1.setX( 1 );
396 odir1.setY( -( point3.x() - point2.x() ) / ( point3.y() - point2.y() ) );
397 odir1.setZ( ( der2X + der3X ) / 2 + ( der2Y + der3Y ) / 2 * odir1.getY() ); //take the linear interpolated cross-derivative
398 }
399 else
400 {
401 odir1.setY( 1 );
402 odir1.setX( -( point3.y() - point2.y() ) / ( point3.x() - point2.x() ) );
403 odir1.setZ( ( der2X + der3X ) / 2 * odir1.getX() + ( der2Y + der3Y ) / 2 );
404 }
405 Vector3D midpoint1cp13( 0, 0, 0 );
406 MathUtils::derVec( &cp9cp16, &odir1, &midpoint1cp13, cp13.x() - midpoint1.x(), cp13.y() - midpoint1.y() );
407 cp13.setZ( midpoint1.z() + midpoint1cp13.getZ() );
408
409
410 cp11.setX( ( point3.x() + cp10.x() + point1.x() ) / 3 );
411 cp11.setY( ( point3.y() + cp10.y() + point1.y() ) / 3 );
412 //find the point in the middle of the bezier curve between point3 and point2
413 const QgsPoint midpoint2( ( cp14.x() + cp6.x() ) / 2, ( cp14.y() + cp6.y() ) / 2, ( cp14.z() + cp6.z() ) / 2 );
414 const Vector3D cp14cp6( cp6.x() - cp14.x(), cp6.y() - cp14.y(), cp6.z() - cp14.z() );
415 Vector3D odir2( 0, 0, 0 ); //direction orthogonal to the line from point 1 to point3
416 if ( ( point3.y() - point1.y() ) != 0 ) //avoid division through zero
417 {
418 odir2.setX( 1 );
419 odir2.setY( -( point3.x() - point1.x() ) / ( point3.y() - point1.y() ) );
420 odir2.setZ( ( der3X + der1X ) / 2 + ( der3Y + der1Y ) / 2 * odir2.getY() ); //take the linear interpolated cross-derivative
421 }
422 else
423 {
424 odir2.setY( 1 );
425 odir2.setX( -( point3.y() - point1.y() ) / ( point3.x() - point1.x() ) );
426 odir2.setZ( ( der3X + der1X ) / 2 * odir2.getX() + ( der3Y + der1Y ) / 2 );
427 }
428 Vector3D midpoint2cp11( 0, 0, 0 );
429 MathUtils::derVec( &cp14cp6, &odir2, &midpoint2cp11, cp11.x() - midpoint2.x(), cp11.y() - midpoint2.y() );
430 cp11.setZ( midpoint2.z() + midpoint2cp11.getZ() );
431
432
433 cp7.setX( cp10.x() + ( point1.x() - cp10.x() ) / 3 );
434 cp7.setY( cp10.y() + ( point1.y() - cp10.y() ) / 3 );
435 //cp7 has to be in the same plane as cp4, cp3 and cp11
436 const Vector3D cp4cp3( cp3.x() - cp4.x(), cp3.y() - cp4.y(), cp3.z() - cp4.z() );
437 const Vector3D cp4cp11( cp11.x() - cp4.x(), cp11.y() - cp4.y(), cp11.z() - cp4.z() );
438 Vector3D cp4cp7( 0, 0, 0 );
439 MathUtils::derVec( &cp4cp3, &cp4cp11, &cp4cp7, cp7.x() - cp4.x(), cp7.y() - cp4.y() );
440 cp7.setZ( cp4.z() + cp4cp7.getZ() );
441
442 cp8.setX( cp10.x() + ( point2.x() - cp10.x() ) / 3 );
443 cp8.setY( cp10.y() + ( point2.y() - cp10.y() ) / 3 );
444 //cp8 has to be in the same plane as cp4, cp5 and cp13
445 const Vector3D cp4cp5( cp5.x() - cp4.x(), cp5.y() - cp4.y(), cp5.z() - cp4.z() );
446 const Vector3D cp4cp13( cp13.x() - cp4.x(), cp13.y() - cp4.y(), cp13.z() - cp4.z() );
447 Vector3D cp4cp8( 0, 0, 0 );
448 MathUtils::derVec( &cp4cp5, &cp4cp13, &cp4cp8, cp8.x() - cp4.x(), cp8.y() - cp4.y() );
449 cp8.setZ( cp4.z() + cp4cp8.getZ() );
450
451 cp12.setX( cp10.x() + ( point3.x() - cp10.x() ) / 3 );
452 cp12.setY( cp10.y() + ( point3.y() - cp10.y() ) / 3 );
453 //cp12 has to be in the same plane as cp13, cp15 and cp11
454 const Vector3D cp13cp11( cp11.x() - cp13.x(), cp11.y() - cp13.y(), cp11.z() - cp13.z() );
455 const Vector3D cp13cp15( cp15.x() - cp13.x(), cp15.y() - cp13.y(), cp15.z() - cp13.z() );
456 Vector3D cp13cp12( 0, 0, 0 );
457 MathUtils::derVec( &cp13cp11, &cp13cp15, &cp13cp12, cp12.x() - cp13.x(), cp12.y() - cp13.y() );
458 cp12.setZ( cp13.z() + cp13cp12.getZ() );
459
460 //cp10 has to be in the same plane as cp7, cp8 and cp12
461 const Vector3D cp7cp8( cp8.x() - cp7.x(), cp8.y() - cp7.y(), cp8.z() - cp7.z() );
462 const Vector3D cp7cp12( cp12.x() - cp7.x(), cp12.y() - cp7.y(), cp12.z() - cp7.z() );
463 Vector3D cp7cp10( 0, 0, 0 );
464 MathUtils::derVec( &cp7cp8, &cp7cp12, &cp7cp10, cp10.x() - cp7.x(), cp10.y() - cp7.y() );
465 cp10.setZ( cp7.z() + cp7cp10.getZ() );
466
467 lpoint1 = point1;
468 lpoint2 = point2;
469 lpoint3 = point3;
470 }
471 else
472 {
473 QgsDebugError( u"warning, null pointer"_s );
474 }
475}
476
477
478#if 0
479void CloughTocherInterpolator::init( double x, double y )//version which has unintended breaklines similar to the Coons interpolator
480{
481 Vector3D v1, v2, v3;//normals at the three data points
482 int ptn1, ptn2, ptn3;//numbers of the vertex points
483 NormVecDecorator::PointState state1, state2, state3;//states of the vertex points (Normal, BreakLine, EndPoint possible)
484
485 if ( mTIN )
486 {
487 mTIN->getTriangle( x, y, &point1, &ptn1, &v1, &state1, &point2, &ptn2, &v2, &state2, &point3, &ptn3, &v3, &state3 );
488
489 if ( point1 == lpoint1 && point2 == lpoint2 && point3 == lpoint3 )//if we are in the same triangle as at the last run, we can leave 'init'
490 {
491 return;
492 }
493
494 //calculate the partial derivatives at the data points
495 der1X = -v1.getX() / v1.getZ();
496 der1Y = -v1.getY() / v1.getZ();
497 der2X = -v2.getX() / v2.getZ();
498 der2Y = -v2.getY() / v2.getZ();
499 der3X = -v3.getX() / v3.getZ();
500 der3Y = -v3.getY() / v3.getZ();
501
502 //calculate the control points
503 cp1.setX( point1.getX() + ( point2.getX() - point1.getX() ) / 3 );
504 cp1.setY( point1.getY() + ( point2.getY() - point1.getY() ) / 3 );
505 cp1.setZ( point1.getZ() + ( cp1.getX() - point1.getX() )*der1X + ( cp1.getY() - point1.getY() )*der1Y );
506
507 cp2.setX( point2.getX() + ( point1.getX() - point2.getX() ) / 3 );
508 cp2.setY( point2.getY() + ( point1.getY() - point2.getY() ) / 3 );
509 cp2.setZ( point2.getZ() + ( cp2.getX() - point2.getX() )*der2X + ( cp2.getY() - point2.getY() )*der2Y );
510
511 cp9.setX( point2.getX() + ( point3.getX() - point2.getX() ) / 3 );
512 cp9.setY( point2.getY() + ( point3.getY() - point2.getY() ) / 3 );
513 cp9.setZ( point2.getZ() + ( cp9.getX() - point2.getX() )*der2X + ( cp9.getY() - point2.getY() )*der2Y );
514
515 cp16.setX( point3.getX() + ( point2.getX() - point3.getX() ) / 3 );
516 cp16.setY( point3.getY() + ( point2.getY() - point3.getY() ) / 3 );
517 cp16.setZ( point3.getZ() + ( cp16.getX() - point3.getX() )*der3X + ( cp16.getY() - point3.getY() )*der3Y );
518
519 cp14.setX( point3.getX() + ( point1.getX() - point3.getX() ) / 3 );
520 cp14.setY( point3.getY() + ( point1.getY() - point3.getY() ) / 3 );
521 cp14.setZ( point3.getZ() + ( cp14.getX() - point3.getX() )*der3X + ( cp14.getY() - point3.getY() )*der3Y );
522
523 cp6.setX( point1.getX() + ( point3.getX() - point1.getX() ) / 3 );
524 cp6.setY( point1.getY() + ( point3.getY() - point1.getY() ) / 3 );
525 cp6.setZ( point1.getZ() + ( cp6.getX() - point1.getX() )*der1X + ( cp6.getY() - point1.getY() )*der1Y );
526
527 //set x- and y-coordinates of the central point
528 cp10.setX( ( point1.getX() + point2.getX() + point3.getX() ) / 3 );
529 cp10.setY( ( point1.getY() + point2.getY() + point3.getY() ) / 3 );
530
531 //do the necessary adjustments in case of breaklines--------------------------------------------------------------------
532
533 //temporary normals and derivatives
534 double tmpx = 0;
535 double tmpy = 0;
536 Vector3D tmp( 0, 0, 0 );
537
538 //point1
539 if ( state1 == NormVecDecorator::BreakLine )
540 {
541 if ( mTIN->calcNormalForPoint( x, y, ptn1, &tmp ) )
542 {
543 tmpx = -tmp.getX() / tmp.getZ();
544 tmpy = -tmp.getY() / tmp.getZ();
545 if ( state3 == NormVecDecorator::Normal )
546 {
547 cp6.setZ( point1.getZ() + ( ( point3.getX() - point1.getX() ) / 3 )*tmpx + ( ( point3.getY() - point1.getY() ) / 3 )*tmpy );
548 }
549 if ( state2 == NormVecDecorator::Normal )
550 {
551 cp1.setZ( point1.getZ() + ( ( point2.getX() - point1.getX() ) / 3 )*tmpx + ( ( point2.getY() - point1.getY() ) / 3 )*tmpy );
552 }
553 }
554 }
555
556 if ( state2 == NormVecDecorator::Breakline )
557 {
558 if ( mTIN->calcNormalForPoint( x, y, ptn2, &tmp ) )
559 {
560 tmpx = -tmp.getX() / tmp.getZ();
561 tmpy = -tmp.getY() / tmp.getZ();
562 if ( state1 == NormVecDecorator::Normal )
563 {
564 cp2.setZ( point2.getZ() + ( ( point1.getX() - point2.getX() ) / 3 )*tmpx + ( ( point1.getY() - point2.getY() ) / 3 )*tmpy );
565 }
566 if ( state3 == NormVecDecorator::Normal )
567 {
568 cp9.setZ( point2.getZ() + ( ( point3.getX() - point2.getX() ) / 3 )*tmpx + ( ( point3.getY() - point2.getY() ) / 3 )*tmpy );
569 }
570 }
571 }
572
573 if ( state3 == NormVecDecorator::Breakline )
574 {
575 if ( mTIN->calcNormalForPoint( x, y, ptn3, &tmp ) )
576 {
577 tmpx = -tmp.getX() / tmp.getZ();
578 tmpy = -tmp.getY() / tmp.getZ();
579 if ( state1 == NormVecDecorator::Normal )
580 {
581 cp14.setZ( point3.getZ() + ( ( point1.getX() - point3.getX() ) / 3 )*tmpx + ( ( point1.getY() - point3.getY() ) / 3 )*tmpy );
582 }
583 if ( state2 == NormVecDecorator::Normal )
584 {
585 cp16.setZ( point3.getZ() + ( ( point2.getX() - point3.getX() ) / 3 )*tmpx + ( ( point2.getY() - point3.getY() ) / 3 )*tmpy );
586 }
587 }
588 }
589
590 //calculate cp3, cp 5 and cp15
591 Vector3D normal;//temporary normal for the vertices on breaklines
592
593 cp3.setX( point1.getX() + ( cp10.getX() - point1.getX() ) / 3 );
594 cp3.setY( point1.getY() + ( cp10.getY() - point1.getY() ) / 3 );
595 if ( state1 == NormVecDecorator::Breakline )
596 {
597 MathUtils::normalFromPoints( &point1, &cp1, &cp6, &normal );
598 //recalculate der1X and der1Y
599 der1X = -normal.getX() / normal.getZ();
600 der1Y = -normal.getY() / normal.getZ();
601 }
602
603 cp3.setZ( point1.getZ() + ( cp3.getX() - point1.getX() )*der1X + ( cp3.getY() - point1.getY() )*der1Y );
604
605
606
607 cp5.setX( point2.getX() + ( cp10.getX() - point2.getX() ) / 3 );
608 cp5.setY( point2.getY() + ( cp10.getY() - point2.getY() ) / 3 );
609 if ( state2 == NormVecDecorator::Breakline )
610 {
611 MathUtils::normalFromPoints( &point2, &cp9, &cp2, &normal );
612 //recalculate der2X and der2Y
613 der2X = -normal.getX() / normal.getZ();
614 der2Y = -normal.getY() / normal.getZ();
615 }
616
617 cp5.setZ( point2.getZ() + ( cp5.getX() - point2.getX() )*der2X + ( cp5.getY() - point2.getY() )*der2Y );
618
619
620 cp15.setX( point3.getX() + ( cp10.getX() - point3.getX() ) / 3 );
621 cp15.setY( point3.getY() + ( cp10.getY() - point3.getY() ) / 3 );
622 if ( state3 == NormVecDecorator::Breakline )
623 {
624 MathUtils::normalFromPoints( &point3, &cp14, &cp16, &normal );
625 //recalculate der3X and der3Y
626 der3X = -normal.getX() / normal.getZ();
627 der3Y = -normal.getY() / normal.getZ();
628 }
629
630 cp15.setZ( point3.getZ() + ( cp15.getX() - point3.getX() )*der3X + ( cp15.getY() - point3.getY() )*der3Y );
631
632
633 cp4.setX( ( point1.getX() + cp10.getX() + point2.getX() ) / 3 );
634 cp4.setY( ( point1.getY() + cp10.getY() + point2.getY() ) / 3 );
635
636 QgsPoint midpoint3( ( cp1.getX() + cp2.getX() ) / 2, ( cp1.getY() + cp2.getY() ) / 2, ( cp1.getZ() + cp2.getZ() ) / 2 );
637 Vector3D cp1cp2( cp2.getX() - cp1.getX(), cp2.getY() - cp1.getY(), cp2.getZ() - cp1.getZ() );
638 Vector3D odir3( 0, 0, 0 );//direction orthogonal to the line from point1 to point2
639 if ( ( point2.getY() - point1.getY() ) != 0 ) //avoid division through zero
640 {
641 odir3.setX( 1 );
642 odir3.setY( -( point2.getX() - point1.getX() ) / ( point2.getY() - point1.getY() ) );
643 odir3.setZ( ( der1X + der2X ) / 2 + ( der1Y + der2Y ) / 2 * odir3.getY() ); //take the linear interpolated cross-derivative
644 }
645 else
646 {
647 odir3.setY( 1 );
648 odir3.setX( -( point2.getY() - point1.getY() ) / ( point2.getX() - point1.getX() ) );
649 odir3.setZ( ( der1X + der2X ) / 2 * odir3.getX() + ( der1Y + der2Y ) / 2 );
650 }
651 Vector3D midpoint3cp4( 0, 0, 0 );
652 MathUtils::derVec( &cp1cp2, &odir3, &midpoint3cp4, cp4.getX() - midpoint3.getX(), cp4.getY() - midpoint3.getY() );
653 cp4.setZ( midpoint3.getZ() + midpoint3cp4.getZ() );
654
655 cp13.setX( ( point2.getX() + cp10.getX() + point3.getX() ) / 3 );
656 cp13.setY( ( point2.getY() + cp10.getY() + point3.getY() ) / 3 );
657 //find the point in the middle of the bezier curve between point2 and point3
658 QgsPoint midpoint1( ( cp9.getX() + cp16.getX() ) / 2, ( cp9.getY() + cp16.getY() ) / 2, ( cp9.getZ() + cp16.getZ() ) / 2 );
659 Vector3D cp9cp16( cp16.getX() - cp9.getX(), cp16.getY() - cp9.getY(), cp16.getZ() - cp9.getZ() );
660 Vector3D odir1( 0, 0, 0 );//direction orthogonal to the line from point2 to point3
661 if ( ( point3.getY() - point2.getY() ) != 0 ) //avoid division through zero
662 {
663 odir1.setX( 1 );
664 odir1.setY( -( point3.getX() - point2.getX() ) / ( point3.getY() - point2.getY() ) );
665 odir1.setZ( ( der2X + der3X ) / 2 + ( der2Y + der3Y ) / 2 * odir1.getY() ); //take the linear interpolated cross-derivative
666 }
667 else
668 {
669 odir1.setY( 1 );
670 odir1.setX( -( point3.getY() - point2.getY() ) / ( point3.getX() - point2.getX() ) );
671 odir1.setZ( ( der2X + der3X ) / 2 * odir1.getX() + ( der2Y + der3Y ) / 2 );
672 }
673 Vector3D midpoint1cp13( 0, 0, 0 );
674 MathUtils::derVec( &cp9cp16, &odir1, &midpoint1cp13, cp13.getX() - midpoint1.getX(), cp13.getY() - midpoint1.getY() );
675 cp13.setZ( midpoint1.getZ() + midpoint1cp13.getZ() );
676
677
678 cp11.setX( ( point3.getX() + cp10.getX() + point1.getX() ) / 3 );
679 cp11.setY( ( point3.getY() + cp10.getY() + point1.getY() ) / 3 );
680 //find the point in the middle of the bezier curve between point3 and point2
681 QgsPoint midpoint2( ( cp14.getX() + cp6.getX() ) / 2, ( cp14.getY() + cp6.getY() ) / 2, ( cp14.getZ() + cp6.getZ() ) / 2 );
682 Vector3D cp14cp6( cp6.getX() - cp14.getX(), cp6.getY() - cp14.getY(), cp6.getZ() - cp14.getZ() );
683 Vector3D odir2( 0, 0, 0 );//direction orthogonal to the line from point 1 to point3
684 if ( ( point3.getY() - point1.getY() ) != 0 ) //avoid division through zero
685 {
686 odir2.setX( 1 );
687 odir2.setY( -( point3.getX() - point1.getX() ) / ( point3.getY() - point1.getY() ) );
688 odir2.setZ( ( der3X + der1X ) / 2 + ( der3Y + der1Y ) / 2 * odir2.getY() ); //take the linear interpolated cross-derivative
689 }
690 else
691 {
692 odir2.setY( 1 );
693 odir2.setX( -( point3.getY() - point1.getY() ) / ( point3.getX() - point1.getX() ) );
694 odir2.setZ( ( der3X + der1X ) / 2 * odir2.getX() + ( der3Y + der1Y ) / 2 );
695 }
696 Vector3D midpoint2cp11( 0, 0, 0 );
697 MathUtils::derVec( &cp14cp6, &odir2, &midpoint2cp11, cp11.getX() - midpoint2.getX(), cp11.getY() - midpoint2.getY() );
698 cp11.setZ( midpoint2.getZ() + midpoint2cp11.getZ() );
699
700
701 cp7.setX( cp10.getX() + ( point1.getX() - cp10.getX() ) / 3 );
702 cp7.setY( cp10.getY() + ( point1.getY() - cp10.getY() ) / 3 );
703 //cp7 has to be in the same plane as cp4, cp3 and cp11
704 Vector3D cp4cp3( cp3.getX() - cp4.getX(), cp3.getY() - cp4.getY(), cp3.getZ() - cp4.getZ() );
705 Vector3D cp4cp11( cp11.getX() - cp4.getX(), cp11.getY() - cp4.getY(), cp11.getZ() - cp4.getZ() );
706 Vector3D cp4cp7( 0, 0, 0 );
707 MathUtils::derVec( &cp4cp3, &cp4cp11, &cp4cp7, cp7.getX() - cp4.getX(), cp7.getY() - cp4.getY() );
708 cp7.setZ( cp4.getZ() + cp4cp7.getZ() );
709
710 cp8.setX( cp10.getX() + ( point2.getX() - cp10.getX() ) / 3 );
711 cp8.setY( cp10.getY() + ( point2.getY() - cp10.getY() ) / 3 );
712 //cp8 has to be in the same plane as cp4, cp5 and cp13
713 Vector3D cp4cp5( cp5.getX() - cp4.getX(), cp5.getY() - cp4.getY(), cp5.getZ() - cp4.getZ() );
714 Vector3D cp4cp13( cp13.getX() - cp4.getX(), cp13.getY() - cp4.getY(), cp13.getZ() - cp4.getZ() );
715 Vector3D cp4cp8( 0, 0, 0 );
716 MathUtils::derVec( &cp4cp5, &cp4cp13, &cp4cp8, cp8.getX() - cp4.getX(), cp8.getY() - cp4.getY() );
717 cp8.setZ( cp4.getZ() + cp4cp8.getZ() );
718
719 cp12.setX( cp10.getX() + ( point3.getX() - cp10.getX() ) / 3 );
720 cp12.setY( cp10.getY() + ( point3.getY() - cp10.getY() ) / 3 );
721 //cp12 has to be in the same plane as cp13, cp15 and cp11
722 Vector3D cp13cp11( cp11.getX() - cp13.getX(), cp11.getY() - cp13.getY(), cp11.getZ() - cp13.getZ() );
723 Vector3D cp13cp15( cp15.getX() - cp13.getX(), cp15.getY() - cp13.getY(), cp15.getZ() - cp13.getZ() );
724 Vector3D cp13cp12( 0, 0, 0 );
725 MathUtils::derVec( &cp13cp11, &cp13cp15, &cp13cp12, cp12.getX() - cp13.getX(), cp12.getY() - cp13.getY() );
726 cp12.setZ( cp13.getZ() + cp13cp12.getZ() );
727
728 //cp10 has to be in the same plane as cp7, cp8 and cp12
729 Vector3D cp7cp8( cp8.getX() - cp7.getX(), cp8.getY() - cp7.getY(), cp8.getZ() - cp7.getZ() );
730 Vector3D cp7cp12( cp12.getX() - cp7.getX(), cp12.getY() - cp7.getY(), cp12.getZ() - cp7.getZ() );
731 Vector3D cp7cp10( 0, 0, 0 );
732 MathUtils::derVec( &cp7cp8, &cp7cp12, &cp7cp10, cp10.getX() - cp7.getX(), cp10.getY() - cp7.getY() );
733 cp10.setZ( cp7.getZ() + cp7cp10.getZ() );
734
735 lpoint1 = point1;
736 lpoint2 = point2;
737 lpoint3 = point3;
738 }
739
740 else
741 {
742 QgsDebugError( u"warning, null pointer"_s );
743 }
744}
745#endif
CloughTocherInterpolator::der2X
double der2X
Derivative in x-direction at point2.
Definition CloughTocherInterpolator.h:85
CloughTocherInterpolator::setTriangulation
virtual void setTriangulation(NormVecDecorator *tin)
Definition CloughTocherInterpolator.cpp:32
CloughTocherInterpolator::cp8
QgsPoint cp8
Control point 8.
Definition CloughTocherInterpolator.h:63
CloughTocherInterpolator::cp1
QgsPoint cp1
Control point 1.
Definition CloughTocherInterpolator.h:49
CloughTocherInterpolator::der3X
double der3X
Derivative in x-direction at point3.
Definition CloughTocherInterpolator.h:89
CloughTocherInterpolator::cp2
QgsPoint cp2
Control point 2.
Definition CloughTocherInterpolator.h:51
CloughTocherInterpolator::cp14
QgsPoint cp14
Control point 14.
Definition CloughTocherInterpolator.h:75
CloughTocherInterpolator::cp13
QgsPoint cp13
Control point 13.
Definition CloughTocherInterpolator.h:73
CloughTocherInterpolator::cp9
QgsPoint cp9
Control point 9.
Definition CloughTocherInterpolator.h:65
CloughTocherInterpolator::point3
QgsPoint point3
Third point of the triangle in x-,y-,z-coordinates.
Definition CloughTocherInterpolator.h:47
CloughTocherInterpolator::der2Y
double der2Y
Derivative in y-direction at point2.
Definition CloughTocherInterpolator.h:87
CloughTocherInterpolator::mTIN
NormVecDecorator * mTIN
Association with a triangulation object.
Definition CloughTocherInterpolator.h:39
CloughTocherInterpolator::cp10
QgsPoint cp10
Control point 10.
Definition CloughTocherInterpolator.h:67
CloughTocherInterpolator::der1X
double der1X
Derivative in x-direction at point1.
Definition CloughTocherInterpolator.h:81
CloughTocherInterpolator::mEdgeTolerance
double mEdgeTolerance
Tolerance of the barycentric coordinates at the borders of the triangles (to prevent errors because o...
Definition CloughTocherInterpolator.h:41
CloughTocherInterpolator::lpoint1
QgsPoint lpoint1
Stores point1 of the last run.
Definition CloughTocherInterpolator.h:93
CloughTocherInterpolator::cp5
QgsPoint cp5
Control point 5.
Definition CloughTocherInterpolator.h:57
CloughTocherInterpolator::cp4
QgsPoint cp4
Control point 4.
Definition CloughTocherInterpolator.h:55
CloughTocherInterpolator::cp3
QgsPoint cp3
Control point 3.
Definition CloughTocherInterpolator.h:53
CloughTocherInterpolator::lpoint2
QgsPoint lpoint2
Stores point2 of the last run.
Definition CloughTocherInterpolator.h:95
CloughTocherInterpolator::calcPoint
bool calcPoint(double x, double y, QgsPoint &result) override
Performs a linear interpolation in a triangle and assigns the x-,y- and z-coordinates to point.
Definition CloughTocherInterpolator.cpp:167
CloughTocherInterpolator::cp7
QgsPoint cp7
Control point 7.
Definition CloughTocherInterpolator.h:61
CloughTocherInterpolator::calcBernsteinPoly
double calcBernsteinPoly(int n, int i, int j, int k, double u, double v, double w)
Calculates the Bernsteinpolynomials to calculate the Beziertriangle. 'n' is three in the cubical case...
Definition CloughTocherInterpolator.cpp:37
CloughTocherInterpolator::point2
QgsPoint point2
Second point of the triangle in x-,y-,z-coordinates.
Definition CloughTocherInterpolator.h:45
CloughTocherInterpolator::der3Y
double der3Y
Derivative in y-direction at point3.
Definition CloughTocherInterpolator.h:91
CloughTocherInterpolator::cp12
QgsPoint cp12
Control point 12.
Definition CloughTocherInterpolator.h:71
CloughTocherInterpolator::cp16
QgsPoint cp16
Control point 16.
Definition CloughTocherInterpolator.h:79
CloughTocherInterpolator::init
void init(double x, double y)
Finds out, in which triangle the point with the coordinates x and y is.
Definition CloughTocherInterpolator.cpp:235
CloughTocherInterpolator::calcNormVec
bool calcNormVec(double x, double y, QgsPoint &result) override
Calculates the normal vector and assigns it to vec (not implemented at the moment).
Definition CloughTocherInterpolator.cpp:60
CloughTocherInterpolator::cp6
QgsPoint cp6
Control point 6.
Definition CloughTocherInterpolator.h:59
CloughTocherInterpolator::cp15
QgsPoint cp15
Control point 15.
Definition CloughTocherInterpolator.h:77
CloughTocherInterpolator::CloughTocherInterpolator
CloughTocherInterpolator()=default
CloughTocherInterpolator::lpoint3
QgsPoint lpoint3
Stores point3 of the last run.
Definition CloughTocherInterpolator.h:97
CloughTocherInterpolator::point1
QgsPoint point1
First point of the triangle in x-,y-,z-coordinates.
Definition CloughTocherInterpolator.h:43
CloughTocherInterpolator::der1Y
double der1Y
Derivative in y-direction at point1.
Definition CloughTocherInterpolator.h:83
CloughTocherInterpolator::cp11
QgsPoint cp11
Control point 11.
Definition CloughTocherInterpolator.h:69
NormVecDecorator
Decorator class which adds the functionality of estimating normals at the data points.
Definition NormVecDecorator.h:41
NormVecDecorator::getTriangle
bool getTriangle(double x, double y, QgsPoint &p1, Vector3D *v1, QgsPoint &p2, Vector3D *v2, QgsPoint &p3, Vector3D *v3)
Finds out, in which triangle a point with coordinates x and y is and assigns the triangle points to p...
Definition NormVecDecorator.cpp:247
NormVecDecorator::PointState
PointState
Enumeration for the state of a point. Normal means, that the point is not on a BreakLine,...
Definition NormVecDecorator.h:45
QgsPoint
Point geometry type, with support for z-dimension and m-values.
Definition qgspoint.h:53
QgsPoint::setY
void setY(double y)
Sets the point's y-coordinate.
Definition qgspoint.h:370
QgsPoint::setX
void setX(double x)
Sets the point's x-coordinate.
Definition qgspoint.h:359
QgsPoint::z
double z
Definition qgspoint.h:58
QgsPoint::x
double x
Definition qgspoint.h:56
QgsPoint::isEmpty
bool isEmpty() const override
Returns true if the geometry is empty.
Definition qgspoint.cpp:760
QgsPoint::setZ
void setZ(double z)
Sets the point's z-coordinate.
Definition qgspoint.h:383
QgsPoint::y
double y
Definition qgspoint.h:57
Vector3D
Represents a 3D-Vector, capable of storing x, y and z-coordinates in double values.
Definition Vector3D.h:37
Vector3D::setX
void setX(double x)
Sets the x-component of the vector.
Definition Vector3D.h:106
Vector3D::getY
double getY() const
Returns the y-component of the vector.
Definition Vector3D.h:96
Vector3D::getX
double getX() const
Returns the x-component of the vector.
Definition Vector3D.h:91
Vector3D::setY
void setY(double y)
Sets the y-component of the vector.
Definition Vector3D.h:111
Vector3D::getZ
double getZ() const
Returns the z-component of the vector.
Definition Vector3D.h:101
Vector3D::setZ
void setZ(double z)
Sets the z-component of the vector.
Definition Vector3D.h:116
MathUtils::derVec
bool ANALYSIS_EXPORT derVec(const Vector3D *v1, const Vector3D *v2, Vector3D *result, double x, double y)
Calculates the z-component of a vector with coordinates 'x' and 'y'which is in the same tangent plane...
Definition MathUtils.cpp:440
MathUtils::faculty
int ANALYSIS_EXPORT faculty(int n)
Faculty function.
Definition MathUtils.cpp:177
MathUtils::BarycentricToXY
bool ANALYSIS_EXPORT BarycentricToXY(double u, double v, double w, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result)
Definition MathUtils.cpp:56
MathUtils::calcBarycentricCoordinates
bool ANALYSIS_EXPORT calcBarycentricCoordinates(double x, double y, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result)
Calculates the barycentric coordinates of a point (x,y) with respect to p1, p2, p3 and stores the thr...
Definition MathUtils.cpp:28
MathUtils::normalFromPoints
void ANALYSIS_EXPORT normalFromPoints(QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, Vector3D *vec)
Calculates the normal vector of the plane through the points p1, p2 and p3 and assigns the result to ...
Definition MathUtils.cpp:565
QgsDebugError
#define QgsDebugError(str)
Definition qgslogger.h:59
Generated on for QGIS API Documentation by doxygen 1.15.0