[Insight-users] how to use fasrmarching in itkshapedetectionlevelsetimagefilter

[xujf]

Hi,everyone:
   I want to employ itkShapeDetectionLevelSetImageFilter to realize 3D Level set Segmentation in Python.
   The ShapeDetectionLevelSetImageFilter expects two inputs, the first being an initial Level Set in the form of an itk::Image, and the second being a feature image.I use an FastMarchingImageFilter to produce the initial level set as the distance function to a set of user-provided seeds.The code referred to FastMarchingImageFilter  is described as below:

itkfastmarchingfilter=itk.itkFastMarchingImageFilterF3F3_New()
seeds=itk.itkNodeContainerF3_New()

seedPosition=itk.itkIndex3()
seedPosition.SetElement(0,80)
seedPosition.SetElement(1,100)
seedPosition.SetElement(2,80)

node=itk.itkLevelSetNodeF3()
node.SetValue(0)
node.SetIndex(seedPosition)
seeds.Initialize()
seeds.InsertElement(0,node)

itkfastmarchingfilter.SetTrialPoints(seeds.GetPointer())
itkfastmarchingfilter.SetSpeedConstant(1.0)
itkfastmarchingfilter.SetStoppingValue(100.0)

itkfastmarchingfilter.SetOutputSize(itkgradfilter.GetOutput().GetBufferedRegion().GetSize())  ##  line 15

itkcurfilter.SetInput(itkimport.GetOutput())   #itkCurvatureAnisotropicDiffusionImageFilter
itkgradfilter.SetInput(itkimport.GetOutput()) #itkGradientMagnitudeRecursiveGaussianImageFilter
itksigfilter.SetInput(itkgradfilter.GetOutput())  #itkSigmoidImageFilter
itkshapedetectionfilter.SetInput(itkfastmarchingfilter.GetOutput()) #itkShapeDetectionLevelSetImageFilter
itkshapedetectionfilter.SetFeatureImage(itksigfilter.GetOutput())

..............

when I ran the code, an error appeared:
ERROR: In F:\\VTK\\Common\\vtkImageData.cxx, line 541
vtkImageData (06B79008): Requesting a point from an empty image.

After I deleted the line 15:itkfastmarchingfilter.SetOutputSize(),this error disappeared.
however,the segmented result is not correct. the extention of the segmented image is (0,15,0,15,0,15)
,which is different from the origin image\'s extention (0,180,0,180,0,180)

It makes me puzzled. what is the matter with my code?  the whole code is in the attachment.(not long)

Best Regards!

Jianfeng
 
-------------- next part --------------
import sys,os
import paths
import InsightToolkit as itk
from vtkMINCReader  import *
import ConnectVTKITK as connect

reader=vtkMINCReader()
vtkexport=vtkImageExport()
reader.SetFileName("F:/atamai/2003-01-27/atamai/examplesPmw/sbrain.mnc")

cast=vtkImageCast()
cast.SetOutputScalarTypeToFloat()
cast.SetInput(reader.GetOutput())

vtkexport.SetInput(cast.GetOutput())

itkimport=itk.itkVTKImageImportF3_New()

connect.ConnectVTKToITKF3(vtkexport,itkimport.GetPointer())

#itkedgefilter=itk.itkEdgePotentialImageFilterF3F3_New()
itkcurfilter=itk.itkCurvatureAnisotropicDiffusionImageFilterF3F3_New()
itkbinfilter=itk.itkBinaryThresholdImageFilterF3US3_New()
itkgradfilter=itk.itkGradientMagnitudeRecursiveGaussianImageFilterF3F3_New()
itksigfilter=itk.itkSigmoidImageFilterF3F3_New()
itkfastmarchingfilter=itk.itkFastMarchingImageFilterF3F3_New()
itkshapedetectionfilter=itk.itkShapeDetectionLevelSetImageFilterF3F3_New()

itkcurfilter.SetTimeStep( 0.0625);
itkcurfilter.SetNumberOfIterations(  5 );
itkcurfilter.SetConductanceParameter( 9.0 );

itksigfilter.SetOutputMinimum(0.0)
itksigfilter.SetOutputMaximum(1.0)
itksigfilter.SetAlpha(-1)
itksigfilter.SetBeta(150)

itkgradfilter.SetSigma(1.0)  
seeds=itk.itkNodeContainerF3_New()

seedPosition=itk.itkIndex3()
seedPosition.SetElement(0,80)
seedPosition.SetElement(1,100)
seedPosition.SetElement(2,80)

node=itk.itkLevelSetNodeF3()
node.SetValue(0)
node.SetIndex(seedPosition)
seeds.Initialize()
seeds.InsertElement(0,node)
itkfastmarchingfilter.SetTrialPoints(seeds.GetPointer())
itkfastmarchingfilter.SetSpeedConstant(1.0)

itkfastmarchingfilter.SetStoppingValue(100.0)

itkbinfilter.SetLowerThreshold(0.0)
itkbinfilter.SetUpperThreshold(85.0 )
itkbinfilter.SetOutsideValue( 0  )
itkbinfilter.SetInsideValue(65535)
print itkimport.GetOutput()  
itkcurfilter.SetInput(itkimport.GetOutput())
itkgradfilter.SetInput(itkimport.GetOutput())
itksigfilter.SetInput(itkgradfilter.GetOutput())
itkshapedetectionfilter.SetInput(itkfastmarchingfilter.GetOutput())
itkshapedetectionfilter.SetFeatureImage(itksigfilter.GetOutput())
itkshapedetectionfilter.SetPropagationScaling(1.0)
itkshapedetectionfilter.SetCurvatureScaling(1.0)
itkshapedetectionfilter.SetMaximumRMSError(0.02)
itkshapedetectionfilter.SetNumberOfIterations(800)
#itkedgefilter.SetInput(itkgradfilter.GetOutput())
#itkfastmarchingfilter.SetInput(itksigfilter.GetOutput())
itkbinfilter.SetInput(itkfastmarchingfilter.GetOutput())

itkfastmarchingfilter.SetOutputSize(itkgradfilter.GetOutput().GetBufferedRegion().GetSize())

vtkimport=vtkImageImport()
itkexport=itk.itkVTKImageExportF3_New()
itkexport.SetInput(itkshapedetectionfilter.GetOutput())

connect.ConnectITKF3ToVTK(itkexport.GetPointer(),vtkimport)
cast1=vtkImageCast()
cast1.SetOutputScalarTypeToFloat()
cast1.SetInput(vtkimport.GetOutput())

lut=vtkLookupTable()
lut.SetTableRange(0,255)
lut.SetSaturationRange(0,0)
lut.SetHueRange(0,0)
lut.SetValueRange(0,255)

data=vtkimport.GetOutput()

vtkimport.Update()

filter=vtkImageDataGeometryFilter()
filter.SetInput(data)
filter.SetExtent(0,180,100,100,0,180)
print data.GetExtent()
mapper=vtkPolyDataMapper()
mapper.SetInput(filter.GetOutput())
mapper.SetLookupTable(lut)
mapper.ScalarVisibilityOn()
mapper.SetScalarRange(0,255)

actor=vtkActor()
actor.SetMapper(mapper)


ren = vtkRenderer()
renWin = vtkRenderWindow()
renWin.AddRenderer(ren)
iren = vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
ren.AddActor(actor)

iren.Initialize()
renWin.Render()
iren.Start()