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 f848958..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_HANDLE(HYDROData_Bathymetry, HYDROData_IAltitudeObject)
 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
 
+inline double sqr( double x )
+{
+  return x*x;
+}
+
+HYDROData_Bathymetry::AltitudePoint::AltitudePoint( double x, double y, double z )
+{
+  X=x; Y=y; Z=z;
+}
+
+double HYDROData_Bathymetry::AltitudePoint::SquareDistance( const HYDROData_Bathymetry::AltitudePoint& p ) const
+{
+  double d = 0;
+  d += sqr( X - p.X );
+  d += sqr( Y - p.Y );
+  //d += sqr( Z - p.Z );
+  return d;
+}
 
 HYDROData_Bathymetry::HYDROData_Bathymetry()
 : HYDROData_IAltitudeObject()
@@ -111,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 )
@@ -160,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 );
@@ -321,7 +421,7 @@ bool interpolZtriangle(const gp_XY& point, vtkPolyData* delaunay2D, vtkIdList* t
   int nbPts = triangle->GetNumberOfIds();
   if (nbPts != 3)
     {
-      DEBTRACE("not a triangle ?");
+      //DEBTRACE("not a triangle ?");
       return false;
     }
   vtkIdType s[3];
@@ -341,7 +441,7 @@ bool interpolZtriangle(const gp_XY& point, vtkPolyData* delaunay2D, vtkIdList* t
   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]);
   if (det == 0)
     {
-      DEBTRACE("flat triangle ?");
+      //DEBTRACE("flat triangle ?");
       return false;
     }
 
@@ -372,9 +472,10 @@ double HYDROData_Bathymetry::GetAltitudeForPoint(const gp_XY& thePoint, int theM
   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");
       return aResAltitude;
@@ -427,7 +528,10 @@ double HYDROData_Bathymetry::GetAltitudeForPoint(const gp_XY& thePoint, int theM
               break;
             }
         }
-      if (!isInside) DEBTRACE("  point outside triangles, nearest z kept");
+      if (!isInside)
+      {
+          DEBTRACE("  point outside triangles, nearest z kept");
+      }
     }
   #endif
   return aResAltitude;
@@ -435,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();
@@ -515,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 )
@@ -552,6 +658,11 @@ bool HYDROData_Bathymetry::IsAltitudesInverted() const
   return aRes;
 }
 
+bool HYDROData_Bathymetry::ImportFromFile( const QString& theFileName )
+{
+  return ImportFromFiles(QStringList(theFileName));
+}
+
 bool HYDROData_Bathymetry::ImportFromFiles( const QStringList& theFileNames )
 {
   AltitudePoints AllPoints;
@@ -584,7 +695,7 @@ bool HYDROData_Bathymetry::ImportFromFiles( const QStringList& theFileNames )
   }
 
   // Convert from global to local CS
-  Handle_HYDROData_Document aDoc = HYDROData_Document::Document( myLab );
+  Handle(HYDROData_Document) aDoc = HYDROData_Document::Document( myLab );
   HYDROData_Bathymetry::AltitudePoints::iterator anIter = AllPoints.begin(), aLast = AllPoints.end();
   for ( ; anIter!=aLast; ++anIter )
   {
@@ -621,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 );
@@ -642,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 ) ||
@@ -760,11 +873,11 @@ bool HYDROData_Bathymetry::importFromASCFile( QFile&          theFile,
 }
 
 
-Handle_HYDROData_PolylineXY HYDROData_Bathymetry::CreateBoundaryPolyline() const
+Handle(HYDROData_PolylineXY) HYDROData_Bathymetry::CreateBoundaryPolyline() const
 {
   Handle(HYDROData_Document) aDocument = HYDROData_Document::Document( myLab );
-  Handle_HYDROData_PolylineXY aResult = 
-    Handle_HYDROData_PolylineXY::DownCast( aDocument->CreateObject( KIND_POLYLINEXY ) );
+  Handle(HYDROData_PolylineXY) aResult = 
+    Handle(HYDROData_PolylineXY)::DownCast( aDocument->CreateObject( KIND_POLYLINEXY ) );
 
   if( aResult.IsNull() )
     return aResult;