From 462f800503df93910ce11e55eac17a4539c6e6a0 Mon Sep 17 00:00:00 2001
From: Mathieu Westphal <mathieu.westphal@kitware.com>
Date: Wed, 19 Feb 2025 17:14:25 +0100
Subject: [PATCH] vtkProbeFilter: Add support for probing oriented images
---
Filters/Core/vtkProbeFilter.cxx | 146 +++++++++++++++-----------------
Filters/Core/vtkProbeFilter.h | 3 +-
2 files changed, 68 insertions(+), 81 deletions(-)
diff --git a/Filters/Core/vtkProbeFilter.cxx b/Filters/Core/vtkProbeFilter.cxx
index 2ca2135ff9f..3249c931fb6 100644
--- a/Filters/Core/vtkProbeFilter.cxx
+++ b/Filters/Core/vtkProbeFilter.cxx
@@ -755,53 +755,72 @@ void vtkProbeFilter::ProbeEmptyPoints(
}
//------------------------------------------------------------------------------
-static void GetPointIdsInRange(double rangeMin, double rangeMax, double start, double stepsize,
- int numSteps, int& minid, int& maxid)
+void vtkProbeFilter::ProbeImagePointsInCell(vtkGenericCell* cell, vtkIdType cellId,
+ vtkDataSet* source, int srcBlockId, vtkImageData* input, vtkPointData* outPD, char* maskArray,
+ double* wtsBuff)
{
- if (stepsize == 0)
- {
- minid = maxid = 0;
- return;
- }
+ vtkPointData* pd = source->GetPointData();
+ // 1. get coordinates of sampling grids
- minid = vtkMath::Ceil((rangeMin - start) / stepsize);
- if (minid < 0)
+ // Recover cell bounds
+ double cellBounds[6];
+ source->GetCellBounds(cellId, cellBounds);
+ vtkBoundingBox cellBB(cellBounds);
+
+ // Recover input bounds, we already know they intersect, but reduce cellbounds
+ // to input bounds to ensure ComputeStructuredCoordinates works as expected on
+ // the edges of the input bounds
+ double inBounds[6];
+ input->GetBounds(inBounds);
+ cellBB.IntersectBox(inBounds);
+
+ // If ComputeTolerance is set, compute a tolerance proportional to the
+ // cell length. Otherwise, use the user specified absolute tolerance.
+ double tol2;
+ if (this->ComputeTolerance)
{
- minid = 0;
+ const vtkBoundingBox bbox(cellBounds);
+ tol2 = CELL_TOLERANCE_FACTOR_SQR * bbox.GetDiagonalLength2();
}
-
- maxid = vtkMath::Floor((rangeMax - start) / stepsize);
- if (maxid > numSteps - 1)
+ else
{
- maxid = numSteps - 1;
+ tol2 = this->Tolerance * this->Tolerance;
}
-}
-//------------------------------------------------------------------------------
-void vtkProbeFilter::ProbeImagePointsInCell(vtkGenericCell* cell, vtkIdType cellId,
- vtkDataSet* source, int srcBlockId, const double start[3], const double spacing[3],
- const int dim[3], vtkPointData* outPD, char* maskArray, double* wtsBuff)
-{
- vtkPointData* pd = source->GetPointData();
-
- // get coordinates of sampling grids
- double cellBounds[6];
- source->GetCellBounds(cellId, cellBounds);
+ // Iterate on each point of the bounding box
+ // to identify the covered grid
+ int* inputExtent = input->GetExtent();
+ int idxBounds[6] = { std::numeric_limits<int>::max(), std::numeric_limits<int>::min(),
+ std::numeric_limits<int>::max(), std::numeric_limits<int>::min(),
+ std::numeric_limits<int>::max(), std::numeric_limits<int>::min() };
+ double corner[3];
+ int ijk[3];
+ double pCoords[3];
+ bool found = false;
+ for (int i = 0; i < 8; i++)
+ {
+ cellBB.GetCorner(i, corner);
+ if (input->ComputeStructuredCoordinates(corner, ijk, pCoords, tol2))
+ {
+ found = true;
+ for (int j = 0; j < 3; j++)
+ {
+ idxBounds[2 * j] = std::min(ijk[j], idxBounds[2 * j]);
- int idxBounds[6];
- GetPointIdsInRange(
- cellBounds[0], cellBounds[1], start[0], spacing[0], dim[0], idxBounds[0], idxBounds[1]);
- GetPointIdsInRange(
- cellBounds[2], cellBounds[3], start[1], spacing[1], dim[1], idxBounds[2], idxBounds[3]);
- GetPointIdsInRange(
- cellBounds[4], cellBounds[5], start[2], spacing[2], dim[2], idxBounds[4], idxBounds[5]);
+ // Force a +1 here to ensure we get even the edge of the input
+ // In some case that can be outside of the input image
+ int externalExtent = std::min(ijk[j] + 1, inputExtent[2 * j + 1]);
+ idxBounds[2 * j + 1] = std::max(externalExtent, idxBounds[2 * j + 1]);
+ }
+ }
+ }
- if ((idxBounds[1] - idxBounds[0]) < 0 || (idxBounds[3] - idxBounds[2]) < 0 ||
- (idxBounds[5] - idxBounds[4]) < 0)
+ if (!found)
{
return;
}
+ // 2. Recover the cell to process
source->GetCell(cellId, cell);
double cpbuf[3];
@@ -814,35 +833,23 @@ void vtkProbeFilter::ProbeImagePointsInCell(vtkGenericCell* cell, vtkIdType cell
closestPoint = nullptr;
}
- // If ComputeTolerance is set, compute a tolerance proportional to the
- // cell length. Otherwise, use the user specified absolute tolerance.
- double tol2;
- if (this->ComputeTolerance)
+ // 3. Check each indices from the grid
+ for (ijk[2] = idxBounds[4]; ijk[2] <= idxBounds[5]; ijk[2]++)
{
- const vtkBoundingBox bbox(cellBounds);
- tol2 = CELL_TOLERANCE_FACTOR_SQR * bbox.GetDiagonalLength2();
- }
- else
- {
- tol2 = this->Tolerance * this->Tolerance;
- }
- for (vtkIdType iz = idxBounds[4]; iz <= idxBounds[5]; iz++)
- {
- double p[3];
- p[2] = start[2] + iz * spacing[2];
- for (vtkIdType iy = idxBounds[2]; iy <= idxBounds[3]; iy++)
+ for (ijk[1] = idxBounds[2]; ijk[1] <= idxBounds[3]; ijk[1]++)
{
- p[1] = start[1] + iy * spacing[1];
- for (vtkIdType ix = idxBounds[0]; ix <= idxBounds[1]; ix++)
+ for (ijk[0] = idxBounds[0]; ijk[0] <= idxBounds[1]; ijk[0]++)
{
// skip processed points
- const vtkIdType ptId = ix + dim[0] * (iy + dim[1] * iz);
+ const vtkIdType ptId = input->ComputePointId(ijk);
if (maskArray[ptId] == 1)
{
continue;
}
- // For each grid point within the cell bound, interpolate values
- p[0] = start[0] + ix * spacing[0];
+
+ // record point coordinates
+ double p[3];
+ input->GetPoint(ptId, p);
double pcoords[3];
int subId;
@@ -876,14 +883,11 @@ class vtkProbeFilter::ProbeImageDataWorklet
{
public:
ProbeImageDataWorklet(vtkProbeFilter* probeFilter, vtkDataSet* source, int srcBlockId,
- const double start[3], const double spacing[3], const int dim[3], vtkPointData* outPD,
- char* maskArray, int maxCellSize)
+ vtkImageData* input, vtkPointData* outPD, char* maskArray, int maxCellSize)
: ProbeFilter(probeFilter)
, Source(source)
+ , Input(input)
, SrcBlockId(srcBlockId)
- , Start(start)
- , Spacing(spacing)
- , Dim(dim)
, OutPointData(outPD)
, MaskArray(maskArray)
, MaxCellSize(maxCellSize)
@@ -924,7 +928,7 @@ public:
}
this->ProbeFilter->ProbeImagePointsInCell(cell, cellId, this->Source, this->SrcBlockId,
- this->Start, this->Spacing, this->Dim, this->OutPointData, this->MaskArray, weights);
+ this->Input, this->OutPointData, this->MaskArray, weights);
}
}
@@ -933,10 +937,8 @@ public:
private:
vtkProbeFilter* ProbeFilter;
vtkDataSet* Source;
+ vtkImageData* Input;
int SrcBlockId;
- const double* Start;
- const double* Spacing;
- const int* Dim;
vtkPointData* OutPointData;
char* MaskArray;
int MaxCellSize;
@@ -952,25 +954,11 @@ void vtkProbeFilter::ProbePointsImageData(
vtkPointData* outPD = output->GetPointData();
char* maskArray = this->MaskPoints->GetPointer(0);
- //----------------------------------------
- double spacing[3];
- input->GetSpacing(spacing);
- int extent[6];
- input->GetExtent(extent);
- int dim[3];
- input->GetDimensions(dim);
- double start[3];
- input->GetOrigin(start);
- start[0] += static_cast<double>(extent[0]) * spacing[0];
- start[1] += static_cast<double>(extent[2]) * spacing[1];
- start[2] += static_cast<double>(extent[4]) * spacing[2];
-
vtkIdType numSrcCells = source->GetNumberOfCells();
-
if (numSrcCells > 0)
{
ProbeImageDataWorklet worklet(
- this, source, srcIdx, start, spacing, dim, outPD, maskArray, source->GetMaxCellSize());
+ this, source, srcIdx, input, outPD, maskArray, source->GetMaxCellSize());
vtkSMPTools::For(0, numSrcCells, worklet);
}
diff --git a/Filters/Core/vtkProbeFilter.h b/Filters/Core/vtkProbeFilter.h
index 14ea5550889..ba11ab73ee8 100644
--- a/Filters/Core/vtkProbeFilter.h
+++ b/Filters/Core/vtkProbeFilter.h
@@ -300,8 +300,7 @@ private:
void ProbePointsImageData(
vtkImageData* input, int srcIdx, vtkDataSet* source, vtkImageData* output);
void ProbeImagePointsInCell(vtkGenericCell* cell, vtkIdType cellId, vtkDataSet* source,
- int srcBlockId, const double start[3], const double spacing[3], const int dim[3],
- vtkPointData* outPD, char* maskArray, double* wtsBuff);
+ int srcBlockId, vtkImageData* input, vtkPointData* outPD, char* maskArray, double* wtsBuff);
class ProbeImageDataWorklet;
--
GitLab