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