#include "HYDROData_Tool.h"
#include "HYDROData_PolylineXY.h"
-#include <boost/math/special_functions/fpclassify.hpp>
-
#include <gp_XY.hxx>
#include <gp_XYZ.hxx>
#include <QPolygonF>
#include <QStringList>
+#include <vtkPoints.h>
+#include <vtkDelaunay2D.h>
+#include <vtkPolyData.h>
+#include <vtkSmartPointer.h>
+#include <vtkIdList.h>
+
+#include <iostream>
+
#include <math.h>
// #define _TIMER
//HYDROData_QuadtreeNode* HYDROData_Bathymetry::myQuadtree = 0;
std::map<int, HYDROData_QuadtreeNode*> HYDROData_Bathymetry::myQuadtrees;
+std::map<int, vtkPolyData*> HYDROData_Bathymetry::myDelaunay2D;
+
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,
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;
if (aLabel.IsNull())
return 0;
int labkey = myLab.Tag();
- int altkey = aLabel.Tag();
+ //int altkey = aLabel.Tag();
//DEBTRACE("GetQuadtreeNodes this labkey altkey "<<this<<" "<<labkey<<" "<<altkey);
-// if (myQuadtree->isEmpty() )
+ // if (myQuadtree->isEmpty() )
if (myQuadtrees.find(labkey) == myQuadtrees.end())
{
DEBTRACE("GetQuadtreeNodes init " << this << " " << labkey);
Nodes_3D* aListOfNodes = new Nodes_3D();
+ int index =0;
for (int i = aCoordsArray->Lower(), n = aCoordsArray->Upper(); i <= n;)
{
if (i + 3 > n + 1)
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");
return myQuadtrees[labkey];
}
+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 "<<this<<" "<<labkey<<" "<<altkey);
+ if (myDelaunay2D.find(labkey) == myDelaunay2D.end())
+ {
+ DEBTRACE("GetVtkDelaunay2D init " << this << " " << labkey);
+
+ TDF_Label aLabel = myLab.FindChild(DataTag_AltitudePoints, false);
+ if (aLabel.IsNull())
+ return 0;
+ Handle(TDataStd_RealArray) aCoordsArray;
+ if (!aLabel.FindAttribute(TDataStd_RealArray::GetID(), aCoordsArray))
+ return 0;
+
+ vtkPoints *points = vtkPoints::New();
+ points->Allocate(aCoordsArray->Upper() +1);
+ for (int i = aCoordsArray->Lower(), n = aCoordsArray->Upper(); i <= n;)
+ {
+ 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];
+
+}
+
void HYDROData_Bathymetry::RemoveAltitudePoints()
{
TDF_Label aLabel = myLab.FindChild( DataTag_AltitudePoints, false );
theResPoint.SetZ( aResVal );
}
-double HYDROData_Bathymetry::GetAltitudeForPoint(const gp_XY& thePoint) const
+bool interpolZtriangle(const gp_XY& point, vtkPolyData* delaunay2D, vtkIdList* triangle, double& z)
{
- DEBTRACE("GetAltitudeForPoint p(" << thePoint.X() << ", " << thePoint.Y() << ")");
+
+ 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;
+}
+
+double HYDROData_Bathymetry::GetAltitudeForPoint(const gp_XY& thePoint, int theMethod) const
+{
+ 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)
{
return aResAltitude;
}
- std::map<double, const gp_XYZ*> dist2nodes;
+ std::map<double, const gpi_XYZ*> dist2nodes;
aQuadtree->NodesAround(thePoint, dist2nodes, aQuadtree->getPrecision());
while (dist2nodes.size() == 0)
{
DEBTRACE("adjust precision to: " << aQuadtree->getPrecision());
aQuadtree->NodesAround(thePoint, dist2nodes, aQuadtree->getPrecision());
}
- aQuadtree->NodesAround(thePoint, dist2nodes, 5.0);
- if (dist2nodes.size())
+ std::map<double, const gpi_XYZ*>::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;
+ if (isBathyInterpolRequired)
{
- std::map<double, const gp_XYZ*>::const_iterator it = dist2nodes.begin();
- aResAltitude = it->second->Z();
- DEBTRACE(" number of points found: " << dist2nodes.size() << " nearest z: " << aResAltitude);
- }
- else
- {
- 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; i<nbCells; i++)
+ {
+ aDelaunay2D->GetCellPoints(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");
}
-
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 )
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
