[Insight-developers] A fix for the bug
Paul Yushkevich
pauly at cognitica . com
Wed, 28 May 2003 14:37:08 -0400
Hi
I think I know the reason I was getting that bug (DiscreteGaussian
filter crashing when using a single-slice RequestedRegion)
I think the problem is in the implementation of the GenerateData()
method of the DG filter. For a 3D image, this implementation applies
NeighborhoodOperatorImageFilter 3 times. First time it's applied to the
input image, the second time it's applied to the output of the first
attempt, the third time it's applied to the result of the second
attempt. However, since there is only one filter object, I think (?)
that the requested region is not being propagated from the output to the
input correctly. That is, the BufferedRegion of the output of the first
application of the filter (call it BR1) ends up being equal to the
RequestedRegion of the output of the entire DiscreteGaussian filter
(call it RR). That can't be right, because BR1 should include some
padding, (i.e., BR1 = BR2 + padding, BR2 = BR3 + padding, BR3 = RR)
Here is a reimplementation of the DiscreteGaussian filter that
constructs a whole mini-pipeline consisting of N filters with N
operators. This mini-pipeline ensures that the correct requested
regions are propagated from the end of the pipeline to the input. I am
not sure if this introduces some inefficiencies into the code, and if
there is a better way to fix the problem.
template< class TInputImage, class TOutputImage >
void
DiscreteGaussianImageFilter<TInputImage, TOutputImage>
::GenerateData()
{
typename TOutputImage::Pointer output = this->GetOutput();
output->SetBufferedRegion(output->GetRequestedRegion());
output->Allocate();
// Type definition for the internal neighborhood filter
typedef NeighborhoodOperatorImageFilter<
InputImageType, OutputImageType> InternalFilterType;
typedef typename InternalFilterType::Pointer InternalFilterPointer;
// Create a chain of filters
InternalFilterPointer *filter =
new InternalFilterPointer[m_FilterDimensionality];
// Create a series of operators
GaussianOperator<OutputPixelType, ImageDimension> *oper =
new
GaussianOperator<OutputPixelType,ImageDimension>[m_FilterDimensionality];
// Initialize and connect the filters
for (unsigned int i = 0; i < m_FilterDimensionality; ++i)
{
// Set up the operator for this dimension
oper[i].SetDirection(i);
if (m_UseImageSpacing == true)
{
if (this->GetInput()->GetSpacing()[i] == 0.0)
{
itkExceptionMacro(<< "Pixel spacing cannot be zero");
}
else
{
oper[i].SetVariance(m_Variance[i] /
this->GetInput()->GetSpacing()[i]);
}
}
else
{
oper[i].SetVariance(m_Variance[i]);
}
oper[i].SetMaximumKernelWidth(m_MaximumKernelWidth);
oper[i].SetMaximumError(m_MaximumError[i]);
oper[i].CreateDirectional();
// Set up the filter
filter[i] = InternalFilterType::New();
filter[i]->SetOperator(oper[i]);
filter[i]->ReleaseDataFlagOn();
// Connect the filter to the input or to the previous filter
filter[i]->SetInput(i == 0 ? this->GetInput() :
filter[i-1]->GetOutput());
}
// Graft this filters output onto the mini-pipeline so that the
mini-pipeline
// has the correct region ivars and will write to this filters bulk data
// output.
filter[m_FilterDimensionality-1]->GraftOutput( output );
// Update the last filter in the chain
filter[m_FilterDimensionality-1]->Update();
// Graft the last output of the mini-pipeline onto this filters output so
// the final output has the correct region ivars and a handle to the final
// bulk data
this->GraftOutput(output);
// Clean up
delete[] oper;
delete[] filter;
}
--
--------------------------------
Paul A. Yushkevich, Ph.D.
President, Cognitica Corporation
17 Flemington Rd
Chapel Hill, NC 27517
Tel: 1-919-929-7652
--------------------------------