Salome HOME
test Strickler
[modules/hydro.git] / src / HYDROData / HYDROData_LandCoverMap.cxx
index 3cc968d5f179fef2cfce8ac62f39ff34b5434aab..18eaeee8f584d2a19cb444de85aa5a08529895a2 100644 (file)
@@ -23,6 +23,7 @@
 #include <HYDROData_ShapeFile.h>
 #include <HYDROData_Document.h>
 #include <HYDROData_StricklerTable.h>
+#include <HYDROData_ShapesTool.h>
 
 #include <BOPAlgo_BOP.hxx>
 #include <BOPAlgo_Builder.hxx>
 #include <TopExp_Explorer.hxx>
 #include <ShapeUpgrade_UnifySameDomain.hxx>
 #include <TopExp.hxx>
-#include <TDataStd_Real.hxx>
 #include <TopTools_IndexedMapOfShape.hxx>
 #include <ShapeBuild_ReShape.hxx>
 #include <ShapeFix_Shape.hxx>
 #include <BRepCheck_Shell.hxx>
 #include <BRepCheck_ListOfStatus.hxx>
 #include <TopTools_SequenceOfShape.hxx>
-#include <Handle_Geom_Curve.hxx>
-#include <Handle_Geom_Line.hxx>
-#include <Handle_Geom_TrimmedCurve.hxx>
+#include <Geom_Curve.hxx>
+#include <Geom_Line.hxx>
+#include <Geom_TrimmedCurve.hxx>
 #include <Geom_TrimmedCurve.hxx>
 #include <TopTools_DataMapOfShapeListOfShape.hxx>
 #include <NCollection_DoubleMap.hxx>
+#include <HYDROData_LCM_FaceClassifier.h>
+#include <QDir>
 
+#include <stdexcept>
 
 #include <QFile>
 #include <QString>
 #include <QTextStream>
 #include <QFileInfo>
 
+#define _DEVDEBUG_
+#include "HYDRO_trace.hxx"
+
 const char TELEMAC_FORMAT = 'f';
 const int TELEMAC_PRECISION = 3;
 
 
-IMPLEMENT_STANDARD_HANDLE(HYDROData_LandCoverMap, HYDROData_Entity)
 IMPLEMENT_STANDARD_RTTIEXT(HYDROData_LandCoverMap, HYDROData_Entity)
 
 /**
@@ -362,7 +367,8 @@ bool EdgeDiscretization( const TopoDS_Edge& theEdge,
 */
 bool HYDROData_LandCoverMap::ExportTelemac( const QString& theFileName,
                                             double theDeflection,
-                                            const Handle(HYDROData_StricklerTable)& theTable ) const
+                                            const Handle(HYDROData_StricklerTable)& theTable,
+                                            QString& statMessage) const
 {
   TopoDS_Shape aLandCoverMapShape = GetShape();
   TopTools_ListOfShape aListOfFaces;
@@ -413,7 +419,19 @@ bool HYDROData_LandCoverMap::ExportTelemac( const QString& theFileName,
 
   QFile aFile( theFileName );
   if( !aFile.open( QFile::WriteOnly | QFile::Text ) )
-    return false;
+  {
+    QString homeFilePath = QDir::home().absoluteFilePath( theFileName );
+    aFile.setFileName(homeFilePath);
+    if (aFile.open( QFile::WriteOnly | QFile::Text ) )
+      statMessage = "Telemac file have been exported to the home directory: " + homeFilePath;
+    else
+      return false;
+  }
+  else
+  {
+    QString absFilePath = QDir::current().absoluteFilePath( theFileName );
+    statMessage = "Telemac file have been exported to the current directory: " + absFilePath;
+  }
 
   QTextStream aStream( &aFile );
   aStream << "# nodes\n";
@@ -642,6 +660,7 @@ TopoDS_Shape HYDROData_LandCoverMap::MergeFaces( const TopTools_ListOfShape& the
                                                  TopTools_IndexedDataMapOfShapeListOfShape* theShHistory, 
                                                  double theTolerance)
 {
+  //DEBTRACE("MergeFaces");
   int anError;
   TopTools_ListIteratorOfListOfShape anIt;
   BOPCol_ListOfShape aLC;
@@ -826,6 +845,7 @@ bool HYDROData_LandCoverMap::LocalPartition( const TopoDS_Shape& theNewShape, co
     aShapesList.Append( anIt.Face() );
   aShapesList.Append( theNewShape );
 
+  //DEBTRACE("theNewType " << theNewType);
   if( aShapesList.Size()==1 && theNewShape.ShapeType()==TopAbs_FACE )
   {
     aNewFaces.Add( TopoDS::Face( theNewShape ), theNewType );
@@ -863,6 +883,7 @@ bool HYDROData_LandCoverMap::LocalPartition( const TopoDS_Shape& theNewShape, co
   for( ; aMIt.More(); aMIt.Next() )
   {
     //std::cout << "   " << aMIt.Value() << std::endl;
+    //DEBTRACE(aMIt.Value());
     int aKey = (int)(uintptr_t)aMIt.Value().TShape().operator->();
     aShapesFromNewFace.Add( aKey );
   }
@@ -872,11 +893,13 @@ bool HYDROData_LandCoverMap::LocalPartition( const TopoDS_Shape& theNewShape, co
   for( ; anIt.More(); anIt.Next() )
   {
     QString aSType = anIt.StricklerType();
+    //DEBTRACE(anIt.StricklerType() << " " << anIt.Face());
     //std::cout << "from " << anIt.Face() << ": " << anIt.StricklerType() << std::endl;
     TopTools_ListOfShape aModified = aBuilder.Modified( anIt.Face() );
     //
     if( aModified.Extent() == 0 )
       aModified.Append( anIt.Face() );
+    //DEBTRACE(anIt.StricklerType() << " " << anIt.Face());
 
     TopTools_ListIteratorOfListOfShape aMIt( aModified );
     for( ; aMIt.More(); aMIt.Next() )
@@ -889,14 +912,21 @@ bool HYDROData_LandCoverMap::LocalPartition( const TopoDS_Shape& theNewShape, co
       if( isFace && !isAlsoFromNew )
         aNewFaces.Add( TopoDS::Face( aShape ), aSType );
     }
+    //DEBTRACE(anIt.StricklerType() << " " << anIt.Face());
   }
-
   //     c. add the new shape if it is face with its type
   if( theNewShape.ShapeType()==TopAbs_FACE )
     aNewFaces.Add( TopoDS::Face( theNewShape ), theNewType );
+  //DEBTRACE(theNewShape << " " << theNewType);
   
   // convert map of shape to type to compound and list of types
   StoreLandCovers( aNewFaces );
+
+//  anIt.Init( *this );
+//  for( ; anIt.More(); anIt.Next() )
+//  {
+//     DEBTRACE(anIt.StricklerType() << " " << anIt.Face());
+//  }
   return true;
 }
 
@@ -908,6 +938,7 @@ void HYDROData_LandCoverMap::StoreLandCovers( const HYDROData_MapOfFaceToStrickl
 {
   TopTools_ListOfShape aListOfFaces;
 
+  //DEBTRACE("theMap.Extent() " << theMap.Extent());
   for( int i = 1; i <= theMap.Extent(); i++ )
   {
     TopoDS_Face aFace = theMap.FindKey(i);
@@ -917,6 +948,7 @@ void HYDROData_LandCoverMap::StoreLandCovers( const HYDROData_MapOfFaceToStrickl
   }
 
   TopTools_IndexedDataMapOfShapeListOfShape ShHistory;
+  ShHistory.Clear();
 
   TopoDS_Shape aResult;
   if( aListOfFaces.Extent() == 1 )
@@ -926,8 +958,10 @@ void HYDROData_LandCoverMap::StoreLandCovers( const HYDROData_MapOfFaceToStrickl
 
   //remove internal edges
   //if nothing changes => the result shape should be the same
-  //in any other case the history of modifications may be broken
-  aResult = RemoveInternal(aResult);
+  //hence the map will be empty
+
+  NCollection_IndexedDataMap<TopoDS_Face, TopoDS_Face> ShF2FHistory;
+  RemoveInternal(aResult, &ShF2FHistory);
 
   //one face => mark as unchanged
   if( aListOfFaces.Extent() == 1 )
@@ -939,21 +973,57 @@ void HYDROData_LandCoverMap::StoreLandCovers( const HYDROData_MapOfFaceToStrickl
   for( int i = 1; i <= theMap.Extent(); i++ )
   {
     TopoDS_Face aFF = theMap.FindKey(i);
+    //DEBTRACE("  --- " << aFF);
     if( aFF.IsNull() )
       continue;
-    TopTools_ListOfShape aLS = ShHistory.FindFromKey(aFF);
+    //DEBTRACE(ShHistory.IsEmpty());
+    //DEBTRACE(aFF.Checked());
+    TopTools_ListOfShape aLS;
+    try
+    {
+      aLS = ShHistory.FindFromKey(aFF); //TODO: bug to fix. Observed on an incomplete split of a face
+    }
+    catch (...)
+    {
+      DEBTRACE("TODO: bug to fix. Observed on an incomplete split of a face");
+      //continue; // No, keep aLS empty and propagate the type of the original face
+    }
     if (aLS.IsEmpty())
     {
-      QString SType = theMap.FindFromKey(aFF);
-      aChF2ST.Add(aFF, SType);
+      //DEBTRACE("--- aLS.IsEmpty()");
+      QString aSType = theMap.FindFromKey(aFF);
+      //DEBTRACE("  --- " << aSType.toStdString());
+      if (ShF2FHistory.Contains(aFF))
+        {
+          //DEBTRACE("ShF2FHistory.FindFromKey(aFF) " << ShF2FHistory.FindFromKey(aFF));
+          aChF2ST.Add(ShF2FHistory.FindFromKey(aFF), aSType);
+        }
+      else
+        {
+          //DEBTRACE("aFF " << aFF);
+          aChF2ST.Add(aFF, aSType);
+        }
     }
     else
     {
+      //DEBTRACE("--- !aLS.IsEmpty()");
       TopTools_ListIteratorOfListOfShape anIt(aLS);
       for (; anIt.More(); anIt.Next())
       {
         QString aSType = theMap.FindFromKey(aFF);
-        aChF2ST.Add(TopoDS::Face(anIt.Value()), aSType);
+        //DEBTRACE("  --- " << aSType.toStdString());
+        const TopoDS_Face& aMF = TopoDS::Face(anIt.Value()); 
+        //if (ShF2FHistory.Contains(aFF))
+        if (ShF2FHistory.Contains(aMF))
+          {
+            //DEBTRACE("ShF2FHistory.FindFromKey(aMF) " << ShF2FHistory.FindFromKey(aFF));
+            aChF2ST.Add(ShF2FHistory.FindFromKey(aMF), aSType);
+          }
+        else
+          {
+            //DEBTRACE("aMF " << aMF);
+            aChF2ST.Add(aMF, aSType);
+          }
       }
     }
   }
@@ -968,6 +1038,7 @@ void HYDROData_LandCoverMap::StoreLandCovers( const HYDROData_MapOfFaceToStrickl
     QString aST = ""; 
     if (aChF2ST.Contains(aFace))
       aST = aChF2ST.FindFromKey(aFace);
+    //DEBTRACE("aFace "  << aFace << " aST " << aST.toStdString());
     aSTypes << aST;
   }
   
@@ -975,6 +1046,7 @@ void HYDROData_LandCoverMap::StoreLandCovers( const HYDROData_MapOfFaceToStrickl
   int k = 0;
   foreach (QString aST, aSTypes)
   {
+    //DEBTRACE("aST " << aST.toStdString());
     aTypes->SetValue( k, HYDROData_Tool::toExtString( aST ) );
     k++;
   }
@@ -1012,7 +1084,7 @@ TopoDS_Face HYDROData_LandCoverMap::FindByPoint( const gp_Pnt2d& thePoint, QStri
 
 void Dump( const QString& theName, const QStringList& theList, QStringList& theLines )
 {
-  theLines.append( QString( "%1 = QStringList()" ).arg( theName ) );
+  theLines.append( QString( "%1 = []" ).arg( theName ) );
   foreach( QString anItem, theList )
     theLines.append( QString( "%1.append( u\"%2\" )" ).arg( theName ).arg( anItem ) );
 }
@@ -1039,66 +1111,79 @@ QStringList HYDROData_LandCoverMap::DumpToPython( const QString&       thePyScri
   HYDROData_Document::Document( myLab )->CollectQGISValues( anAttr, anAttrValues, aTypes );
   ExportDBF( aDbfFileName, anAttr, anAttrValues, aTypes );
 
-  aResList << QString( "%1.ImportSHP( '%2' )" ).
+  aResList << QString( "if not(%1.ImportSHP( '%2' )):" ).
     arg( aName ).arg( QFileInfo( aShpFileName ).fileName() );
+  aResList << QString( "  raise ValueError('problem while loading LandCoverMap shape')" );
 
   Dump( "attr_values", anAttrValues, aResList );
   Dump( "types", aTypes, aResList );
-  aResList << QString( "%1.ImportDBF( '%2', '%3', attr_values, types )" ).
+  aResList << QString( "if (%1.ImportDBF( '%2', '%3', attr_values, types ) != %1.DBFStatus_OK):" ).
     arg( aName ).arg( QFileInfo( aDbfFileName ).fileName() ).arg( anAttr );
+  aResList << QString( "  raise ValueError('problem while loading LandCoverMap data base')" );
 
   return aResList;
 }
 
-TopoDS_Shape HYDROData_LandCoverMap::RemoveInternal(const TopoDS_Shape& InSh)
+void HYDROData_LandCoverMap::RemoveInternal(TopoDS_Shape& ShToRebuild, NCollection_IndexedDataMap<TopoDS_Face, TopoDS_Face>* aF2FReplace)
 {
+  //DEBTRACE("RemoveInternal");
   //Shape must be topologically correct
-  TopExp_Explorer anExp(InSh, TopAbs_EDGE);
-  TopTools_ListOfShape anEdgesToRemove;
-
-  for(; anExp.More(); anExp.Next() )
+  TopExp_Explorer anExpF(ShToRebuild, TopAbs_FACE);
+  //
+  for(; anExpF.More(); anExpF.Next() )
   {
-    TopoDS_Edge CurEdge = TopoDS::Edge(anExp.Current());
-    if (CurEdge.Orientation() == TopAbs_INTERNAL)
-      anEdgesToRemove.Append(CurEdge);
+    TopoDS_Face CurFace = TopoDS::Face(anExpF.Current());
+    //
+    TopExp_Explorer anExp(CurFace, TopAbs_EDGE);
+    TopTools_ListOfShape anEdgesToRemove;
+    //
+    for(; anExp.More(); anExp.Next() )
+    {
+      TopoDS_Edge CurEdge = TopoDS::Edge(anExp.Current());
+      if (CurEdge.Orientation() == TopAbs_INTERNAL)
+        anEdgesToRemove.Append(CurEdge);
+    }
+    //
+    if (!anEdgesToRemove.IsEmpty())
+    {
+      Handle(ShapeBuild_ReShape) aReshape = new ShapeBuild_ReShape();
+      TopoDS_Shape OutF = aReshape->Apply(CurFace);
+      TopTools_ListIteratorOfListOfShape aIt(anEdgesToRemove);
+      for (; aIt.More(); aIt.Next()) 
+        aReshape->Remove(aIt.Value());
+      OutF = aReshape->Apply(CurFace);
+
+      Handle(ShapeFix_Shape) sfs = new ShapeFix_Shape;
+      sfs->Init(OutF);
+      sfs->Perform();
+      OutF = sfs->Shape();
+      aF2FReplace->Add(CurFace, TopoDS::Face(OutF));
+    }
   }
-  if (!anEdgesToRemove.IsEmpty())
+  //
+  Handle(ShapeBuild_ReShape) anExtReshape = new ShapeBuild_ReShape();
+  for (int i = 1; i <= aF2FReplace->Extent(); i++)
   {
-    Handle_ShapeBuild_ReShape aReshape = new ShapeBuild_ReShape();
-    TopoDS_Shape OutSh = aReshape->Apply(InSh);
-    TopTools_ListIteratorOfListOfShape aIt(anEdgesToRemove);
-    for (; aIt.More(); aIt.Next()) 
-      aReshape->Remove(aIt.Value());
-    OutSh = aReshape->Apply(InSh);
-
-    Handle(ShapeFix_Shape) sfs = new ShapeFix_Shape;
-    sfs->Init(OutSh);
-    sfs->Perform();
-    return sfs->Shape();
+    TopoDS_Face aFK = aF2FReplace->FindKey(i);
+    TopoDS_Face aFV = aF2FReplace->FindFromIndex(i);
+    anExtReshape->Replace(aFK, aFV);
+    ShToRebuild = anExtReshape->Apply(ShToRebuild);
   }
-  else
-    return InSh;
+//  for (int i = 1; i <= aF2FReplace->Extent(); i++)
+//    {
+//      DEBTRACE("aF2FReplace key,value " << aF2FReplace->FindKey(i) << " " << aF2FReplace->FindFromIndex(i));
+//    }
+
 }
 
 void HYDROData_LandCoverMap::SetTransparency( double theTransparency )
 {
-  Handle(TDataStd_Real) anAttr;
-  TDF_Label aLabel = myLab.FindChild( DataTag_Transparency );
-  if( !aLabel.FindAttribute( TDataStd_Real::GetID(), anAttr ) )
-    aLabel.AddAttribute( anAttr = new TDataStd_Real() );
-  anAttr->Set( theTransparency );
+  SetDouble( DataTag_Transparency, theTransparency );
 }
 
 double HYDROData_LandCoverMap::GetTransparency() const
 {
-  Handle(TDataStd_Real) anAttr;
-  TDF_Label aLabel = myLab.FindChild( DataTag_Transparency );
-  if( !aLabel.FindAttribute( TDataStd_Real::GetID(), anAttr ) )
-    return 0.5;
-
-  return anAttr->Get();
-  
+  return GetDouble( DataTag_Transparency, 0.25 );
 }
 
 bool HYDROData_LandCoverMap::ImportSHP( const QString& theSHPFileName,
@@ -1108,7 +1193,7 @@ bool HYDROData_LandCoverMap::ImportSHP( const QString& theSHPFileName,
   QStringList aPolyList; 
   TopTools_SequenceOfShape aFaces;
   int aSHapeType = -1;
-  int Stat = anImporter.ImportPolygons(theSHPFileName, aPolyList, aFaces, aSHapeType);
+  int Stat = anImporter.ImportPolygons(HYDROData_Document::Document(1), theSHPFileName, aPolyList, aFaces, aSHapeType);
   //
   if (Stat != 1)
     return false;
@@ -1149,7 +1234,7 @@ bool HYDROData_LandCoverMap::ExportSHP( const QString& theSHPFileName, bool bUse
 {
   HYDROData_ShapeFile anExporter;
   QStringList aList;
-  anExporter.Export(theSHPFileName, this, aList, bUseDiscr, theDefl );
+  anExporter.Export(HYDROData_Document::Document(1), theSHPFileName, this, aList, bUseDiscr, theDefl );
   if (aList.empty())
     return true;
   else 
@@ -1164,7 +1249,7 @@ bool HYDROData_LandCoverMap::CheckLinear()
   {
     TopoDS_Edge E = TopoDS::Edge(anEdgeEx.Current());
     double aFP, aLP;
-    Handle_Geom_Curve aCur = BRep_Tool::Curve(E, aFP, aLP);
+    Handle(Geom_Curve) aCur = BRep_Tool::Curve(E, aFP, aLP);
     Handle(Geom_Line) aLine = Handle(Geom_Line)::DownCast(aCur);
     if (aLine.IsNull())
     {
@@ -1181,3 +1266,69 @@ bool HYDROData_LandCoverMap::CheckLinear()
   }
   return true;
 }
+
+void HYDROData_LandCoverMap::UpdateLocalCS( double theDx, double theDy )
+{
+  TopoDS_Shape aShape = GetShape();
+  TopoDS_Shape aLocatedShape = HYDROData_ShapesTool::Translated( aShape, theDx, theDy, 0 );
+  SetShape( aLocatedShape );
+}
+
+void HYDROData_LandCoverMap::ClassifyPoints( const std::vector<gp_XY>& thePoints,
+                                             std::vector<std::set <QString> >& theTypes ) const
+{
+  HYDROData_LCM_FaceClassifier FC(this);
+  FC.Classify(thePoints, theTypes, NULL);
+}
+
+void HYDROData_LandCoverMap::ClassifyPoints( const std::vector<gp_XY>& thePoints,
+                                             const Handle(HYDROData_StricklerTable)& theTable,
+                                             std::vector<int>& theTypes ) const
+{
+  std::vector<std::set <QString> > types;
+  HYDROData_LCM_FaceClassifier FC(this);
+  FC.Classify(thePoints, types, NULL);
+
+  size_t n = types.size();
+  theTypes.resize( n );
+  for( size_t i=0; i<n; i++ )
+  {
+    const std::set<QString>& sub_types = types[i];
+    int aType = 0;
+    if( !sub_types.empty() )
+    {
+      const QString& aSType = *sub_types.begin();
+      aType = theTable->GetAttrValue( aSType ).toInt();
+    }
+    theTypes[i] = aType;
+  }
+}
+
+void HYDROData_LandCoverMap::ClassifyPoints( const std::vector<gp_XY>& thePoints,  
+  const Handle(HYDROData_StricklerTable)& theTable, 
+  std::vector<double>& theCoeffs, double DefValue, bool UseMax ) const
+{
+  std::vector<std::set <QString> > Types;
+  HYDROData_LCM_FaceClassifier FC(this);
+  FC.Classify(thePoints, Types, NULL);
+  theCoeffs.resize(thePoints.size());
+  for (size_t i = 0; i < Types.size(); i++)
+  {
+    const std::set<QString>& SStr = Types[i];
+    if (SStr.empty())
+      theCoeffs[i] = DefValue;
+    else
+    {
+      std::set<QString>::const_iterator it;
+      std::vector<double> C1;
+      for (it = SStr.begin(); it != SStr.end(); ++it)
+        C1.push_back(theTable->Get( *it, DefValue ));
+      double Val;
+      if (UseMax)
+        Val = *(std::max_element( C1.begin(), C1.end() ) );
+      else
+        Val = *(std::min_element( C1.begin(), C1.end() ) );
+      theCoeffs[i] = Val;
+    }
+  }
+}