mesh scalar renderer optimization (#40821)

src/core/mesh/qgsmeshlayerinterpolator.cpp

@@ -116,7 +116,7 @@ QgsRasterBlock *QgsMeshLayerInterpolator::block( int, const QgsRectangle &extent
     const QgsMeshFace &face = mTriangularMesh.triangles()[triangleIndex];
 
     const int v1 = face[0], v2 = face[1], v3 = face[2];
-    const QgsPoint p1 = vertices[v1], p2 = vertices[v2], p3 = vertices[v3];
+    const QgsPointXY &p1 = vertices[v1], &p2 = vertices[v2], &p3 = vertices[v3];
 
     const int nativeFaceIndex = mTriangularMesh.trianglesToNativeFaces()[triangleIndex];
     const bool isActive = mActiveFaceFlagValues.active( nativeFaceIndex );
@@ -131,6 +131,18 @@ QgsRasterBlock *QgsMeshLayerInterpolator::block( int, const QgsRectangle &extent
     int topLim, bottomLim, leftLim, rightLim;
     QgsMeshLayerUtils::boundingBoxToScreenRectangle( mContext.mapToPixel(), mOutputSize, bbox, leftLim, rightLim, topLim, bottomLim );
 
+    double value( 0 ), value1( 0 ), value2( 0 ), value3( 0 );
+    const int faceIdx = mTriangularMesh.trianglesToNativeFaces()[triangleIndex];
+
+    if ( mDataType == QgsMeshDatasetGroupMetadata::DataType::DataOnVertices )
+    {
+      value1 = mDatasetValues[v1];
+      value2 = mDatasetValues[v2];
+      value3 = mDatasetValues[v3];
+    }
+    else
+      value = mDatasetValues[faceIdx];
+
     // interpolate in the bounding box of the face
     for ( int j = topLim; j <= bottomLim; j++ )
     {
@@ -144,18 +156,17 @@ QgsRasterBlock *QgsMeshLayerInterpolator::block( int, const QgsRectangle &extent
                   p1,
                   p2,
                   p3,
-                  mDatasetValues[v1],
-                  mDatasetValues[v2],
-                  mDatasetValues[v3],
+                  value1,
+                  value2,
+                  value3,
                   p );
         else
         {
-          const int faceIdx = mTriangularMesh.trianglesToNativeFaces()[triangleIndex];
           val = QgsMeshLayerUtils::interpolateFromFacesData(
                   p1,
                   p2,
                   p3,
-                  mDatasetValues[faceIdx],
+                  value,
                   p
                 );
         }

src/core/mesh/qgsmeshlayerutils.cpp

@@ -245,23 +245,28 @@ static void lamTol( double &lam )
 static bool E3T_physicalToBarycentric( const QgsPointXY &pA, const QgsPointXY &pB, const QgsPointXY &pC, const QgsPointXY &pP,
                                        double &lam1, double &lam2, double &lam3 )
 {
-  if ( pA == pB || pA == pC || pB == pC )
-    return false; // this is not a valid triangle!
-
   // Compute vectors
-  QgsVector v0( pC - pA );
-  QgsVector v1( pB - pA );
-  QgsVector v2( pP - pA );
+  double xa = pA.x();
+  double ya = pA.y();
+  double v0x = pC.x() - xa ;
+  double v0y = pC.y() - ya ;
+  double v1x = pB.x() - xa ;
+  double v1y = pB.y() - ya ;
+  double v2x = pP.x() - xa ;
+  double v2y = pP.y() - ya ;
 
   // Compute dot products
-  double dot00 = v0 * v0;
-  double dot01 = v0 * v1;
-  double dot02 = v0 * v2;
-  double dot11 = v1 * v1;
-  double dot12 = v1 * v2;
+  double dot00 = v0x * v0x + v0y * v0y;
+  double dot01 = v0x * v1x + v0y * v1y;
+  double dot02 = v0x * v2x + v0y * v2y;
+  double dot11 = v1x * v1x + v1y * v1y;
+  double dot12 = v1x * v2x + v1y * v2y;
 
   // Compute barycentric coordinates
-  double invDenom = 1.0 / ( dot00 * dot11 - dot01 * dot01 );
+  double invDenom =  dot00 * dot11 - dot01 * dot01;
+  if ( invDenom == 0 )
+    return false;
+  invDenom = 1.0 / invDenom;
   lam1 = ( dot11 * dot02 - dot01 * dot12 ) * invDenom;
   lam2 = ( dot00 * dot12 - dot01 * dot02 ) * invDenom;
   lam3 = 1.0 - lam1 - lam2;