[Insight-developers] A fix for the bug

[Paul Yushkevich]

James,

Thank you for your detailed reply.  I have a couple of questions:

Could you explain why ReleaseDataFlag should not be activated for the 
first filter in the mini-pipeline, but only for filters 2..N?  I was 
mistakingly calling it for each fitler, but the code did not generate 
errors.

Should I check my fix into CVS?  If so, could you tell me where I can 
get the login/password for doing so, please?  Would you mind glancing 
over the code after I check it in just to make sure I didn't make a 
mistake somewhere?  I'm still learning the ropes of ITK filter internals.

Thanks again,

Paul.

Miller, James V (Research) wrote:

>I see Paul's point.  The GenerateInputRequestedRegion() for DiscreteGaussian
>is asking for the proper amount of the original input.  The problem is that
>the output requested region for the internal filter cannot set simple as the
>output requested region.  So the GenerateInputRequestedRegion() is correct,
>is the managing of the output requested region for the internal
>mini-pipeline that is problematic.  Paul's solution of actually establishing
>a complete mini-pipeline should work (I haven't looked at that code yet but
>it should have the standard GraftOutput() calls and set the ReleaseData flag
>on 2nd through Nth stage of the mini-pipeline.
>
>Another solution is to manage the output requested region for the
>mini-pipeline.  The original output of the DiscreteGaussian is grafted onto
>the mini-pipeline.  This sets the region ivars for the mini-pipeline as well
>as sets up the mini-pipeline to use the bulk data container of the
>DiscreteGaussian output.  
>
>What needs to be done is that each call to the mini-pipeline needs to have a
>different output requested region.  The first call will be for the original
>output requested region padded in N-1 of the dimensions (the mini-pipeline
>will pad the last dimension for the input requested region).  The second
>call will be for the original output requested region padded in N-2 of the
>dimensions, etc. The last call for the original output requested region.
>
>There are other places in the toolkit where we probably need to fix this.
>Pretty much anywhere we are using a single filter in a loop where the filter
>is NOT set up to require all the input and produce all the output.  
>
>Setting up a full mini-pipeline as Paul suggested will also address this
>problem. I think if the release data mechanism is used on the 2nd thru Nth
>stage of the mini-pipeline, the two approaches will have identical memory
>utilization.
>
>
>
>-----Original Message-----
>From: Paul Yushkevich [mailto:pauly at cognitica . com]
>Sent: Wednesday, May 28, 2003 2:37 PM
>To: insight-dev-list
>Subject: [Insight-developers] A fix for the bug
>
>
>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
--------------------------------