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