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