[Insight-users] Problem with itkN3MRIBiasFieldCorrectionImageFilter

[Ricardo Ferrari]

Dear itk users,

I am trying to run the itk N3 algorithm implemented by Nicholas J. Tustison
and James C. Gee  and I am getting a "segmentation fault" error.

Above I am sending the main part of my code and the error I am getting.

-------------------------------------
    /// m_inputImage is a minc image - you can download it from the
following site:
http://www.4shared.com/file/140891877/c11d6f0c/20090808_0_DUAL.html

     int m_shrinkFactor = 2;

    /// m_maskImage is purposely equal to zero so that it can be computed
using Otsu threshold

    ///
    /// Use Otsu threshold to compute a simple mask, if one is not provided
    ///
    if( ! m_maskImage )
    {
        typedef itk::OtsuThresholdImageFilter< TInputImageType,
TMaskImageType > ThresholderType;
        typename ThresholderType::Pointer otsu = ThresholderType::New();
        otsu->SetInput( m_inputImage );
        otsu->SetNumberOfHistogramBins( 200 );
        otsu->SetInsideValue( 0 );
        otsu->SetOutsideValue( 1 );
        otsu->Update();

        m_maskImage = otsu->GetOutput();
    }

    ///
    /// DownSampling the input image to reduce computation time
    ///
    typedef itk::ShrinkImageFilter< TInputImageType, TRealImageType >
ShrinkerImageType;
    typename ShrinkerImageType::Pointer shrinker1 =
ShrinkerImageType::New();
    shrinker1->SetInput( m_inputImage );
    shrinker1->SetShrinkFactors( 1 );

    typedef itk::ShrinkImageFilter< TMaskImageType, TMaskImageType >
ShrinkerMaskType;
    typename ShrinkerMaskType::Pointer shrinker2 = ShrinkerMaskType::New();
    shrinker2->SetInput( m_maskImage );
    shrinker2->SetShrinkFactors( 1 );

    if( m_shrinkFactor > 1 )
    {
        shrinker1->SetShrinkFactors( m_shrinkFactor );
        shrinker2->SetShrinkFactors( m_shrinkFactor );
        shrinker1->Update();
        shrinker2->Update();
    }

    typename TRealImageType::Pointer realInputImage =
shrinker1->GetOutput();
    typename TMaskImageType::Pointer maskShrinked = shrinker2->GetOutput();


    ///
    /// Run bias field correction
    ///
    typedef itk::N3MRIBiasFieldCorrectionImageFilter< TRealImageType,
TMaskImageType, TRealImageType > CorrectorType;
      typename CorrectorType::Pointer corrector = CorrectorType::New();
      corrector->SetInput( realInputImage );
      corrector->SetMaskImage( maskShrinked );
     corrector->SetMaskLabel( NumericTraits< TMaskImageType::PixelType
>::One );
    corrector->SetMaximumNumberOfIterations( m_maxNumIter );
    corrector->SetNumberOfFittingLevels( m_numFittingLevels );

    try
    {
        corrector->Update();
    }
    catch ( itk::ExceptionObject & err )
    {
        std::cout << "Caught an exception: " << std::endl;
        std::cout << err << " " << __FILE__ << " " << __LINE__ << std::endl;
        throw err;
    }
    catch (... )
    {
        std::cout << "Error while applying bias field correction" <<
std::endl;
        throw;
    }


OUTPUT
--------------------------
Starting
[-116.98, -103.715, -34.2648]
[0.857143, 0.857143, 6]
[280, 280, 25]
0.999208 -0.00996255 -0.0385217
0.0108005 0.999708 0.0216062
0.0382952 -0.0220051 0.999024

[3, 3, 3]
3
[4, 4, 4]
Segmentation fault
--------------------------


After spending some time looking for the problem I found out that the it is
happening in the following part of the
itkN3MRIBiasFieldCorrectionImageFilter.txx file.  The problem does not occur
if I use the original image without downsampling it.

I really appreciate any help on this !!!

Thanks,
Ricardo


cout << "Starting" << endl;

  cout << fieldEstimate->GetOrigin() << endl;
  cout << fieldEstimate->GetSpacing() << endl;
  cout << fieldEstimate->GetLargestPossibleRegion().GetSize() << endl;
  cout << fieldEstimate->GetDirection() << endl;
  cout << this->m_NumberOfFittingLevels << endl;
  cout << this->m_SplineOrder << endl;
  cout << this->m_NumberOfControlPoints << endl;

  typename BSplineFilterType::Pointer bspliner = BSplineFilterType::New();
  bspliner->SetOrigin( fieldEstimate->GetOrigin() );
  bspliner->SetSpacing( fieldEstimate->GetSpacing() );
  bspliner->SetSize( fieldEstimate->GetLargestPossibleRegion().GetSize() );
  bspliner->SetDirection( fieldEstimate->GetDirection() );
  bspliner->SetGenerateOutputImage( true );
  bspliner->SetNumberOfLevels( this->m_NumberOfFittingLevels );
  bspliner->SetSplineOrder( this->m_SplineOrder );
  bspliner->SetNumberOfControlPoints( this->m_NumberOfControlPoints );
  bspliner->SetInput( fieldPoints );
  bspliner->Update();

cout << "Finished" << endl;
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20091014/c72a209f/attachment.htm>