src/HYDROData/HYDROData_Tool.cxx

   1 // Copyright (C) 2014-2015  EDF-R&D
   2 // This library is free software; you can redistribute it and/or
   3 // modify it under the terms of the GNU Lesser General Public
   4 // License as published by the Free Software Foundation; either
   5 // version 2.1 of the License, or (at your option) any later version.
   6 //
   7 // This library is distributed in the hope that it will be useful,
   8 // but WITHOUT ANY WARRANTY; without even the implied warranty of
   9 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  10 // Lesser General Public License for more details.
  11 //
  12 // You should have received a copy of the GNU Lesser General Public
  13 // License along with this library; if not, write to the Free Software
  14 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  15 //
  16 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  17 //
  18
  19 #include <HYDROData_Tool.h>
  20 #include <HYDROData_ArtificialObject.h>
  21 #include <HYDROData_Document.h>
  22 #include <HYDROData_Entity.h>
  23 #include <HYDROData_Iterator.h>
  24 #include <HYDROData_NaturalObject.h>
  25 #include <HYDROData_ShapesGroup.h>
  26 #include <HYDROData_PolylineXY.h>
  27 #include <QColor>
  28 #include <QFile>
  29 #include <QStringList>
  30 #include <QTextStream>
  31 #include <QSet>
  32 #include <BRep_Tool.hxx>
  33 #include <BRepAdaptor_Surface.hxx>
  34 #include <BRepTopAdaptor_FClass2d.hxx>
  35 #include <ElSLib.hxx>
  36 #include <Geom_Curve.hxx>
  37 #include <gp_Pln.hxx>
  38 #include <Quantity_Color.hxx>
  39 #include <TopExp_Explorer.hxx>
  40 #include <TopoDS.hxx>
  41 #include <TopoDS_Face.hxx>
  42 #include <TopoDS_Shape.hxx>
  43 #include <TopoDS_Wire.hxx>
  44 #include <limits>
  45 #include <math.h>
  46
  47
  48 #include <BRepTools.hxx>
  49 #include <NCollection_Map.hxx>
  50 #include <TopTools_ShapeMapHasher.hxx>
  51 #include <BRep_Builder.hxx>
  52 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
  53 #include <TopExp.hxx>
  54 #include <NCollection_List.hxx>
  55 #include <TopTools_ListIteratorOfListOfShape.hxx>
  56
  57 #include <BRepBuilderAPI_MakeFace.hxx>
  58
  59 #include <TopExp_Explorer.hxx>
  60 #include <TopTools_ListIteratorOfListOfShape.hxx>
  61 #include <TopTools_HSequenceOfShape.hxx>
  62
  63 #include <BRep_Builder.hxx>
  64 #include <BRepAlgo_FaceRestrictor.hxx>
  65 #include <BRepCheck_Analyzer.hxx>
  66
  67 #include <ShapeAnalysis.hxx>
  68 #include <ShapeAnalysis_FreeBounds.hxx>
  69
  70 #include <HYDROData_PolylineXY.h>
  71 #include <HYDROData_Polyline3D.h>
  72 #include <HYDROData_Bathymetry.h>
  73
  74 #include <QFileInfo>
  75 #define BLOCK_SIZE 10000
  76
  77 HYDROData_Tool::ExecStatus HYDROData_Tool::myTriangulationStatus = ExecStatus::None;
  78
  79 static int aMaxNameId = INT_MAX;
  80 static int aMaxColorNb = 92000;
  81 void HYDROData_Tool::WriteStringsToFile( QFile&             theFile,
  82                                          const QStringList& theStrings,
  83                                          const QString&     theSep )
  84 {
  85   if ( !theFile.isOpen() || theStrings.isEmpty() )
  86     return;
  87
  88   QString aWriteStr = theStrings.join( theSep );
  89   if ( aWriteStr.isEmpty() )
  90     return;
  91
  92   QTextStream anOutStream( &theFile );
  93   anOutStream << aWriteStr.toUtf8() << theSep << theSep;
  94 }
  95
  96 bool HYDROData_Tool::ExtractGeneratedObjectName(const QString& theName, int& outValue, QString& thePrefName)
  97 {
  98   QStringList aLs = theName.split('_');
  99   bool ok = false;
 100   QString last = aLs.last();
 101   outValue = last.toInt(&ok);
 102   if (!ok)
 103     return false;
 104   int last_len = last.length();
 105   int total_len = theName.length();
 106   thePrefName = theName.left(total_len-last_len-1);
 107   return true;
 108 }
 109
 110 QString HYDROData_Tool::GenerateObjectName( const Handle(HYDROData_Document)& theDoc,
 111                                             const QString&                    thePrefix,
 112                                             const QStringList&                theUsedNames,
 113                                             const bool                        theIsTryToUsePurePrefix )
 114 {
 115   QStringList aNamesList( theUsedNames );
 116
 117   // Collect all used names in the document
 118   HYDROData_Iterator anIter( theDoc );
 119   for( ; anIter.More(); anIter.Next() )
 120   {
 121     Handle(HYDROData_Entity) anObject = anIter.Current();
 122     if( anObject.IsNull() )
 123       continue;
 124
 125     QString anObjName = anObject->GetName();
 126     if ( anObjName.isEmpty() )
 127       continue;
 128
 129     aNamesList.append( anObjName );
 130   }
 131
 132   QString aName;
 133
 134   if ( theIsTryToUsePurePrefix && !aNamesList.contains( thePrefix ) ) {
 135     aName = thePrefix;
 136   } else {
 137     int anId = 1;
 138     while( anId < aMaxNameId )
 139     {
 140       aName = QString( "%1_%2" ).arg( thePrefix ).arg( QString::number( anId++ ) );
 141
 142       // check that there are no other objects with the same name in the document
 143       if ( !aNamesList.contains( aName ) )
 144         break;
 145     }
 146   }
 147
 148   return aName;
 149 }
 150
 151 bool HYDROData_Tool::IsGeometryObject( const Handle(HYDROData_Entity)& theObject )
 152 {
 153   if ( theObject.IsNull() )
 154     return false;
 155
 156   return theObject->IsKind( STANDARD_TYPE(HYDROData_ArtificialObject) ) ||
 157          theObject->IsKind( STANDARD_TYPE(HYDROData_NaturalObject) );
 158 }
 159
 160 void HYDROData_Tool::UpdateChildObjectName( const QString&                  theOldStr,
 161                                             const QString&                  theNewStr,
 162                                             const Handle(HYDROData_Entity)& theObject )
 163 {
 164   if ( theObject.IsNull() )
 165     return;
 166
 167   QString anObjName = theObject->GetName();
 168   if ( theOldStr.isEmpty() )
 169   {
 170     while ( anObjName.startsWith( '_' ) )
 171       anObjName.remove( 0, 1 );
 172
 173     anObjName.prepend( theNewStr + "_" );
 174   }
 175   else if ( anObjName.startsWith( theOldStr ) )
 176   {
 177     anObjName.replace( 0, theOldStr.length(), theNewStr );
 178   }
 179   else
 180     return;
 181
 182   theObject->SetName( anObjName );
 183 }
 184
 185 QString HYDROData_Tool::GenerateNameForPython( const MapOfTreatedObjects& theTreatedObjects,
 186                                                const QString&             thePrefix )
 187 {
 188   QString aName = thePrefix;
 189   if ( !theTreatedObjects.contains( aName ) )
 190     return aName;
 191
 192   int anId = 1;
 193   while( anId < aMaxNameId )
 194   {
 195     aName = QString( "%1_%2" ).arg( thePrefix ).arg( QString::number( anId++ ) );
 196
 197     // check that there are no other objects with the same name
 198     if ( !theTreatedObjects.contains( aName ) )
 199       break;
 200   }
 201
 202   return aName;
 203 }
 204 //======================================================================================================
 205 TopAbs_State HYDROData_Tool::ComputePointState( const gp_XY& theXY, const TopoDS_Face& theFace )
 206 {
 207   TopAbs_State aState(TopAbs_UNKNOWN);
 208   if(theFace.IsNull()) return aState;
 209   Standard_Real aTol = BRep_Tool::Tolerance(theFace);
 210   BRepAdaptor_Surface Ads ( theFace, Standard_False );
 211   Standard_Real toluv = Min ( Ads.UResolution(aTol), Ads.VResolution(aTol) );
 212   if (toluv < 0.01)
 213     toluv = 0.01; // there is no need to be more precise than 1cm a any case !
 214                   // another solution could be to compute a tolerance related to the distance between the border nodes
 215   const gp_Pln& aPlane = Ads.Surface().Plane();
 216   gp_Pnt aPnt(theXY.X(), theXY.Y(), 0.);
 217   Standard_Real aU1, aV1;
 218   ElSLib::Parameters(aPlane,aPnt, aU1, aV1);
 219   BRepTopAdaptor_FClass2d aClassifier( theFace, toluv );
 220   aState = aClassifier.Perform( gp_Pnt2d(aU1, aV1), Standard_False );
 221   return aState;
 222 }
 223
 224 double HYDROData_Tool::GetAltitudeForEdge( const TopoDS_Edge& theEdge,
 225                                            const gp_XY& thePoint,
 226                                            double theParameterTolerance,
 227                                            double theSquareDistanceTolerance,
 228                                            double theInvalidAltitude )
 229 {
 230   double aFirst, aLast;
 231   Handle(Geom_Curve) aCurve = BRep_Tool::Curve( theEdge, aFirst, aLast );
 232   if( aCurve.IsNull() )
 233     return theInvalidAltitude;
 234
 235   gp_Pnt aFirstPnt, aLastPnt;
 236
 237   aCurve->D0( aFirst, aFirstPnt );
 238   aCurve->D0( aLast, aLastPnt );
 239
 240   gp_Pnt2d aFirstPnt2d( aFirstPnt.X(), aFirstPnt.Y() );
 241   gp_Pnt2d aLastPnt2d( aLastPnt.X(), aLastPnt.Y() );
 242
 243   double aFirstDist = 0;
 244   double aLastDist = aFirstPnt2d.SquareDistance( aLastPnt2d );
 245   double aNecDist = aFirstPnt2d.SquareDistance( thePoint );
 246
 247   while( fabs( aLast - aFirst ) > theParameterTolerance )
 248   {
 249     double aMid = ( aFirst + aLast ) / 2;
 250     gp_Pnt aMidPnt;
 251     aCurve->D0( aMid, aMidPnt );
 252     double aDist = aFirstPnt2d.SquareDistance( gp_Pnt2d( aMidPnt.X(), aMidPnt.Y() ) );
 253
 254     if( aDist < aNecDist )
 255       aFirst = aMid;
 256     else
 257       aLast = aMid;
 258   }
 259
 260   double aMid = ( aFirst + aLast ) / 2;
 261   gp_Pnt aMidPnt;
 262   aCurve->D0( aMid, aMidPnt );
 263
 264   gp_Pnt2d aMidPnt2d( aMidPnt.X(), aMidPnt.Y() );
 265   if( aMidPnt2d.SquareDistance( thePoint ) < theSquareDistanceTolerance )
 266     return aMidPnt.Z();
 267   else
 268     return theInvalidAltitude;
 269 }
 270
 271 double HYDROData_Tool::GetAltitudeForWire( const TopoDS_Wire& theWire,
 272                                            const gp_XY& thePoint,
 273                                            double theParameterTolerance,
 274                                            double theSquareDistanceTolerance,
 275                                            double theInvalidAltitude )
 276 {
 277   TopExp_Explorer anExp( theWire, TopAbs_EDGE );
 278   for( ; anExp.More(); anExp.Next() )
 279   {
 280     double anAltitude = GetAltitudeForEdge( TopoDS::Edge( anExp.Current() ), thePoint,
 281       theParameterTolerance, theSquareDistanceTolerance, theInvalidAltitude );
 282     if( anAltitude != theInvalidAltitude )
 283       return anAltitude;
 284   }
 285   return theInvalidAltitude;
 286 }
 287
 288 TopoDS_Shape HYDROData_Tool::getFirstShapeFromGroup( const HYDROData_SequenceOfObjects& theGroups,
 289                                                      const int                          theGroupId )
 290 {
 291   TopoDS_Shape aResShape;
 292   if ( theGroupId < 1 || theGroupId > theGroups.Length() )
 293     return aResShape;
 294
 295   Handle(HYDROData_ShapesGroup) aGroup =
 296     Handle(HYDROData_ShapesGroup)::DownCast( theGroups.Value( theGroupId ) );
 297   if ( aGroup.IsNull() )
 298     return aResShape;
 299
 300   TopTools_SequenceOfShape aGroupShapes;
 301   aGroup->GetShapes( aGroupShapes );
 302
 303   if ( !aGroupShapes.IsEmpty() )
 304     aResShape = aGroupShapes.First();
 305
 306   return aResShape;
 307 }
 308
 309 TCollection_ExtendedString HYDROData_Tool::toExtString( const QString& theStr )
 310 {
 311   TCollection_ExtendedString aRes;
 312   if( !theStr.isEmpty() )
 313   {
 314           Standard_ExtString extStr = new Standard_ExtCharacter[ ( theStr.length() + 1 ) * 2 ];
 315           memcpy( (void*)extStr, theStr.unicode(), theStr.length() * 2 );
 316           ((short*)extStr)[theStr.length()] = '\0';
 317     aRes = TCollection_ExtendedString( extStr );
 318           delete [] extStr;
 319   }
 320   return aRes;
 321 }
 322
 323 QString HYDROData_Tool::toQString( const TCollection_ExtendedString& theStr )
 324 {
 325   return QString( (QChar*)theStr.ToExtString(), theStr.Length() );
 326 }
 327
 328 Quantity_Color HYDROData_Tool::toOccColor( const QColor& theColor )
 329 {
 330   double r = theColor.red() / 255.0;
 331   double g = theColor.green() / 255.0;
 332   double b = theColor.blue() / 255.0;
 333
 334   return Quantity_Color( r, g, b, Quantity_TOC_RGB );
 335 }
 336
 337 QColor HYDROData_Tool::toQtColor( const Quantity_Color& theColor )
 338 {
 339   int r = 255 * theColor.Red();
 340   int g = 255 * theColor.Green();
 341   int b = 255 * theColor.Blue();
 342   return QColor( r, g, b );
 343 }
 344
 345 QColor HYDROData_Tool::GenerateRandColor()
 346 {
 347   float aHue = ( rand()%1000 ) * 0.001f;
 348
 349   QColor aColor;
 350   aColor.setHsl( (int)(aHue*255.), 200, 128 );
 351   int r = aColor.red();
 352   int g = aColor.green();
 353   int b = aColor.blue();
 354   return ( aColor.isValid() ? aColor : Qt::darkBlue );
 355 }
 356
 357 void HYDROData_Tool::GenerateRepeatableRandColors(int nbColorsToGen, QVector<QColor>& theColors)
 358 {
 359   for (int i = 1; i <= nbColorsToGen; i++)
 360     theColors.append(HYDROData_Tool::GenerateRandColor());
 361 }
 362
 363 bool HYDROData_Tool::GenerateNonRepeatableRandColors(int nbColorsToGen, QVector<QColor>& theColors)
 364 {
 365   if (nbColorsToGen > aMaxColorNb)
 366     return false;
 367   QSet<int> Codes;
 368   float aHue;
 369   theColors.clear();
 370   do
 371   {
 372     QColor aColor;
 373     int H = rand()%360;
 374     int S = rand()%256;
 375     int Code = H*256+S;
 376     if (Codes.contains(Code))
 377       continue;
 378     aColor.setHsl( H, S, 128 );
 379     if (aColor.isValid())
 380     {
 381       theColors.append(aColor);
 382       Codes.insert(Code);
 383     }
 384   } while (theColors.size() <= nbColorsToGen);
 385   return true;
 386 }
 387
 388
 389 bool HYDROData_Tool::IsNan( double theValue )
 390 {
 391 #ifdef WIN32
 392   return _isnan( theValue );
 393 #else
 394   return isnan( theValue );
 395 #endif
 396 }
 397
 398 bool HYDROData_Tool::IsInf( double theValue )
 399 {
 400 #ifdef WIN32
 401   return (!_finite( theValue  ) );
 402 #else
 403   return isinf( theValue );
 404 #endif
 405 }
 406
 407 static void MakeShellG(const NCollection_Map<TopoDS_Face, TopTools_ShapeMapHasher>& FG,
 408   TopoDS_Shape& outSh)
 409 {
 410   BRep_Builder bb;
 411   NCollection_Map<TopoDS_Face, TopTools_ShapeMapHasher>::Iterator itFG(FG);
 412   if (FG.Extent() > 1)
 413   {
 414     //face nb > 1 => make shell
 415     TopoDS_Shell outShell;
 416     bb.MakeShell(outShell);
 417     for (;itFG.More();itFG.Next())
 418       bb.Add(outShell, itFG.Value());
 419     outSh = outShell;
 420   }
 421   else if (FG.Extent() == 1)
 422   {
 423     outSh = itFG.Value(); //one face
 424   }
 425 }
 426
 427 TopoDS_Shape HYDROData_Tool::RebuildCmp(const TopoDS_Shape& in)
 428 {
 429   TopTools_IndexedDataMapOfShapeListOfShape mE2LF;
 430   TopExp::MapShapesAndAncestors(in, TopAbs_EDGE, TopAbs_FACE, mE2LF);
 431   if (mE2LF.IsEmpty())
 432     return TopoDS_Shape();
 433   NCollection_Map<TopoDS_Face, TopTools_ShapeMapHasher> dfm;
 434   //TopExp::MapShapes(aFuseShape, TopAbs_FACE, dfm);
 435   TopExp_Explorer expf(in, TopAbs_FACE);
 436   for (;expf.More(); expf.Next())
 437     dfm.Add(TopoDS::Face(expf.Current()));
 438
 439   int nbF = dfm.Extent();
 440   TopExp_Explorer exp_f(in, TopAbs_FACE);
 441   const TopoDS_Face& FF = TopoDS::Face(exp_f.Current());
 442   NCollection_List<TopoDS_Face> CurrFS;
 443   NCollection_List<TopoDS_Face> NeighFS;
 444   NCollection_Map<TopoDS_Face, TopTools_ShapeMapHasher> PrF;
 445   CurrFS.Append(FF);
 446   NCollection_List<NCollection_Map<TopoDS_Face, TopTools_ShapeMapHasher>> GL_F;
 447   NCollection_Map<TopoDS_Face, TopTools_ShapeMapHasher> OneGr;
 448   bool end = false;
 449   while (!end)
 450   {
 451     NCollection_List<TopoDS_Face>::Iterator it_currfs(CurrFS);
 452     NeighFS.Clear();
 453     for (;it_currfs.More();it_currfs.Next())
 454     {
 455       const TopoDS_Face& CF = it_currfs.Value();
 456       TopExp_Explorer exp_edge(CF, TopAbs_EDGE);
 457       for (;exp_edge.More();exp_edge.Next())
 458       {
 459         const TopoDS_Shape& CE = exp_edge.Current();
 460         const TopTools_ListOfShape& lsf = mE2LF.FindFromKey(CE);
 461         TopTools_ListIteratorOfListOfShape ls_it(lsf); //always one face (since all faces are planar)
 462         for (;ls_it.More();ls_it.Next())
 463         {
 464           const TopoDS_Face& F = TopoDS::Face(ls_it.Value());
 465           if (F.IsSame(CF))
 466             continue;
 467           if (!PrF.Contains(F))
 468           {
 469             OneGr.Add(F);
 470             NeighFS.Append(F);
 471             PrF.Add(F);
 472           }
 473         }
 474       }
 475       OneGr.Add(CF);
 476       PrF.Add(CF);
 477     }
 478     if (NeighFS.IsEmpty())
 479     {
 480       GL_F.Append(OneGr);
 481       OneGr.Clear();
 482       dfm.Subtract(PrF);
 483       if (dfm.IsEmpty())
 484         end = true;
 485       else
 486       {
 487         NCollection_Map<TopoDS_Face, TopTools_ShapeMapHasher>::Iterator itDm(dfm);
 488         const TopoDS_Face& nsh = itDm.Key();
 489         NeighFS.Append(nsh);
 490       }
 491     }
 492     CurrFS = NeighFS;
 493   }
 494
 495   TopoDS_Shape sh;
 496
 497   if (GL_F.Extent() > 1)
 498   {
 499     TopoDS_Compound cmp;
 500     NCollection_List<NCollection_Map<TopoDS_Face, TopTools_ShapeMapHasher>>::Iterator itGL_F(GL_F);
 501     BRep_Builder bb;
 502     bb.MakeCompound(cmp);
 503     for (;itGL_F.More();itGL_F.Next())
 504     {
 505       MakeShellG(itGL_F.Value(), sh);
 506       if (!sh.IsNull())
 507         bb.Add(cmp, sh);
 508     }
 509     return  cmp;
 510   }
 511   else if (GL_F.Extent() == 1)
 512   {
 513     MakeShellG(GL_F.First(), sh);
 514     return sh;
 515   }
 516
 517 }
 518
 519 TopoDS_Shape HYDROData_Tool::PolyXY2Face( const Handle(HYDROData_PolylineXY)& aPolyline )
 520 {
 521   //DEBTRACE("generateTopShape");
 522   TopoDS_Face aResultFace = TopoDS_Face(); // --- result: default = null face
 523
 524   if (!aPolyline.IsNull())
 525     {
 526       TopoDS_Shape aPolylineShape = aPolyline->GetShape();
 527 #ifdef DEB_IMZ
 528       std::string brepName = "imz.brep";
 529       BRepTools::Write(aPolylineShape, brepName.c_str());
 530 #endif
 531       TopTools_ListOfShape aWiresList;
 532
 533       if (!aPolylineShape.IsNull() && aPolylineShape.ShapeType() == TopAbs_WIRE)
 534         {
 535           // --- only one wire: try to make a face
 536           //DEBTRACE("one wire: try to build a face");
 537           const TopoDS_Wire& aPolylineWire = TopoDS::Wire(aPolylineShape);
 538           if (!aPolylineWire.IsNull())
 539             {
 540               BRepBuilderAPI_MakeFace aMakeFace(aPolylineWire, Standard_True);
 541               aMakeFace.Build();
 542               if (aMakeFace.IsDone())
 543                 {
 544                   //DEBTRACE(" a face with the only wire given");
 545                   aResultFace = aMakeFace.Face();
 546                 }
 547             }
 548         }
 549       else
 550         {
 551           // --- a list of wires ? inventory of wires and edges
 552           Handle(TopTools_HSequenceOfShape) aSeqWires = new TopTools_HSequenceOfShape;
 553           Handle(TopTools_HSequenceOfShape) aSeqEdges = new TopTools_HSequenceOfShape;
 554           TopExp_Explorer anExp(aPolylineShape, TopAbs_WIRE);
 555           //DEBTRACE("list of wires ?");
 556           for (; anExp.More(); anExp.Next())
 557             {
 558               if (!anExp.Current().IsNull())
 559                 {
 560                   const TopoDS_Wire& aWire = TopoDS::Wire(anExp.Current());
 561                   aWiresList.Append(aWire);
 562                   //DEBTRACE("  append wire");
 563                   TopExp_Explorer it2(aWire, TopAbs_EDGE);
 564                   for (; it2.More(); it2.Next())
 565                     aSeqEdges->Append(it2.Current());
 566                 }
 567             }
 568           if (aWiresList.IsEmpty())
 569             return aResultFace; // --- no wires: null result
 570
 571           if (aSeqEdges->Length() > 1)
 572             {
 573               //DEBTRACE("try to connect all the edges together, build a unique wire and a face");
 574               // --- try to create one wire by connecting edges with a distance tolerance (no necessity of sharing points)
 575               ShapeAnalysis_FreeBounds::ConnectEdgesToWires(aSeqEdges, 1E-5, Standard_False, aSeqWires);
 576
 577               if (aSeqWires->Length() == 1)
 578                 {
 579                   // --- one wire: try to make a face
 580                   const TopoDS_Wire& aPolylineWire = TopoDS::Wire(aSeqWires->Value(1));
 581                   if (!aPolylineWire.IsNull())
 582                     {
 583                       BRepBuilderAPI_MakeFace aMakeFace(aPolylineWire, Standard_True);
 584                       aMakeFace.Build();
 585                       if (aMakeFace.IsDone())
 586                         {
 587                           //DEBTRACE("  a face from all the wires connected");
 588                           aResultFace = aMakeFace.Face();
 589                         }
 590                     }
 591                 }
 592             }
 593
 594           if (aResultFace.IsNull())
 595             {
 596               //DEBTRACE("try to make a face with the first wire of the list and other wires as restrictions");
 597               // --- try to make a face with the first wire of the list and other wires as restrictions
 598               BRepAlgo_FaceRestrictor aFR;
 599               TopoDS_Face aRefFace;
 600               TopoDS_Shape aS = aWiresList.First();
 601               BRepBuilderAPI_MakeFace aMakeFace(TopoDS::Wire(aWiresList.First()), Standard_True);
 602               aMakeFace.Build();
 603               if (aMakeFace.IsDone())
 604                 {
 605                   //DEBTRACE("  a face with first wire");
 606                   aRefFace = aMakeFace.Face();
 607                 }
 608               if (!aRefFace.IsNull())
 609                 {
 610                   aFR.Init(aRefFace, Standard_False, Standard_True);
 611                   TopTools_ListIteratorOfListOfShape anIt(aWiresList);
 612                   for (; anIt.More(); anIt.Next())
 613                     {
 614                       TopoDS_Wire& aWire = TopoDS::Wire(anIt.Value());
 615                       if (aWire.IsNull())
 616                         continue;
 617                       aFR.Add(aWire);
 618                     }
 619                   aFR.Perform();
 620                   if (aFR.IsDone())
 621                     {
 622                       for (; aFR.More(); aFR.Next())
 623                         {
 624                           //DEBTRACE("  a restricted face");
 625                           aResultFace = aFR.Current();
 626                           break;
 627                         }
 628                     }
 629                 }
 630             }
 631         }
 632     }
 633
 634   if (aResultFace.IsNull())
 635     return aResultFace;
 636
 637   //DEBTRACE("check the face");
 638   BRepCheck_Analyzer anAnalyzer(aResultFace);
 639   if (anAnalyzer.IsValid() && aResultFace.ShapeType() == TopAbs_FACE)
 640   {
 641     //DEBTRACE("face OK");
 642     return aResultFace;
 643   }
 644   else
 645   {
 646     //DEBTRACE("bad face");
 647   }
 648   return TopoDS_Face();
 649 }
 650
 651 void HYDROData_Tool::SetSIProgress(const Handle(Message_ProgressIndicator)& thePI)
 652 {
 653   StricklerInterpolationProgress() = thePI;
 654 }
 655
 656 const Handle(Message_ProgressIndicator)& HYDROData_Tool::GetSIProgress()
 657 {
 658   return StricklerInterpolationProgress();
 659 }
 660
 661 Handle(Message_ProgressIndicator)& HYDROData_Tool::StricklerInterpolationProgress()
 662 {
 663   static Handle(Message_ProgressIndicator) aPI = NULL;
 664   return aPI;
 665 }
 666
 667 void HYDROData_Tool::SetZIProgress(const Handle(Message_ProgressIndicator)& thePI)
 668 {
 669   BathymetryInterpolationProgress() = thePI;
 670 }
 671
 672 const Handle(Message_ProgressIndicator)& HYDROData_Tool::GetZIProgress()
 673 {
 674   return BathymetryInterpolationProgress();
 675 }
 676
 677 Handle(Message_ProgressIndicator)& HYDROData_Tool::BathymetryInterpolationProgress()
 678 {
 679   static Handle(Message_ProgressIndicator) aPI = NULL;
 680   return aPI;
 681 }
 682
 683 void HYDROData_Tool::SetTriangulationStatus(const ExecStatus& theStatus)
 684 {
 685   myTriangulationStatus = theStatus;
 686 }
 687
 688 const HYDROData_Tool::ExecStatus& HYDROData_Tool::GetTriangulationStatus()
 689 {
 690   return myTriangulationStatus;
 691 }
 692
 693 static bool AddXYZ(bool bImportXY,
 694                    double x,
 695                    double y,
 696                    double z,
 697                    std::vector<gp_XYZ>& thePointsXYZ,
 698                    std::vector<gp_XY>& thePointsXY)
 699 {
 700   if (!bImportXY)
 701   {
 702     if ( HYDROData_Tool::IsNan( x ) || HYDROData_Tool::IsInf( x ) ||
 703       HYDROData_Tool::IsNan( y ) || HYDROData_Tool::IsInf( y ) ||
 704       HYDROData_Tool::IsNan( z ) || HYDROData_Tool::IsInf( z ) )
 705       return false;
 706
 707     if( thePointsXYZ.size()>=thePointsXYZ.capacity() )
 708       thePointsXYZ.reserve( thePointsXYZ.size()+BLOCK_SIZE );
 709
 710     thePointsXYZ.push_back(gp_XYZ(x,y,z));
 711   }
 712   else
 713   {
 714     if ( HYDROData_Tool::IsNan( x ) || HYDROData_Tool::IsInf( x ) ||
 715       HYDROData_Tool::IsNan( y ) || HYDROData_Tool::IsInf( y ) )
 716       return false;
 717
 718     if( thePointsXY.size()>=thePointsXY.capacity() )
 719       thePointsXY.reserve( thePointsXY.size()+BLOCK_SIZE );
 720
 721     thePointsXY.push_back(gp_XY(x,y));
 722   }
 723   return true;
 724 }
 725
 726 bool HYDROData_Tool::importFromXYZ( QString& theFileName,
 727                                     bool bImportXY,
 728                                     std::vector<gp_XYZ>& thePointsXYZ,
 729                                     std::vector<gp_XY>& thePointsXY)
 730 {
 731   QFile aFile( theFileName );
 732   if ( !aFile.exists() || !aFile.open( QIODevice::ReadOnly ) )
 733     return false;
 734
 735   QString aFileSuf = QFileInfo( aFile ).suffix().toLower();
 736
 737   double x,y;
 738   if ( aFileSuf == "xyz" )
 739   {
 740     double z;
 741     while ( !aFile.atEnd() )
 742     {
 743       std::string aLine = aFile.readLine().simplified().toStdString();
 744       if ( aLine.empty() )
 745         continue;
 746
 747       x = 0;
 748       y = 0;
 749       z = 0;
 750
 751       if( sscanf( aLine.c_str(), "%lf %lf %lf", &x, &y, &z )!=3 )
 752       {
 753         aFile.close();
 754         return false;
 755       }
 756
 757       if (!AddXYZ(bImportXY, x, y, z, thePointsXYZ, thePointsXY ))
 758       {
 759         aFile.close();
 760         return false;
 761       }
 762     }
 763   }
 764   else if (aFileSuf == "xy" )
 765   {
 766     while ( !aFile.atEnd() )
 767     {
 768       std::string aLine = aFile.readLine().simplified().toStdString();
 769       if ( aLine.empty() )
 770         continue;
 771
 772       x = 0;
 773       y = 0;
 774
 775       if( sscanf( aLine.c_str(), "%lf %lf", &x, &y )!=2 )
 776       {
 777         aFile.close();
 778         return false;
 779       }
 780
 781       if (!AddXYZ(true, x, y, 0, thePointsXYZ, thePointsXY ))
 782       {
 783         aFile.close();
 784         return false;
 785       }
 786     }
 787   }
 788
 789   aFile.close();
 790
 791   return true;
 792 }
 793
 794 bool HYDROData_Tool::importPolylineFromXYZ(QString aFileName, Handle(HYDROData_Document) theDocument,
 795   bool importXY, NCollection_Sequence<Handle(HYDROData_Entity)>& importedEntities)
 796 {
 797   if (importXY)
 798   {
 799     std::vector<gp_XY> aPoints2d;
 800     std::vector<gp_XYZ> aDPoints3d;
 801
 802     if (HYDROData_Tool::importFromXYZ(aFileName, importXY, aDPoints3d, aPoints2d))
 803     {
 804       QString basename = QFileInfo( aFileName ).baseName();
 805
 806       Handle(HYDROData_PolylineXY) aPolylineXY = Handle(HYDROData_PolylineXY)::DownCast( theDocument->CreateObject( KIND_POLYLINEXY ) );
 807       HYDROData_PolylineXY::SectionType aSectType = HYDROData_PolylineXY::SECTION_POLYLINE;
 808       bool IsClosed = false;
 809       if ((aPoints2d.front()-aPoints2d.back()).Modulus()<Precision::Confusion())
 810       {
 811         IsClosed = true;
 812         aPolylineXY->AddSection( TCollection_AsciiString("poly_section"), aSectType, true);
 813       }
 814       else
 815         aPolylineXY->AddSection( TCollection_AsciiString("poly_section"), aSectType, false);
 816
 817       int n = aPoints2d.size();
 818       if (IsClosed)
 819         n--;
 820
 821       for ( int i = 0; i < n; i++ )
 822       {
 823         gp_XY aSectPoint = aPoints2d[i];
 824         theDocument->Transform(aSectPoint, true);
 825         aPolylineXY->AddPoint( 0, aSectPoint );
 826       }
 827
 828       aPolylineXY->SetWireColor( HYDROData_PolylineXY::DefaultWireColor() );
 829       aPolylineXY->SetName( basename + "_PolyXY_" );
 830       aPolylineXY->Update();
 831       importedEntities.Append(aPolylineXY);
 832       return true;
 833     }
 834     else
 835       return false;
 836   }
 837   else //xyz
 838   {
 839     std::vector<gp_XY> aDPoints2d;
 840     std::vector<gp_XYZ> aPoints3d;
 841     if (HYDROData_Tool::importFromXYZ(aFileName, false, aPoints3d, aDPoints2d))
 842     {
 843       QString basename = QFileInfo( aFileName ).baseName();
 844       Handle(HYDROData_PolylineXY) aPolylineXY = Handle(HYDROData_PolylineXY)::DownCast( theDocument->CreateObject( KIND_POLYLINEXY ) );
 845       Handle(HYDROData_Polyline3D) aPolylineObj = Handle(HYDROData_Polyline3D)::DownCast( theDocument->CreateObject( KIND_POLYLINE ) );
 846       Handle(HYDROData_Bathymetry) aBath = Handle(HYDROData_Bathymetry)::DownCast( theDocument->CreateObject( KIND_BATHYMETRY ) );
 847       HYDROData_Bathymetry::AltitudePoints aAPoints;
 848       HYDROData_PolylineXY::SectionType aSectType = HYDROData_PolylineXY::SECTION_POLYLINE;
 849       bool IsClosed = false;
 850       if ((aPoints3d.front()-aPoints3d.back()).Modulus()<Precision::Confusion())
 851       {
 852         IsClosed = true;
 853         aPolylineXY->AddSection( TCollection_AsciiString("poly_section"), aSectType, true);
 854       }
 855       else
 856         aPolylineXY->AddSection( TCollection_AsciiString("poly_section"), aSectType, false);
 857
 858       int n = aPoints3d.size();
 859       if (IsClosed)
 860         n--;
 861
 862       for ( int i = 0; i < n; i++ )
 863       {
 864         gp_XY aSectPoint(aPoints3d[i].X(), aPoints3d[i].Y());
 865         theDocument->Transform(aSectPoint, true);
 866         aPolylineXY->AddPoint( 0, aSectPoint );
 867         HYDROData_Bathymetry::AltitudePoint p;
 868         p.X = aSectPoint.X();
 869         p.Y = aSectPoint.Y();
 870         p.Z = aPoints3d[i].Z();
 871         aAPoints.push_back(p);
 872       }
 873
 874       QString aBathName = basename + "_bath_";
 875       QString aPolyXYName = basename + "_polyXY_";
 876       QString aPoly3DName = basename + "_poly3D_";
 877
 878       aPolylineXY->SetName( aPolyXYName );
 879       aPolylineXY->SetWireColor(HYDROData_PolylineXY::DefaultWireColor());
 880       aPolylineXY->Update();
 881
 882       aBath->SetAltitudePoints(aAPoints);
 883       aBath->SetName( aBathName );
 884
 885       aPolylineObj->SetPolylineXY (aPolylineXY, false);
 886       aPolylineObj->SetAltitudeObject(aBath);
 887
 888       aPolylineObj->SetBorderColor( aPolylineObj->DefaultBorderColor() );
 889       aPolylineObj->SetName( aPoly3DName );
 890
 891       aPolylineObj->Update();
 892       importedEntities.Append(aPolylineXY);
 893       importedEntities.Append(aPolylineObj);
 894       return true;
 895     }
 896     else
 897       return false;
 898   }
 899 }
 900
 901
 902 std::ostream& operator<<( std::ostream& theStream, const QString& theText )
 903 {
 904   theStream << theText.toStdString();
 905   return theStream;
 906 }
 907
 908 std::ostream& operator<<( std::ostream& theStream, const QColor& theColor )
 909 {
 910   theStream << "[" << theColor.red() << ", " << theColor.green() << ", " << theColor.blue() << "]";
 911   return theStream;
 912 }
 913
 914 std::ostream& operator<<( std::ostream& theStream, const TopoDS_Shape& theShape )
 915 {
 916   theStream << "[" << theShape.TShape().operator->() << "]";
 917   return theStream;
 918 }
 919
 920 std::ostream& operator<<( std::ostream& theStream, const TopoDS_Face& theFace )
 921 {
 922   theStream << "[" << theFace.TShape().operator->() << "]";
 923   return theStream;
 924 }
 925
 926 std::ostream& operator<<( std::ostream& theStream, const gp_XY& theXY )
 927 {
 928   theStream << "(" << theXY.X() << "; " << theXY.Y() << ")";
 929   return theStream;
 930 }
 931
 932 bool operator == ( const gp_XY& thePoint1, const gp_XY& thePoint2 )
 933 {
 934   const double EPS = 1E-3;
 935   return
 936     fabs( thePoint1.X() - thePoint2.X() ) < EPS &&
 937     fabs( thePoint1.Y() - thePoint2.Y() ) < EPS;
 938
 939 }