Salome HOME
Merge branch 'BR_H2018_2' of https://codev-tuleap.cea.fr/plugins/git/salome/hydro...
[modules/hydro.git] / src / HYDROData / HYDROData_Bathymetry.cxx
old mode 100644 (file)
new mode 100755 (executable)
index 2e6ac1c..c274ff5
 #include <OSD_Timer.hxx>
 #endif
 
-#define _DEVDEBUG_
+//#define _DEVDEBUG_
 #include "HYDRO_trace.hxx"
 
-const int BLOCK_SIZE = 1000;
+const int BLOCK_SIZE = 10000;
 
 IMPLEMENT_STANDARD_RTTIEXT(HYDROData_Bathymetry, HYDROData_IAltitudeObject)
 
-//HYDROData_QuadtreeNode* HYDROData_Bathymetry::myQuadtree = 0;
-
+//int HYDROData_Bathymetry::myQuadTreeNumber = 0;
 std::map<int, HYDROData_QuadtreeNode*> HYDROData_Bathymetry::myQuadtrees;
 
 #ifndef LIGHT_MODE
+//int HYDROData_Bathymetry::myDelaunayNumber = 0;
 std::map<int, vtkPolyData*> HYDROData_Bathymetry::myDelaunay2D;
 #endif
 
@@ -128,6 +128,7 @@ void HYDROData_Bathymetry::SetAltitudePoints( const HYDROData_Bathymetry::Altitu
   // Save coordinates
   Handle(TDataStd_RealArray) aCoordsArray = 
     TDataStd_RealArray::Set( myLab.FindChild( DataTag_AltitudePoints ), 0, thePoints.size() * 3 - 1 );
+  aCoordsArray->SetID(TDataStd_RealArray::GetID());
 
   HYDROData_Bathymetry::AltitudePoints::const_iterator anIter = thePoints.begin(), aLast = thePoints.end();
   for ( int i = 0 ; anIter!=aLast; ++i, ++anIter )
@@ -177,98 +178,180 @@ HYDROData_Bathymetry::AltitudePoints HYDROData_Bathymetry::GetAltitudePoints(boo
 
 HYDROData_QuadtreeNode* HYDROData_Bathymetry::GetQuadtreeNodes() const
 {
-  TDF_Label aLabel = myLab.FindChild(DataTag_AltitudePoints, false);
-  if (aLabel.IsNull())
-    return 0;
-  int labkey = myLab.Tag();
-  //int altkey = aLabel.Tag();
-  //DEBTRACE("GetQuadtreeNodes this labkey altkey "<<this<<" "<<labkey<<" "<<altkey);
-  // if (myQuadtree->isEmpty() )
-  if (myQuadtrees.find(labkey) == myQuadtrees.end())
+  TDF_Label aLabel2 = myLab.FindChild(DataTag_Quadtree, false);
+  if (aLabel2.IsNull())
     {
-      //DEBTRACE("GetQuadtreeNodes init " << this << " " << labkey);
-      HYDROData_QuadtreeNode* aQuadtree = new HYDROData_QuadtreeNode(0, 30, 5, 0.);
-      myQuadtrees[labkey] = aQuadtree;
       TDF_Label aLabel = myLab.FindChild(DataTag_AltitudePoints, false);
       if (aLabel.IsNull())
         return 0;
 
-      Handle(TDataStd_RealArray) aCoordsArray;
-      if (!aLabel.FindAttribute(TDataStd_RealArray::GetID(), aCoordsArray))
-        return 0;
-
-      Nodes_3D* aListOfNodes = new Nodes_3D();
-
-      int index =0;
-      for (int i = aCoordsArray->Lower(), n = aCoordsArray->Upper(); i <= n;)
+      int aQuadTreeNumber = 0;
+      Handle(HYDROData_Document) aDocument = HYDROData_Document::Document( myLab );
+      if ( ! aDocument.IsNull() )
         {
-          if (i + 3 > n + 1)
-            break;
-
-          double x = aCoordsArray->Value(i++);
-          double y = aCoordsArray->Value(i++);
-          double z = aCoordsArray->Value(i++);
-          gpi_XYZ* aPoint = new gpi_XYZ(x, y, z, index);
-          index++;
-          aListOfNodes->push_back(aPoint);
+          aQuadTreeNumber = aDocument->GetCountQuadtree();
+          DEBTRACE("aQuadTreeNumber " << aQuadTreeNumber);
+          aQuadTreeNumber++;
+          aDocument->SetCountQuadtree(aQuadTreeNumber);
         }
-      //DEBTRACE("  GetQuadtreeNodes call setNodesAndCompute");
-      aQuadtree->setNodesAndCompute(aListOfNodes);
+      else
+          DEBTRACE("document.IsNull()");
+      DEBTRACE("compute Quadtree "<< aQuadTreeNumber);
+      HYDROData_QuadtreeNode* aQuadtree = ComputeQuadtreeNodes(aQuadTreeNumber);
       return aQuadtree;
     }
   else
-    return myQuadtrees[labkey];
+    {
+      Handle(TDataStd_Integer) aQuadtreeNum;
+      if ( aLabel2.FindAttribute( TDataStd_Integer::GetID(), aQuadtreeNum ) )
+        {
+          if (myQuadtrees.find(aQuadtreeNum->Get()) != myQuadtrees.end())
+            return myQuadtrees[aQuadtreeNum->Get()];
+          else
+            {
+              DEBTRACE("recompute Quadtree "<< aQuadtreeNum->Get());
+              HYDROData_QuadtreeNode* aQuadtree = ComputeQuadtreeNodes(aQuadtreeNum->Get());
+              return aQuadtree;
+            }
+        }
+      else DEBTRACE("no attribute TDataStd_Integer");
+    }
+  return 0;
 }
 
-#ifndef LIGHT_MODE
-vtkPolyData* HYDROData_Bathymetry::GetVtkDelaunay2D() const
+HYDROData_QuadtreeNode* HYDROData_Bathymetry::ComputeQuadtreeNodes( int key) const
 {
   TDF_Label aLabel = myLab.FindChild(DataTag_AltitudePoints, false);
   if (aLabel.IsNull())
     return 0;
-  int labkey = myLab.Tag();
-  //int altkey = aLabel.Tag();
-  //DEBTRACE("GetVtkDelaunay2D this labkey altkey "<<this<<" "<<labkey<<" "<<altkey);
-  if (myDelaunay2D.find(labkey) == myDelaunay2D.end())
+
+  Handle(TDataStd_RealArray) aCoordsArray;
+  if (!aLabel.FindAttribute(TDataStd_RealArray::GetID(), aCoordsArray))
+    return 0;
+
+  Handle(TDataStd_Integer) anAttr = TDataStd_Integer::Set( myLab.FindChild( DataTag_Quadtree ), key );
+  anAttr->SetID(TDataStd_Integer::GetID());
+  DEBTRACE("GetQuadtreeNodes init " << this << " " << key);
+  HYDROData_QuadtreeNode* aQuadtree = new HYDROData_QuadtreeNode(0, 30, 5, 0.);
+
+  Nodes_3D* aListOfNodes = new Nodes_3D();
+
+  int index =0;
+  for (int i = aCoordsArray->Lower(), n = aCoordsArray->Upper(); i <= n;)
     {
-      //DEBTRACE("GetVtkDelaunay2D init " << this << " " << labkey);
+      if (i + 3 > n + 1)
+        break;
+
+      double x = aCoordsArray->Value(i++);
+      double y = aCoordsArray->Value(i++);
+      double z = aCoordsArray->Value(i++);
+      gpi_XYZ* aPoint = new gpi_XYZ(x, y, z, index);
+      index++;
+      aListOfNodes->push_back(aPoint);
+    }
+  DEBTRACE("  GetQuadtreeNodes call setNodesAndCompute");
+  aQuadtree->setNodesAndCompute(aListOfNodes);
+
+  Handle(Message_ProgressIndicator) aZIProgress = HYDROData_Tool::GetZIProgress();
+  if ( aZIProgress && aZIProgress->UserBreak() ) {
+    return 0;
+  }
+
+  myQuadtrees[key] = aQuadtree;
+
+  return aQuadtree;
+}
 
+#ifndef LIGHT_MODE
+vtkPolyData* HYDROData_Bathymetry::GetVtkDelaunay2D() const
+{
+  TDF_Label aLabel2 = myLab.FindChild(DataTag_Delaunay, false);
+  if (aLabel2.IsNull())
+    {
       TDF_Label aLabel = myLab.FindChild(DataTag_AltitudePoints, false);
       if (aLabel.IsNull())
         return 0;
-      Handle(TDataStd_RealArray) aCoordsArray;
-      if (!aLabel.FindAttribute(TDataStd_RealArray::GetID(), aCoordsArray))
-        return 0;
 
-      vtkPoints *points = vtkPoints::New();
-      points->Allocate(aCoordsArray->Upper() +1);
-      for (int i = aCoordsArray->Lower(), n = aCoordsArray->Upper(); i <= n;)
+      int aDelaunayNumber = 0;
+      Handle(HYDROData_Document) aDocument = HYDROData_Document::Document( myLab );
+      if ( ! aDocument.IsNull() )
         {
-          if (i + 3 > n + 1)
-            break;
-          double x = aCoordsArray->Value(i++);
-          double y = aCoordsArray->Value(i++);
-          double z = aCoordsArray->Value(i++);
-          vtkIdType index = points->InsertNextPoint(x, y, z); // same index than in GetQuadtreeNodes
-          //DEBTRACE("  " << index);
+          aDelaunayNumber = aDocument->GetCountDelaunay();
+          DEBTRACE("aDelaunayNumber " << aDelaunayNumber);
+          aDelaunayNumber++;
+          aDocument->SetCountDelaunay(aDelaunayNumber);
         }
-      vtkPolyData* profile = vtkPolyData::New();
-      profile->SetPoints(points);
-      //DEBTRACE("Number of Points: "<< points->GetNumberOfPoints());
-
-      vtkDelaunay2D* delaunay2D = vtkDelaunay2D::New();
-      delaunay2D->SetInputData(profile);
-      delaunay2D->Update();
-      vtkPolyData* data = delaunay2D->GetOutput();
-      data->BuildLinks();
-      myDelaunay2D[labkey] = data;
+      else
+          DEBTRACE("document.IsNull()");
+      DEBTRACE("compute Delaunay "<< aDelaunayNumber);
+      vtkPolyData* data = ComputeVtkDelaunay2D(aDelaunayNumber);
       return data;
     }
   else
-    return myDelaunay2D[labkey];
+    {
+      Handle(TDataStd_Integer) aDelaunayNum;
+      if ( aLabel2.FindAttribute( TDataStd_Integer::GetID(), aDelaunayNum ) )
+        {
+          if (myDelaunay2D.find(aDelaunayNum->Get()) != myDelaunay2D.end())
+            return myDelaunay2D[aDelaunayNum->Get()];
+          else
+            {
+              DEBTRACE("recompute Delaunay "<< aDelaunayNum->Get());
+              vtkPolyData* data = ComputeVtkDelaunay2D(aDelaunayNum->Get());
+              return data;
+            }
+        }
+      else DEBTRACE("no attribute TDataStd_Integer");
+    }
+  return 0;
+}
+
+vtkPolyData* HYDROData_Bathymetry::ComputeVtkDelaunay2D(int key) const
+{
+  TDF_Label aLabel = myLab.FindChild(DataTag_AltitudePoints, false);
+  if (aLabel.IsNull())
+    return 0;
+
+  Handle(TDataStd_RealArray) aCoordsArray;
+  if (!aLabel.FindAttribute(TDataStd_RealArray::GetID(), aCoordsArray))
+    return 0;
 
+  HYDROData_Tool::SetTriangulationStatus(HYDROData_Tool::Running);
+
+  Handle(TDataStd_Integer) anAttr = TDataStd_Integer::Set( myLab.FindChild( DataTag_Delaunay ), key );
+  anAttr->SetID(TDataStd_Integer::GetID());
+  DEBTRACE("GetVtkDelaunay2D init " << this << " " << key);
+  vtkPoints *points = vtkPoints::New();
+  points->Allocate(aCoordsArray->Upper() +1);
+  for (int i = aCoordsArray->Lower(), n = aCoordsArray->Upper(); i <= n;)
+    {
+      if (i + 3 > n + 1)
+        break;
+      double x = aCoordsArray->Value(i++);
+      double y = aCoordsArray->Value(i++);
+      double z = aCoordsArray->Value(i++);
+      vtkIdType index = points->InsertNextPoint(x, y, z); // same index than in GetQuadtreeNodes
+      //DEBTRACE("  " << index);
+    }
+  vtkPolyData* profile = vtkPolyData::New();
+  profile->SetPoints(points);
+  DEBTRACE("Number of Points: "<< points->GetNumberOfPoints());
+
+  vtkDelaunay2D* delaunay2D = vtkDelaunay2D::New();
+  delaunay2D->SetInputData(profile);
+  delaunay2D->Update();
+  vtkPolyData* data = delaunay2D->GetOutput();
+  data->BuildLinks();
+  myDelaunay2D[key] = data;
+
+  HYDROData_Tool::SetTriangulationStatus(HYDROData_Tool::Finished);
+
+  return data;
 }
+
 #endif
+
+
 void HYDROData_Bathymetry::RemoveAltitudePoints()
 {
   TDF_Label aLabel = myLab.FindChild( DataTag_AltitudePoints, false );
@@ -384,16 +467,17 @@ bool interpolZtriangle(const gp_XY& point, vtkPolyData* delaunay2D, vtkIdList* t
 
 double HYDROData_Bathymetry::GetAltitudeForPoint(const gp_XY& thePoint, int theMethod) const
 {
-  //DEBTRACE("GetAltitudeForPoint p(" << thePoint.X() << ", " << thePoint.Y() << "), interpolation method: " << theMethod);
+  DEBTRACE("GetAltitudeForPoint p(" << thePoint.X() << ", " << thePoint.Y() << "), interpolation method: " << theMethod);
   double anInvalidAltitude = GetInvalidAltitude();
   double aResAltitude = anInvalidAltitude;
 
   // --- find the nearest point in the bathymetry cloud, with quadtree
+  Handle(Message_ProgressIndicator) aZIProgress = HYDROData_Tool::GetZIProgress();
 
   HYDROData_QuadtreeNode* aQuadtree = GetQuadtreeNodes();
-  if (!aQuadtree)
+  if (!aQuadtree || (aZIProgress && aZIProgress->UserBreak()))
     {
-      //DEBTRACE("  no Quadtree");
+      DEBTRACE("  no Quadtree");
       return aResAltitude;
     }
 
@@ -402,13 +486,13 @@ double HYDROData_Bathymetry::GetAltitudeForPoint(const gp_XY& thePoint, int theM
   while (dist2nodes.size() == 0)
     {
       aQuadtree->setPrecision(aQuadtree->getPrecision() *2);
-      //DEBTRACE("adjust precision to: " << aQuadtree->getPrecision());
+      DEBTRACE("adjust precision to: " << aQuadtree->getPrecision());
       aQuadtree->NodesAround(thePoint, dist2nodes, aQuadtree->getPrecision());
     }
   std::map<double, const gpi_XYZ*>::const_iterator it = dist2nodes.begin();
   aResAltitude = it->second->Z();
   int nodeIndex = it->second->getIndex();
-  //DEBTRACE("  number of points found: " << dist2nodes.size() << " nearest z: " << aResAltitude << " point index: " << nodeIndex);
+  DEBTRACE("  number of points found: " << dist2nodes.size() << " nearest z: " << aResAltitude << " point index: " << nodeIndex);
 
   // --- for coarse bathymetry clouds (when the TELEMAC mesh is more refined than the bathymetry cloud)
   //     interpolation is required.
@@ -430,7 +514,7 @@ double HYDROData_Bathymetry::GetAltitudeForPoint(const gp_XY& thePoint, int theM
       points->Allocate(64);
       aDelaunay2D->GetPointCells(nodeIndex, cells);
       vtkIdType nbCells = cells->GetNumberOfIds();
-      //DEBTRACE("  triangles on nearest point: " << nbCells);
+      DEBTRACE("  triangles on nearest point: " << nbCells);
       bool isInside = false;
       for (int i=0; i<nbCells; i++)
         {
@@ -440,13 +524,13 @@ double HYDROData_Bathymetry::GetAltitudeForPoint(const gp_XY& thePoint, int theM
           if (isInside)
             {
               aResAltitude = z;
-              //DEBTRACE("  interpolated z: " << z);
+              DEBTRACE("  interpolated z: " << z);
               break;
             }
         }
       if (!isInside)
       {
-      //    DEBTRACE("  point outside triangles, nearest z kept");
+          DEBTRACE("  point outside triangles, nearest z kept");
       }
     }
   #endif
@@ -455,13 +539,15 @@ double HYDROData_Bathymetry::GetAltitudeForPoint(const gp_XY& thePoint, int theM
 
 void HYDROData_Bathymetry::SetFilePath( const TCollection_AsciiString& theFilePath )
 {
-  TDataStd_AsciiString::Set( myLab.FindChild( DataTag_FilePath ), theFilePath );
+  Handle(TDataStd_AsciiString) anAttr = TDataStd_AsciiString::Set( myLab.FindChild( DataTag_FilePath ), theFilePath );
+  anAttr->SetID(TDataStd_AsciiString::GetID());
 }
 
 void HYDROData_Bathymetry::SetFilePaths( const QStringList& theFilePaths )
 {
   int i = 1;
   Handle_TDataStd_ExtStringArray TExtStrArr = TDataStd_ExtStringArray::Set( myLab.FindChild( DataTag_FilePaths ), 1, theFilePaths.size() );
+  TExtStrArr->SetID(TDataStd_ExtStringArray::GetID());
   foreach (QString filepath, theFilePaths)
   {
     std::string sstr = filepath.toStdString();
@@ -535,8 +621,8 @@ void HYDROData_Bathymetry::SetAltitudesInverted( const bool theIsInverted,
   if ( anIsAltitudesInverted == theIsInverted )
     return;
 
-  TDataStd_Integer::Set( myLab.FindChild( DataTag_AltitudesInverted ), (Standard_Integer)theIsInverted );
-
+  Handle(TDataStd_Integer) anAttr = TDataStd_Integer::Set( myLab.FindChild( DataTag_AltitudesInverted ), (Standard_Integer)theIsInverted );
+  anAttr->SetID(TDataStd_Integer::GetID());
   Changed( Geom_Z );
 
   if ( !theIsUpdate )
@@ -646,15 +732,17 @@ bool HYDROData_Bathymetry::importFromXYZFile( QFile&          theFile,
   bool anIsAltitudesInverted = IsAltitudesInverted();
   while ( !theFile.atEnd() )
   {
-    QString aLine = theFile.readLine().simplified();
-    if ( aLine.isEmpty() )
+    std::string aLine = theFile.readLine().simplified().toStdString();
+    if ( aLine.empty() )
       continue;
 
-    QStringList aValues = aLine.split( ' ', QString::SkipEmptyParts );
-    if ( aValues.length() < 3 )
+    HYDROData_Bathymetry::AltitudePoint aPoint;
+    if( sscanf( aLine.c_str(), "%lf %lf %lf", &aPoint.X, &aPoint.Y, &aPoint.Z )!=3 )
       return false;
 
-    HYDROData_Bathymetry::AltitudePoint aPoint;
+    /*QStringList aValues = aLine.split( ' ', QString::SkipEmptyParts );
+    if ( aValues.length() < 3 )
+      return false;
     
     QString anX = aValues.value( 0 );
     QString anY = aValues.value( 1 );
@@ -667,7 +755,7 @@ bool HYDROData_Bathymetry::importFromXYZFile( QFile&          theFile,
     aPoint.Z = aZ.toDouble( &isZOk );
 
     if ( !isXOk || !isYOk || !isZOk )
-      return false;
+      return false;*/
 
     if ( HYDROData_Tool::IsNan( aPoint.X ) || HYDROData_Tool::IsInf( aPoint.X ) ||
          HYDROData_Tool::IsNan( aPoint.Y ) || HYDROData_Tool::IsInf( aPoint.Y ) ||