Commit 1c287e4f authored by David Gobbi's avatar David Gobbi Committed by Kitware Robot
Browse files

Merge topic 'imagedata-findcell-tolerance'

01c23323 BUG 11550: Add tolerance check to FindCell and FindAndGetCell
parents 51eff545 01c23323
......@@ -16,6 +16,7 @@ CREATE_TEST_SOURCELIST(Tests ${KIT}CxxTests.cxx
TestAMRBox.cxx
TestInterpolationFunctions.cxx
TestInterpolationDerivs.cxx
TestImageDataFindCell.cxx
TestImageIterator.cxx
TestGenericCell.cxx
TestGraph.cxx
......
/*=========================================================================
Program: Visualization Toolkit
Module: TestImageIterator.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.
=========================================================================*/
// .NAME Test FindCell methods for image data
// .SECTION Description
// This program tests the FindCell methods for vtkImageData to
// ensure that they give correct results near the boundaries and
// to ensure that tolerance is handled properly. Even when the
// tolerance is zero, points on the boundary must be considered
// to be inside the dataset.
#include "vtkDebugLeaks.h"
#include "vtkImageData.h"
#include "vtkSmartPointer.h"
inline int DoTest(
int extent[6], double origin[3], double spacing[3])
{
vtkSmartPointer<vtkImageData> image =
vtkSmartPointer<vtkImageData>::New();
image->SetExtent(extent);
image->SetOrigin(origin);
image->SetSpacing(spacing);
image->AllocateScalars();
double bounds[6];
image->GetBounds(bounds);
int subId = 0;
double pcoords[3];
double weights[8];
double x[3];
vtkIdType cellId;
double tol = 1e-4;
for (int i = 0; i < 3; i++)
{
x[0] = 0.5*(bounds[0] + bounds[1]);
x[1] = 0.5*(bounds[2] + bounds[3]);
x[2] = 0.5*(bounds[4] + bounds[5]);
for (int j = 0; j < 2; j++)
{
// test point right on the boundary with zero tolerance
x[i] = bounds[2*i+j];
cellId = image->FindCell(
x, 0, 0, 0.0, subId, pcoords, weights);
if (cellId < 0)
{
cerr << "point (" << x[0] << ", " << x[1] << ", " << x[2] << ")"
<< " should be in bounds (" << bounds[0] << ", " << bounds[1]
<< ", " << bounds[2] << ", " << bounds[3] << ", "
<< bounds[4] << ", " << bounds[5] << ") with tol 0.0\n";
return 1;
}
// test point just outside boundary with zero tolerance
double offset = ((j == 0) ? (-tol*0.5) : (tol*0.5));
x[i] = bounds[2*i+j] + offset;
cellId = image->FindCell(
x, 0, 0, 0.0, subId, pcoords, weights);
if (cellId >= 0)
{
cerr << "point (" << x[0] << ", " << x[1] << ", " << x[2] << ")"
<< " should be out of bounds (" << bounds[0] << ", " << bounds[1]
<< ", " << bounds[2] << ", " << bounds[3] << ", "
<< bounds[4] << ", " << bounds[5] << ") with tol 0.0\n";
return 1;
}
// test point just outside boundary with nonzero tolerance
x[i] = bounds[2*i+j] + ((j == 0) ? (-tol*0.5) : (tol*0.5));
cellId = image->FindCell(
x, 0, 0, tol*tol, subId, pcoords, weights);
if (cellId < 0)
{
cerr << "point (" << x[0] << ", " << x[1] << ", " << x[2] << ")"
<< " should be inside bounds (" << bounds[0] << ", " << bounds[1]
<< ", " << bounds[2] << ", " << bounds[3] << ", "
<< bounds[4] << ", " << bounds[5] << ") with tol " << tol << "\n";
return 1;
}
// check pcoords at boundaries
int isUpperBound = ((j == 1) ^ (spacing[i] < 0));
int isOnePixelThick = (extent[2*i] == extent[2*i+1]);
if (isUpperBound && !isOnePixelThick)
{
if (pcoords[i] != 1.0)
{
cerr << "at upper bounds, pcoord should be 1, but is "
<< pcoords[i] << "\n";
return 1;
}
}
else
{
if (pcoords[i] != 0.0)
{
cerr << "at lower bounds and for 0,1,2D cells, pcoord should be 0, "
<< "but is " << pcoords[i] << " " << extent[2*i] << ", " << extent[2*i+1] << ", " << j << "\n";
return 1;
}
}
// validate the cellId
x[i] = bounds[2*i+j];
double pcoords2[3];
int idx[3];
if (image->ComputeStructuredCoordinates(x, idx, pcoords2) == 0)
{
cerr << "ComputeStructuredCoordinates failed for "
<< "point (" << x[0] << ", " << x[1] << ", " << x[2] << ")"
<< " and bounds (" << bounds[0] << ", " << bounds[1]
<< ", " << bounds[2] << ", " << bounds[3] << ", "
<< bounds[4] << ", " << bounds[5] << ")\n";
return 1;
}
if (image->ComputeCellId(idx) != cellId)
{
cerr << "cellId = " << cellId << ", should be "
<< image->ComputeCellId(idx) << "\n";
return 1;
}
// validate the pcoords, allow a tolerance
double diff = (pcoords[i] - pcoords2[i]);
if (diff*diff > (1e-15)*(1e-15))
{
cerr << "pcoords[" << i << "] = " << pcoords[i] << ", should be "
<< pcoords2[i] << "\n";
return 1;
}
}
}
return 0;
}
int TestImageDataFindCell(int,char *[])
{
// test 0D, 1D, 2D, 3D data with various extents, spacings, origins
static int dims[4][3] = {
{ 1, 1, 1 }, { 3, 1, 1 }, { 3, 3, 1 }, { 3, 3, 3 } };
static int starts[4][3] = {
{ 0, 0, 0 }, { -1, 0, -1 }, { 2, 3, 6 }, { -10, 0, 5 } };
static double spacings[4][3] = {
{ 1, 1, 1 }, { 1.0/7, 1, 1 }, { 1, -1, 1 }, { -1, 1, -1/13.0 } };
static double origins[4][3] = {
{ 0, 0, 0 }, { 1.0/13, 0, 0 }, { 0, -1, 0 }, { -1, 0, -1/7.0 } };
int extent[6];
double *spacing;
double *origin;
int failed = 0;
for (int i = 0; i < 4; i++)
{
for (int j = 0; j < 4; j++)
{
for (int k = 0; k < 4; k++)
{
for (int l = 0; l < 4; l++)
{
for (int ii = 0; ii < 3; ii++)
{
extent[2*ii] = starts[i][ii];
extent[2*ii+1] = starts[i][ii] + dims[j][ii] - 1;
}
spacing = spacings[k];
origin = origins[l];
if (DoTest(extent, origin, spacing))
{
failed = 1;
}
}
}
}
}
return failed;
}
......@@ -741,155 +741,103 @@ vtkIdType vtkImageData::FindPoint(double x[3])
vtkIdType vtkImageData::FindCell(double x[3], vtkCell *vtkNotUsed(cell),
vtkGenericCell *vtkNotUsed(gencell),
vtkIdType vtkNotUsed(cellId),
double vtkNotUsed(tol2),
int& subId, double pcoords[3],
double *weights)
double tol2,
int& subId, double pcoords[3],
double *weights)
{
return
this->FindCell( x, NULL, 0, 0.0, subId, pcoords, weights );
this->FindCell( x, NULL, 0, tol2, subId, pcoords, weights );
}
//----------------------------------------------------------------------------
vtkIdType vtkImageData::FindCell(double x[3], vtkCell *vtkNotUsed(cell),
vtkIdType vtkNotUsed(cellId),
double vtkNotUsed(tol2),
double tol2,
int& subId, double pcoords[3], double *weights)
{
int loc[3];
int idx[3];
if ( this->ComputeStructuredCoordinates(x, loc, pcoords) == 0 )
// Compute the voxel index
if ( this->ComputeStructuredCoordinates(x, idx, pcoords) == 0 )
{
return -1;
// If voxel index is out of bounds, check point "x" against the
// bounds to see if within tolerance of the bounds.
const int* extent = this->Extent;
const double* origin = this->Origin;
const double* spacing = this->Spacing;
// Compute squared distance of point x from the boundary
double dist2 = 0.0;
for (int i=0; i<3; i++)
{
int minIdx = extent[i*2];
int maxIdx = extent[i*2+1];
double minBound = minIdx*spacing[i] + origin[i];
double maxBound = maxIdx*spacing[i] + origin[i];
if ( idx[i] < minIdx )
{
idx[i] = minIdx;
pcoords[i] = 0.0;
double dist = x[i] - minBound;
dist2 += dist*dist;
}
else if ( idx[i] >= maxIdx )
{
if (maxIdx == minIdx)
{
idx[i] = minIdx;
pcoords[i] = 0.0;
}
else
{
idx[i] = maxIdx-1;
pcoords[i] = 1.0;
}
double dist = x[i] - maxBound;
dist2 += dist*dist;
}
}
// Check squared distance against the tolerance
if (dist2 > tol2)
{
return -1;
}
}
// NOTE: Do not use the Voxel ivar for this. That ivar may be NULL
// if the dimensionality of the image data is less than 3.
vtkVoxel::InterpolationFunctions(pcoords,weights);
if (weights)
{
vtkVoxel::InterpolationFunctions(pcoords,weights);
}
//
// From this location get the cell id
//
subId = 0;
return this->ComputeCellId(loc);
return this->ComputeCellId(idx);
}
//----------------------------------------------------------------------------
vtkCell *vtkImageData::FindAndGetCell(double x[3],
vtkCell *vtkNotUsed(cell),
vtkIdType vtkNotUsed(cellId),
double vtkNotUsed(tol2), int& subId,
double tol2, int& subId,
double pcoords[3], double *weights)
{
int i, j, k, ijk[3], loc[3];
vtkIdType npts, idx;
double xOut[3];
int iMax = 0;
int jMax = 0;
int kMax = 0;;
vtkCell *cell = NULL;
const double *origin = this->Origin;
const double *spacing = this->Spacing;
const int* extent = this->Extent;
vtkIdType dims[3];
dims[0] = extent[1] - extent[0] + 1;
dims[1] = extent[3] - extent[2] + 1;
dims[2] = extent[5] - extent[4] + 1;
vtkIdType cellId = this->FindCell(x, 0, 0, tol2, subId, pcoords, 0);
if ( this->ComputeStructuredCoordinates(x, loc, pcoords) == 0 )
if (cellId < 0)
{
return NULL;
}
//
// Get the parametric coordinates and weights for interpolation
//
switch (this->DataDescription)
{
case VTK_EMPTY:
return NULL;
case VTK_SINGLE_POINT: // cellId can only be = 0
iMax = loc[0];
jMax = loc[1];
kMax = loc[2];
cell = this->Vertex;
break;
case VTK_X_LINE:
iMax = loc[0] + 1;
jMax = loc[1];
kMax = loc[2];
cell = this->Line;
break;
case VTK_Y_LINE:
iMax = loc[0];
jMax = loc[1] + 1;
kMax = loc[2];
cell = this->Line;
break;
case VTK_Z_LINE:
iMax = loc[0];
jMax = loc[1];
kMax = loc[2] + 1;
cell = this->Line;
break;
case VTK_XY_PLANE:
iMax = loc[0] + 1;
jMax = loc[1] + 1;
kMax = loc[2];
cell = this->Pixel;
break;
case VTK_YZ_PLANE:
iMax = loc[0];
jMax = loc[1] + 1;
kMax = loc[2] + 1;
cell = this->Pixel;
break;
case VTK_XZ_PLANE:
iMax = loc[0] + 1;
jMax = loc[1];
kMax = loc[2] + 1;
cell = this->Pixel;
break;
case VTK_XYZ_GRID:
iMax = loc[0] + 1;
jMax = loc[1] + 1;
kMax = loc[2] + 1;
cell = this->Voxel;
break;
}
vtkCell *cell = this->GetCell(cellId);
cell->InterpolateFunctions(pcoords, weights);
npts = 0;
for (k = loc[2]; k <= kMax; k++)
{
ijk[2] = k;
xOut[2] = origin[2] + k * spacing[2];
for (j = loc[1]; j <= jMax; j++)
{
ijk[1] = j;
xOut[1] = origin[1] + j * spacing[1];
for (i = loc[0]; i <= iMax; i++)
{
ijk[0] = i;
xOut[0] = origin[0] + i * spacing[0];
idx = this->ComputePointId(ijk);
cell->PointIds->SetId(npts,idx);
cell->Points->SetPoint(npts++,xOut);
idx++;
}
}
}
subId = 0;
return cell;
}
......@@ -1107,44 +1055,62 @@ int vtkImageData::ComputeStructuredCoordinates(double x[3], int ijk[3],
const double *spacing = this->Spacing;
const int* extent = this->Extent;
vtkIdType dims[3];
dims[0] = extent[1] - extent[0] + 1;
dims[1] = extent[3] - extent[2] + 1;
dims[2] = extent[5] - extent[4] + 1;
//
// Compute the ijk location
//
int isInBounds = 1;
for (i=0; i<3; i++)
{
d = x[i] - origin[i];
doubleLoc = d / spacing[i];
// Floor for negative indexes.
ijk[i] = vtkMath::Floor(doubleLoc);
if ( ijk[i] >= extent[i*2] && ijk[i] < extent[i*2 + 1] )
{
pcoords[i] = doubleLoc - static_cast<double>(ijk[i]);
}
else if ( ijk[i] < extent[i*2] || ijk[i] > extent[i*2+1] )
// low boundary check
if ( ijk[i] < extent[i*2])
{
return 0;
if ( x[i] == origin[i] + spacing[i]*extent[i*2] )
{
pcoords[i] = 0.0;
ijk[i] = extent[i*2];
}
else
{
isInBounds = 0;
}
}
else //if ( ijk[i] == extent[i*2+1] )
// high boundary check
else if ( ijk[i] >= extent[i*2 + 1] )
{
if (dims[i] == 1)
if ( x[i] == origin[i] + spacing[i]*extent[i*2 + 1] )
{
pcoords[i] = 0.0;
// make sure index is within the allowed cell index range
if (extent[i*2] < extent[i*2 + 1])
{
pcoords[i] = 1.0;
ijk[i] = extent[i*2 + 1] - 1;
}
else
{
pcoords[i] = 0.0;
ijk[i] = extent[i*2];
}
}
else
{
ijk[i] -= 1;
pcoords[i] = 1.0;
isInBounds = 0;
}
}
// else index is definitely within bounds
else
{
pcoords[i] = doubleLoc - static_cast<double>(ijk[i]);
}
}
return 1;
return isInBounds;
}
......
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