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