26
26
#include " qgsvector.h"
27
27
#include " qgspoint.h"
28
28
#include " qgspointxy.h"
29
-
30
- static void bbox2rect (
31
- const QgsMapToPixel &mtp,
32
- const QSize &outputSize,
33
- const QgsRectangle &bbox,
34
- int &leftLim, int &rightLim, int &topLim, int &bottomLim )
35
- {
36
- QgsPointXY ll = mtp.transform ( bbox.xMinimum (), bbox.yMinimum () );
37
- QgsPointXY ur = mtp.transform ( bbox.xMaximum (), bbox.yMaximum () );
38
- topLim = std::max ( int ( ur.y () ), 0 );
39
- bottomLim = std::min ( int ( ll.y () ), outputSize.height () - 1 );
40
- leftLim = std::max ( int ( ll.x () ), 0 );
41
- rightLim = std::min ( int ( ur.x () ), outputSize.width () - 1 );
42
- }
43
-
44
- static void lamTol ( double &lam )
45
- {
46
- const static double eps = 1e-6 ;
47
- if ( ( lam < 0.0 ) && ( lam > -eps ) )
48
- {
49
- lam = 0.0 ;
50
- }
51
- }
52
-
53
- static bool E3T_physicalToBarycentric ( const QgsPointXY &pA, const QgsPointXY &pB, const QgsPointXY &pC, const QgsPointXY &pP,
54
- double &lam1, double &lam2, double &lam3 )
55
- {
56
- if ( pA == pB || pA == pC || pB == pC )
57
- return false ; // this is not a valid triangle!
58
-
59
- // Compute vectors
60
- QgsVector v0 ( pC - pA );
61
- QgsVector v1 ( pB - pA );
62
- QgsVector v2 ( pP - pA );
63
-
64
- // Compute dot products
65
- double dot00 = v0 * v0;
66
- double dot01 = v0 * v1;
67
- double dot02 = v0 * v2;
68
- double dot11 = v1 * v1;
69
- double dot12 = v1 * v2;
70
-
71
- // Compute barycentric coordinates
72
- double invDenom = 1.0 / ( dot00 * dot11 - dot01 * dot01 );
73
- lam1 = ( dot11 * dot02 - dot01 * dot12 ) * invDenom;
74
- lam2 = ( dot00 * dot12 - dot01 * dot02 ) * invDenom;
75
- lam3 = 1.0 - lam1 - lam2;
76
-
77
- // Apply some tolerance to lam so we can detect correctly border points
78
- lamTol ( lam1 );
79
- lamTol ( lam2 );
80
- lamTol ( lam3 );
81
-
82
- // Return if POI is outside triangle
83
- if ( ( lam1 < 0 ) || ( lam2 < 0 ) || ( lam3 < 0 ) )
84
- {
85
- return false ;
86
- }
87
-
88
- return true ;
89
- }
90
-
91
- double QgsMeshLayerInterpolator::interpolateFromVerticesData ( const QgsPointXY &p1, const QgsPointXY &p2, const QgsPointXY &p3,
92
- double val1, double val2, double val3, const QgsPointXY &pt )
93
- {
94
- double lam1, lam2, lam3;
95
- if ( !E3T_physicalToBarycentric ( p1, p2, p3, pt, lam1, lam2, lam3 ) )
96
- return std::numeric_limits<double >::quiet_NaN ();
97
-
98
- return lam1 * val3 + lam2 * val2 + lam3 * val1;
99
- }
100
-
101
- double interpolateFromFacesData ( const QgsPointXY &p1, const QgsPointXY &p2, const QgsPointXY &p3, double val, const QgsPointXY &pt )
102
- {
103
- double lam1, lam2, lam3;
104
- if ( !E3T_physicalToBarycentric ( p1, p2, p3, pt, lam1, lam2, lam3 ) )
105
- return std::numeric_limits<double >::quiet_NaN ();
106
-
107
- return val;
108
- }
29
+ #include " qgsmeshlayerutils.h"
109
30
110
31
QgsMeshLayerInterpolator::QgsMeshLayerInterpolator ( const QgsTriangularMesh &m,
111
32
const QVector<double > &datasetValues, const QVector<bool > &activeFaceFlagValues,
@@ -167,16 +88,13 @@ QgsRasterBlock *QgsMeshLayerInterpolator::block( int, const QgsRectangle &extent
167
88
if ( !isActive )
168
89
continue ;
169
90
170
- QgsRectangle bbox;
171
- bbox.combineExtentWith ( p1.x (), p1.y () );
172
- bbox.combineExtentWith ( p2.x (), p2.y () );
173
- bbox.combineExtentWith ( p3.x (), p3.y () );
91
+ QgsRectangle bbox = QgsMeshLayerUtils::triangleBoundingBox ( p1, p2, p3 );
174
92
if ( !extent.intersects ( bbox ) )
175
93
continue ;
176
94
177
95
// Get the BBox of the element in pixels
178
96
int topLim, bottomLim, leftLim, rightLim;
179
- bbox2rect ( mContext .mapToPixel (), mOutputSize , bbox, leftLim, rightLim, topLim, bottomLim );
97
+ QgsMeshLayerUtils::boundingBoxToScreenRectangle ( mContext .mapToPixel (), mOutputSize , bbox, leftLim, rightLim, topLim, bottomLim );
180
98
181
99
// interpolate in the bounding box of the face
182
100
for ( int j = topLim; j <= bottomLim; j++ )
@@ -187,7 +105,7 @@ QgsRasterBlock *QgsMeshLayerInterpolator::block( int, const QgsRectangle &extent
187
105
double val;
188
106
const QgsPointXY p = mContext .mapToPixel ().toMapCoordinates ( k, j );
189
107
if ( mDataOnVertices )
190
- val = interpolateFromVerticesData (
108
+ val = QgsMeshLayerUtils:: interpolateFromVerticesData (
191
109
p1,
192
110
p2,
193
111
p3,
@@ -198,7 +116,7 @@ QgsRasterBlock *QgsMeshLayerInterpolator::block( int, const QgsRectangle &extent
198
116
else
199
117
{
200
118
int face = mTriangularMesh .trianglesToNativeFaces ()[i];
201
- val = interpolateFromFacesData (
119
+ val = QgsMeshLayerUtils:: interpolateFromFacesData (
202
120
p1,
203
121
p2,
204
122
p3,
0 commit comments