QGIS API Documentation  3.2.0-Bonn (bc43194)
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( "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( "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( "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( "warning, null pointer" );
732  }
733 }
734 #endif
735 
736 
NormVecDecorator
Decorator class which adds the functionality of estimating normals at the data points.
Definition: NormVecDecorator.h:37
Vector3D::setX
void setX(double x)
Sets the x-component of the vector.
Definition: Vector3D.h:103
QgsPoint::y
double y
Definition: qgspoint.h:42
CloughTocherInterpolator::lpoint3
QgsPoint lpoint3
Stores point3 of the last run.
Definition: CloughTocherInterpolator.h:95
QgsPoint::setZ
void setZ(double z)
Sets the point&#39;s z-coordinate.
Definition: qgspoint.h:237
CloughTocherInterpolator::der2Y
double der2Y
Derivative in y-direction at point2.
Definition: CloughTocherInterpolator.h:85
QgsDebugMsg
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
CloughTocherInterpolator::cp9
QgsPoint cp9
Control point 9.
Definition: CloughTocherInterpolator.h:63
Vector3D::standardise
void standardise()
Standardises the vector.
Definition: Vector3D.cpp:24
CloughTocherInterpolator::cp6
QgsPoint cp6
Control point 6.
Definition: CloughTocherInterpolator.h:57
CloughTocherInterpolator::setTriangulation
virtual void setTriangulation(NormVecDecorator *tin)
Definition: CloughTocherInterpolator.cpp:28
CloughTocherInterpolator::cp2
QgsPoint cp2
Control point 2.
Definition: CloughTocherInterpolator.h:49
MathUtils::BarycentricToXY
bool ANALYSIS_EXPORT BarycentricToXY(double u, double v, double w, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result)
Definition: MathUtils.cpp:50
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::calcPoint
bool calcPoint(double x, double y, QgsPoint &result SIP_OUT) override
Performs a linear interpolation in a triangle and assigns the x-,y- and z-coordinates to point...
Definition: CloughTocherInterpolator.cpp:154
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
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
CloughTocherInterpolator::cp14
QgsPoint cp14
Control point 14.
Definition: CloughTocherInterpolator.h:73
CloughTocherInterpolator::calcNormVec
bool calcNormVec(double x, double y, Vector3D *result SIP_OUT) override
Calculates the normal vector and assigns it to vec (not implemented at the moment) ...
Definition: CloughTocherInterpolator.cpp:45
CloughTocherInterpolator::cp8
QgsPoint cp8
Control point 8.
Definition: CloughTocherInterpolator.h:61
Vector3D::getZ
double getZ() const
Returns the z-component of the vector.
Definition: Vector3D.h:98
CloughTocherInterpolator::point3
QgsPoint point3
Third point of the triangle in x-,y-,z-coordinates.
Definition: CloughTocherInterpolator.h:45
CloughTocherInterpolator::cp10
QgsPoint cp10
Control point 10.
Definition: CloughTocherInterpolator.h:65
CloughTocherInterpolator::der3X
double der3X
Derivative in x-direction at point3.
Definition: CloughTocherInterpolator.h:87
CloughTocherInterpolator::cp13
QgsPoint cp13
Control point 13.
Definition: CloughTocherInterpolator.h:71
CloughTocherInterpolator::cp15
QgsPoint cp15
Control point 15.
Definition: CloughTocherInterpolator.h:75
Vector3D
Class Vector3D represents a 3D-Vector, capable to store x-,y- and z-coordinates in double values...
Definition: Vector3D.h:33
Vector3D::setZ
void setZ(double z)
Sets the z-component of the vector.
Definition: Vector3D.h:113
CloughTocherInterpolator::der1X
double der1X
Derivative in x-direction at point1.
Definition: CloughTocherInterpolator.h:79
CloughTocherInterpolator::cp1
QgsPoint cp1
Control point 1.
Definition: CloughTocherInterpolator.h:47
CloughTocherInterpolator::cp5
QgsPoint cp5
Control point 5.
Definition: CloughTocherInterpolator.h:55
CloughTocherInterpolator::cp12
QgsPoint cp12
Control point 12.
Definition: CloughTocherInterpolator.h:69
CloughTocherInterpolator::lpoint1
QgsPoint lpoint1
Stores point1 of the last run.
Definition: CloughTocherInterpolator.h:91
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. &#39;n&#39; is three in the cubical case...
Definition: CloughTocherInterpolator.cpp:33
QgsPoint
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
QgsPoint::setX
void setX(double x)
Sets the point&#39;s x-coordinate.
Definition: qgspoint.h:213
QgsPoint::setY
void setY(double y)
Sets the point&#39;s y-coordinate.
Definition: qgspoint.h:224
CloughTocherInterpolator::cp4
QgsPoint cp4
Control point 4.
Definition: CloughTocherInterpolator.h:53
NormVecDecorator::PointState
PointState
Enumeration for the state of a point. Normal means, that the point is not on a BreakLine, BreakLine means that the point is on a breakline (but not an end point of it) and EndPoint means, that it is an endpoint of a breakline.
Definition: NormVecDecorator.h:41
CloughTocherInterpolator::cp7
QgsPoint cp7
Control point 7.
Definition: CloughTocherInterpolator.h:59
Vector3D::getY
double getY() const
Returns the y-component of the vector.
Definition: Vector3D.h:93
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::der3Y
double der3Y
Derivative in y-direction at point3.
Definition: CloughTocherInterpolator.h:89
Vector3D::setY
void setY(double y)
Sets the y-component of the vector.
Definition: Vector3D.h:108
CloughTocherInterpolator::CloughTocherInterpolator
CloughTocherInterpolator()=default
Standard constructor.
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
QgsPoint::z
double z
Definition: qgspoint.h:43
Vector3D::getX
double getX() const
Returns the x-component of the vector.
Definition: Vector3D.h:88
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::der1Y
double der1Y
Derivative in y-direction at point1.
Definition: CloughTocherInterpolator.h:81
CloughTocherInterpolator::mEdgeTolerance
double mEdgeTolerance
Tolerance of the barycentric coordinates at the borders of the triangles (to prevent errors because o...
Definition: CloughTocherInterpolator.h:39
CloughTocherInterpolator::der2X
double der2X
Derivative in x-direction at point2.
Definition: CloughTocherInterpolator.h:83
CloughTocherInterpolator::cp11
QgsPoint cp11
Control point 11.
Definition: CloughTocherInterpolator.h:67
CloughTocherInterpolator::point2
QgsPoint point2
Second point of the triangle in x-,y-,z-coordinates.
Definition: CloughTocherInterpolator.h:43
CloughTocherInterpolator::cp16
QgsPoint cp16
Control point 16.
Definition: CloughTocherInterpolator.h:77
CloughTocherInterpolator::point1
QgsPoint point1
First point of the triangle in x-,y-,z-coordinates.
Definition: CloughTocherInterpolator.h:41
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 &#39;x&#39; and &#39;y&#39;which is in the same tangent plane...
Definition: MathUtils.cpp:505
MathUtils::faculty
int ANALYSIS_EXPORT faculty(int n)
Faculty function.
Definition: MathUtils.cpp:221
CloughTocherInterpolator::mTIN
NormVecDecorator * mTIN
Association with a triangulation object.
Definition: CloughTocherInterpolator.h:37
QgsPoint::x
double x
Definition: qgspoint.h:41
Generated on Sat Jun 23 2018 20:58:15 for QGIS API Documentation by   doxygen 1.8.13