Commit 71bffccf authored by Yohann Bearzi's avatar Yohann Bearzi Committed by Kitware Robot

Merge topic 'cell-bounding-sphere'

0faa45e8 Optimization of overlapping cell detector
b90f56fa New bounding sphere computation for vtkCell
77a1f63f now vtkMath::Distance2BeetweenPoints range aware
Acked-by: Kitware Robot's avatarKitware Robot <kwrobot@kitware.com>
Acked-by: Utkarsh Ayachit's avatarUtkarsh Ayachit <utkarsh.ayachit@kitware.com>
Merge-request: !6909
parents 9531f58b 0faa45e8
......@@ -521,6 +521,27 @@ public:
static bool ProjectVector2D(const double a[2], const double b[2], double projection[2]);
//@}
/**
* Compute distance squared between two points p1 and p2.
* This version allows for custom range and iterator types to be used. These
* types must implement `operator[]`, and at least one of them must have
* a `value_type` typedef.
*
* The first template parameter `ReturnTypeT` sets the return type of this method.
* By default, it is set to `double`, but it can be overridden.
*
* The `EnableT` template parameter is used to make sure that this version
* doesn't capture the `float*`, `float[]`, `double*`, `double[]` as
* those should go to the other `Distance2BetweenPoints` functions.
*
* @warning This method assumes that both parameters have 3 components.
*/
template <typename ReturnTypeT = double, typename TupleRangeT1, typename TupleRangeT2,
typename EnableT = typename std::conditional<!std::is_pointer<TupleRangeT1>::value &&
!std::is_array<TupleRangeT1>::value,
TupleRangeT1, TupleRangeT2>::type::value_type>
static ReturnTypeT Distance2BetweenPoints(const TupleRangeT1& p1, const TupleRangeT2& p2);
/**
* Compute distance squared between two points p1 and p2.
* (float version).
......@@ -1436,6 +1457,14 @@ inline double vtkMath::Distance2BetweenPoints(const double p1[3], const double p
(p1[2] - p2[2]) * (p1[2] - p2[2]));
}
//----------------------------------------------------------------------------
template <typename ReturnTypeT, typename TupleRangeT1, typename TupleRangeT2, typename EnableT>
inline ReturnTypeT vtkMath::Distance2BetweenPoints(const TupleRangeT1& p1, const TupleRangeT2& p2)
{
return ((p1[0] - p2[0]) * (p1[0] - p2[0]) + (p1[1] - p2[1]) * (p1[1] - p2[1]) +
(p1[2] - p2[2]) * (p1[2] - p2[2]));
}
//----------------------------------------------------------------------------
// Cross product of two 3-vectors. Result (a x b) is stored in c[3].
inline void vtkMath::Cross(const float a[3], const float b[3], float c[3])
......
......@@ -20,6 +20,8 @@
#include "vtkPoints.h"
#include "vtkPolygon.h"
#include "vtkSmartPointer.h"
#include "vtkTetra.h"
#include "vtkTriangle.h"
#include <vector>
......@@ -160,6 +162,117 @@ void vtkCell::DeepCopy(vtkCell* c)
this->PointIds->DeepCopy(c->PointIds);
}
//------------------------------------------------------------------------------
double vtkCell::ComputeBoundingSphere(double center[3]) const
{
// We do easy cases first for number of points <= 4
switch (this->Points->GetNumberOfPoints())
{
case 0:
center[0] = std::numeric_limits<double>::quiet_NaN();
center[1] = std::numeric_limits<double>::quiet_NaN();
center[2] = std::numeric_limits<double>::quiet_NaN();
return std::numeric_limits<double>::quiet_NaN();
case 1:
this->Points->GetPoint(0, center);
return 0.0;
case 2:
{
auto points = vtk::DataArrayTupleRange(this->Points->GetData());
auto p0 = points[0], p1 = points[1];
center[0] = 0.5 * (p0[0] + p1[0]);
center[1] = 0.5 * (p0[1] + p1[1]);
center[2] = 0.5 * (p0[2] + p1[2]);
return vtkMath::Distance2BetweenPoints(center, p0);
}
case 3:
{
vtkTriangle::ComputeCentroid(this->Points, nullptr, center);
auto points = vtk::DataArrayTupleRange(this->Points->GetData());
return vtkMath::Distance2BetweenPoints(center, points[0]);
}
case 4:
{
vtkTetra::ComputeCentroid(this->Points, nullptr, center);
auto points = vtk::DataArrayTupleRange(this->Points->GetData());
return vtkMath::Distance2BetweenPoints(center, points[0]);
}
default:
break;
}
// For more complex cells, we follow Ritter's bounding sphere algorithm
// 1. Pick a point x (first point in our case) in the cell, and look for
// a point y the furthest from x.
// 2. Look for a point z the furthest from y.
// 3. Create a sphere centered at [z,y] with appropriate radius
// 4. Until all points are not in the sphere, take a point outside the sphere,
// and update the sphere to include former sphere + this point
auto points = vtk::DataArrayTupleRange(this->Points->GetData());
using ConstRefType = typename decltype(points)::ConstTupleReferenceType;
ConstRefType x = points[0];
vtkIdType yid = 1, zid = 0;
double dist2 = 0.0;
for (vtkIdType id = 1; id < points.size(); ++id)
{
double tmpdist2 = vtkMath::Distance2BetweenPoints(points[id], x);
if (tmpdist2 > dist2)
{
dist2 = tmpdist2;
yid = id;
}
}
ConstRefType y = points[yid];
dist2 = 0.0;
for (vtkIdType id = 0; id < points.size(); ++id)
{
double tmpdist2 = vtkMath::Distance2BetweenPoints(points[id], y);
if (tmpdist2 > dist2)
{
dist2 = tmpdist2;
zid = id;
}
}
ConstRefType z = points[zid];
center[0] = 0.5 * (y[0] + z[0]);
center[1] = 0.5 * (y[1] + z[1]);
center[2] = 0.5 * (y[2] + z[2]);
dist2 = vtkMath::Distance2BetweenPoints(y, center);
double v[3];
vtkIdType pointId;
do
{
for (pointId = 0; pointId < points.size(); ++pointId)
{
if (vtkMath::Distance2BetweenPoints(points[pointId], center) > dist2)
{
break;
}
}
if (pointId != points.size())
{
auto p = points[pointId];
v[0] = p[0] - center[0];
v[1] = p[1] - center[1];
v[2] = p[2] - center[2];
double d = 0.5 * (vtkMath::Norm(v) - std::sqrt(dist2));
vtkMath::Normalize(v);
center[0] += d * v[0];
center[1] += d * v[0];
center[2] += d * v[2];
dist2 = vtkMath::Distance2BetweenPoints(p, center);
}
} while (pointId != points.size());
return dist2;
}
//------------------------------------------------------------------------------
void vtkCell::Inflate(double dist)
{
......
......@@ -247,6 +247,16 @@ public:
*/
virtual void Inflate(double dist);
/**
* Computes the bounding sphere of the cell. If the number of points in the cell is lower
* or equal to 4, an exact bounding sphere is computed. If not, Ritter's algorithm
* is followed. If the input sphere has zero points, then each coordinate of
* center is set to NaN, as well as the returned distance.
*
* This method computes the center of the sphere, and returns its squared radius.
*/
virtual double ComputeBoundingSphere(double center[3]) const;
/**
* 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.
......
......@@ -237,6 +237,17 @@ void vtkPixel::Inflate(double dist)
}
}
//------------------------------------------------------------------------------
double vtkPixel::ComputeBoundingSphere(double center[3]) const
{
auto points = vtk::DataArrayTupleRange(this->Points->GetData());
auto p0 = points[0], p3 = points[3];
center[0] = 0.5 * (p0[0] + p3[0]);
center[1] = 0.5 * (p0[1] + p3[1]);
center[2] = 0.5 * (p0[2] + p3[2]);
return vtkMath::Distance2BetweenPoints(center, p0);
}
//----------------------------------------------------------------------------
int vtkPixel::CellBoundary(int vtkNotUsed(subId), const double pcoords[3], vtkIdList* pts)
{
......
......@@ -61,6 +61,11 @@ public:
void Inflate(double dist) override;
//@}
/**
* Computes exact bounding sphere of this pixel.
*/
double ComputeBoundingSphere(double center[3]) const override;
/**
* Return the center of the triangle in parametric coordinates.
*/
......
......@@ -105,6 +105,17 @@ bool vtkVoxel::IsInsideOut()
return (pt2[0] - pt1[0]) * (pt2[1] - pt1[1]) * (pt2[2] - pt1[2]) < 0.0;
}
//------------------------------------------------------------------------------
double vtkVoxel::ComputeBoundingSphere(double center[3]) const
{
auto points = vtk::DataArrayTupleRange(this->Points->GetData());
auto p0 = points[0], p7 = points[7];
center[0] = 0.5 * (p0[0] + p7[0]);
center[1] = 0.5 * (p0[1] + p7[1]);
center[2] = 0.5 * (p0[2] + p7[2]);
return vtkMath::Distance2BetweenPoints(center, p0);
}
//------------------------------------------------------------------------------
int vtkVoxel::EvaluatePosition(const double x[3], double closestPoint[3], int& subId,
double pcoords[3], double& dist2, double weights[])
......
......@@ -64,6 +64,11 @@ public:
bool IsInsideOut() override;
//@}
/**
* Computes exact bounding sphere of this voxel.
*/
double ComputeBoundingSphere(double center[3]) const override;
/**
* static constexpr handle on the number of points.
*/
......
......@@ -49,8 +49,6 @@
#include <unordered_map>
#include <vector>
#include "vtkXMLDataObjectWriter.h"
// clang-format off
#include "vtk_diy2.h"
#include VTK_DIY2(diy/mpi.hpp)
......@@ -103,12 +101,12 @@ vtkPointSet* ConvertCellsToBoundingSpheres(vtkDataSet* ds, std::vector<vtkBoundi
for (vtkIdType id = 0; id < size; ++id)
{
bboxes.emplace_back(ds->GetCell(id)->GetBounds());
vtkBoundingBox& bbox = bboxes[id];
vtkCell* cell = ds->GetCell(id);
bboxes.emplace_back(ds->GetBounds());
double center[3];
bbox.GetCenter(center);
double radius2 = cell->ComputeBoundingSphere(center);
SphereRadiusArray->SetValue(id, std::sqrt(radius2));
points->SetPoint(id, center);
SphereRadiusArray->SetValue(id, bbox.GetDiagonalLength() * 0.5);
}
pointCloud->GetPointData()->AddArray(SphereRadiusArray);
......@@ -440,7 +438,7 @@ int vtkOverlappingCellsDetector::ExposeOverlappingCellsAmongBlocks(vtkPointSet*
}
// We share overlapping candidates with neighbor blocks.
master.foreach ([&output, &overlappingCellCandidatesDataSets](
master.foreach ([&overlappingCellCandidatesDataSets](
void*, const diy::Master::ProxyWithLink& cp) {
// enqueue
for (int i = 0; i < static_cast<int>(cp.link()->size()); ++i)
......@@ -450,8 +448,7 @@ int vtkOverlappingCellsDetector::ExposeOverlappingCellsAmongBlocks(vtkPointSet*
}
});
master.exchange();
master.foreach ([&output, &overlappingCellCandidatesDataSets](
detail::Block* block, const diy::Master::ProxyWithLink& cp) {
master.foreach ([](detail::Block* block, const diy::Master::ProxyWithLink& cp) {
// dequeue
std::vector<int> incoming;
cp.incoming(incoming);
......@@ -527,7 +524,7 @@ int vtkOverlappingCellsDetector::ExposeOverlappingCellsAmongBlocks(vtkPointSet*
// We need to send back collision information to the original block, so they
// can add the collisions they couldn't detect.
master.foreach ([&output, &collisionListMapList](void*, const diy::Master::ProxyWithLink& cp) {
master.foreach ([&collisionListMapList](void*, const diy::Master::ProxyWithLink& cp) {
// enqueue
for (int i = 0; i < static_cast<int>(cp.link()->size()); ++i)
{
......@@ -536,7 +533,7 @@ int vtkOverlappingCellsDetector::ExposeOverlappingCellsAmongBlocks(vtkPointSet*
}
});
master.exchange();
master.foreach ([&output](detail::Block* block, const diy::Master::ProxyWithLink& cp) {
master.foreach ([](detail::Block* block, const diy::Master::ProxyWithLink& cp) {
// dequeue
std::vector<int> incoming;
cp.incoming(incoming);
......
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