From d3e944f8eb3794a42e8228c9c40d777bf9ae0409 Mon Sep 17 00:00:00 2001 From: Paul RASCLE Date: Fri, 29 Apr 2016 11:24:57 +0200 Subject: [PATCH] z interpolation on Delaunay triangulation of bathymetry: prototype --- src/HYDROData/HYDROData_Bathymetry.cxx | 357 +++++++++++++++---------- src/HYDROData/HYDROData_Bathymetry.h | 20 ++ 2 files changed, 232 insertions(+), 145 deletions(-) diff --git a/src/HYDROData/HYDROData_Bathymetry.cxx b/src/HYDROData/HYDROData_Bathymetry.cxx index bc7e75b4..68d14daf 100644 --- a/src/HYDROData/HYDROData_Bathymetry.cxx +++ b/src/HYDROData/HYDROData_Bathymetry.cxx @@ -35,6 +35,16 @@ #include #include +#include +#include +#include +#include +#include +#include +#include + +#include + #include // #define _TIMER @@ -50,6 +60,37 @@ IMPLEMENT_STANDARD_RTTIEXT(HYDROData_Bathymetry, HYDROData_IAltitudeObject) //HYDROData_QuadtreeNode* HYDROData_Bathymetry::myQuadtree = 0; std::map HYDROData_Bathymetry::myQuadtrees; +std::map HYDROData_Bathymetry::myDelaunay2D; + + +hydroDelaunay2D* hydroDelaunay2D::New() +{ + DEBTRACE("hydroDelaunay2D::New()"); + return new hydroDelaunay2D(); +} +hydroDelaunay2D::hydroDelaunay2D() : + vtkDelaunay2D() +{ + DEBTRACE("hydroDelaunay2D"); +} + +hydroDelaunay2D::~hydroDelaunay2D() +{ + DEBTRACE("~hydroDelaunay2D"); +} + +int hydroDelaunay2D::RequestData(vtkInformation *request , vtkInformationVector **inv, vtkInformationVector *outv) +{ + DEBTRACE("hydroDelaunay2D::RequestData"); + vtkDelaunay2D::RequestData(request, inv, outv); +} + +int hydroDelaunay2D::ProcessRequest(vtkInformation *request , vtkInformationVector **inv, vtkInformationVector *outv) +{ + DEBTRACE("hydroDelaunay2D::ProcessRequest"); + vtkDelaunay2D::ProcessRequest(request, inv, outv); +} + HYDROData_Bathymetry::HYDROData_Bathymetry() : HYDROData_IAltitudeObject() @@ -152,9 +193,9 @@ 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); @@ -170,6 +211,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) @@ -178,7 +220,8 @@ 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"); @@ -189,6 +232,83 @@ HYDROData_QuadtreeNode* HYDROData_Bathymetry::GetQuadtreeNodes() const 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 "<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()); + + hydroDelaunay2D* delaunay2D = hydroDelaunay2D::New(); + delaunay2D->GetInputPortInformation(0)->Print(std::cerr); + DEBTRACE("set the input data"); + delaunay2D->SetInputData(profile); + delaunay2D->GetInputPortInformation(0)->Print(std::cerr); + delaunay2D->GetOutputPortInformation(0)->Print(std::cerr); + DEBTRACE("---"); + delaunay2D->GetOutputInformation(0)->Print(std::cerr); + delaunay2D->GetExecutive()->GetOutputInformation(0)->Print(std::cerr); + DEBTRACE("---"); + vtkInformationVector** inv = delaunay2D->GetExecutive()->GetInputInformation(); + vtkInformationVector* outv = delaunay2D->GetExecutive()->GetOutputInformation(); + delaunay2D->GetExecutive()->GetInputInformation(0,0)->Print(std::cerr); + delaunay2D->GetExecutive()->GetOutputInformation(0)->Print(std::cerr); + vtkInformation* request = vtkInformation::New(); + + //delaunay2D->ProcessRequest(request, inv,outv); + + delaunay2D->GetOutputPortInformation(0)->Print(std::cerr); + delaunay2D->GetExecutive()->GetInputInformation(0,0)->Print(std::cerr); + delaunay2D->GetExecutive()->GetOutputInformation(0)->Print(std::cerr); + + DEBTRACE("---"); + delaunay2D->Update(); + + DEBTRACE("---"); + delaunay2D->UpdateDataObject(); + + vtkPolyData* data = delaunay2D->GetOutput(); + + data->BuildLinks(); + myDelaunay2D[labkey] = data; + DEBTRACE("---"); + return data; + } + else + return myDelaunay2D[labkey]; + +} + void HYDROData_Bathymetry::RemoveAltitudePoints() { TDF_Label aLabel = myLab.FindChild( DataTag_AltitudePoints, false ); @@ -252,12 +372,62 @@ void interpolateAltitudeForPoints( const gp_XY& th theResPoint.SetZ( aResVal ); } +bool interpolZtriangle(const gp_XY& point, vtkPolyData* delaunay2D, vtkIdList* triangle, double& z) +{ + + 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) and (l0<=1) and (l1>=0) and (l1<=1) and (l2>=0) and (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) const { DEBTRACE("GetAltitudeForPoint p(" << thePoint.X() << ", " << thePoint.Y() << ")"); double anInvalidAltitude = GetInvalidAltitude(); double aResAltitude = anInvalidAltitude; + // --- find the nearest point in the bathymetry cloud, with quadtree + HYDROData_QuadtreeNode* aQuadtree = GetQuadtreeNodes(); if (!aQuadtree) { @@ -265,7 +435,7 @@ double HYDROData_Bathymetry::GetAltitudeForPoint(const gp_XY& thePoint) const return aResAltitude; } - std::map dist2nodes; + std::map dist2nodes; aQuadtree->NodesAround(thePoint, dist2nodes, aQuadtree->getPrecision()); while (dist2nodes.size() == 0) { @@ -273,150 +443,47 @@ double HYDROData_Bathymetry::GetAltitudeForPoint(const gp_XY& thePoint) const 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 + //aQuadtree->NodesAround(thePoint, dist2nodes, 5.0); + 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 = true; + 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); + DEBTRACE("---"); + aDelaunay2D->GetPointCells(nodeIndex, cells); + vtkIdType nbCells = cells->GetNumberOfIds(); + DEBTRACE(nbCells); + for (int i=0; iGetCellPoints(cells->GetId(i), points); + double z = 0; + if (interpolZtriangle(thePoint, aDelaunay2D, points, z)) + { + aResAltitude = z; + DEBTRACE("interpolated z: " << z); + break; + } + else + { + 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 ) diff --git a/src/HYDROData/HYDROData_Bathymetry.h b/src/HYDROData/HYDROData_Bathymetry.h index 0987afbe..e2828c85 100644 --- a/src/HYDROData/HYDROData_Bathymetry.h +++ b/src/HYDROData/HYDROData_Bathymetry.h @@ -21,11 +21,27 @@ #include "HYDROData_IAltitudeObject.h" #include "HYDROData_QuadtreeNode.hxx" +#include +#include + +#include +class hydroDelaunay2D: public vtkDelaunay2D +{ +public: + static hydroDelaunay2D *New(); + virtual int ProcessRequest(vtkInformation*, vtkInformationVector**, vtkInformationVector*); +protected: + hydroDelaunay2D(); + ~hydroDelaunay2D(); + virtual int RequestData(vtkInformation *, vtkInformationVector **, vtkInformationVector *); +}; class QFile; class gp_XYZ; +class gp_XY; class Handle_HYDROData_PolylineXY; + DEFINE_STANDARD_HANDLE(HYDROData_Bathymetry, HYDROData_IAltitudeObject) @@ -86,6 +102,9 @@ public: HYDRODATA_EXPORT virtual AltitudePoints GetAltitudePoints(bool IsConvertToGlobal = false) const; HYDRODATA_EXPORT virtual HYDROData_QuadtreeNode* GetQuadtreeNodes() const; + HYDRODATA_EXPORT virtual vtkPolyData* GetVtkDelaunay2D() const; + //HYDRODATA_EXPORT bool interpolZtriangle(const gp_XY& point, vtkIdList* triangle, double* z); + /** * Remove all altitude points. */ @@ -145,6 +164,7 @@ private: bool importFromXYZFile( QFile& theFile, AltitudePoints& thePoints ) const; static std::map myQuadtrees; + static std::map myDelaunay2D; bool importFromASCFile( QFile& theFile, AltitudePoints& thePoints ) const; -- 2.39.2