Salome HOME
Minor change.
[modules/hydro.git] / src / HYDROData / HYDROData_Bathymetry.cxx
index 31bd8cc7435de7c709beda555bbf18e51073eec5..60e2d9bf300fd069819fa1c8188df4803b4f14ac 100644 (file)
@@ -1,20 +1,36 @@
 
 #include "HYDROData_Bathymetry.h"
+#include "HYDROData_Document.h"
+#include "HYDROData_Tool.h"
 
 #include <gp_XY.hxx>
 #include <gp_XYZ.hxx>
 
 #include <TDataStd_RealArray.hxx>
+#include <TDataStd_AsciiString.hxx>
+#include <TDataStd_Integer.hxx>
 
 #include <QFile>
 #include <QFileInfo>
 #include <QPointF>
+#include <QPolygonF>
 #include <QStringList>
 
-IMPLEMENT_STANDARD_HANDLE(HYDROData_Bathymetry, HYDROData_Object)
-IMPLEMENT_STANDARD_RTTIEXT(HYDROData_Bathymetry, HYDROData_Object)
+#include <math.h>
+
+#define _TIMER
+#ifdef _TIMER
+#include <OSD_Timer.hxx>
+#endif
+
+#define PYTHON_BATHYMETRY_ID "KIND_BATHYMETRY"
+
+
+IMPLEMENT_STANDARD_HANDLE(HYDROData_Bathymetry, HYDROData_IAltitudeObject)
+IMPLEMENT_STANDARD_RTTIEXT(HYDROData_Bathymetry, HYDROData_IAltitudeObject)
 
 HYDROData_Bathymetry::HYDROData_Bathymetry()
+: HYDROData_IAltitudeObject()
 {
 }
 
@@ -22,6 +38,39 @@ HYDROData_Bathymetry::~HYDROData_Bathymetry()
 {
 }
 
+QStringList HYDROData_Bathymetry::DumpToPython( MapOfTreatedObjects& theTreatedObjects ) const
+{
+  QStringList aResList;
+
+  Handle(HYDROData_Document) aDocument = HYDROData_Document::Document( myLab );
+  if ( aDocument.IsNull() )
+    return aResList;
+                             
+  QString aDocName = aDocument->GetDocPyName();
+  QString aBathymetryName = GetName();
+
+  aResList << QString( "%1 = %2.CreateObject( %3 );" )
+              .arg( aBathymetryName ).arg( aDocName ).arg( PYTHON_BATHYMETRY_ID );
+  aResList << QString( "%1.SetName( \"%2\" );" )
+              .arg( aBathymetryName ).arg( aBathymetryName );
+
+  aResList << QString( "%1.SetAltitudesInverted( %2 );" )
+              .arg( aBathymetryName ).arg( IsAltitudesInverted() );
+
+  QString aFilePath = GetFilePath();
+  if ( !aFilePath.isEmpty() )
+  {
+    aResList << QString( "%1.ImportFromFile( \"%2\" );" )
+                .arg( aBathymetryName ).arg( aFilePath );
+  }
+  else
+  {
+    // TODO : bathymetry is composed from other bathymetry(ies)
+  }
+
+  return aResList;
+}
+
 void HYDROData_Bathymetry::SetAltitudePoints( const AltitudePoints& thePoints )
 {
   RemoveAltitudePoints();
@@ -43,21 +92,25 @@ void HYDROData_Bathymetry::SetAltitudePoints( const AltitudePoints& thePoints )
     aCoordsArray->SetValue( i * 3 + 1, aPoint.Y() );
     aCoordsArray->SetValue( i * 3 + 2, aPoint.Z() );
   }
+
+  SetToUpdate( true );
 }
 
 HYDROData_Bathymetry::AltitudePoints HYDROData_Bathymetry::GetAltitudePoints() const
 {
   AltitudePoints aPoints;
 
+  TDF_Label aLabel = myLab.FindChild( DataTag_AltitudePoints, false );
+  if ( aLabel.IsNull() )
+    return aPoints;
+
   Handle(TDataStd_RealArray) aCoordsArray;
-  if ( !myLab.FindChild( DataTag_AltitudePoints ).FindAttribute( TDataStd_RealArray::GetID(), aCoordsArray ) )
+  if ( !aLabel.FindAttribute( TDataStd_RealArray::GetID(), aCoordsArray ) )
     return aPoints;
 
-  int aLowerIdx = aCoordsArray->Lower();
-  int anUpperIdx = aCoordsArray->Upper();
   for ( int i = aCoordsArray->Lower(), n = aCoordsArray->Upper(); i <= n; )
   {
-    if ( i + 3 > n )
+    if ( i + 3 > n + 1 )
       break;
 
     AltitudePoint aPoint;
@@ -73,27 +126,270 @@ HYDROData_Bathymetry::AltitudePoints HYDROData_Bathymetry::GetAltitudePoints() c
 
 void HYDROData_Bathymetry::RemoveAltitudePoints()
 {
-  TDF_Label aLab = myLab.FindChild( DataTag_AltitudePoints );
-  aLab.ForgetAllAttributes();
+  TDF_Label aLabel = myLab.FindChild( DataTag_AltitudePoints, false );
+  if ( !aLabel.IsNull() )
+  {
+    aLabel.ForgetAllAttributes();
+    SetToUpdate( true );
+  }
 }
 
-double HYDROData_Bathymetry::GetAltitudeForPoint( const QPointF& thePoint ) const
+void interpolateAltitudeForPoints( const gp_XY&                               thePoint,
+                                   const HYDROData_Bathymetry::AltitudePoint& theFirstPoint,
+                                   const HYDROData_Bathymetry::AltitudePoint& theSecPoint,
+                                   HYDROData_Bathymetry::AltitudePoint&       theResPoint,
+                                   const bool&                                theIsVertical )
 {
-  gp_XY aGpPoint( thePoint.x(), thePoint.y() );
-  return GetAltitudeForPoint( aGpPoint );
+  double aCoordX = thePoint.X();
+  double aCoordY = thePoint.Y();
+
+  if ( theIsVertical )
+  {
+    aCoordX = theFirstPoint.X();
+
+    if ( !ValuesEquals( theFirstPoint.X(), theSecPoint.X() ) )
+    {
+      // Recalculate X coordinate by equation of line from two points
+      aCoordX = ( ( ( thePoint.Y() - theFirstPoint.Y() ) * ( theSecPoint.X() - theFirstPoint.X() ) ) /
+                  ( theSecPoint.Y() - theFirstPoint.Y() ) ) + theFirstPoint.X();
+    }
+  }
+  else
+  {
+    aCoordY = theFirstPoint.Y();
+
+    if ( !ValuesEquals( theFirstPoint.Y(), theSecPoint.Y() ) )
+    {
+      // Recalculate y by equation of line from two points
+      aCoordY = ( ( ( thePoint.X() - theFirstPoint.X() ) * ( theSecPoint.Y() - theFirstPoint.Y() ) ) /
+                  ( theSecPoint.X() - theFirstPoint.X() ) ) + theFirstPoint.Y();
+    }
+  }
+
+  theResPoint.SetX( aCoordX );
+  theResPoint.SetY( aCoordY );
+
+  // Calculate coefficient for interpolation
+  double aLength = Sqrt( Pow( theSecPoint.Y() - theFirstPoint.Y(), 2 ) +
+                         Pow( theSecPoint.X() - theFirstPoint.X(), 2 ) );
+
+  double aInterCoeff = 0;
+  if ( aLength != 0 )
+   aInterCoeff = ( theSecPoint.Z() - theFirstPoint.Z() ) / aLength;
+
+
+  double aNewLength = Sqrt( Pow( theResPoint.Y() - theFirstPoint.Y(), 2 ) +
+                            Pow( theResPoint.X() - theFirstPoint.X(), 2 ) );
+
+  // Calculate interpolated value
+  double aResVal = theFirstPoint.Z() + aInterCoeff * aNewLength;
+
+  theResPoint.SetZ( aResVal );
 }
 
 double HYDROData_Bathymetry::GetAltitudeForPoint( const gp_XY& thePoint ) const
 {
-  double aResAltitude = -9999.90;
+  double anInvalidAltitude = GetInvalidAltitude();
+  double aResAltitude = anInvalidAltitude;
   
   AltitudePoints anAltitudePoints = GetAltitudePoints();
+  if ( anAltitudePoints.isEmpty() )
+    return aResAltitude;
+
+  QPolygonF aBoundingRect;
+
+  // Boundary plane
+  // [ 0 (top-left) ]          [ 1 (top-right) ]
+  //                  thePoint
+  // [ 2 (bot-left) ]          [ 3 (bot-right) ] 
+  AltitudePoint aBounds[ 4 ] = { AltitudePoint( -DBL_MAX, -DBL_MAX, anInvalidAltitude ),
+                                 AltitudePoint(  DBL_MAX, -DBL_MAX, anInvalidAltitude ),
+                                 AltitudePoint( -DBL_MAX,  DBL_MAX, anInvalidAltitude ),
+                                 AltitudePoint(  DBL_MAX,  DBL_MAX, anInvalidAltitude ) }; 
+
+  AltitudePoints::const_iterator aListItBeg = anAltitudePoints.constBegin();
+  AltitudePoints::const_iterator aListItEnd = anAltitudePoints.constEnd();
+  for ( ; aListItBeg != aListItEnd; ++aListItBeg )
+  {
+    const AltitudePoint& aPoint = *aListItBeg;
+
+    double aDeltaX = Abs( aPoint.X() ) - Abs( thePoint.X() );
+    double aDeltaY = Abs( aPoint.Y() ) - Abs( thePoint.Y() );
+
+    if ( ValuesEquals( aDeltaX, 0.0 ) ) // Both left and right sides
+    {
+      if ( ValuesEquals( aDeltaY, 0.0 ) ) // Both top and bottom sides
+      {
+        aResAltitude = aPoint.Z();
+        return aResAltitude;
+      }
+      else if ( aDeltaY < 0 ) // top side
+      {
+        // top border
+        if ( ValuesMoreEquals( aPoint.X(), aBounds[ 0 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 0 ].Y() ) )
+          aBounds[ 0 ] = aPoint;
+        if ( ValuesLessEquals( aPoint.X(), aBounds[ 1 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 1 ].Y() ) )
+          aBounds[ 1 ] = aPoint;
+      }
+      else
+      {
+        // bottom border
+        if ( ValuesMoreEquals( aPoint.X(), aBounds[ 2 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 2 ].Y() ) )
+          aBounds[ 2 ] = aPoint;
+        if ( ValuesLessEquals( aPoint.X(), aBounds[ 3 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 3 ].Y() ) )
+          aBounds[ 3 ] = aPoint;
+      }
+    }
+    else if ( aDeltaX < 0 ) // left side
+    {
+      if ( ValuesEquals( aDeltaY, 0.0 ) )
+      {
+        // Left border
+        if ( ValuesMoreEquals( aPoint.X(), aBounds[ 0 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 0 ].Y() ) )
+          aBounds[ 0 ] = aPoint;
+        if ( ValuesMoreEquals( aPoint.X(), aBounds[ 2 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 2 ].Y() ) )
+          aBounds[ 2 ] = aPoint;
+      }
+      else if ( aDeltaY < 0 )
+      {
+        // top left corner
+        if ( ValuesMoreEquals( aPoint.X(), aBounds[ 0 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 0 ].Y() ) )
+          aBounds[ 0 ] = aPoint;
+      }
+      else
+      {
+        // bottom left corner
+        if ( ValuesMoreEquals( aPoint.X(), aBounds[ 2 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 2 ].Y() ) )
+          aBounds[ 2 ] = aPoint;
+      }
+    }
+    else // right side
+    {
+      if ( ValuesEquals( aDeltaY, 0.0 ) )
+      {
+        // Right border
+        if ( ValuesLessEquals( aPoint.X(), aBounds[ 1 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 1 ].Y() ) )
+          aBounds[ 1 ] = aPoint;
+        if ( ValuesLessEquals( aPoint.X(), aBounds[ 3 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 3 ].Y() ) )
+          aBounds[ 3 ] = aPoint;
+      }
+      else if ( aDeltaY < 0 )
+      {
+        // top right corner
+        if ( ValuesLessEquals( aPoint.X(), aBounds[ 1 ].X() ) && ValuesMoreEquals( aPoint.Y(), aBounds[ 1 ].Y() ) )
+          aBounds[ 1 ] = aPoint;
+      }
+      else
+      {
+        // bottom right corner
+        if ( ValuesLessEquals( aPoint.X(), aBounds[ 3 ].X() ) && ValuesLessEquals( aPoint.Y(), aBounds[ 3 ].Y() ) )
+          aBounds[ 3 ] = aPoint;
+      }
+    }
+
+    // Update bounding rectangle of our global grid
+    aBoundingRect << QPointF( aPoint.X(), aPoint.Y() );
+  }
+
+  const double LIMIT = 1E300;
+  if( fabs( aBounds[ 0 ].X() ) > LIMIT || fabs( aBounds[ 0 ].Y() ) > LIMIT ||
+      fabs( aBounds[ 1 ].X() ) > LIMIT || fabs( aBounds[ 1 ].Y() ) > LIMIT ||
+      fabs( aBounds[ 2 ].X() ) > LIMIT || fabs( aBounds[ 2 ].Y() ) > LIMIT ||
+      fabs( aBounds[ 3 ].X() ) > LIMIT || fabs( aBounds[ 3 ].Y() ) > LIMIT )
+    return anInvalidAltitude;
+
+
+  // Check if requested point is inside of our bounding rectangle
+  if ( !aBoundingRect.boundingRect().contains( thePoint.X(), thePoint.Y() ) )
+    return aResAltitude;
+
+  // Calculate result altitude for point
+  AltitudePoint aFirstPoint( aBounds[ 0 ] ), aSecPoint( aBounds[ 1 ] );
 
-  // TODO : implement
+  // At first we merge top and bottom borders
+  if ( aBounds[ 0 ].Y() != aBounds[ 2 ].Y() || aBounds[ 0 ].X() != aBounds[ 2 ].X() )
+    interpolateAltitudeForPoints( thePoint, aBounds[ 0 ], aBounds[ 2 ], aFirstPoint, true );
+
+  if ( aBounds[ 1 ].Y() != aBounds[ 3 ].Y() || aBounds[ 1 ].X() != aBounds[ 3 ].X() )
+    interpolateAltitudeForPoints( thePoint, aBounds[ 1 ], aBounds[ 3 ], aSecPoint, true );
+
+  AltitudePoint aResPoint( aFirstPoint );
+
+  // At last we merge left and right borders
+  if ( aFirstPoint.Y() != aSecPoint.Y() || aFirstPoint.X() != aSecPoint.X() )
+    interpolateAltitudeForPoints( thePoint, aFirstPoint, aSecPoint, aResPoint, false );
+    
+  aResAltitude = aResPoint.Z();
 
   return aResAltitude;
 }
 
+void HYDROData_Bathymetry::SetFilePath(const QString& theFilePath)
+{
+  TCollection_AsciiString anAsciiStr( theFilePath.toStdString().c_str() );
+  TDataStd_AsciiString::Set( myLab.FindChild( DataTag_FilePath ), anAsciiStr );
+}
+
+QString HYDROData_Bathymetry::GetFilePath() const
+{
+  QString aRes;
+
+  TDF_Label aLabel = myLab.FindChild( DataTag_FilePath, false );
+  if ( !aLabel.IsNull() )
+  {
+    Handle(TDataStd_AsciiString) anAsciiStr;
+    if ( aLabel.FindAttribute( TDataStd_AsciiString::GetID(), anAsciiStr ) )
+      aRes = QString( anAsciiStr->Get().ToCString() );
+  }
+
+  return aRes;
+}
+
+void HYDROData_Bathymetry::SetAltitudesInverted( const bool theIsInverted,
+                                                 const bool theIsUpdate )
+{
+  bool anIsAltitudesInverted = IsAltitudesInverted();
+  if ( anIsAltitudesInverted == theIsInverted )
+    return;
+
+  TDataStd_Integer::Set( myLab.FindChild( DataTag_AltitudesInverted ), (Standard_Integer)theIsInverted );
+
+  SetToUpdate( true );
+
+  if ( !theIsUpdate )
+    return;
+
+  // Update altitude points
+  AltitudePoints anAltitudePoints = GetAltitudePoints();
+  if ( anAltitudePoints.isEmpty() )
+    return;
+
+  AltitudePoints::iterator aListItBeg = anAltitudePoints.begin();
+  AltitudePoints::iterator aListItEnd = anAltitudePoints.end();
+  for ( ; aListItBeg != aListItEnd; ++aListItBeg )
+  {
+    AltitudePoint& aPoint = *aListItBeg;
+    aPoint.SetZ( aPoint.Z() * -1 );
+  }
+
+  SetAltitudePoints( anAltitudePoints );
+}
+
+bool HYDROData_Bathymetry::IsAltitudesInverted() const
+{
+  bool aRes = false;
+
+  TDF_Label aLabel = myLab.FindChild( DataTag_AltitudesInverted, false );
+  if ( !aLabel.IsNull() )
+  {
+    Handle(TDataStd_Integer) anIntVal;
+    if ( aLabel.FindAttribute( TDataStd_Integer::GetID(), anIntVal ) )
+      aRes = (bool)anIntVal->Get();
+  }
+
+  return aRes;
+}
+
 bool HYDROData_Bathymetry::ImportFromFile( const QString& theFileName )
 {
   // Try to open the file
@@ -116,7 +412,8 @@ bool HYDROData_Bathymetry::ImportFromFile( const QString& theFileName )
 
   if ( aRes )
   {
-    // Update altitude points of this Bathymetry
+    // Update file path and altitude points of this Bathymetry
+    SetFilePath( theFileName );
     SetAltitudePoints( aPoints );
   }
 
@@ -124,7 +421,7 @@ bool HYDROData_Bathymetry::ImportFromFile( const QString& theFileName )
 }
 
 bool HYDROData_Bathymetry::importFromXYZFile( QFile&          theFile,
-                                              AltitudePoints& thePoints )
+                                              AltitudePoints& thePoints ) const
 {
   if ( !theFile.isOpen() )
     return false;
@@ -134,13 +431,19 @@ bool HYDROData_Bathymetry::importFromXYZFile( QFile&          theFile,
   //  2. X(float) Y(float) Z(float)
   //  ...
 
+#ifdef _TIMER
+  OSD_Timer aTimer;
+  aTimer.Start();
+#endif
+
+  bool anIsAltitudesInverted = IsAltitudesInverted();
   while ( !theFile.atEnd() )
   {
-    QString aLine = theFile.readLine();
+    QString aLine = theFile.readLine().simplified();
     if ( aLine.isEmpty() )
       continue;
 
-    QStringList aValues = aLine.split( QRegExp( "\\s+" ), QString::SkipEmptyParts );
+    QStringList aValues = aLine.split( ' ', QString::SkipEmptyParts );
     if ( aValues.length() < 3 )
       return false;
 
@@ -154,14 +457,24 @@ bool HYDROData_Bathymetry::importFromXYZFile( QFile&          theFile,
 
     aPoint.SetX( anX.toDouble( &isXOk ) );
     aPoint.SetY( anY.toDouble( &isYOk ) );
-    aPoint.SetZ(  aZ.toDouble( &isZOk ) );
+    aPoint.SetZ( aZ.toDouble( &isZOk ) );
 
     if ( !isXOk || !isYOk || !isZOk )
       return false;
 
+    // Invert the z value if requested
+    if ( anIsAltitudesInverted )
+      aPoint.SetZ( -aPoint.Z() );
+
     thePoints << aPoint;
   }
 
+#ifdef _TIMER
+  aTimer.Stop();
+  std::ofstream stream( "W:/HYDRO/WORK/log.txt", std::ofstream::out );
+  aTimer.Show( stream );
+#endif
+
   return true;
 }