z interpolation on Delaunay triangulation: better trace

diff --git a/src/HYDROData/HYDROData_Bathymetry.cxx b/src/HYDROData/HYDROData_Bathymetry.cxx
index 68d14dafc10ae3e52d273e3724677c6d90775fa0..d3a6d8447ef41dcc1673082220297765626b5119 100644 (file)
--- a/src/HYDROData/HYDROData_Bathymetry.cxx
+++ b/src/HYDROData/HYDROData_Bathymetry.cxx
@@ -392,6 +392,7 @@ bool interpolZtriangle(const gp_XY& point, vtkPolyData* delaunay2D, vtkIdList* t
   //DEBTRACE("triangle node 0: " << v[0][0]  << " " << v[0][1] << " " << v[0][2]);
   //DEBTRACE("triangle node 1: " << v[1][0]  << " " << v[1][1] << " " << v[1][2]);
   //DEBTRACE("triangle node 2: " << v[2][0]  << " " << v[2][1] << " " << v[2][2]);
+
   // compute barycentric coordinates (https://en.wikipedia.org/wiki/Barycentric_coordinate_system)
   //     det = (y2-y3)(x1-x3)+(x3-x2)(y1-y3)
   double det = (v[1][1]-v[2][1])*(v[0][0]-v[2][0]) + (v[2][0]-v[1][0])*(v[0][1]-v[2][1]);
@@ -443,7 +444,6 @@ double HYDROData_Bathymetry::GetAltitudeForPoint(const gp_XY& thePoint) const
       DEBTRACE("adjust precision to: " << aQuadtree->getPrecision());
       aQuadtree->NodesAround(thePoint, dist2nodes, aQuadtree->getPrecision());
     }
-  //aQuadtree->NodesAround(thePoint, dist2nodes, 5.0);
   std::map<double, const gpi_XYZ*>::const_iterator it = dist2nodes.begin();
   aResAltitude = it->second->Z();
   int nodeIndex = it->second->getIndex();
@@ -463,25 +463,23 @@ double HYDROData_Bathymetry::GetAltitudeForPoint(const gp_XY& thePoint) const
       cells->Allocate(64);
       vtkIdList* points= vtkIdList::New();
       points->Allocate(64);
-      DEBTRACE("---");
       aDelaunay2D->GetPointCells(nodeIndex, cells);
       vtkIdType nbCells = cells->GetNumberOfIds();
-      DEBTRACE(nbCells);
+      DEBTRACE("  triangles on nearest point: " << nbCells);
+      bool isInside = false;
       for (int i=0; i<nbCells; i++)
         {
           aDelaunay2D->GetCellPoints(cells->GetId(i), points);
           double z = 0;
-          if (interpolZtriangle(thePoint, aDelaunay2D, points, z))
+          isInside = interpolZtriangle(thePoint, aDelaunay2D, points, z);
+          if (isInside)
             {
               aResAltitude = z;
-              DEBTRACE("interpolated z: " << z);
+              DEBTRACE("  interpolated z: " << z);
               break;
             }
-          else
-            {
-              DEBTRACE("point outside triangles, nearest z kept");
-            }
         }
+      if (!isInside) DEBTRACE("  point outside triangles, nearest z kept");
     }
   return aResAltitude;
 }