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

Amardeep Singh amar.singh at gmx.de
Fri Mar 6 13:15:09 EST 2009


Dear Baoyun

I would suggest, you have a look at the following posting by Luis Ibanez:
http://www.nabble.com/units-of-SetNoiseLevel-in-itkGradientVectorFlowImageFilter-td21957034.html

So, you can leave out the SetTimeStep(). I found that a small number of 
iterations may be sufficient. At least, it was for me.
The exact amount of the noise level depends highly on your data. On a 
synthetic example, a value as high as 20000 worked for me.
On real data, I used values between 20 and 200 mostly, but I have not 
finished evaluating the best parameters yet. So, you should just give it 
a try and see what works best for you.

Best regards
Amardeep

Baoyun Li wrote:
> Dear Amardeep:
>  
> Finally, my code can work now.
>  
> I have to add geodesicActiveContour->GenerateSpeedImage(). Without it, 
> the program always shows segmentation error.
>  
>  
> Can you please share some experience about the selection of GVF 
> parameters for 3D images"
>  
> > >  m_GVFFilter->SetNoiseLevel(500);
> > >  m_GVFFilter->SetTimeStep(0.05);
> > >  m_GVFFilter->SetIterationNum(5);
>  
> Thanks for your help.
>  
> Baoyun
>  
>
>
>  
>
> ------------------------------------------------------------------------
> *From:* Amardeep Singh <amar.singh at gmx.de>
> *To:* Baoyun Li <baoyun_li123 at yahoo.com>; insight-users at itk.org
> *Sent:* Friday, March 6, 2009 10:44:54 AM
> *Subject:* Re: Segmenation Error--after set advection image by GVF flow
>
> Dear Baoyun
>
> I am not sure, I understand what you mean by saying "Since you just 
> recently fixed the GVF as advection term".
> I had hacked the code, because I did not know of the flag 
> "SetAutoGenerateSpeedAdvection(...)". In the end, I ended up using the 
> same code as you.
> You might want to have a look at the full posting:
> http://www.nabble.com/Use-of-advection-image-to-guide-geodesic-level-set-td22083962.html#a22159804
>
> Best regards
> Amardeep
>
> Baoyun Li wrote:
> > Dear Amardeep:
> >  Thank you for your help. Since you just recently fixed the GVF as 
> advection term, can you please send me a copy of your test code.
> >  I saw your original code (GVFtest.cc <http://gvftest.cc/>) pasted 
> on the web. But that is not the working version.
> >  Thanks.
> >  Baoyun
> >
> > ------------------------------------------------------------------------
> > *From:* Amardeep Singh <amar.singh at gmx.de <mailto:amar.singh at gmx.de>>
> > *To:* Baoyun Li <baoyun_li123 at yahoo.com <mailto:baoyun_li123 at yahoo.com>>
> > *Cc:* insight-users at itk.org <mailto:insight-users at itk.org>
> > *Sent:* Friday, March 6, 2009 4:36:35 AM
> > *Subject:* Re: Segmenation Error--after set advection image by GVF flow
> >
> > 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 
> <mailto:baoyun_li123 at yahoo.com> <mailto:baoyun_li123 at yahoo.com 
> <mailto:baoyun_li123 at yahoo.com>>>
> > > *To:* Amardeep Singh <amar.singh at gmx.de <mailto:amar.singh at gmx.de> 
> <mailto:amar.singh at gmx.de <mailto:amar.singh at gmx.de>>>
> > > *Cc:* insight-users at itk.org <mailto:insight-users at itk.org> 
> <mailto:insight-users at itk.org <mailto: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 
> <mailto:baoyun_li123 at yahoo.com> <mailto:baoyun_li123 at yahoo.com 
> <mailto:baoyun_li123 at yahoo.com>>>
> > > *To:* Amardeep Singh <amar.singh at gmx.de <mailto:amar.singh at gmx.de> 
> <mailto:amar.singh at gmx.de <mailto:amar.singh at gmx.de>>>
> > > *Cc:* insight-users at itk.org <mailto:insight-users at itk.org> 
> <mailto:insight-users at itk.org <mailto: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();
> > >
> > >
> > >
> > >
> >
>


More information about the Insight-users mailing list