[Insight-users] Segmenation Error--after set advection image by GVF flow

[Amardeep Singh]

Dear Baoyun

Unfortunately, I have not yet encountered your error message, so I am 
not sure, if I can help you.
I would try to save the gradient vector field as a *.vtk file, not in 
Analyze format, as the gradient vector flow output is a vector field. 
Can Analyze deal with that? (I don't know.)
Furthermore, I would suggest that you encapsulate every Update() call in 
a try-catch block:
try
        {
           ...->Update();
        }
        catch( itk::ExceptionObject & excep )
        {
          std::cerr << "Exception caught !" << std::endl;
          std::cerr << excep << std::endl;
        }

This might give you more meaningful error messages.

Best regards
Amardeep

Baoyun Li wrote:
> Dear All:
>  
> The error message is: Program *received* signal *SIGFPE*, *Arithmetic 
> exception*.
>  
> It happens in itkImageHelp.h
>  
> index[NLoop] = static_cast<IndexValueType>(offset / offsetTable[NLoop]);
>  
>  
> inline static void ComputeIndexInner(const IndexType &bufferedRegionIndex,
>                                        OffsetValueType &offset,
>                                        const OffsetValueType 
> offsetTable[],
>                                        IndexType &index,
>                                        const UniqueTypeBoolFalse& )
>     {
>     index[NLoop] = static_cast<IndexValueType>(offset / 
> offsetTable[NLoop]);
>     offset = offset - (index[NLoop] * offsetTable[NLoop]);
>     index[NLoop] = index[NLoop] + bufferedRegionIndex[NLoop];
>     ImageHelper<NImageDimension, NLoop-1>::
>       ComputeIndexInner(bufferedRegionIndex,
>                         offset,
>                         offsetTable,
>                         index,
>                         Concept::Detail::UniqueType_bool<(NLoop==1)>());
>     }
>  
>  
> Does this mean some overflow, or my datatype has some problme.
>  
> Thanks
>  
> Baoyun
>
>  
>
> ------------------------------------------------------------------------
> *From:* Baoyun Li <baoyun_li123 at yahoo.com>
> *To:* Amardeep Singh <amar.singh at gmx.de>
> *Cc:* insight-users at itk.org
> *Sent:* Thursday, March 5, 2009 4:51:41 PM
> *Subject:* Re: Segmenation Error--after set advection image by GVF flow
>
> Dear Amardeep:
>  
> After I called
>
>  geodesicActiveContour->SetAdvectionImage(castflowfilter->GetOutput());
>  
> Do I need to call GenerateAdvectioniImage() since Atumatic Generate is 
> set off.
>  
> After I add : geodesicActiveContour->GenerateAdvectionImage();
>  
> I still got error, please give me some suggestions.
>  
> Baoyun
> ------------------------------------------------------------------------
> *From:* Baoyun Li <baoyun_li123 at yahoo.com>
> *To:* Amardeep Singh <amar.singh at gmx.de>
> *Cc:* insight-users at itk.org
> *Sent:* Thursday, March 5, 2009 3:25:12 PM
> *Subject:* Segmenation Error--after set advection image by GVF flow
>
> Dear Amardeep and All:
>  
> Now I can build the program after seting advection image by GVF flow. 
> However, when I run the program, the program shows segmeation error 
> when geodesicActiveContour->Update();
> The program stop at: itkSegmentationLevelSetFunction_txx  at  return ( 
> m_AdvectionImage->GetPixel(idx) );
>  
> template <class TImageType, class TFeatureImageType>
> typename SegmentationLevelSetFunction<TImageType, 
> TFeatureImageType>::VectorType
> SegmentationLevelSetFunction<TImageType, TFeatureImageType>
> ::AdvectionField(const NeighborhoodType &neighborhood,
>                  const FloatOffsetType &offset, GlobalDataStruct *)  const
> {
>   IndexType idx = neighborhood.GetIndex();
>   ContinuousIndexType cdx;
>   for (unsigned i = 0; i < ImageDimension; ++i)
>     {
>     cdx[i] = static_cast<double>(idx[i]) - offset[i];
>     }
>   if ( m_VectorInterpolator->IsInsideBuffer(cdx) )
>     {
>     return ( 
> m_VectorCast(m_VectorInterpolator->EvaluateAtContinuousIndex(cdx)));
>     }
>   //Just return the default else
>     return ( m_AdvectionImage->GetPixel(idx) );
>  
> }
>  
> Clearly, something is still wrong with SetAdvectionImage. Below is 
> related code.
>  
> Can you please help me figure out which setting is wrong?
>  
> Thanks
>  
> Baoyun
>
>  
>   typedef itk::CovariantVector<double, Dim> myGradientType;
>   typedef itk::Image<myGradientType, Dim>   myGradientImageType;
>   typedef itk::GradientVectorFlowImageFilter<myGradientImageType, 
> myGradientImageType>
>                                               myGVFFilterType;
>   typedef itk::GradientImageFilter<DoublelImageType, double,double>
>                                               myGFilterType;
>   typedef itk::GradientToMagnitudeImageFilter<myGradientImageType, 
> DoublelImageType>
>                                               myGToMFilterType;
>
>   typedef itk::GradientRecursiveGaussianImageFilter<
>                                             DoublelImageType,
>                                             myGradientImageType
>                                                   >  myFilterType;
>   typedef typename GeodesicActiveContourFilterType::VectorImageType   
> AdvectionImageType;
>   typedef itk::CastImageFilter<myGradientImageType,AdvectionImageType>
>    CastFlowFilterType;
>   typedef itk::LaplacianImageFilter<DoublelImageType, 
> DoublelImageType> myLaplacianFilterType;
>   typename myFilterType::Pointer filter = myFilterType::New();
>   typename myGFilterType::Pointer gfilter = myGFilterType::New();
>   typename myGToMFilterType::Pointer gtomfilter = myGToMFilterType::New();
>   typename myLaplacianFilterType::Pointer m_LFilter = 
> myLaplacianFilterType::New();
>   typename myGVFFilterType::Pointer m_GVFFilter = myGVFFilterType::New();
>   typename CastFlowFilterType::Pointer castflowfilter = 
> CastFlowFilterType::New();
>   typedef itk::ImageFileWriter< myGradientImageType >  WriterType2;
>   typename WriterType2::Pointer writer21=WriterType2::New();
>  
>
>   typename WriterType1::Pointer writer4=WriterType1::New();
>   writer4->SetInput(caster11->GetOutput());
>   writer4->SetFileName("../data/intial_contour.hdr");
>   writer4->Update();
>
>   sigmoid->SetOutputMinimum(  0.0  );
>   sigmoid->SetOutputMaximum(  1.0  );
>   caster12->SetInput(inputfilter);
>   smoothing->SetInput(caster12->GetOutput() );
>   gradientMagnitude->SetInput( smoothing->GetOutput() );
>   sigmoid->SetInput( gradientMagnitude->GetOutput() );
>   caster31->SetInput(inputfilter);
>   filter->SetInput( caster31->GetOutput() ); //caster short to double
>   filter->SetSigma( 1.0);
>   filter->Update();
>   gtomfilter->SetInput(filter->GetOutput());
>   gtomfilter->Update();
>   gfilter->SetInput(gtomfilter->GetOutput());
>   gfilter->Update();
>   m_GVFFilter->SetInput(gfilter->GetOutput());
>   m_GVFFilter->SetLaplacianFilter(m_LFilter);
>   m_GVFFilter->SetNoiseLevel(500);
>   m_GVFFilter->SetTimeStep(0.05);
>   m_GVFFilter->SetIterationNum(5);
>   m_GVFFilter->Update();
>   writer21->SetInput(m_GVFFilter->GetOutput());
>   writer21->SetFileName("../data/GVF_flow.hdr");
>   writer21->Update();
> //to get gradient manitude of GVF low
>  
>   geodesicActiveContour->SetAutoGenerateSpeedAdvection(false);
>   //typename myGradientImageType::Pointer m_GVFField;
>   //m_GVFField = m_GVFFilter->GetOutput();
>   castflowfilter->SetInput(m_GVFFilter->GetOutput());
>   geodesicActiveContour->SetAdvectionImage(castflowfilter->GetOutput());
>  // geodesicActiveContour->GenerateSpeedImage();
>   //geodesicActiveContour->SetAdvectionImage(m_GVFFilter->GetOutput());
>  
>  
>   
>   
>  
>   geodesicActiveContour->SetInput( caster11->GetOutput() );
>   geodesicActiveContour->SetFeatureImage( sigmoid->GetOutput() );
>   geodesicActiveContour->GenerateSpeedImage();
>   geodesicActiveContour->Update();
>   thresholder->SetInput( geodesicActiveContour->GetOutput() );
>   thresholder->Update();
>   outputfilter=thresholder->GetOutput();
>
>
>
>