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