QGIS API Documentation 3.41.0-Master (af5edcb665c)
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
26
31
32double CloughTocherInterpolator::calcBernsteinPoly( int n, int i, int j, int k, double u, double v, double w )
33{
34 if ( i < 0 || j < 0 || k < 0 )
35 {
36 QgsDebugError( QStringLiteral( "Invalid parameters for Bernstein poly calculation!" ) );
37 return 0;
38 }
39
40 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 ) );
41 return result;
42}
43
44static void normalize( QgsPoint &point )
45{
46 const double length = sqrt( pow( point.x(), 2 ) + pow( point.y(), 2 ) + pow( point.z(), 2 ) );
47 if ( length > 0 )
48 {
49 point.setX( point.x() / length );
50 point.setY( point.y() / length );
51 point.setZ( point.z() / length );
52 }
53}
54
55bool CloughTocherInterpolator::calcNormVec( double x, double y, QgsPoint &result )
56{
57 if ( !result.isEmpty() )
58 {
59 init( x, y );
60 QgsPoint barycoord( 0, 0, 0 ); //barycentric coordinates of (x,y) with respect to the triangle
61 QgsPoint endpointUXY( 0, 0, 0 ); //endpoint of the derivative in u-direction (in xy-coordinates)
62 const QgsPoint endpointV( 0, 0, 0 ); //endpoint of the derivative in v-direction (in barycentric coordinates)
63 QgsPoint endpointVXY( 0, 0, 0 ); //endpoint of the derivative in v-direction (in xy-coordinates)
64
65 //find out, in which subtriangle the point (x,y) is
66 //is the point in the first subtriangle (point1,point2,cp10)?
67 MathUtils::calcBarycentricCoordinates( x, y, &point1, &point2, &cp10, &barycoord );
68 if ( barycoord.x() >= -mEdgeTolerance && barycoord.x() <= ( 1 + mEdgeTolerance ) && barycoord.y() >= -mEdgeTolerance && barycoord.y() <= ( 1 + mEdgeTolerance ) )
69 {
70 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() );
71 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() );
72 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() );
73 MathUtils::BarycentricToXY( barycoord.x() + 1, barycoord.y() - 1, barycoord.z(), &point1, &point2, &cp10, &endpointUXY );
74 endpointUXY.setZ( 3 * ( zu - zv ) );
75 MathUtils::BarycentricToXY( barycoord.x(), barycoord.y() + 1, barycoord.z() - 1, &point1, &point2, &cp10, &endpointVXY );
76 endpointVXY.setZ( 3 * ( zv - zw ) );
77 const Vector3D v1( endpointUXY.x() - x, endpointUXY.y() - y, endpointUXY.z() );
78 const Vector3D v2( endpointVXY.x() - x, endpointVXY.y() - y, endpointVXY.z() );
79 result.setX( v1.getY() * v2.getZ() - v1.getZ() * v2.getY() );
80 result.setY( v1.getZ() * v2.getX() - v1.getX() * v2.getZ() );
81 result.setZ( v1.getX() * v2.getY() - v1.getY() * v2.getX() );
82 normalize( result );
83 return true;
84 }
85 //is the point in the second subtriangle (point2,point3,cp10)?
86 MathUtils::calcBarycentricCoordinates( x, y, &point2, &point3, &cp10, &barycoord );
87 if ( barycoord.x() >= -mEdgeTolerance && barycoord.x() <= ( 1 + mEdgeTolerance ) && barycoord.y() >= -mEdgeTolerance && barycoord.y() <= ( 1 + mEdgeTolerance ) )
88 {
89 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() );
90 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() );
91 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() );
92 MathUtils::BarycentricToXY( barycoord.x() + 1, barycoord.y() - 1, barycoord.z(), &point2, &point3, &cp10, &endpointUXY );
93 endpointUXY.setZ( 3 * ( zu - zv ) );
94 MathUtils::BarycentricToXY( barycoord.x(), barycoord.y() + 1, barycoord.z() - 1, &point2, &point3, &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 third subtriangle (point3,point1,cp10)?
105 MathUtils::calcBarycentricCoordinates( x, y, &point3, &point1, &cp10, &barycoord );
106 if ( barycoord.x() >= -mEdgeTolerance && barycoord.x() <= ( 1 + mEdgeTolerance ) && barycoord.y() >= -mEdgeTolerance && barycoord.y() <= ( 1 + mEdgeTolerance ) )
107 {
108 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() );
109 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() );
110 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() );
111 MathUtils::BarycentricToXY( barycoord.x() + 1, barycoord.y() - 1, barycoord.z(), &point3, &point1, &cp10, &endpointUXY );
112 endpointUXY.setZ( 3 * ( zu - zv ) );
113 MathUtils::BarycentricToXY( barycoord.x(), barycoord.y() + 1, barycoord.z() - 1, &point3, &point1, &cp10, &endpointVXY );
114 endpointVXY.setZ( 3 * ( zv - zw ) );
115 const Vector3D v1( endpointUXY.x() - x, endpointUXY.y() - y, endpointUXY.z() );
116 const Vector3D v2( endpointVXY.x() - x, endpointVXY.y() - y, endpointVXY.z() );
117 result.setX( v1.getY() * v2.getZ() - v1.getZ() * v2.getY() );
118 result.setY( v1.getZ() * v2.getX() - v1.getX() * v2.getZ() );
119 result.setZ( v1.getX() * v2.getY() - v1.getY() * v2.getX() );
120 normalize( result );
121 return true;
122 }
123
124 //the point is in none of the subtriangles, test if point has the same position as one of the vertices
125 if ( x == point1.x() && y == point1.y() )
126 {
127 result.setX( -der1X );
128 result.setY( -der1Y );
129 result.setZ( 1 );
130 normalize( result );
131 return true;
132 }
133 else if ( x == point2.x() && y == point2.y() )
134 {
135 result.setX( -der2X );
136 result.setY( -der2Y );
137 result.setZ( 1 );
138 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 );
146 normalize( result );
147 return true;
148 }
149
150 result.setX( 0 ); //return a vertical normal if failed
151 result.setY( 0 );
152 result.setZ( 1 );
153 return false;
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 else
467 {
468 QgsDebugError( QStringLiteral( "warning, null pointer" ) );
469 }
470}
471
472
473#if 0
474void CloughTocherInterpolator::init( double x, double y )//version which has unintended breaklines similar to the Coons interpolator
475{
476 Vector3D v1, v2, v3;//normals at the three data points
477 int ptn1, ptn2, ptn3;//numbers of the vertex points
478 NormVecDecorator::PointState state1, state2, state3;//states of the vertex points (Normal, BreakLine, EndPoint possible)
479
480 if ( mTIN )
481 {
482 mTIN->getTriangle( x, y, &point1, &ptn1, &v1, &state1, &point2, &ptn2, &v2, &state2, &point3, &ptn3, &v3, &state3 );
483
484 if ( point1 == lpoint1 && point2 == lpoint2 && point3 == lpoint3 )//if we are in the same triangle as at the last run, we can leave 'init'
485 {
486 return;
487 }
488
489 //calculate the partial derivatives at the data points
490 der1X = -v1.getX() / v1.getZ();
491 der1Y = -v1.getY() / v1.getZ();
492 der2X = -v2.getX() / v2.getZ();
493 der2Y = -v2.getY() / v2.getZ();
494 der3X = -v3.getX() / v3.getZ();
495 der3Y = -v3.getY() / v3.getZ();
496
497 //calculate the control points
498 cp1.setX( point1.getX() + ( point2.getX() - point1.getX() ) / 3 );
499 cp1.setY( point1.getY() + ( point2.getY() - point1.getY() ) / 3 );
500 cp1.setZ( point1.getZ() + ( cp1.getX() - point1.getX() )*der1X + ( cp1.getY() - point1.getY() )*der1Y );
501
502 cp2.setX( point2.getX() + ( point1.getX() - point2.getX() ) / 3 );
503 cp2.setY( point2.getY() + ( point1.getY() - point2.getY() ) / 3 );
504 cp2.setZ( point2.getZ() + ( cp2.getX() - point2.getX() )*der2X + ( cp2.getY() - point2.getY() )*der2Y );
505
506 cp9.setX( point2.getX() + ( point3.getX() - point2.getX() ) / 3 );
507 cp9.setY( point2.getY() + ( point3.getY() - point2.getY() ) / 3 );
508 cp9.setZ( point2.getZ() + ( cp9.getX() - point2.getX() )*der2X + ( cp9.getY() - point2.getY() )*der2Y );
509
510 cp16.setX( point3.getX() + ( point2.getX() - point3.getX() ) / 3 );
511 cp16.setY( point3.getY() + ( point2.getY() - point3.getY() ) / 3 );
512 cp16.setZ( point3.getZ() + ( cp16.getX() - point3.getX() )*der3X + ( cp16.getY() - point3.getY() )*der3Y );
513
514 cp14.setX( point3.getX() + ( point1.getX() - point3.getX() ) / 3 );
515 cp14.setY( point3.getY() + ( point1.getY() - point3.getY() ) / 3 );
516 cp14.setZ( point3.getZ() + ( cp14.getX() - point3.getX() )*der3X + ( cp14.getY() - point3.getY() )*der3Y );
517
518 cp6.setX( point1.getX() + ( point3.getX() - point1.getX() ) / 3 );
519 cp6.setY( point1.getY() + ( point3.getY() - point1.getY() ) / 3 );
520 cp6.setZ( point1.getZ() + ( cp6.getX() - point1.getX() )*der1X + ( cp6.getY() - point1.getY() )*der1Y );
521
522 //set x- and y-coordinates of the central point
523 cp10.setX( ( point1.getX() + point2.getX() + point3.getX() ) / 3 );
524 cp10.setY( ( point1.getY() + point2.getY() + point3.getY() ) / 3 );
525
526 //do the necessary adjustments in case of breaklines--------------------------------------------------------------------
527
528 //temporary normals and derivatives
529 double tmpx = 0;
530 double tmpy = 0;
531 Vector3D tmp( 0, 0, 0 );
532
533 //point1
534 if ( state1 == NormVecDecorator::BreakLine )
535 {
536 if ( mTIN->calcNormalForPoint( x, y, ptn1, &tmp ) )
537 {
538 tmpx = -tmp.getX() / tmp.getZ();
539 tmpy = -tmp.getY() / tmp.getZ();
540 if ( state3 == NormVecDecorator::Normal )
541 {
542 cp6.setZ( point1.getZ() + ( ( point3.getX() - point1.getX() ) / 3 )*tmpx + ( ( point3.getY() - point1.getY() ) / 3 )*tmpy );
543 }
544 if ( state2 == NormVecDecorator::Normal )
545 {
546 cp1.setZ( point1.getZ() + ( ( point2.getX() - point1.getX() ) / 3 )*tmpx + ( ( point2.getY() - point1.getY() ) / 3 )*tmpy );
547 }
548 }
549 }
550
551 if ( state2 == NormVecDecorator::Breakline )
552 {
553 if ( mTIN->calcNormalForPoint( x, y, ptn2, &tmp ) )
554 {
555 tmpx = -tmp.getX() / tmp.getZ();
556 tmpy = -tmp.getY() / tmp.getZ();
557 if ( state1 == NormVecDecorator::Normal )
558 {
559 cp2.setZ( point2.getZ() + ( ( point1.getX() - point2.getX() ) / 3 )*tmpx + ( ( point1.getY() - point2.getY() ) / 3 )*tmpy );
560 }
561 if ( state3 == NormVecDecorator::Normal )
562 {
563 cp9.setZ( point2.getZ() + ( ( point3.getX() - point2.getX() ) / 3 )*tmpx + ( ( point3.getY() - point2.getY() ) / 3 )*tmpy );
564 }
565 }
566 }
567
568 if ( state3 == NormVecDecorator::Breakline )
569 {
570 if ( mTIN->calcNormalForPoint( x, y, ptn3, &tmp ) )
571 {
572 tmpx = -tmp.getX() / tmp.getZ();
573 tmpy = -tmp.getY() / tmp.getZ();
574 if ( state1 == NormVecDecorator::Normal )
575 {
576 cp14.setZ( point3.getZ() + ( ( point1.getX() - point3.getX() ) / 3 )*tmpx + ( ( point1.getY() - point3.getY() ) / 3 )*tmpy );
577 }
578 if ( state2 == NormVecDecorator::Normal )
579 {
580 cp16.setZ( point3.getZ() + ( ( point2.getX() - point3.getX() ) / 3 )*tmpx + ( ( point2.getY() - point3.getY() ) / 3 )*tmpy );
581 }
582 }
583 }
584
585 //calculate cp3, cp 5 and cp15
586 Vector3D normal;//temporary normal for the vertices on breaklines
587
588 cp3.setX( point1.getX() + ( cp10.getX() - point1.getX() ) / 3 );
589 cp3.setY( point1.getY() + ( cp10.getY() - point1.getY() ) / 3 );
590 if ( state1 == NormVecDecorator::Breakline )
591 {
592 MathUtils::normalFromPoints( &point1, &cp1, &cp6, &normal );
593 //recalculate der1X and der1Y
594 der1X = -normal.getX() / normal.getZ();
595 der1Y = -normal.getY() / normal.getZ();
596 }
597
598 cp3.setZ( point1.getZ() + ( cp3.getX() - point1.getX() )*der1X + ( cp3.getY() - point1.getY() )*der1Y );
599
600
601
602 cp5.setX( point2.getX() + ( cp10.getX() - point2.getX() ) / 3 );
603 cp5.setY( point2.getY() + ( cp10.getY() - point2.getY() ) / 3 );
604 if ( state2 == NormVecDecorator::Breakline )
605 {
606 MathUtils::normalFromPoints( &point2, &cp9, &cp2, &normal );
607 //recalculate der2X and der2Y
608 der2X = -normal.getX() / normal.getZ();
609 der2Y = -normal.getY() / normal.getZ();
610 }
611
612 cp5.setZ( point2.getZ() + ( cp5.getX() - point2.getX() )*der2X + ( cp5.getY() - point2.getY() )*der2Y );
613
614
615 cp15.setX( point3.getX() + ( cp10.getX() - point3.getX() ) / 3 );
616 cp15.setY( point3.getY() + ( cp10.getY() - point3.getY() ) / 3 );
617 if ( state3 == NormVecDecorator::Breakline )
618 {
619 MathUtils::normalFromPoints( &point3, &cp14, &cp16, &normal );
620 //recalculate der3X and der3Y
621 der3X = -normal.getX() / normal.getZ();
622 der3Y = -normal.getY() / normal.getZ();
623 }
624
625 cp15.setZ( point3.getZ() + ( cp15.getX() - point3.getX() )*der3X + ( cp15.getY() - point3.getY() )*der3Y );
626
627
628 cp4.setX( ( point1.getX() + cp10.getX() + point2.getX() ) / 3 );
629 cp4.setY( ( point1.getY() + cp10.getY() + point2.getY() ) / 3 );
630
631 QgsPoint midpoint3( ( cp1.getX() + cp2.getX() ) / 2, ( cp1.getY() + cp2.getY() ) / 2, ( cp1.getZ() + cp2.getZ() ) / 2 );
632 Vector3D cp1cp2( cp2.getX() - cp1.getX(), cp2.getY() - cp1.getY(), cp2.getZ() - cp1.getZ() );
633 Vector3D odir3( 0, 0, 0 );//direction orthogonal to the line from point1 to point2
634 if ( ( point2.getY() - point1.getY() ) != 0 ) //avoid division through zero
635 {
636 odir3.setX( 1 );
637 odir3.setY( -( point2.getX() - point1.getX() ) / ( point2.getY() - point1.getY() ) );
638 odir3.setZ( ( der1X + der2X ) / 2 + ( der1Y + der2Y ) / 2 * odir3.getY() ); //take the linear interpolated cross-derivative
639 }
640 else
641 {
642 odir3.setY( 1 );
643 odir3.setX( -( point2.getY() - point1.getY() ) / ( point2.getX() - point1.getX() ) );
644 odir3.setZ( ( der1X + der2X ) / 2 * odir3.getX() + ( der1Y + der2Y ) / 2 );
645 }
646 Vector3D midpoint3cp4( 0, 0, 0 );
647 MathUtils::derVec( &cp1cp2, &odir3, &midpoint3cp4, cp4.getX() - midpoint3.getX(), cp4.getY() - midpoint3.getY() );
648 cp4.setZ( midpoint3.getZ() + midpoint3cp4.getZ() );
649
650 cp13.setX( ( point2.getX() + cp10.getX() + point3.getX() ) / 3 );
651 cp13.setY( ( point2.getY() + cp10.getY() + point3.getY() ) / 3 );
652 //find the point in the middle of the bezier curve between point2 and point3
653 QgsPoint midpoint1( ( cp9.getX() + cp16.getX() ) / 2, ( cp9.getY() + cp16.getY() ) / 2, ( cp9.getZ() + cp16.getZ() ) / 2 );
654 Vector3D cp9cp16( cp16.getX() - cp9.getX(), cp16.getY() - cp9.getY(), cp16.getZ() - cp9.getZ() );
655 Vector3D odir1( 0, 0, 0 );//direction orthogonal to the line from point2 to point3
656 if ( ( point3.getY() - point2.getY() ) != 0 ) //avoid division through zero
657 {
658 odir1.setX( 1 );
659 odir1.setY( -( point3.getX() - point2.getX() ) / ( point3.getY() - point2.getY() ) );
660 odir1.setZ( ( der2X + der3X ) / 2 + ( der2Y + der3Y ) / 2 * odir1.getY() ); //take the linear interpolated cross-derivative
661 }
662 else
663 {
664 odir1.setY( 1 );
665 odir1.setX( -( point3.getY() - point2.getY() ) / ( point3.getX() - point2.getX() ) );
666 odir1.setZ( ( der2X + der3X ) / 2 * odir1.getX() + ( der2Y + der3Y ) / 2 );
667 }
668 Vector3D midpoint1cp13( 0, 0, 0 );
669 MathUtils::derVec( &cp9cp16, &odir1, &midpoint1cp13, cp13.getX() - midpoint1.getX(), cp13.getY() - midpoint1.getY() );
670 cp13.setZ( midpoint1.getZ() + midpoint1cp13.getZ() );
671
672
673 cp11.setX( ( point3.getX() + cp10.getX() + point1.getX() ) / 3 );
674 cp11.setY( ( point3.getY() + cp10.getY() + point1.getY() ) / 3 );
675 //find the point in the middle of the bezier curve between point3 and point2
676 QgsPoint midpoint2( ( cp14.getX() + cp6.getX() ) / 2, ( cp14.getY() + cp6.getY() ) / 2, ( cp14.getZ() + cp6.getZ() ) / 2 );
677 Vector3D cp14cp6( cp6.getX() - cp14.getX(), cp6.getY() - cp14.getY(), cp6.getZ() - cp14.getZ() );
678 Vector3D odir2( 0, 0, 0 );//direction orthogonal to the line from point 1 to point3
679 if ( ( point3.getY() - point1.getY() ) != 0 ) //avoid division through zero
680 {
681 odir2.setX( 1 );
682 odir2.setY( -( point3.getX() - point1.getX() ) / ( point3.getY() - point1.getY() ) );
683 odir2.setZ( ( der3X + der1X ) / 2 + ( der3Y + der1Y ) / 2 * odir2.getY() ); //take the linear interpolated cross-derivative
684 }
685 else
686 {
687 odir2.setY( 1 );
688 odir2.setX( -( point3.getY() - point1.getY() ) / ( point3.getX() - point1.getX() ) );
689 odir2.setZ( ( der3X + der1X ) / 2 * odir2.getX() + ( der3Y + der1Y ) / 2 );
690 }
691 Vector3D midpoint2cp11( 0, 0, 0 );
692 MathUtils::derVec( &cp14cp6, &odir2, &midpoint2cp11, cp11.getX() - midpoint2.getX(), cp11.getY() - midpoint2.getY() );
693 cp11.setZ( midpoint2.getZ() + midpoint2cp11.getZ() );
694
695
696 cp7.setX( cp10.getX() + ( point1.getX() - cp10.getX() ) / 3 );
697 cp7.setY( cp10.getY() + ( point1.getY() - cp10.getY() ) / 3 );
698 //cp7 has to be in the same plane as cp4, cp3 and cp11
699 Vector3D cp4cp3( cp3.getX() - cp4.getX(), cp3.getY() - cp4.getY(), cp3.getZ() - cp4.getZ() );
700 Vector3D cp4cp11( cp11.getX() - cp4.getX(), cp11.getY() - cp4.getY(), cp11.getZ() - cp4.getZ() );
701 Vector3D cp4cp7( 0, 0, 0 );
702 MathUtils::derVec( &cp4cp3, &cp4cp11, &cp4cp7, cp7.getX() - cp4.getX(), cp7.getY() - cp4.getY() );
703 cp7.setZ( cp4.getZ() + cp4cp7.getZ() );
704
705 cp8.setX( cp10.getX() + ( point2.getX() - cp10.getX() ) / 3 );
706 cp8.setY( cp10.getY() + ( point2.getY() - cp10.getY() ) / 3 );
707 //cp8 has to be in the same plane as cp4, cp5 and cp13
708 Vector3D cp4cp5( cp5.getX() - cp4.getX(), cp5.getY() - cp4.getY(), cp5.getZ() - cp4.getZ() );
709 Vector3D cp4cp13( cp13.getX() - cp4.getX(), cp13.getY() - cp4.getY(), cp13.getZ() - cp4.getZ() );
710 Vector3D cp4cp8( 0, 0, 0 );
711 MathUtils::derVec( &cp4cp5, &cp4cp13, &cp4cp8, cp8.getX() - cp4.getX(), cp8.getY() - cp4.getY() );
712 cp8.setZ( cp4.getZ() + cp4cp8.getZ() );
713
714 cp12.setX( cp10.getX() + ( point3.getX() - cp10.getX() ) / 3 );
715 cp12.setY( cp10.getY() + ( point3.getY() - cp10.getY() ) / 3 );
716 //cp12 has to be in the same plane as cp13, cp15 and cp11
717 Vector3D cp13cp11( cp11.getX() - cp13.getX(), cp11.getY() - cp13.getY(), cp11.getZ() - cp13.getZ() );
718 Vector3D cp13cp15( cp15.getX() - cp13.getX(), cp15.getY() - cp13.getY(), cp15.getZ() - cp13.getZ() );
719 Vector3D cp13cp12( 0, 0, 0 );
720 MathUtils::derVec( &cp13cp11, &cp13cp15, &cp13cp12, cp12.getX() - cp13.getX(), cp12.getY() - cp13.getY() );
721 cp12.setZ( cp13.getZ() + cp13cp12.getZ() );
722
723 //cp10 has to be in the same plane as cp7, cp8 and cp12
724 Vector3D cp7cp8( cp8.getX() - cp7.getX(), cp8.getY() - cp7.getY(), cp8.getZ() - cp7.getZ() );
725 Vector3D cp7cp12( cp12.getX() - cp7.getX(), cp12.getY() - cp7.getY(), cp12.getZ() - cp7.getZ() );
726 Vector3D cp7cp10( 0, 0, 0 );
727 MathUtils::derVec( &cp7cp8, &cp7cp12, &cp7cp10, cp10.getX() - cp7.getX(), cp10.getY() - cp7.getY() );
728 cp10.setZ( cp7.getZ() + cp7cp10.getZ() );
729
730 lpoint1 = point1;
731 lpoint2 = point2;
732 lpoint3 = point3;
733 }
734
735 else
736 {
737 QgsDebugError( QStringLiteral( "warning, null pointer" ) );
738 }
739}
740#endif
CloughTocherInterpolator::der2X
double der2X
Derivative in x-direction at point2.
Definition CloughTocherInterpolator.h:85
CloughTocherInterpolator::setTriangulation
virtual void setTriangulation(NormVecDecorator *tin)
Definition CloughTocherInterpolator.cpp:27
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:32
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:55
CloughTocherInterpolator::cp6
QgsPoint cp6
Control point 6.
Definition CloughTocherInterpolator.h:59
CloughTocherInterpolator::cp15
QgsPoint cp15
Control point 15.
Definition CloughTocherInterpolator.h:77
CloughTocherInterpolator::CloughTocherInterpolator
CloughTocherInterpolator()=default
CloughTocherInterpolator::lpoint3
QgsPoint lpoint3
Stores point3 of the last run.
Definition CloughTocherInterpolator.h:97
CloughTocherInterpolator::point1
QgsPoint point1
First point of the triangle in x-,y-,z-coordinates.
Definition CloughTocherInterpolator.h:43
CloughTocherInterpolator::der1Y
double der1Y
Derivative in y-direction at point1.
Definition CloughTocherInterpolator.h:83
CloughTocherInterpolator::cp11
QgsPoint cp11
Control point 11.
Definition CloughTocherInterpolator.h:69
NormVecDecorator
Decorator class which adds the functionality of estimating normals at the data points.
Definition NormVecDecorator.h: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:245
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:41
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:738
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: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: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:560
QgsDebugError
#define QgsDebugError(str)
Definition qgslogger.h:38
Generated on Thu Dec 26 2024 02:19:01 for QGIS API Documentation by doxygen 1.9.8