-
-
Notifications
You must be signed in to change notification settings - Fork 3k
/
LinTriangleInterpolator.cpp
139 lines (109 loc) · 4.3 KB
/
LinTriangleInterpolator.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
/***************************************************************************
LinTriangleInterpolator.cc - description
-------------------
copyright : (C) 2004 by Marco Hugentobler
email : mhugent@geo.unizh.ch
***************************************************************************/
/***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************/
#include "LinTriangleInterpolator.h"
#include "qgslogger.h"
bool LinTriangleInterpolator::calcFirstDerX( double x, double y, Vector3D *vec )
{
if ( vec && mTIN )
{
QgsPoint pt1( 0, 0, 0 );
QgsPoint pt2( 0, 0, 0 );
QgsPoint pt3( 0, 0, 0 );
if ( !mTIN->getTriangle( x, y, &pt1, &pt2, &pt3 ) )
{
return false;//point outside the convex hull or numerical problems
}
vec->setX( 1.0 );
vec->setY( 0.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() ) ) );
return true;
}
else
{
QgsDebugMsg( "warning, null pointer" );
return false;
}
}
bool LinTriangleInterpolator::calcFirstDerY( double x, double y, Vector3D *vec )
{
if ( vec && mTIN )
{
QgsPoint pt1( 0, 0, 0 );
QgsPoint pt2( 0, 0, 0 );
QgsPoint pt3( 0, 0, 0 );
if ( !mTIN->getTriangle( x, y, &pt1, &pt2, &pt3 ) )
{
return false;
}
vec->setX( 0 );
vec->setY( 1.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() ) ) );
return true;
}
else
{
QgsDebugMsg( "warning, null pointer" );
return false;
}
}
bool LinTriangleInterpolator::calcNormVec( double x, double y, Vector3D *vec )
{
//calculate vector product of the two derivative vectors in x- and y-direction and set the length to 1
if ( vec && mTIN )
{
Vector3D vec1;
Vector3D vec2;
if ( !calcFirstDerX( x, y, &vec1 ) )
{return false;}
if ( !calcFirstDerY( x, y, &vec2 ) )
{return false;}
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
double absvec3 = std::sqrt( vec3.getX() * vec3.getX() + vec3.getY() * vec3.getY() + vec3.getZ() * vec3.getZ() );//length of vec3
vec->setX( vec3.getX() / absvec3 );//standardize vec3 and assign it to vec
vec->setY( vec3.getY() / absvec3 );
vec->setZ( vec3.getZ() / absvec3 );
return true;
}
else
{
QgsDebugMsg( "warning, null pointer" );
return false;
}
}
bool LinTriangleInterpolator::calcPoint( double x, double y, QgsPoint *point )
{
if ( point && mTIN )
{
QgsPoint pt1( 0, 0, 0 );
QgsPoint pt2( 0, 0, 0 );
QgsPoint pt3( 0, 0, 0 );
if ( !mTIN->getTriangle( x, y, &pt1, &pt2, &pt3 ) )
{
return false;//point is outside the convex hull or numerical problems
}
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() ) );
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() ) );
double c = pt1.z() - a * pt1.x() - b * pt1.y();
point->setX( x );
point->setY( y );
point->setZ( a * x + b * y + c );
return true;
}
else
{
QgsDebugMsg( "warning, null pointer" );
return false;
}
}