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