diff --git a/Documentation/release/dev/add-oriented-image-support-probe-filter.md b/Documentation/release/dev/add-oriented-image-support-probe-filter.md
new file mode 100644
index 0000000000000000000000000000000000000000..740ee57053f1f616c49b5f1d0ff10768d9fbe92f
--- /dev/null
+++ b/Documentation/release/dev/add-oriented-image-support-probe-filter.md
@@ -0,0 +1,5 @@
+## Add support for oriented images in vtkProbeFilter
+
+vtkProbeFilter has been improved to support oriented images.
+The implementation now rely on the native support of orientation
+in the API of vtkImageData.
diff --git a/Filters/Core/Testing/Cxx/TestProbeFilter.cxx b/Filters/Core/Testing/Cxx/TestProbeFilter.cxx
index 0cfbc0560a6e1e39942d7b64d111745ec18d9d04..d1b72049755a28add9bb6d347e6aab4aedf36b3b 100644
--- a/Filters/Core/Testing/Cxx/TestProbeFilter.cxx
+++ b/Filters/Core/Testing/Cxx/TestProbeFilter.cxx
@@ -4,18 +4,20 @@
#include "vtkArrayCalculator.h"
#include "vtkDataArray.h"
#include "vtkDataSet.h"
+#include "vtkImageData.h"
#include "vtkLineSource.h"
#include "vtkNew.h"
#include "vtkPointData.h"
#include "vtkProbeFilter.h"
+#include "vtkRTAnalyticSource.h"
// Gets the number of points the probe filter counted as valid.
// The parameter should be the output of the probe filter
-int GetNumberOfValidPoints(vtkDataSet* pd)
+vtkIdType GetNumberOfValidPoints(vtkDataSet* pd)
{
vtkDataArray* data = pd->GetPointData()->GetScalars("vtkValidPointMask");
- int numValid = 0;
- for (int i = 0; i < data->GetNumberOfTuples(); ++i)
+ vtkIdType numValid = 0;
+ for (vtkIdType i = 0; i < data->GetNumberOfTuples(); ++i)
{
if (data->GetVariantValue(i).ToDouble() == 1)
{
@@ -25,8 +27,27 @@ int GetNumberOfValidPoints(vtkDataSet* pd)
return numValid;
}
+bool TestProbeFilterWithProvidedData(
+ vtkDataSet* input, vtkDataSet* source, vtkIdType expectedNValidPoints)
+{
+ vtkNew<vtkProbeFilter> probe;
+ probe->SetInputData(input);
+ probe->SetSourceData(source);
+ probe->Update();
+
+ vtkIdType nValidPoints = GetNumberOfValidPoints(probe->GetOutput());
+ if (nValidPoints != expectedNValidPoints)
+ {
+ std::cerr << "Unexpected number of valid points, got " << nValidPoints << " instead of "
+ << expectedNValidPoints << std::endl;
+ return false;
+ }
+
+ return true;
+}
+
// Tests the CompteThreshold and Threshold parameters on the vtkProbeFilter
-int TestProbeFilterThreshold()
+bool TestProbeFilterThreshold()
{
vtkNew<vtkLineSource> line1;
line1->SetPoint1(-1, 0, 0);
@@ -51,7 +72,7 @@ int TestProbeFilterThreshold()
int validDefault = GetNumberOfValidPoints(probe->GetOutput());
if (validDefault != 2)
{
- return 1;
+ return false;
}
// turn off computing tolerance and set it to 11 times what is was.
// 11 is magic number to get all the points within line1 selected.
@@ -63,19 +84,56 @@ int TestProbeFilterThreshold()
if (validNext != 11)
{
- return 1;
+ return false;
}
// threshold is still set high, but we tell it to ignore it
probe->SetComputeTolerance(true);
probe->Update();
int validIgnore = GetNumberOfValidPoints(probe->GetOutput());
- return (validIgnore == 2) ? 0 : 1;
+ return (validIgnore == 2) ? true : false;
+}
+
+// Test probing one image into another
+bool TestProbeFilterWithImages()
+{
+ // Create Pipeline
+ vtkNew<vtkRTAnalyticSource> wavelet;
+ wavelet->SetWholeExtent(0, 16, 0, 16, 0, 16);
+ wavelet->SetCenter(8, 8, 8);
+ wavelet->Update();
+
+ vtkNew<vtkImageData> img;
+ img->SetExtent(1, 15, 1, 15, 1, 15);
+ img->SetOrigin(1, 1, 1);
+
+ return TestProbeFilterWithProvidedData(img, wavelet->GetOutput(), 3375);
+}
+
+// Test probing one image into an oriented one
+bool TestProbeFilterWithOrientedImages()
+{
+ // Create Pipeline
+ vtkNew<vtkRTAnalyticSource> wavelet;
+ wavelet->SetWholeExtent(0, 16, 0, 16, 0, 16);
+ wavelet->SetCenter(8, 8, 8);
+ wavelet->Update();
+
+ vtkNew<vtkImageData> img;
+ img->SetExtent(1, 15, 1, 15, 1, 15);
+ img->SetOrigin(1, 1, 1);
+ img->SetDirectionMatrix(0.7, -0.7, 0, 0.7, 0.7, 0, 0, 0, 1);
+
+ return TestProbeFilterWithProvidedData(img, wavelet->GetOutput(), 1575);
}
// Currently only tests the ComputeThreshold and Threshold. Other tests
// should be added
int TestProbeFilter(int, char*[])
{
- return TestProbeFilterThreshold();
+ bool status = true;
+ status &= TestProbeFilterThreshold();
+ status &= TestProbeFilterWithImages();
+ status &= TestProbeFilterWithOrientedImages();
+ return status ? EXIT_SUCCESS : EXIT_FAILURE;
}
diff --git a/Filters/Core/Testing/Cxx/TestResampleToImage.cxx b/Filters/Core/Testing/Cxx/TestResampleToImage.cxx
index 106ef390f22e901daa3285001adaf0dbe44877b4..1016b9e27a115a0042b7ee5b9f9d326732f6cbba 100644
--- a/Filters/Core/Testing/Cxx/TestResampleToImage.cxx
+++ b/Filters/Core/Testing/Cxx/TestResampleToImage.cxx
@@ -87,6 +87,7 @@ int TestResampleToImage(int, char*[])
}
// Test for ParaView issue #19856
+ // Please note that #19856 is not fully fixed yet, but the test below still has value
vtkNew<vtkCellTypeSource> cellTypeSource;
cellTypeSource->SetCellOrder(1);
@@ -116,9 +117,9 @@ int TestResampleToImage(int, char*[])
{
zeros += !val;
}
- if (zeros != 744732)
+ if (zeros != 744725)
{
- std::cout << "Caught " << zeros << " invalid points, it should have been 744732" << std::endl;
+ std::cout << "Caught " << zeros << " invalid points, it should have been 744725" << std::endl;
status = 1;
}
diff --git a/Filters/Core/Testing/Cxx/TestResampleToImage2D.cxx b/Filters/Core/Testing/Cxx/TestResampleToImage2D.cxx
index 6209ff48977e00d4fe5335b39e2b63277e046b07..c18b0c2beb50faa1748d3f7996714cc6aade46de 100644
--- a/Filters/Core/Testing/Cxx/TestResampleToImage2D.cxx
+++ b/Filters/Core/Testing/Cxx/TestResampleToImage2D.cxx
@@ -31,6 +31,7 @@ int TestResampleToImage2D(int argc, char* argv[])
if (range[1] - range[0] < 0.01)
{
cerr << "Error resampling along X" << endl;
+ return EXIT_FAILURE;
}
// test on Y
@@ -41,6 +42,7 @@ int TestResampleToImage2D(int argc, char* argv[])
if (range[1] - range[0] < 0.01)
{
cerr << "Error resampling along Y" << endl;
+ return EXIT_FAILURE;
}
// test on Z
@@ -51,6 +53,7 @@ int TestResampleToImage2D(int argc, char* argv[])
if (range[1] - range[0] < 0.01)
{
cerr << "Error resampling along Z" << endl;
+ return EXIT_FAILURE;
}
return EXIT_SUCCESS;
diff --git a/Filters/Core/vtkProbeFilter.cxx b/Filters/Core/vtkProbeFilter.cxx
index 2ca2135ff9f8b82dc4fd3fb032e750e9efc8814c..3249c931fb6accde3c36451e3cc6ff7e29f6457a 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 14ea55508899ddbbc31561bc5d2005cabfa7e91b..ba11ab73ee8fcde895353d99794d4119267b7494 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;