How to evaluate the min distance between a point and a Poly_Triangulation

Hi All,

Any hint is welcome. Thanks in advance.

Evgeny Lodyzhehsky's picture

Dear Cauchy Ding.
An object of type Poly_Triangulation contains nodes, and triangles.
Each node is just 3D point.
Each triangle contain indices of 3 nodes above.

- Try to explode Poly_Triangulation on nodes and find the min distance between your point and node point(i) .

- To make more precisable you can find the min distance between a your point and triangle .
___
PS.
You will face with the following insuperable obstacles:
The formula for distance between two 3D poins is top secret thing.
Nobody know how to implement the algorithm that calculates the distance between 3D point and triangle.

Cauchy Ding's picture

Dear Lodyzhensky Evgeny Nicolaich,

Thank you for your reply firstly. In fact the definition of Distance(Point,Mesh) is somewhat fuzzy. I have searched several Computer Geometry book, but no found. But I think the distance could be defined as follows:

Def 1. The min distance of between point and 3 Nodes of triangle. dis1 = Min( Dis(Point, N1), Dis(Point, N2), Dis(Point, N3))

Def 2. Project point to 3 edge of triangle and got three projected points, they may locate outside of the range of edge. If it's outside, abandon the projected point.
dis2 = Min( bOutside1?MAX:Dis(Point,E1), bOutside2?MAX:Dis(Point,E2), bOutside3?MAX:Dis(Point,E3));

Def 3. Project point to the triangle, if the projected point is outside of the triangle region, abandon it.
dis3 = Min(bOutside?MAX:Dis(Point, T));

so, the min distance between a point and a triangle is dis4 = min(dis1, dis2, dis3).
And the final min distance between a point and a mesh if min{dis4}.

BTW: To evaluate the min distance between a point and a mesh is just because I want to evaluete the overlapped region between two neighbour meshes.

Thanks a lot.

Evgeny Lodyzhehsky's picture

Dear Cauchy Ding.

Another approach can be the following:

1. Find nearest node Nn using binary search algo. for e.g. To sort Bnd_Box(-es) you can use the class NCollection_UBTreeFiller. So you will have the value of min distance Dmin as a first approximation

2. Find all triangles (T1, T2...Tn) that have Nn as shared node. Use TColStd_DataMapOfIntegerInteger as auxiliary tool.
3. For each triangle Ti (i=1,2...n) {
3.1. Obtain the plane of Ti PL
3.2 Project your point (P) on PL; //you will have projection point Px
3.3 Classify Px for triangle Ti ; // you will have state IN,ON,OUT
3.4 In case when te state=(ON or IN) {
obtain the distance between P and Px // Di
}
3.5 Find minimum distance if(Di

Cauchy Ding's picture

Dear Lodyzhensky Evgeny Nicolaich,

Thank you so much for your detailed explanation.
I find another method referred in , chapter 10.3.2 "POINT TO TRIANGLE".

-Cauchy Ding

Hugues Delorme's picture

Hello Cauchy Ding,

I have taken a look at the article you found (http://www.geometrictools.com/Documentation/Documentation.html). Did you succeed to implement it ?

I'm a bit confused with how a triangle is defined ( T(s,t) = B + sE0 + tE1 ).

Cauchy Ding's picture

Hello Hugues,

I implement it, but I haven't tested it using much samples. If you find any error, please tell me freely.

double SolidGeomMain::Mesh_SquarDis_PointToTriangle(const gp_Pnt& p, const gp_Pnt& v0, const gp_Pnt& v1, const gp_Pnt& v2, gp_Pnt& prjPnt, bool *pInside)
{
gp_XYZ e0, e1, D;
double a, b, c, d, e, f, det, s, t;
int region;

e0 = v1.XYZ()-v0.XYZ();
e1 = v2.XYZ()-v0.XYZ();
D = v0.XYZ()-p.XYZ();

a = e0.Dot(e0);
b = e0.Dot(e1);
c = e1.Dot(e1);
d = e0.Dot(D);
e = e1.Dot(D);
f = D.Dot(D);

det = a*c-b*b;
s = b*e - c*d;
t = b*d - a*e;

if(pInside)
*pInside = false;

region = 0;
if(s+t<=det)
{
if(s<0)
{
if(t<0)
region = 4;
else
region = 3;
}
else if(t<0)
region = 5;
}
else
{
if(s<0)
region = 2;
else if(t<0)
region = 6;
else
region = 1;
}

if(region==0)
{
double invDet = 1/det;
s*=invDet;
t*=invDet;
}
else if(region==1)
{
double numer = c+e-b-d;
if(numer<=0)
s = 0;
else
{
double denom = a - 2*b +c;
s = (numer>=denom?1:numer/denom);
}
t = 1-s;
}
else if(region==3)
{
s = 0;
t = (e>=0?0:(-e>=c?1:-e/c));
}
else if(region==5)
{
t = 0;
s = (d>=0?0:(-d>=a?1:-d/a));
}
else if(region==2) //notice: 2, 4, 6 is not same with the method refered in Geometric Tools For Computer Graphics
{
s = 0;
t = 1;
}
else if(region==4)
{
s = 0;
t = 0;
}
else if(region==6)
{
s = 1;
t = 0;
}

if(pInside && region==0)
*pInside = true;

prjPnt = v0.XYZ()+e0*s+e1*t;
return prjPnt.SquareDistance(p);
}

Hugues Delorme's picture

Thanks for the posting. I am going to exercise with NCollection_UBTree to find the triangles being candidates. And use your code to finally get the minimal distance.

Hugues Delorme's picture

I have implemented the binary search with NCollection_UBTree as Nicolaich suggested it. It works pretty well. The "object type" of my UB-tree is std::pair.
The first item of the pair is an index of a node that is valid in the triangulation hold in the second item of the pair.
This way I can fill the UB-tree with many triangulations (needed when you want to compute the distance between a point and a set of faces).

The "bound type" is simply Bnd_Box, for each node a Bnd_Box is created and associated.

The search of the interesting node is done with a specialized selector : I give the code here, hope it will helpful :

typedef std::pair
NodeIndexInTriangulation_t;
typedef NCollection_UBTree
UBTreeOfNodeIndices_t;
typedef NCollection_UBTreeFiller
UBTreeOfNodeIndicesFiller_t;

class NodeBndBoxSelector : public UBTreeOfNodeIndices_t::Selector
{
public:
NodeBndBoxSelector(const gp_Pnt& pntToProject) :
_pntToProject(pntToProject),
_currMinDist(std::numeric_limits::max()),
_currMinDistNodeId(-1, Handle_Poly_Triangulation())
{
}

virtual Standard_Boolean Reject(const Bnd_Box& bb) const
{
const bool result = bb.IsOpenXmin() == Standard_True ||
bb.IsOpenXmax() == Standard_True ||
bb.IsOpenYmin() == Standard_True ||
bb.IsOpenYmax() == Standard_True ||
bb.IsOpenZmin() == Standard_True ||
bb.IsOpenZmax() == Standard_True ||
bb.IsWhole() == Standard_True ||
bb.IsVoid() == Standard_True;
return result ? Standard_True : Standard_False;
}

virtual Standard_Boolean Accept(const NodeIndexInTriangulation_t& nodeId)
{
if (nodeId.second.IsNull())
return Standard_False;
const TColgp_Array1OfPnt& nodes = nodeId.second->Nodes();
if (!(nodes.Lower() <= nodeId.first && nodeId.first <= nodes.Upper()))
return Standard_False;
const gp_Pnt& pnt = nodes(nodeId.first);
const double dist = this->_pntToProject.SquareDistance(pnt);
if (dist < this->_currMinDist || this->_currMinDistNodeId.first == -1)
{
this->_currMinDistNodeId = nodeId;
this->_currMinDist = dist;
return Standard_True;
}
return Standard_False;
}

const double minDistance() const
{
return std::sqrt(this->_currMinDist);
}

const NodeIndexInTriangulation_t& minDistanceNodeIndex() const
{
return this->_currMinDistNodeId;
}

private:
const gp_Pnt _pntToProject;
double _currMinDist;
NodeIndexInTriangulation_t _currMinDistNodeId;
}; // class NodeBndBoxSelector