[Insight-developers] A fix for the bug

[Paul Yushkevich]

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
--------------------------------