Commit 86f7895b authored by David C. Lonie's avatar David C. Lonie
Browse files

Preserve extent ranges when extracting unsampled VOIs.

Previously, extent ranges would always start at 0. This makes sense for
downsampled extractions, but for dimensions where SampleRate was 1, it's
more intuitive to keep the existing sample ranges.
parent f2e626cb
......@@ -18,6 +18,7 @@
#include "vtkBoundingBox.h"
#include "vtkCellData.h"
#include "vtkIdList.h"
#include "vtkMath.h"
#include "vtkNew.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
......@@ -49,6 +50,8 @@ namespace vtk
namespace detail
{
// Index mapping works as:
// inputExtent = Mapping[dim][outputExtent - this->OutputWholeExtent[2*dim]]
struct vtkIndexMap
{
std::vector<int> Mapping[3];
......@@ -104,6 +107,11 @@ void vtkExtractStructuredGridHelper::Invalidate()
this->OutputWholeExtent[3]=-1;
this->OutputWholeExtent[4]= 0;
this->OutputWholeExtent[5]=-1;
for(int i=0; i < 3; ++i)
{
this->IndexMap->Mapping[ i ].clear();
}
}
//-----------------------------------------------------------------------------
......@@ -126,6 +134,17 @@ void vtkExtractStructuredGridHelper::Initialize(
return;
}
// Is the VOI valid?
if (voi[1] < voi[0] || voi[3] < voi[2] || voi[5] < voi[4])
{
this->Invalidate();
vtkWarningMacro("Invalid volume of interest: ["
<< " [ " << voi[0] << ", " << voi[1] << " ], "
<< " [ " << voi[2] << ", " << voi[3] << " ], "
<< " [ " << voi[4] << ", " << voi[5] << " ] ]");
return;
}
// Save the input parameters so we'll know when the map is out of date
std::copy(voi, voi + 6, this->VOI);
std::copy(wholeExtent, wholeExtent + 6, this->InputWholeExtent);
......@@ -138,10 +157,6 @@ void vtkExtractStructuredGridHelper::Initialize(
if(!wExtB.Intersects(voiB))
{
for(int i=0; i < 3; ++i)
{
this->IndexMap->Mapping[ i ].clear();
}
this->Invalidate();
vtkDebugMacro(<< "Extent ["
<< wholeExtent[0] << ", " << wholeExtent[1] << ", "
......@@ -160,6 +175,8 @@ void vtkExtractStructuredGridHelper::Initialize(
// Compute the output whole extent in the process.
for(int dim=0; dim < 3; ++dim)
{
// +2: +1 to include start/end points, +1 in case we need to append an
// extra point for includeBoundary edge cases.
this->IndexMap->Mapping[dim].resize(voi[2*dim+1]-voi[2*dim]+2);
int outIdx = 0;
......@@ -179,9 +196,13 @@ void vtkExtractStructuredGridHelper::Initialize(
}
this->IndexMap->Mapping[dim].resize(outIdx);
// Preserve the extent range when sample rate is 1, otherwise extents start
// at 0 if downsampling.
int offset = this->SampleRate[dim] == 1 ? voi[2*dim] : 0;
// Update output whole extent
this->OutputWholeExtent[2*dim] = 0;
this->OutputWholeExtent[2*dim+1] =
this->OutputWholeExtent[2*dim] = offset;
this->OutputWholeExtent[2*dim+1] = offset +
static_cast<int>( this->IndexMap->Mapping[dim].size()-1 );
} // END for all dimensions
}
......@@ -613,25 +634,37 @@ void vtkExtractStructuredGridHelper::GetPartitionedOutputExtent(
// cleaned up by the parallel filter using vtkStructuredImplicitConnectivity.
for (int dim = 0; dim < 3; ++dim)
{
// Ex: 0 | 4
EMIN(partitionedOutputExtent, dim) =
(EMIN(partitionedVOI, dim) - EMIN(globalVOI, dim)) / sampleRate[dim];
if (includeBoundary && EMAX(partitionedVOI, dim) == EMAX(globalVOI, dim))
if (sampleRate[dim] == 1)
{
int length = EMAX(partitionedVOI, dim) - EMIN(globalVOI, dim);
EMAX(partitionedOutputExtent, dim) = length / sampleRate[dim];
EMAX(partitionedOutputExtent, dim) +=
((length % sampleRate[dim]) == 0) ? 0 : 1;
}
else {
// Ex: 3 | 7
EMAX(partitionedOutputExtent, dim) =
(EMAX(partitionedVOI, dim) - EMIN(globalVOI, dim)) / sampleRate[dim];
// If we're not downsampling, just return the partitioned VOI:
EMIN(partitionedOutputExtent, dim) = EMIN(partitionedVOI, dim);
EMAX(partitionedOutputExtent, dim) = EMAX(partitionedVOI, dim);
}
else
{
// If we downsample, the global output VOI will be offset to start at 0,
// so we'll adjust the minimum
// Ex: 0 | 4
EMIN(partitionedOutputExtent, dim) =
(EMIN(partitionedVOI, dim) - EMIN(globalVOI, dim)) / sampleRate[dim];
// Account for any offsets in the OutputWholeExtent:
EMIN(partitionedOutputExtent, dim) += EMIN(outputWholeExtent, dim);
EMAX(partitionedOutputExtent, dim) += EMIN(outputWholeExtent, dim);
if (includeBoundary && EMAX(partitionedVOI, dim) == EMAX(globalVOI, dim))
{
int length = EMAX(partitionedVOI, dim) - EMIN(globalVOI, dim);
EMAX(partitionedOutputExtent, dim) = length / sampleRate[dim];
EMAX(partitionedOutputExtent, dim) +=
((length % sampleRate[dim]) == 0) ? 0 : 1;
}
else {
// Ex: 3 | 7
EMAX(partitionedOutputExtent, dim) =
(EMAX(partitionedVOI, dim) - EMIN(globalVOI, dim)) /
sampleRate[dim];
}
// Account for any offsets in the OutputWholeExtent:
EMIN(partitionedOutputExtent, dim) += EMIN(outputWholeExtent, dim);
EMAX(partitionedOutputExtent, dim) += EMIN(outputWholeExtent, dim);
}
}
}
......@@ -1017,6 +1017,12 @@ void vtkRectilinearGrid::Crop(const int* updateExtent)
{
return;
}
// Invalid extents would lead to unpleasant results:
else if (ext[1] < ext[0] || ext[3] < ext[2] || ext[5] < ext[4] ||
uExt[1] < uExt[0] || uExt[3] < uExt[2] || uExt[5] < uExt[4])
{
return;
}
else
{
vtkRectilinearGrid *newGrid;
......
......@@ -45,6 +45,7 @@ vtkExtractGrid::~vtkExtractGrid()
}
}
//------------------------------------------------------------------------------
int vtkExtractGrid::RequestInformation(
vtkInformation *vtkNotUsed(request),
vtkInformationVector **inputVector,
......@@ -60,6 +61,13 @@ int vtkExtractGrid::RequestInformation(
this->Internal->Initialize(
this->VOI,wholeExtent,this->SampleRate,(this->IncludeBoundary==1));
if (!this->Internal->IsValid())
{
vtkWarningMacro("Error while initializing filter.");
return 0;
}
this->Internal->GetOutputWholeExtent(outWholeExt);
outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),
......@@ -72,6 +80,11 @@ int vtkExtractGrid::RequestUpdateExtent(
vtkInformationVector **inputVector,
vtkInformationVector *outputVector)
{
if (!this->Internal->IsValid())
{
return 0;
}
int i;
// get the info objects
......@@ -97,16 +110,18 @@ int vtkExtractGrid::RequestUpdateExtent(
int oUExt[6];
outputVector->GetInformationObject(0)->Get(
vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), oUExt);
int oWExt[6]; // For parallel parititon this will be different.
this->Internal->GetOutputWholeExtent(oWExt);
for (i=0; i<3; i++)
{
int idx = oUExt[2*i];
if (idx < 0 || oUExt[2*i] >= (int)this->Internal->GetSize(i))
int idx = oUExt[2*i] - oWExt[2*i]; // Extent value to index
if (idx < 0 || idx >= (int)this->Internal->GetSize(i))
{
vtkWarningMacro("Requested extent outside whole extent.")
idx = 0;
}
uExt[2*i] = this->Internal->GetMappedExtentValueFromIndex(i, idx);
int jdx = oUExt[2*i+1];
int jdx = oUExt[2*i+1] - oWExt[2*i]; // Extent value to index
if (jdx < idx || jdx >= (int)this->Internal->GetSize(i))
{
vtkWarningMacro("Requested extent outside whole extent.")
......@@ -128,12 +143,22 @@ int vtkExtractGrid::RequestData(
vtkInformationVector **inputVector,
vtkInformationVector *outputVector)
{
return this->RequestDataImpl(this->VOI, inputVector, outputVector) ? 1 : 0;
if (!this->Internal->IsValid())
{
return 0;
}
// Set the output extent -- this is how RequestDataImpl knows what to copy.
vtkInformation *outInfo = outputVector->GetInformationObject(0);
vtkStructuredGrid *output = vtkStructuredGrid::SafeDownCast(
outInfo->Get(vtkDataObject::DATA_OBJECT()));
output->SetExtent(this->Internal->GetOutputWholeExtent());
return this->RequestDataImpl(inputVector, outputVector) ? 1 : 0;
}
//------------------------------------------------------------------------------
bool vtkExtractGrid::RequestDataImpl(int vtkNotUsed(voi)[6],
vtkInformationVector **inputVector,
bool vtkExtractGrid::RequestDataImpl(vtkInformationVector **inputVector,
vtkInformationVector *outputVector)
{
if( (this->SampleRate[0] < 1) ||
......@@ -154,31 +179,23 @@ bool vtkExtractGrid::RequestDataImpl(int vtkNotUsed(voi)[6],
vtkStructuredGrid *output = vtkStructuredGrid::SafeDownCast(
outInfo->Get(vtkDataObject::DATA_OBJECT()));
vtkPointData *pd=input->GetPointData();
vtkCellData *cd=input->GetCellData();
vtkPointData *outPD=output->GetPointData();
vtkCellData *outCD=output->GetCellData();
int *inExt;
int *outExt;
vtkPoints *newPts, *inPts;
vtkDebugMacro(<< "Extracting Grid");
if (input->GetNumberOfPoints() == 0)
{
return true;
}
inPts = input->GetPoints();
inExt = input->GetExtent();
vtkPointData *pd=input->GetPointData();
vtkCellData *cd=input->GetCellData();
vtkPointData *outPD=output->GetPointData();
vtkCellData *outCD=output->GetCellData();
vtkPoints *inPts = input->GetPoints();
int *inExt = input->GetExtent();
int begin[3];
int end[3];
this->Internal->ComputeBeginAndEnd(inExt,this->VOI,begin,end);
output->SetExtent(begin[0], end[0], begin[1], end[1], begin[2], end[2]);
vtkPoints *newPts = inPts->NewInstance();
int *outExt = output->GetExtent();
newPts = inPts->NewInstance();
outExt = output->GetExtent();
vtkDebugMacro(<< "Extracting Grid");
this->Internal->CopyPointsAndPointData(inExt,outExt,pd,inPts,outPD,newPts);
output->SetPoints(newPts);
......
......@@ -88,9 +88,9 @@ protected:
// Description:
// Implementation for RequestData using a specified VOI. This is because the
// parallel filter needs to muck around with the VOI to get spacing and
// partitioning to play nice.
bool RequestDataImpl(int voi[6],
vtkInformationVector **inputVector,
// partitioning to play nice. The VOI is calculated from the output
// data object's extents in this implementation.
bool RequestDataImpl(vtkInformationVector **inputVector,
vtkInformationVector *outputVector);
......
......@@ -54,6 +54,11 @@ int vtkExtractRectilinearGrid::RequestUpdateExtent(
vtkInformationVector **inputVector,
vtkInformationVector *outputVector)
{
if (!this->Internal->IsValid())
{
return 0;
}
int i;
// get the info objects
vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
......@@ -78,16 +83,18 @@ int vtkExtractRectilinearGrid::RequestUpdateExtent(
int oUExt[6];
outputVector->GetInformationObject(0)->Get(
vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), oUExt);
int oWExt[6]; // For parallel parititon this will be different.
this->Internal->GetOutputWholeExtent(oWExt);
for (i=0; i<3; i++)
{
int idx = oUExt[2*i];
if (idx < 0 || oUExt[2*i] >= (int)this->Internal->GetSize(i))
int idx = oUExt[2*i] - oWExt[2*i]; // Extent value to index
if (idx < 0 || idx >= (int)this->Internal->GetSize(i))
{
vtkWarningMacro("Requested extent outside whole extent.")
idx = 0;
}
uExt[2*i] = this->Internal->GetMappedExtentValueFromIndex(i, idx);
int jdx = oUExt[2*i+1];
int jdx = oUExt[2*i+1] - oWExt[2*i]; // Extent value to index
if (jdx < idx || jdx >= (int)this->Internal->GetSize(i))
{
vtkWarningMacro("Requested extent outside whole extent.")
......@@ -121,6 +128,12 @@ int vtkExtractRectilinearGrid::RequestInformation(
this->VOI,wholeExtent,this->SampleRate,(this->IncludeBoundary==1));
this->Internal->GetOutputWholeExtent(outWholeExt);
if (!this->Internal->IsValid())
{
vtkWarningMacro("Error while initializing filter.");
return 0;
}
outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),
outWholeExt, 6);
return 1;
......@@ -132,12 +145,23 @@ int vtkExtractRectilinearGrid::RequestData(
vtkInformationVector **inputVector,
vtkInformationVector *outputVector)
{
return this->RequestDataImpl(this->VOI, inputVector, outputVector) ? 1 : 0;
if (!this->Internal->IsValid())
{
return 0;
}
// Set the output extent -- this is how RequestDataImpl knows what to copy.
vtkInformation *outInfo = outputVector->GetInformationObject(0);
vtkRectilinearGrid *output = vtkRectilinearGrid::SafeDownCast(
outInfo->Get(vtkDataObject::DATA_OBJECT()));
output->SetExtent(this->Internal->GetOutputWholeExtent());
return this->RequestDataImpl(inputVector, outputVector) ? 1 : 0;
}
//----------------------------------------------------------------------------
bool vtkExtractRectilinearGrid::RequestDataImpl(
int voi[6], vtkInformationVector **inputVector,
vtkInformationVector **inputVector,
vtkInformationVector *outputVector)
{
if( (this->SampleRate[0] < 1) ||
......@@ -158,42 +182,24 @@ bool vtkExtractRectilinearGrid::RequestDataImpl(
vtkRectilinearGrid *output = vtkRectilinearGrid::SafeDownCast(
outInfo->Get(vtkDataObject::DATA_OBJECT()));
vtkPointData *pd=input->GetPointData();
vtkCellData *cd=input->GetCellData();
vtkPointData *outPD=output->GetPointData();
vtkCellData *outCD=output->GetCellData();
int *inExt;
int *outExt;
vtkDebugMacro(<< "Extracting Grid");
if (input->GetNumberOfPoints() == 0)
{
return true;
}
inExt = input->GetExtent();
this->Internal->Initialize(voi, inExt, this->SampleRate,
(this->IncludeBoundary == 1));
if (!this->Internal->IsValid())
{
vtkErrorMacro("Error initializing index map.");
return false;
}
int begin[3];
int end[3];
this->Internal->ComputeBeginAndEnd(inExt, voi, begin, end);
output->SetExtent(begin[0], end[0], begin[1], end[1], begin[2], end[2]);
vtkPointData *pd=input->GetPointData();
vtkCellData *cd=input->GetCellData();
vtkPointData *outPD=output->GetPointData();
vtkCellData *outCD=output->GetCellData();
outExt = output->GetExtent();
int *inExt = input->GetExtent();
int *outExt = output->GetExtent();
int outDims[3];
vtkStructuredData::GetDimensionsFromExtent(outExt,outDims);
vtkDebugMacro(<< "Extracting Grid");
this->Internal->CopyPointsAndPointData(inExt,outExt,pd,NULL,outPD,NULL);
this->Internal->CopyCellData(inExt,outExt,cd,outCD);
// copy coordinates
......@@ -211,10 +217,11 @@ bool vtkExtractRectilinearGrid::RequestDataImpl(
vtkDataArray::CreateDataArray( in_coords[ dim ]->GetDataType() );
out_coords[ dim ]->SetNumberOfTuples( outDims[ dim ] );
for (int outIdx = begin[dim]; outIdx <= end[dim]; ++outIdx)
for (int oExtVal = outExt[2*dim]; oExtVal <= outExt[2*dim+1]; ++oExtVal)
{
int inIdx = this->Internal->GetMappedIndex(dim, outIdx);
out_coords[dim]->SetTuple(outIdx, inIdx, in_coords[dim]);
int outExtIdx = oExtVal - outExt[2*dim];
int inExtIdx = this->Internal->GetMappedIndex(dim, outExtIdx);
out_coords[dim]->SetTuple(outExtIdx, inExtIdx, in_coords[dim]);
} // END for all points along this dimension in the output
} // END for all dimensions
......
......@@ -77,9 +77,9 @@ protected:
// Description:
// Implementation for RequestData using a specified VOI. This is because the
// parallel filter needs to muck around with the VOI to get spacing and
// partitioning to play nice.
bool RequestDataImpl(int voi[6],
vtkInformationVector **inputVector,
// partitioning to play nice. The VOI is calculated from the output
// data object's extents in this implementation.
bool RequestDataImpl(vtkInformationVector **inputVector,
vtkInformationVector *outputVector);
int VOI[6];
......
......@@ -106,6 +106,11 @@ int vtkPExtractGrid::RequestData(
return this->Superclass::RequestData(request, inputVector, outputVector);
}
if (!this->Internal->IsValid())
{
return 0;
}
// Collect information:
vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
vtkInformation *outInfo = outputVector->GetInformationObject(0);
......@@ -155,11 +160,11 @@ int vtkPExtractGrid::RequestData(
// 1) Compute ParititonedVOI that will allow the base class to produce as much
// of the output data set as possible from the partitioned piece.
//
// 2) Extract PartitionedVOI using the base class's implementation.
//
// 3) Update the output dataset's extents to match PartitionedOutputExtent (it
// 2) Update the output dataset's extents to match PartitionedOutputExtent (it
// will be [0, L] in each dimension by default).
//
// 3) Extract PartitionedVOI using the base class's implementation.
//
// 4) Close gaps using vtkStructuredImplicitConnectivity (e.g. [3, 4] in the
// above example).
......@@ -171,6 +176,7 @@ int vtkPExtractGrid::RequestData(
}
DEBUG_EXTENT("InputWholeExtent", inputWholeExtent);
DEBUG_EXTENT("OutputWholeExtent", outputWholeExtent);
DEBUG_EXTENT("GlobalVOI", globalVOI);
DEBUG_EXTENT("InputPartitionedExtent", inputExtent);
......@@ -185,29 +191,31 @@ int vtkPExtractGrid::RequestData(
vtkExtractStructuredGridHelper::GetPartitionedVOI(
globalVOI, inputExtent, this->SampleRate, this->IncludeBoundary != 0,
partitionedVOI);
}
DEBUG_EXTENT("PartitionedVOI", partitionedVOI);
////////////////////////////////////////////////////////////
// 2) Extract actual VOI using superclass implementation: //
////////////////////////////////////////////////////////////
if (!this->Superclass::RequestDataImpl(partitionedVOI, inputVector,
outputVector))
{
return 0;
}
if (partitionContainsVOI)
{
////////////////////////////////////////////////////////////////
// 3) Compute and update the output dataset's actual extents. //
// 2) Compute and update the output dataset's actual extents. //
////////////////////////////////////////////////////////////////
vtkExtractStructuredGridHelper::GetPartitionedOutputExtent(
globalVOI, partitionedVOI, outputWholeExtent, this->SampleRate,
this->IncludeBoundary != 0, partitionedOutputExtent);
output->SetExtent(partitionedOutputExtent);
}
DEBUG_EXTENT("PartitionedVOI", partitionedVOI);
DEBUG_EXTENT("PartitionedOutputExtent", partitionedOutputExtent);
DEBUG_EXTENT("OutputWholeExtent", outputWholeExtent);
if (partitionContainsVOI)
{
////////////////////////////////////////////////////////////
// 3) Extract actual VOI using superclass implementation: //
////////////////////////////////////////////////////////////
if (!this->Superclass::RequestDataImpl(inputVector, outputVector))
{
return 0;
}
}
//////////////////////////////
// 4: Detect & resolve gaps //
......
......@@ -105,6 +105,11 @@ int vtkPExtractRectilinearGrid::RequestData(
return this->Superclass::RequestData(request, inputVector, outputVector);
}
if (!this->Internal->IsValid())
{
return 0;
}
// Collect information:
vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
vtkInformation *outInfo = outputVector->GetInformationObject(0);
......@@ -154,11 +159,11 @@ int vtkPExtractRectilinearGrid::RequestData(
// 1) Compute ParititonedVOI that will allow the base class to produce as much
// of the output data set as possible from the partitioned piece.
//
// 2) Extract PartitionedVOI using the base class's implementation.
//
// 3) Update the output dataset's extents to match PartitionedOutputExtent (it
// 2) Update the output dataset's extents to match PartitionedOutputExtent (it
// will be [0, L] in each dimension by default).
//
// 3) Extract PartitionedVOI using the base class's implementation.
//
// 4) Close gaps using vtkStructuredImplicitConnectivity (e.g. [3, 4] in the
// above example).
......@@ -170,6 +175,7 @@ int vtkPExtractRectilinearGrid::RequestData(
}
DEBUG_EXTENT("InputWholeExtent", inputWholeExtent);
DEBUG_EXTENT("OutputWholeExtent", outputWholeExtent);
DEBUG_EXTENT("GlobalVOI", globalVOI);
DEBUG_EXTENT("InputPartitionedExtent", inputExtent);
......@@ -184,30 +190,31 @@ int vtkPExtractRectilinearGrid::RequestData(
vtkExtractStructuredGridHelper::GetPartitionedVOI(
globalVOI, inputExtent, this->SampleRate, this->IncludeBoundary != 0,
partitionedVOI);
}
DEBUG_EXTENT("PartitionedVOI", partitionedVOI);
////////////////////////////////////////////////////////////
// 2) Extract actual VOI using superclass implementation: //
////////////////////////////////////////////////////////////
if (!this->Superclass::RequestDataImpl(partitionedVOI, inputVector,
outputVector))
{
return 0;
}
if (partitionContainsVOI)
{
////////////////////////////////////////////////////////////////
// 3) Compute and update the output dataset's actual extents. //
// 2) Compute and update the output dataset's actual extents. //
////////////////////////////////////////////////////////////////
vtkExtractStructuredGridHelper::GetPartitionedOutputExtent(
globalVOI, partitionedVOI, outputWholeExtent, this->SampleRate,
this->IncludeBoundary != 0, partitionedOutputExtent);
output->SetExtent(partitionedOutputExtent);
}
DEBUG_EXTENT("PartitionedVOI", partitionedVOI);
DEBUG_EXTENT("PartitionedOutputExtent", partitionedOutputExtent);
DEBUG_EXTENT("OutputWholeExtent", outputWholeExtent);
if (partitionContainsVOI)
{
////////////////////////////////////////////////////////////
// 3) Extract actual VOI using superclass implementation: //
////////////////////////////////////////////////////////////
if (!this->Superclass::RequestDataImpl(inputVector, outputVector))
{
return 0;
}
}