X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FHYDROData%2FHYDROData_Bathymetry.cxx;h=2e6ac1cc2638a741eeb9f8880fe14b1c1be83cb7;hb=f86cf7ecf17dbae2a1d84e0ebbab07c732208c2f;hp=9fabcb648d6175233c79c51efa9170d050e5782a;hpb=4d422c4f2dc05c85f0924db65ab87a8d5deacea8;p=modules%2Fhydro.git diff --git a/src/HYDROData/HYDROData_Bathymetry.cxx b/src/HYDROData/HYDROData_Bathymetry.cxx index 9fabcb64..2e6ac1cc 100644 --- a/src/HYDROData/HYDROData_Bathymetry.cxx +++ b/src/HYDROData/HYDROData_Bathymetry.cxx @@ -20,8 +20,7 @@ #include "HYDROData_Document.h" #include "HYDROData_Tool.h" #include "HYDROData_PolylineXY.h" - -#include +#include "HYDROData_QuadtreeNode.hxx" #include #include @@ -29,6 +28,7 @@ #include #include #include +#include #include #include @@ -36,6 +36,17 @@ #include #include #include +#include + +#ifndef LIGHT_MODE +#include +#include +#include +#include +#include +#endif + +#include #include @@ -47,30 +58,44 @@ #define _DEVDEBUG_ #include "HYDRO_trace.hxx" -IMPLEMENT_STANDARD_HANDLE(HYDROData_Bathymetry, HYDROData_IAltitudeObject) +const int BLOCK_SIZE = 1000; + IMPLEMENT_STANDARD_RTTIEXT(HYDROData_Bathymetry, HYDROData_IAltitudeObject) //HYDROData_QuadtreeNode* HYDROData_Bathymetry::myQuadtree = 0; + std::map HYDROData_Bathymetry::myQuadtrees; +#ifndef LIGHT_MODE +std::map 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() { - //DEBTRACE("HYDROData_Bathymetry constructor start " << this); -// if (! myQuadtree) -// myQuadtree = new HYDROData_QuadtreeNode(0, 30, 5, 0.); - //DEBTRACE("HYDROData_Bathymetry constructor end " << this); } HYDROData_Bathymetry::~HYDROData_Bathymetry() { - //DEBTRACE("HYDROData_Bathymetry destructor start " << this); -// if (myQuadtree) -// delete myQuadtree; -// Nodes_3D::iterator it = myListOfNodes.begin(); -// for( ; it != myListOfNodes.end(); ++it) -// delete *it; -// myListOfNodes.clear(); } QStringList HYDROData_Bathymetry::DumpToPython( const QString& thePyScriptPath, @@ -79,39 +104,39 @@ QStringList HYDROData_Bathymetry::DumpToPython( const QString& thePyScriptPath, QStringList aResList = dumpObjectCreation( theTreatedObjects ); QString aBathymetryName = GetObjPyName(); - aResList << QString( "%1.SetAltitudesInverted( %2 );" ) + aResList << QString( "%1.SetAltitudesInverted( %2 )" ) .arg( aBathymetryName ).arg( IsAltitudesInverted() ); TCollection_AsciiString aFilePath = GetFilePath(); - aResList << QString( "%1.ImportFromFile( \"%2\" );" ) + aResList << QString( "if not(%1.ImportFromFile( \"%2\" )):" ) .arg( aBathymetryName ).arg( aFilePath.ToCString() ); - + aResList << QString( " raise ValueError('problem while loading bathymetry')" ); aResList << QString( "" ); - aResList << QString( "%1.Update();" ).arg( aBathymetryName ); + aResList << QString( "%1.Update()" ).arg( aBathymetryName ); aResList << QString( "" ); return aResList; } -void HYDROData_Bathymetry::SetAltitudePoints( const AltitudePoints& thePoints ) +void HYDROData_Bathymetry::SetAltitudePoints( const HYDROData_Bathymetry::AltitudePoints& thePoints ) { RemoveAltitudePoints(); - if ( thePoints.IsEmpty() ) + if ( thePoints.empty() ) return; // Save coordinates Handle(TDataStd_RealArray) aCoordsArray = - TDataStd_RealArray::Set( myLab.FindChild( DataTag_AltitudePoints ), 0, thePoints.Length() * 3 - 1 ); + TDataStd_RealArray::Set( myLab.FindChild( DataTag_AltitudePoints ), 0, thePoints.size() * 3 - 1 ); - AltitudePoints::Iterator anIter( thePoints ); - for ( int i = 0 ; anIter.More(); ++i, anIter.Next() ) + HYDROData_Bathymetry::AltitudePoints::const_iterator anIter = thePoints.begin(), aLast = thePoints.end(); + for ( int i = 0 ; anIter!=aLast; ++i, ++anIter ) { - const AltitudePoint& aPoint = anIter.Value(); + const HYDROData_Bathymetry::AltitudePoint& aPoint = *anIter; - aCoordsArray->SetValue( i * 3, aPoint.X() ); - aCoordsArray->SetValue( i * 3 + 1, aPoint.Y() ); - aCoordsArray->SetValue( i * 3 + 2, aPoint.Z() ); + aCoordsArray->SetValue( i * 3, aPoint.X ); + aCoordsArray->SetValue( i * 3 + 1, aPoint.Y ); + aCoordsArray->SetValue( i * 3 + 2, aPoint.Z ); } Changed( Geom_Z ); @@ -119,7 +144,7 @@ void HYDROData_Bathymetry::SetAltitudePoints( const AltitudePoints& thePoints ) HYDROData_Bathymetry::AltitudePoints HYDROData_Bathymetry::GetAltitudePoints(bool IsConvertToGlobal) const { - AltitudePoints aPoints; + HYDROData_Bathymetry::AltitudePoints aPoints; TDF_Label aLabel = myLab.FindChild( DataTag_AltitudePoints, false ); if ( aLabel.IsNull() ) @@ -130,19 +155,21 @@ HYDROData_Bathymetry::AltitudePoints HYDROData_Bathymetry::GetAltitudePoints(boo return aPoints; Handle(HYDROData_Document) aDoc = HYDROData_Document::Document( myLab ); + int q = ( aCoordsArray->Upper() - aCoordsArray->Lower() + 1 ) / 3; + aPoints.reserve( q ); for ( int i = aCoordsArray->Lower(), n = aCoordsArray->Upper(); i <= n; ) { if ( i + 3 > n + 1 ) break; - AltitudePoint aPoint; - aPoint.SetX( aCoordsArray->Value( i++ ) ); - aPoint.SetY( aCoordsArray->Value( i++ ) ); - aPoint.SetZ( aCoordsArray->Value( i++ ) ); + HYDROData_Bathymetry::AltitudePoint aPoint; + aPoint.X = aCoordsArray->Value( i++ ); + aPoint.Y = aCoordsArray->Value( i++ ); + aPoint.Z = aCoordsArray->Value( i++ ); if( IsConvertToGlobal ) - aDoc->Transform( aPoint, false ); - aPoints.Append( aPoint ); + aDoc->Transform( aPoint.X, aPoint.Y, aPoint.Z, false ); + aPoints.push_back( aPoint ); } return aPoints; @@ -154,12 +181,12 @@ HYDROData_QuadtreeNode* HYDROData_Bathymetry::GetQuadtreeNodes() const if (aLabel.IsNull()) return 0; int labkey = myLab.Tag(); - int altkey = aLabel.Tag(); + //int altkey = aLabel.Tag(); //DEBTRACE("GetQuadtreeNodes this labkey altkey "<isEmpty() ) + // if (myQuadtree->isEmpty() ) if (myQuadtrees.find(labkey) == myQuadtrees.end()) { - DEBTRACE("GetQuadtreeNodes init " << this << " " << labkey); + //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); @@ -172,6 +199,7 @@ HYDROData_QuadtreeNode* HYDROData_Bathymetry::GetQuadtreeNodes() const Nodes_3D* aListOfNodes = new Nodes_3D(); + int index =0; for (int i = aCoordsArray->Lower(), n = aCoordsArray->Upper(); i <= n;) { if (i + 3 > n + 1) @@ -180,10 +208,11 @@ HYDROData_QuadtreeNode* HYDROData_Bathymetry::GetQuadtreeNodes() const double x = aCoordsArray->Value(i++); double y = aCoordsArray->Value(i++); double z = aCoordsArray->Value(i++); - gp_XYZ* aPoint = new gp_XYZ(x, y, z); + gpi_XYZ* aPoint = new gpi_XYZ(x, y, z, index); + index++; aListOfNodes->push_back(aPoint); } - DEBTRACE(" GetQuadtreeNodes call setNodesAndCompute"); + //DEBTRACE(" GetQuadtreeNodes call setNodesAndCompute"); aQuadtree->setNodesAndCompute(aListOfNodes); return aQuadtree; } @@ -191,6 +220,55 @@ HYDROData_QuadtreeNode* HYDROData_Bathymetry::GetQuadtreeNodes() const return myQuadtrees[labkey]; } +#ifndef LIGHT_MODE +vtkPolyData* HYDROData_Bathymetry::GetVtkDelaunay2D() 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 "<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[labkey] = data; + return data; + } + else + return myDelaunay2D[labkey]; + +} +#endif void HYDROData_Bathymetry::RemoveAltitudePoints() { TDF_Label aLabel = myLab.FindChild( DataTag_AltitudePoints, false ); @@ -201,224 +279,178 @@ void HYDROData_Bathymetry::RemoveAltitudePoints() } } -void interpolateAltitudeForPoints( const gp_XY& thePoint, +void interpolateAltitudeForPoints( const gp_XY& thePoint, const HYDROData_Bathymetry::AltitudePoint& theFirstPoint, const HYDROData_Bathymetry::AltitudePoint& theSecPoint, HYDROData_Bathymetry::AltitudePoint& theResPoint, - const bool& theIsVertical ) + const bool& theIsVertical ) { double aCoordX = thePoint.X(); double aCoordY = thePoint.Y(); if ( theIsVertical ) { - aCoordX = theFirstPoint.X(); + aCoordX = theFirstPoint.X; - if ( !ValuesEquals( theFirstPoint.X(), theSecPoint.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(); + aCoordX = ( ( ( thePoint.Y() - theFirstPoint.Y ) * ( theSecPoint.X - theFirstPoint.X ) ) / + ( theSecPoint.Y - theFirstPoint.Y ) ) + theFirstPoint.X; } } else { - aCoordY = theFirstPoint.Y(); + aCoordY = theFirstPoint.Y; - if ( !ValuesEquals( theFirstPoint.Y(), theSecPoint.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(); + aCoordY = ( ( ( thePoint.X() - theFirstPoint.X ) * ( theSecPoint.Y - theFirstPoint.Y ) ) / + ( theSecPoint.X - theFirstPoint.X ) ) + theFirstPoint.Y; } } - theResPoint.SetX( aCoordX ); - theResPoint.SetY( aCoordY ); + theResPoint.X = aCoordX; + theResPoint.Y = aCoordY; // Calculate coefficient for interpolation - double aLength = Sqrt( Pow( theSecPoint.Y() - theFirstPoint.Y(), 2 ) + - Pow( theSecPoint.X() - theFirstPoint.X(), 2 ) ); + 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; + aInterCoeff = ( theSecPoint.Z - theFirstPoint.Z ) / aLength; - double aNewLength = Sqrt( Pow( theResPoint.Y() - theFirstPoint.Y(), 2 ) + - Pow( theResPoint.X() - theFirstPoint.X(), 2 ) ); + double aNewLength = Sqrt( Pow( theResPoint.Y - theFirstPoint.Y, 2 ) + + Pow( theResPoint.X - theFirstPoint.X, 2 ) ); // Calculate interpolated value - double aResVal = theFirstPoint.Z() + aInterCoeff * aNewLength; + double aResVal = theFirstPoint.Z + aInterCoeff * aNewLength; - theResPoint.SetZ( aResVal ); + theResPoint.Z = aResVal; } +#ifndef LIGHT_MODE +bool interpolZtriangle(const gp_XY& point, vtkPolyData* delaunay2D, vtkIdList* triangle, double& z) +{ -double HYDROData_Bathymetry::GetAltitudeForPoint(const gp_XY& thePoint) const + int nbPts = triangle->GetNumberOfIds(); + if (nbPts != 3) + { + //DEBTRACE("not a triangle ?"); + return false; + } + vtkIdType s[3]; + double v[3][3]; // v[i][j] = j coordinate of node i + for (int i=0; i<3; i++) + { + s[i] = triangle->GetId(i); + delaunay2D->GetPoint(s[i],v[i]); + } + //DEBTRACE("triangle node id: " << s[0] << " " << s[1] << " " << s[2]); + //DEBTRACE("triangle node 0: " << v[0][0] << " " << v[0][1] << " " << v[0][2]); + //DEBTRACE("triangle node 1: " << v[1][0] << " " << v[1][1] << " " << v[1][2]); + //DEBTRACE("triangle node 2: " << v[2][0] << " " << v[2][1] << " " << v[2][2]); + + // compute barycentric coordinates (https://en.wikipedia.org/wiki/Barycentric_coordinate_system) + // det = (y2-y3)(x1-x3)+(x3-x2)(y1-y3) + 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 ?"); + return false; + } + + // l0 = ((y2-y3)(x -x3)+(x3-x2)(y -y3))/det + double l0 = (v[1][1]-v[2][1])*(point.X()-v[2][0]) + (v[2][0]-v[1][0])*(point.Y()-v[2][1]); + l0 = l0/det; + + // l1 = ((y3-y1)(x -x3)+(x1-x3)(y -y3))/det + double l1 = (v[2][1]-v[0][1])*(point.X()-v[2][0]) + (v[0][0]-v[2][0])*(point.Y()-v[2][1]); + l1 = l1/det; + + double l2 = 1 -l0 -l1; + //DEBTRACE("l0, l1, l2: " << l0 << " " << l1 << " " << l2); + + if ((l0>=0) && (l0<=1) && (l1>=0) && (l1<=1) && (l2>=0) && (l2<=1)) + { + z = l0*v[0][2] + l1*v[1][2] + l2*v[2][2]; + return true; + } + return false; +} +#endif + +double HYDROData_Bathymetry::GetAltitudeForPoint(const gp_XY& thePoint, int theMethod) const { - DEBTRACE("GetAltitudeForPoint p(" << thePoint.X() << ", " << thePoint.Y() << ")"); + //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 + HYDROData_QuadtreeNode* aQuadtree = GetQuadtreeNodes(); if (!aQuadtree) { - DEBTRACE(" no Quadtree"); + //DEBTRACE(" no Quadtree"); return aResAltitude; } - std::map dist2nodes; + std::map dist2nodes; aQuadtree->NodesAround(thePoint, dist2nodes, aQuadtree->getPrecision()); 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()); } - aQuadtree->NodesAround(thePoint, dist2nodes, 5.0); - if (dist2nodes.size()) - { - std::map::const_iterator it = dist2nodes.begin(); - aResAltitude = it->second->Z(); - DEBTRACE(" number of points found: " << dist2nodes.size() << " nearest z: " << aResAltitude); - } - else + std::map::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); + + // --- for coarse bathymetry clouds (when the TELEMAC mesh is more refined than the bathymetry cloud) + // interpolation is required. + // - get a Delaunay2D mesh on the bathymetry cloud, + // - get the triangle containing the point in the Delaunay2D mesh, + // - interpolate altitude + + bool isBathyInterpolRequired = false; + if (theMethod) + isBathyInterpolRequired =true; + +#ifndef LIGHT_MODE + if (isBathyInterpolRequired) { - DEBTRACE(" number of points found: 0"); + vtkPolyData* aDelaunay2D = GetVtkDelaunay2D(); + vtkIdList* cells= vtkIdList::New(); + cells->Allocate(64); + vtkIdList* points= vtkIdList::New(); + points->Allocate(64); + aDelaunay2D->GetPointCells(nodeIndex, cells); + vtkIdType nbCells = cells->GetNumberOfIds(); + //DEBTRACE(" triangles on nearest point: " << nbCells); + bool isInside = false; + for (int i=0; iGetCellPoints(cells->GetId(i), points); + double z = 0; + isInside = interpolZtriangle(thePoint, aDelaunay2D, points, z); + if (isInside) + { + aResAltitude = z; + //DEBTRACE(" interpolated z: " << z); + break; + } + } + if (!isInside) + { + // DEBTRACE(" point outside triangles, nearest z kept"); + } } - + #endif return aResAltitude; - - -// 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::Iterator anIter( anAltitudePoints ); -// for ( ; anIter.More(); anIter.Next() ) -// { -// const AltitudePoint& aPoint = anIter.Value(); -// -// 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 ] ); -// -// // 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 TCollection_AsciiString& theFilePath ) @@ -426,6 +458,19 @@ void HYDROData_Bathymetry::SetFilePath( const TCollection_AsciiString& theFilePa TDataStd_AsciiString::Set( myLab.FindChild( DataTag_FilePath ), theFilePath ); } +void HYDROData_Bathymetry::SetFilePaths( const QStringList& theFilePaths ) +{ + int i = 1; + Handle_TDataStd_ExtStringArray TExtStrArr = TDataStd_ExtStringArray::Set( myLab.FindChild( DataTag_FilePaths ), 1, theFilePaths.size() ); + foreach (QString filepath, theFilePaths) + { + std::string sstr = filepath.toStdString(); + const char* Val = sstr.c_str(); + TExtStrArr->SetValue(i, TCollection_ExtendedString(Val)); + i++; + } +} + TCollection_AsciiString HYDROData_Bathymetry::GetFilePath() const { TCollection_AsciiString aRes; @@ -437,10 +482,52 @@ TCollection_AsciiString HYDROData_Bathymetry::GetFilePath() const if ( aLabel.FindAttribute( TDataStd_AsciiString::GetID(), anAsciiStr ) ) aRes = anAsciiStr->Get(); } + else + { + aLabel = myLab.FindChild( DataTag_FilePaths, false ); + if ( !aLabel.IsNull() ) + { + Handle(TDataStd_ExtStringArray) anExtStrArr; + if ( aLabel.FindAttribute( TDataStd_ExtStringArray::GetID(), anExtStrArr ) ) + aRes = anExtStrArr->Value(1); //try take the first; convert extstring to asciistring + } + } return aRes; } +QStringList HYDROData_Bathymetry::GetFilePaths() const +{ + QStringList aResL; + + TDF_Label aLabel = myLab.FindChild( DataTag_FilePaths, false ); + if ( !aLabel.IsNull() ) + { + Handle(TDataStd_ExtStringArray) anExtStrArr; + if ( aLabel.FindAttribute( TDataStd_ExtStringArray::GetID(), anExtStrArr ) ) + { + for (int i = anExtStrArr->Lower(); i <= anExtStrArr->Upper(); i++ ) + { + Standard_ExtString str = anExtStrArr->Value(i).ToExtString(); + TCollection_AsciiString aText (str); + aResL << QString(aText.ToCString()); + } + } + } + else //backward compatibility + { + TDF_Label anOldLabel = myLab.FindChild( DataTag_FilePath, false ); + if ( !anOldLabel.IsNull() ) + { + Handle(TDataStd_AsciiString) anAsciiStr; + if ( anOldLabel.FindAttribute( TDataStd_AsciiString::GetID(), anAsciiStr ) ) + aResL << QString(anAsciiStr->Get().ToCString()); + } + } + + return aResL; +} + void HYDROData_Bathymetry::SetAltitudesInverted( const bool theIsInverted, const bool theIsUpdate ) { @@ -456,15 +543,15 @@ void HYDROData_Bathymetry::SetAltitudesInverted( const bool theIsInverted, return; // Update altitude points - AltitudePoints anAltitudePoints = GetAltitudePoints(); - if ( anAltitudePoints.IsEmpty() ) + HYDROData_Bathymetry::AltitudePoints anAltitudePoints = GetAltitudePoints(); + if ( anAltitudePoints.empty() ) return; - AltitudePoints::Iterator anIter( anAltitudePoints ); - for ( ; anIter.More(); anIter.Next() ) + HYDROData_Bathymetry::AltitudePoints::iterator anIter = anAltitudePoints.begin(), aLast = anAltitudePoints.end(); + for ( ; anIter!=aLast; ++anIter ) { - AltitudePoint& aPoint = anIter.ChangeValue(); - aPoint.SetZ( aPoint.Z() * -1 ); + HYDROData_Bathymetry::AltitudePoint& aPoint = *anIter; + aPoint.Z *= -1; } SetAltitudePoints( anAltitudePoints ); @@ -485,50 +572,63 @@ bool HYDROData_Bathymetry::IsAltitudesInverted() const return aRes; } -bool HYDROData_Bathymetry::ImportFromFile( const TCollection_AsciiString& theFileName ) +bool HYDROData_Bathymetry::ImportFromFile( const QString& theFileName ) { - // Try to open the file - QFile aFile( theFileName.ToCString() ); - if ( !aFile.exists() || !aFile.open( QIODevice::ReadOnly ) ) - return false; + return ImportFromFiles(QStringList(theFileName)); +} - bool aRes = false; +bool HYDROData_Bathymetry::ImportFromFiles( const QStringList& theFileNames ) +{ + AltitudePoints AllPoints; + bool Stat = false; - QString aFileSuf = QFileInfo( aFile ).suffix().toLower(); + foreach (QString theFileName, theFileNames) + { + // Try to open the file + QFile aFile( theFileName ); + if ( !aFile.exists() || !aFile.open( QIODevice::ReadOnly ) ) + continue; - AltitudePoints aPoints; + QString aFileSuf = QFileInfo( aFile ).suffix().toLower(); - // Try to import the file - if ( aFileSuf == "xyz" ) - aRes = importFromXYZFile( aFile, aPoints ); - else if ( aFileSuf == "asc" ) - aRes = importFromASCFile( aFile, aPoints ); + HYDROData_Bathymetry::AltitudePoints aPoints; - // Close the file - aFile.close(); - + // Try to import the file + if ( aFileSuf == "xyz" ) + Stat = importFromXYZFile( aFile, aPoints ); + else if ( aFileSuf == "asc" ) + Stat = importFromASCFile( aFile, aPoints ); + + if (!Stat) + continue; //ignore this points + + // Close the file + aFile.close(); + + AllPoints.insert(AllPoints.end(), aPoints.begin(), aPoints.end()); + } // Convert from global to local CS - Handle_HYDROData_Document aDoc = HYDROData_Document::Document( myLab ); - AltitudePoints::Iterator anIter( aPoints ); - for ( ; anIter.More(); anIter.Next() ) + Handle(HYDROData_Document) aDoc = HYDROData_Document::Document( myLab ); + HYDROData_Bathymetry::AltitudePoints::iterator anIter = AllPoints.begin(), aLast = AllPoints.end(); + for ( ; anIter!=aLast; ++anIter ) { - AltitudePoint& aPoint = anIter.ChangeValue(); - aDoc->Transform( aPoint, true ); + HYDROData_Bathymetry::AltitudePoint& aPoint = *anIter; + aDoc->Transform( aPoint.X, aPoint.Y, aPoint.Z, true ); } - if ( aRes ) + if ( Stat ) { // Update file path and altitude points of this Bathymetry - SetFilePath( theFileName ); - SetAltitudePoints( aPoints ); + SetFilePaths (theFileNames ); + SetAltitudePoints( AllPoints ); } - return aRes && !aPoints.IsEmpty(); + return Stat && !AllPoints.empty(); } bool HYDROData_Bathymetry::importFromXYZFile( QFile& theFile, - AltitudePoints& thePoints ) const + HYDROData_Bathymetry::AltitudePoints& thePoints ) const { if ( !theFile.isOpen() ) return false; @@ -554,7 +654,7 @@ bool HYDROData_Bathymetry::importFromXYZFile( QFile& theFile, if ( aValues.length() < 3 ) return false; - AltitudePoint aPoint; + HYDROData_Bathymetry::AltitudePoint aPoint; QString anX = aValues.value( 0 ); QString anY = aValues.value( 1 ); @@ -562,23 +662,26 @@ bool HYDROData_Bathymetry::importFromXYZFile( QFile& theFile, bool isXOk = false, isYOk = false, isZOk = false; - aPoint.SetX( anX.toDouble( &isXOk ) ); - aPoint.SetY( anY.toDouble( &isYOk ) ); - aPoint.SetZ( aZ.toDouble( &isZOk ) ); + aPoint.X = anX.toDouble( &isXOk ); + aPoint.Y = anY.toDouble( &isYOk ); + aPoint.Z = aZ.toDouble( &isZOk ); if ( !isXOk || !isYOk || !isZOk ) return false; - if ( boost::math::isnan( aPoint.X() ) || boost::math::isinf( aPoint.X() ) || - boost::math::isnan( aPoint.Y() ) || boost::math::isinf( aPoint.Y() ) || - boost::math::isnan( aPoint.Z() ) || boost::math::isinf( aPoint.Z() ) ) + if ( HYDROData_Tool::IsNan( aPoint.X ) || HYDROData_Tool::IsInf( aPoint.X ) || + HYDROData_Tool::IsNan( aPoint.Y ) || HYDROData_Tool::IsInf( aPoint.Y ) || + HYDROData_Tool::IsNan( aPoint.Z ) || HYDROData_Tool::IsInf( aPoint.Z ) ) return false; // Invert the z value if requested if ( anIsAltitudesInverted ) - aPoint.SetZ( -aPoint.Z() ); + aPoint.Z = -aPoint.Z; + + if( thePoints.size()>=thePoints.capacity() ) + thePoints.reserve( thePoints.size()+BLOCK_SIZE ); - thePoints.Append( aPoint ); + thePoints.push_back( aPoint ); } #ifdef _TIMER @@ -591,7 +694,7 @@ bool HYDROData_Bathymetry::importFromXYZFile( QFile& theFile, } bool HYDROData_Bathymetry::importFromASCFile( QFile& theFile, - AltitudePoints& thePoints ) const + HYDROData_Bathymetry::AltitudePoints& thePoints ) const { if ( !theFile.isOpen() ) return false; @@ -662,31 +765,31 @@ bool HYDROData_Bathymetry::importFromASCFile( QFile& theFile, { if (aStrList[j].toDouble() != aNoDataValue) { - AltitudePoint aPoint; - aPoint.SetX(anXllCorner + aCellSize*(j + 0.5)); - aPoint.SetY(anYllCorner + aCellSize*(aNRows - i + 0.5)); - aPoint.SetZ(aStrList[j].toDouble()); + HYDROData_Bathymetry::AltitudePoint aPoint; + aPoint.X = anXllCorner + aCellSize*(j + 0.5); + aPoint.Y = anYllCorner + aCellSize*(aNRows - i + 0.5); + aPoint.Z = aStrList[j].toDouble(); if ( anIsAltitudesInverted ) - aPoint.SetZ( -aPoint.Z() ); + aPoint.Z = -aPoint.Z; - thePoints.Append(aPoint); + if( thePoints.size()>=thePoints.capacity() ) + thePoints.reserve( thePoints.size()+BLOCK_SIZE ); + thePoints.push_back(aPoint); } } i++; - } return true; - } -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; @@ -698,14 +801,14 @@ Handle_HYDROData_PolylineXY HYDROData_Bathymetry::CreateBoundaryPolyline() const double Xmin = 0.0, Xmax = 0.0, Ymin = 0.0, Ymax = 0.0; bool isFirst = true; - AltitudePoints aPoints = GetAltitudePoints(); + HYDROData_Bathymetry::AltitudePoints aPoints = GetAltitudePoints(); - AltitudePoints::Iterator anIter( aPoints ); - for ( ; anIter.More(); anIter.Next() ) + HYDROData_Bathymetry::AltitudePoints::const_iterator anIter = aPoints.begin(), aLast = aPoints.end(); + for ( ; anIter!=aLast; ++anIter ) { - const AltitudePoint& aPoint = anIter.Value(); + const HYDROData_Bathymetry::AltitudePoint& aPoint = *anIter; - double x = aPoint.X(), y = aPoint.Y(); + double x = aPoint.X, y = aPoint.Y; if( isFirst || xXmax ) @@ -733,12 +836,14 @@ Handle_HYDROData_PolylineXY HYDROData_Bathymetry::CreateBoundaryPolyline() const void HYDROData_Bathymetry::UpdateLocalCS( double theDx, double theDy ) { gp_XYZ aDelta( theDx, theDy, 0 ); - AltitudePoints aPoints = GetAltitudePoints(); - AltitudePoints::Iterator anIter( aPoints ); - for ( int i = 0 ; anIter.More(); ++i, anIter.Next() ) + HYDROData_Bathymetry::AltitudePoints aPoints = GetAltitudePoints(); + HYDROData_Bathymetry::AltitudePoints::iterator anIter = aPoints.begin(), aLast = aPoints.end(); + for ( int i = 0; anIter!=aLast; ++i, ++anIter ) { - AltitudePoint& aPoint = anIter.ChangeValue(); - aPoint += aDelta; + HYDROData_Bathymetry::AltitudePoint& aPoint = *anIter; + aPoint.X += aDelta.X(); + aPoint.Y += aDelta.Y(); + aPoint.Z += aDelta.Z(); } SetAltitudePoints( aPoints ); }