Commit 0cb54edb authored by Yohann Bearzi's avatar Yohann Bearzi

bug fix in vtkQuad::IntersectWithLine

parent 208a56c6
Pipeline #181296 failed with stage
......@@ -16,6 +16,7 @@
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkDataArrayRange.h"
#include "vtkIncrementalPointLocator.h"
#include "vtkLine.h"
#include "vtkMath.h"
......@@ -523,100 +524,38 @@ vtkCell* vtkQuad::GetEdge(int edgeId)
int vtkQuad::IntersectWithLine(const double p1[3], const double p2[3], double tol, double& t,
double x[3], double pcoords[3], int& subId)
{
int diagonalCase;
double d1 = vtkMath::Distance2BetweenPoints(this->Points->GetPoint(0), this->Points->GetPoint(2));
double d2 = vtkMath::Distance2BetweenPoints(this->Points->GetPoint(1), this->Points->GetPoint(3));
subId = 0;
// Figure out how to uniquely tessellate the quad. Watch out for
// equivalent triangulations (i.e., the triangulation is equivalent
// no matter where the diagonal). In this case use the point ids as
// a tie breaker to ensure unique triangulation across the quad.
//
if (d1 == d2) // rare case; discriminate based on point id
{
int i, id, maxId = 0, maxIdx = 0;
for (i = 0; i < 4; i++) // find the maximum id
{
if ((id = this->PointIds->GetId(i)) > maxId)
{
maxId = id;
maxIdx = i;
}
}
if (maxIdx == 0 || maxIdx == 2)
diagonalCase = 0;
else
diagonalCase = 1;
}
else if (d1 < d2)
{
diagonalCase = 0;
}
else // d2 < d1
{
diagonalCase = 1;
}
// Note: in the following code the parametric coords must be adjusted to
// reflect the use of the triangle parametric coordinate system.
switch (diagonalCase)
this->Triangle->Points->SetPoint(0, this->Points->GetPoint(0));
this->Triangle->Points->SetPoint(1, this->Points->GetPoint(1));
this->Triangle->Points->SetPoint(2, this->Points->GetPoint(2));
if (this->Triangle->IntersectWithLine(p1, p2, tol, t, x, pcoords, subId))
{
case 0:
this->Triangle->Points->SetPoint(0, this->Points->GetPoint(0));
this->Triangle->Points->SetPoint(1, this->Points->GetPoint(1));
this->Triangle->Points->SetPoint(2, this->Points->GetPoint(2));
if (this->Triangle->IntersectWithLine(p1, p2, tol, t, x, pcoords, subId))
{
pcoords[0] = pcoords[0] + pcoords[1];
return 1;
}
this->Triangle->Points->SetPoint(0, this->Points->GetPoint(2));
this->Triangle->Points->SetPoint(1, this->Points->GetPoint(3));
this->Triangle->Points->SetPoint(2, this->Points->GetPoint(0));
if (this->Triangle->IntersectWithLine(p1, p2, tol, t, x, pcoords, subId))
{
pcoords[0] = 1.0 - (pcoords[0] + pcoords[1]);
pcoords[1] = 1.0 - pcoords[1];
return 1;
}
return 0;
case 1:
this->Triangle->Points->SetPoint(0, this->Points->GetPoint(0));
this->Triangle->Points->SetPoint(1, this->Points->GetPoint(1));
this->Triangle->Points->SetPoint(2, this->Points->GetPoint(3));
if (this->Triangle->IntersectWithLine(p1, p2, tol, t, x, pcoords, subId))
{
return 1;
}
this->Triangle->Points->SetPoint(0, this->Points->GetPoint(2));
this->Triangle->Points->SetPoint(1, this->Points->GetPoint(3));
this->Triangle->Points->SetPoint(2, this->Points->GetPoint(1));
if (this->Triangle->IntersectWithLine(p1, p2, tol, t, x, pcoords, subId))
{
pcoords[0] = 1.0 - pcoords[0];
pcoords[1] = 1.0 - pcoords[1];
return 1;
}
return 0;
pcoords[0] = pcoords[0] + pcoords[1];
return 1;
}
this->Triangle->Points->SetPoint(0, this->Points->GetPoint(2));
this->Triangle->Points->SetPoint(1, this->Points->GetPoint(3));
this->Triangle->Points->SetPoint(2, this->Points->GetPoint(0));
if (this->Triangle->IntersectWithLine(p1, p2, tol, t, x, pcoords, subId))
{
pcoords[0] = 1.0 - (pcoords[0] + pcoords[1]);
pcoords[1] = 1.0 - pcoords[1];
return 1;
}
return 0;
}
//------------------------------------------------------------------------------
int vtkQuad::Triangulate(int vtkNotUsed(index), vtkIdList* ptIds, vtkPoints* pts)
{
double d1, d2;
pts->Reset();
ptIds->Reset();
// use minimum diagonal (Delaunay triangles) - assumed convex
d1 = vtkMath::Distance2BetweenPoints(this->Points->GetPoint(0), this->Points->GetPoint(2));
d2 = vtkMath::Distance2BetweenPoints(this->Points->GetPoint(1), this->Points->GetPoint(3));
auto pointRange = vtk::DataArrayTupleRange<3>(this->Points->GetData());
double d1 = vtkMath::Distance2BetweenPoints(pointRange[0], pointRange[2]);
double d2 = vtkMath::Distance2BetweenPoints(pointRange[1], pointRange[3]);
if (d1 <= d2)
{
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment