Commit 73c7971b authored by Yohann Bearzi's avatar Yohann Bearzi

New cell intersection API in `vtkCell`

One can now test if a cell intersects another cell
with method `vtkCell::IntersectCell`.
This method works for linear cells only.
parent e999e0ac
......@@ -14,8 +14,86 @@
=========================================================================*/
#include "vtkCell.h"
#include "vtkDataArrayRange.h"
#include "vtkMath.h"
#include "vtkNew.h"
#include "vtkPoints.h"
#include "vtkPolygon.h"
#include "vtkSmartPointer.h"
#include <vector>
namespace detail
{
//----------------------------------------------------------------------------
// Strategy:
//
// We throw all edges from one cell to an other and look if they intersect.
// In the case of a cell of one point, we just check if it lies inside
// the other cell.
int IntersectWithCell(vtkCell* self, vtkCell* other, double tol)
{
// Strategy:
// We throw edges from on cell to the other to look for intersections, and vice versa.
// If we catch one intersection, bingo.
//
// But first... we handle the case of a one point cell.
if (!other->GetNumberOfPoints() || !self->GetNumberOfPoints())
{
return 0;
}
double x[3], pcoords[3];
if (other->GetNumberOfPoints() == 1)
{
double closestPoint[3];
double* point = other->GetPoints()->GetPoint(0);
int subId;
double dist2, *weights = new double[self->GetNumberOfPoints()];
self->EvaluatePosition(point, closestPoint, subId, pcoords, dist2, weights);
delete[] weights;
return dist2 <= tol * tol;
}
if (self->GetNumberOfPoints() == 1)
{
double closestPoint[3];
double* point = self->GetPoints()->GetPoint(0);
int subId;
double dist2, *weights = new double[other->GetNumberOfPoints()];
other->EvaluatePosition(point, closestPoint, subId, pcoords, dist2, weights);
delete[] weights;
return dist2 <= tol * tol;
}
double p1[3], p2[3];
for (vtkIdType edgeId = 0; edgeId < self->GetNumberOfEdges(); ++edgeId)
{
double t;
int subId;
vtkCell* edge = self->GetEdge(edgeId);
vtkPoints* ends = edge->GetPoints();
ends->GetPoint(0, p1);
ends->GetPoint(1, p2);
if (other->IntersectWithLine(p1, p2, tol, t, x, pcoords, subId))
{
return 1;
}
}
for (vtkIdType edgeId = 0; edgeId < other->GetNumberOfEdges(); ++edgeId)
{
double t;
int subId;
vtkCell* edge = other->GetEdge(edgeId);
vtkPoints* ends = edge->GetPoints();
ends->GetPoint(0, p1);
ends->GetPoint(1, p2);
if (self->IntersectWithLine(p1, p2, tol, t, x, pcoords, subId))
{
std::cout << "second " << std::endl;
return 1;
}
}
return 0;
}
}
//------------------------------------------------------------------------------
// Construct cell.
......@@ -89,6 +167,104 @@ void vtkCell::DeepCopy(vtkCell* c)
}
//------------------------------------------------------------------------------
void vtkCell::Inflate(double dist)
{
if (this->GetNumberOfFaces() != 0)
{
vtkWarningMacro(<< "Base version of vtkCell::Inflate only implements cell inflation"
<< " for linear non 3D cells. Class " << this->GetClassName()
<< " needs to overload this method. Ignoring this cell.");
return;
}
std::vector<double> buf(3 * this->Points->GetNumberOfPoints(), 0.0);
vtkIdType pointId = 0;
auto pointRange = vtk::DataArrayTupleRange(this->Points->GetData());
using TupleRef = typename decltype(pointRange)::TupleReferenceType;
auto stlPoints = vtk::DataArrayTupleRange(this->Points->GetData());
auto postPointIt = stlPoints.begin();
++postPointIt;
auto prePointIt = stlPoints.end();
--prePointIt;
double v[3];
for (auto pointIt = stlPoints.begin(); pointIt != stlPoints.end();
++pointIt, ++postPointIt, ++prePointIt, ++pointId)
{
if (postPointIt == stlPoints.end())
{
postPointIt = stlPoints.begin();
}
if (prePointIt == stlPoints.end())
{
prePointIt = stlPoints.begin();
}
double* vIt = v;
for (auto compIt = pointIt->begin(), postCompIt = postPointIt->begin();
compIt != pointIt->end(); ++compIt, ++postCompIt, ++vIt)
{
*vIt = *compIt - *postCompIt;
}
vtkMath::Normalize(v);
buf[pointId * 3] += v[0] * dist;
buf[pointId * 3 + 1] += v[1] * dist;
buf[pointId * 3 + 2] += v[2] * dist;
vIt = v;
for (auto compIt = pointIt->begin(), preCompIt = prePointIt->begin(); compIt != pointIt->end();
++compIt, ++preCompIt, ++vIt)
{
*vIt = *compIt - *preCompIt;
}
vtkMath::Normalize(v);
buf[pointId * 3] += v[0] * dist;
buf[pointId * 3 + 1] += v[1] * dist;
buf[pointId * 3 + 2] += v[2] * dist;
}
auto it = buf.begin();
using TupleRef = typename decltype(stlPoints)::TupleReferenceType;
using CompRef = typename decltype(stlPoints)::ComponentReferenceType;
for (TupleRef point : stlPoints)
{
for (CompRef comp : point)
{
comp += *(it++);
}
}
}
//----------------------------------------------------------------------------
int vtkCell::IntersectWithCell(vtkCell* other, const vtkBoundingBox& boundingBox,
const vtkBoundingBox& otherBoundingBox, double tol)
{
if (!boundingBox.Intersects(otherBoundingBox))
{
return 0;
}
/**
* Given the strategy of detail::IntersectWithCell,
* the intersection detection is likely to be speeded up
* if exchanging other given this condition.
* The implementation first throws edges from first cell
* to look if it intersects with second cell, then it checks
* the other way.
* Since when one intersection is found, algorithm stops,
* we'd rather check embedded bounding box's cell's edges first.
*/
if (otherBoundingBox.IsSubsetOf(boundingBox))
{
return detail::IntersectWithCell(other, this, tol);
}
return detail::IntersectWithCell(this, other, tol);
}
//----------------------------------------------------------------------------
int vtkCell::IntersectWithCell(vtkCell* other, double tol)
{
return this->IntersectWithCell(
other, vtkBoundingBox(this->GetBounds()), vtkBoundingBox(other->GetBounds()), tol);
}
//----------------------------------------------------------------------------
// Compute cell bounding box (xmin,xmax,ymin,ymax,zmin,zmax). Return pointer
// to array of six double values.
double* vtkCell::GetBounds()
......
......@@ -43,8 +43,9 @@
#include "vtkCommonDataModelModule.h" // For export macro
#include "vtkObject.h"
#include "vtkCellType.h" // Needed to define cell types
#include "vtkIdList.h" // Needed for inline methods
#include "vtkBoundingBox.h" // Needed for IntersectWithCell
#include "vtkCellType.h" // Needed to define cell types
#include "vtkIdList.h" // Needed for inline methods
class vtkCellArray;
class vtkCellData;
......@@ -239,6 +240,13 @@ public:
vtkCellArray* connectivity, vtkPointData* inPd, vtkPointData* outPd, vtkCellData* inCd,
vtkIdType cellId, vtkCellData* outCd, int insideOut) = 0;
/**
* Inflates the cell by dist. Each point is displaced by dist in the direction of every halfedge.
* Input can be negative. If so, the cell will shrink. Artifacts will appear if one shrinks
* with -dist larger than any edge length, probably causing the cell to self-intersect.
*/
virtual void Inflate(double dist);
/**
* Intersect with a ray. Return parametric coordinates (both line and cell)
* and global intersection coordinates, given ray definition p1[3], p2[3] and tolerance tol.
......@@ -250,6 +258,19 @@ public:
virtual int IntersectWithLine(const double p1[3], const double p2[3], double tol, double& t,
double x[3], double pcoords[3], int& subId) = 0;
//@{
/**
* Intersects with an other cell. Returns 1 if cells intersect, 0 otherwise.
* If an exact intersection detection with regards to floating point precision is wanted,
* tol should be disregarded.
* `vtkBoundingBox` are optional parameters which slightly improve the speed of the computation
* if bounding boxes are already available to user.
*/
virtual int IntersectWithCell(vtkCell* other, double tol = 0.0);
virtual int IntersectWithCell(vtkCell* other, const vtkBoundingBox& boudingBox,
const vtkBoundingBox& otherBoundingBox, double tol = 0.0);
//@}
/**
* Generate simplices of proper dimension. If cell is 3D, tetrahedron are
* generated; if 2D triangles; if 1D lines; if 0D points. The form of the
......
......@@ -16,6 +16,7 @@
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkDataArrayRange.h"
#include "vtkDoubleArray.h"
#include "vtkMath.h"
#include "vtkOrderedTriangulator.h"
......@@ -25,6 +26,8 @@
#include "vtkPolygon.h"
#include "vtkTetra.h"
#include <vector>
vtkCell3D::vtkCell3D()
{
this->Triangulator = nullptr;
......@@ -74,6 +77,48 @@ bool vtkCell3D::IsInsideOut()
return signedDistanceToCentroid > 0.0;
}
//----------------------------------------------------------------------------
void vtkCell3D::Inflate(double dist)
{
// The strategy is the following:
// - For each point, get points from its one-ring
// - Displace each of theses points following its halfedge direction,
// by a distance of dist.
std::vector<double> buf(3 * this->Points->GetNumberOfPoints(), 0.0);
double sign = this->IsInsideOut() ? -1.0 : 1.0;
double p[3];
vtkIdType pointId = 0;
auto pointRange = vtk::DataArrayTupleRange(this->Points->GetData());
using ConstTupleRef = typename decltype(pointRange)::ConstTupleReferenceType;
for (ConstTupleRef point : pointRange)
{
const vtkIdType* ringIds;
vtkIdType ringSize = this->GetPointToOneRingPoints(pointId, ringIds);
for (vtkIdType id = 0; id < ringSize; ++id)
{
auto it = point->begin();
this->Points->GetPoint(ringIds[id], p);
double v[3] = { p[0] - *(it++), p[1] - *(it++), p[2] - *(it++) };
vtkMath::Normalize(v);
buf[3 * pointId] -= sign * v[0] * dist;
buf[3 * pointId + 1] -= sign * v[1] * dist;
buf[3 * pointId + 2] -= sign * v[2] * dist;
}
++pointId;
}
auto it = buf.begin();
using TupleRef = typename decltype(pointRange)::TupleReferenceType;
using CompRef = typename decltype(pointRange)::ComponentReferenceType;
for (TupleRef point : pointRange)
{
for (CompRef comp : point)
{
comp += *(it++);
}
}
}
//----------------------------------------------------------------------------
void vtkCell3D::Contour(double value, vtkDataArray* cellScalars,
vtkIncrementalPointLocator* locator, vtkCellArray* verts, vtkCellArray* lines,
vtkCellArray* polys, vtkPointData* inPd, vtkPointData* outPd, vtkCellData* inCd, vtkIdType cellId,
......
......@@ -177,6 +177,13 @@ public:
*/
int GetCellDimension() override { return 3; }
/**
* Inflates the cell by dist. Each point is displaced by dist in the direction of every halfedge.
* Input can be negative. If so, the cell will shrink. Artifacts will appear if one shrinks
* with -dist larger than any edge length, probably causing the cell to self-intersect.
*/
void Inflate(double dist) override;
//@{
/**
* Set the tolerance for merging clip intersection points that are near
......
......@@ -16,6 +16,7 @@
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkDataArrayRange.h"
#include "vtkIncrementalPointLocator.h"
#include "vtkLine.h"
#include "vtkMarchingSquaresLineCases.h"
......@@ -27,6 +28,9 @@
#include "vtkQuad.h"
#include "vtkTriangle.h"
#include <algorithm>
#include <array>
vtkStandardNewMacro(vtkPixel);
//------------------------------------------------------------------------------
......@@ -158,6 +162,82 @@ void vtkPixel::EvaluateLocation(int& subId, const double pcoords[3], double x[3]
}
//------------------------------------------------------------------------------
int vtkPixel::ComputeNormal(double n[3])
{
double p0[3], p1[3], p2[3];
this->Points->GetPoint(0, p0);
this->Points->GetPoint(1, p1);
this->Points->GetPoint(2, p2);
p1[0] -= p0[0];
p1[1] -= p0[1];
p1[2] -= p0[2];
p2[0] -= p0[0];
p2[1] -= p0[1];
p2[2] -= p0[2];
vtkMath::Cross(p1, p2, n);
if (std::abs(n[0]) < VTK_DBL_EPSILON && std::abs(n[1]) < VTK_DBL_EPSILON &&
std::abs(n[2]) < VTK_DBL_EPSILON)
{
return -1;
}
vtkMath::Normalize(n);
return (std::abs(n[1]) > 0.5) + (std::abs(n[2]) > 0.5) * 2;
}
//----------------------------------------------------------------------------
void vtkPixel::Inflate(double dist)
{
auto range = vtk::DataArrayTupleRange(this->Points->GetData());
using TupleRef = typename decltype(range)::TupleReferenceType;
double n[3];
int normalDirection = this->ComputeNormal(n);
int degeneratePixelDirection = -1;
if (normalDirection == -1)
{
auto it = range.cbegin(), jt = range.cbegin();
auto compIt = it->cbegin(), compJt = jt->cbegin();
++jt;
std::array<double, 3> diff{ *(compJt++) - *(compIt++), *(compJt++) - *(compIt++),
*(compJt++) - *(compIt++) };
auto maxIt = std::max_element(diff.cbegin(), diff.cend(),
[](double a, double b) -> bool { return std::abs(a) < std::abs(b); });
degeneratePixelDirection = std::distance(diff.cbegin(), maxIt);
}
int index = 0;
for (TupleRef point : range)
{
auto it = point.begin();
switch (normalDirection)
{
case 0:
*(++it++) += dist * (index % 2 ? 1.0 : -1.0);
*(it++) += dist * (index / 2 ? 1.0 : -1.0);
break;
case 1:
*(it++) += dist * (index % 2 ? 1.0 : -1.0);
*(++it++) += dist * (index / 2 ? 1.0 : -1.0);
break;
case 2:
*(it++) += dist * (index % 2 ? 1.0 : -1.0);
*(it++) += dist * (index / 2 ? 1.0 : -1.0);
break;
case -1:
if (degeneratePixelDirection > 0)
{
++it;
if (degeneratePixelDirection > 1)
{
++it;
}
}
*it += dist * (index % 2 ? 1.0 : -1.0);
break;
}
++index;
}
}
//----------------------------------------------------------------------------
int vtkPixel::CellBoundary(int vtkNotUsed(subId), const double pcoords[3], vtkIdList* pts)
{
double t1 = pcoords[0] - pcoords[1];
......
......@@ -58,6 +58,7 @@ public:
int EvaluatePosition(const double x[3], double closestPoint[3], int& subId, double pcoords[3],
double& dist2, double weights[]) override;
void EvaluateLocation(int& subId, const double pcoords[3], double x[3], double* weights) override;
void Inflate(double dist) override;
//@}
/**
......@@ -89,6 +90,13 @@ public:
}
//@}
/**
* vtkPixel's normal cannot be computed using vtkPolygon::ComputeNormal because
* its points are not sorted such that circulating on them forms the pixel.
* This is a convenient method so one can compute normals on a pixel.
*/
int ComputeNormal(double n[3]);
protected:
vtkPixel();
~vtkPixel() override;
......
......@@ -17,6 +17,7 @@
#include "vtkBox.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkDataArrayRange.h"
#include "vtkIncrementalPointLocator.h"
#include "vtkLine.h"
#include "vtkMarchingCubesTriangleCases.h"
......@@ -188,6 +189,22 @@ void vtkVoxel::EvaluateLocation(
}
//------------------------------------------------------------------------------
void vtkVoxel::Inflate(double dist)
{
int index = 0;
auto range = vtk::DataArrayTupleRange(this->Points->GetData());
using TupleRef = typename decltype(range)::TupleReferenceType;
for (TupleRef point : range)
{
auto it = point.begin();
*(it++) += dist * (index % 2 ? 1.0 : -1.0);
*(it++) += dist * ((index / 2) % 2 ? 1.0 : -1.0);
*(it++) += dist * (index / 4 ? 1.0 : -1.0);
++index;
}
}
//----------------------------------------------------------------------------
//
// Compute Interpolation functions
//
......@@ -420,7 +437,7 @@ static constexpr vtkIdType pointToOneRingPoints[vtkVoxel::NumberOfPoints]
};
}
//------------------------------------------------------------------------------
//----------------------------------------------------------------------------
//
// Marching cubes case table
//
......
......@@ -114,6 +114,7 @@ public:
int Triangulate(int index, vtkIdList* ptIds, vtkPoints* pts) override;
void Derivatives(
int subId, const double pcoords[3], const double* values, int dim, double* derivs) override;
void Inflate(double dist) override;
//@}
static void InterpolationDerivs(const double pcoords[3], double derivs[24]);
......
/*=========================================================================
Program: Visualization Toolkit
Module: TestOverlappingCellsDetector.cxx
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
#if VTK_MODULE_ENABLE_VTK_ParallelMPI
#include "vtkMPIController.h"
#else
#include "vtkDummyController.h"
#endif
#include "vtkCellData.h"
#include "vtkDataArray.h"
#include "vtkDataArrayRange.h"
#include "vtkDataSet.h"
#include "vtkGenerateGlobalIds.h"
#include "vtkMultiProcessController.h"
#include "vtkNew.h"
#include "vtkOverlappingCellsDetector.h"
#include "vtkRedistributeDataSetFilter.h"
#include "vtkTestUtilities.h"
#include "vtkUnstructuredGrid.h"
#include "vtkXMLUnstructuredGridReader.h"
namespace
{
static constexpr vtkIdType Collisions[72] = { 6, 0, 6, 0, 4, 4, 6, 0, 10, 7, 4, 0, 0, 7, 9, 0, 0, 5,
5, 0, 0, 0, 5, 9, 0, 6, 4, 6, 4, 4, 0, 6, 1, 0, 6, 9, 8, 7, 1, 9, 5, 0, 0, 5, 9, 7, 2, 3, 0, 0, 0,
0, 0, 6, 1, 4, 0, 1, 0, 0, 4, 0, 0, 0, 0, 0, 3, 6, 0, 0, 0, 0 };
}
int TestOverlappingCellsDetector(int argc, char* argv[])
{
#if VTK_MODULE_ENABLE_VTK_ParallelMPI
vtkNew<vtkMPIController> contr;
#else
vtkNew<vtkDummyController> contr;
#endif
contr->Initialize(&argc, &argv);
vtkMultiProcessController::SetGlobalController(contr);
int retVal = EXIT_SUCCESS;
int myrank = contr->GetLocalProcessId();
const char* name =
vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/overlapping_tetras.vtu");
vtkNew<vtkGenerateGlobalIds> globalIds;
if (myrank == 0)
{
vtkNew<vtkXMLUnstructuredGridReader> reader;
reader->SetFileName(name);
globalIds->SetInputConnection(reader->GetOutputPort());
globalIds->Update();
}
else
{
vtkNew<vtkUnstructuredGrid> ug;
ug->Initialize();
globalIds->SetInputDataObject(ug);
}
vtkNew<vtkRedistributeDataSetFilter> redistribute;
redistribute->SetInputConnection(globalIds->GetOutputPort());
std::cout << myrank << "coucou " << std::endl;
vtkNew<vtkOverlappingCellsDetector> detector;
detector->SetInputConnection(redistribute->GetOutputPort());
detector->Update();
vtkDataSet* output = vtkDataSet::SafeDownCast(detector->GetOutput(0));
vtkDataArray* data =
output->GetCellData()->GetArray(vtkOverlappingCellsDetector::NumberOfCollisionsPerCell);
vtkDataArray* ids = output->GetCellData()->GetArray("GlobalCellIds");
auto valIt = vtk::DataArrayValueRange<3>(data).cbegin();
auto idIt = vtk::DataArrayValueRange<3>(ids).cbegin();
for (; valIt != vtk::DataArrayValueRange<3>(data).cend(); ++valIt, ++idIt)
{
if (Collisions[static_cast<vtkIdType>(*idIt)] != *valIt)
{
retVal = EXIT_FAILURE;
break;
}
}
vtkMultiProcessController::SetGlobalController(nullptr);
contr->Finalize();
return retVal;
}
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