Branch data Line data Source code
1 : : /***************************************************************************
2 : : LinTriangleInterpolator.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 "LinTriangleInterpolator.h"
18 : : #include "qgslogger.h"
19 : :
20 : 0 : bool LinTriangleInterpolator::calcFirstDerX( double x, double y, Vector3D *vec )
21 : : {
22 : :
23 : 0 : if ( vec && mTIN )
24 : : {
25 : 0 : QgsPoint pt1( 0, 0, 0 );
26 : 0 : QgsPoint pt2( 0, 0, 0 );
27 : 0 : QgsPoint pt3( 0, 0, 0 );
28 : :
29 : 0 : if ( !mTIN->triangleVertices( x, y, pt1, pt2, pt3 ) )
30 : : {
31 : 0 : return false;//point outside the convex hull or numerical problems
32 : : }
33 : :
34 : 0 : vec->setX( 1.0 );
35 : 0 : vec->setY( 0.0 );
36 : 0 : vec->setZ( ( pt1.z() * ( pt2.y() - pt3.y() ) + pt2.z() * ( pt3.y() - pt1.y() ) + pt3.z() * ( pt1.y() - pt2.y() ) ) / ( ( pt1.x() - pt2.x() ) * ( pt2.y() - pt3.y() ) - ( pt2.x() - pt3.x() ) * ( pt1.y() - pt2.y() ) ) );
37 : 0 : return true;
38 : 0 : }
39 : :
40 : : else
41 : : {
42 : 0 : QgsDebugMsg( QStringLiteral( "warning, null pointer" ) );
43 : 0 : return false;
44 : : }
45 : 0 : }
46 : :
47 : 0 : bool LinTriangleInterpolator::calcFirstDerY( double x, double y, Vector3D *vec )
48 : : {
49 : 0 : if ( vec && mTIN )
50 : : {
51 : 0 : QgsPoint pt1( 0, 0, 0 );
52 : 0 : QgsPoint pt2( 0, 0, 0 );
53 : 0 : QgsPoint pt3( 0, 0, 0 );
54 : :
55 : 0 : if ( !mTIN->triangleVertices( x, y, pt1, pt2, pt3 ) )
56 : : {
57 : 0 : return false;
58 : : }
59 : :
60 : 0 : vec->setX( 0 );
61 : 0 : vec->setY( 1.0 );
62 : 0 : vec->setZ( ( pt1.z() * ( pt2.x() - pt3.x() ) + pt2.z() * ( pt3.x() - pt1.x() ) + pt3.z() * ( pt1.x() - pt2.x() ) ) / ( ( pt1.y() - pt2.y() ) * ( pt2.x() - pt3.x() ) - ( pt2.y() - pt3.y() ) * ( pt1.x() - pt2.x() ) ) );
63 : 0 : return true;
64 : 0 : }
65 : :
66 : : else
67 : : {
68 : 0 : QgsDebugMsg( QStringLiteral( "warning, null pointer" ) );
69 : 0 : return false;
70 : : }
71 : 0 : }
72 : :
73 : 0 : bool LinTriangleInterpolator::calcNormVec( double x, double y, QgsPoint &vec )
74 : : {
75 : : //calculate vector product of the two derivative vectors in x- and y-direction and set the length to 1
76 : 0 : if ( mTIN )
77 : : {
78 : 0 : Vector3D vec1;
79 : 0 : Vector3D vec2;
80 : 0 : if ( !calcFirstDerX( x, y, &vec1 ) )
81 : 0 : {return false;}
82 : 0 : if ( !calcFirstDerY( x, y, &vec2 ) )
83 : 0 : {return false;}
84 : 0 : Vector3D vec3( vec1.getY()*vec2.getZ() - vec1.getZ()*vec2.getY(), vec1.getZ()*vec2.getX() - vec1.getX()*vec2.getZ(), vec1.getX()*vec2.getY() - vec1.getY()*vec2.getX() );//calculate vector product
85 : 0 : double absvec3 = std::sqrt( vec3.getX() * vec3.getX() + vec3.getY() * vec3.getY() + vec3.getZ() * vec3.getZ() );//length of vec3
86 : 0 : vec.setX( vec3.getX() / absvec3 );//standardize vec3 and assign it to vec
87 : 0 : vec.setY( vec3.getY() / absvec3 );
88 : 0 : vec.setZ( vec3.getZ() / absvec3 );
89 : 0 : return true;
90 : : }
91 : :
92 : : else
93 : : {
94 : 0 : QgsDebugMsg( QStringLiteral( "warning, null pointer" ) );
95 : 0 : return false;
96 : : }
97 : :
98 : 0 : }
99 : :
100 : 0 : bool LinTriangleInterpolator::calcPoint( double x, double y, QgsPoint &point )
101 : : {
102 : 0 : if ( mTIN )
103 : : {
104 : 0 : QgsPoint pt1( 0, 0, 0 );
105 : 0 : QgsPoint pt2( 0, 0, 0 );
106 : 0 : QgsPoint pt3( 0, 0, 0 );
107 : :
108 : 0 : if ( !mTIN->triangleVertices( x, y, pt1, pt2, pt3 ) )
109 : : {
110 : 0 : return false;//point is outside the convex hull or numerical problems
111 : : }
112 : :
113 : 0 : double a = ( pt1.z() * ( pt2.y() - pt3.y() ) + pt2.z() * ( pt3.y() - pt1.y() ) + pt3.z() * ( pt1.y() - pt2.y() ) ) / ( ( pt1.x() - pt2.x() ) * ( pt2.y() - pt3.y() ) - ( pt2.x() - pt3.x() ) * ( pt1.y() - pt2.y() ) );
114 : 0 : double b = ( pt1.z() * ( pt2.x() - pt3.x() ) + pt2.z() * ( pt3.x() - pt1.x() ) + pt3.z() * ( pt1.x() - pt2.x() ) ) / ( ( pt1.y() - pt2.y() ) * ( pt2.x() - pt3.x() ) - ( pt2.y() - pt3.y() ) * ( pt1.x() - pt2.x() ) );
115 : 0 : double c = pt1.z() - a * pt1.x() - b * pt1.y();
116 : :
117 : 0 : point.setX( x );
118 : 0 : point.setY( y );
119 : 0 : point.setZ( a * x + b * y + c );
120 : 0 : return true;
121 : 0 : }
122 : : else
123 : : {
124 : 0 : QgsDebugMsg( QStringLiteral( "warning, null pointer" ) );
125 : 0 : return false;
126 : : }
127 : :
128 : 0 : }
129 : :
130 : :
131 : :
132 : :
133 : :
134 : :
135 : :
136 : :
137 : :
138 : :
|
