[Insight-users] Problem in 3D segmentation

[Luis Ibanez]

Hi Rodrigo,


Two questions:


1) What is the size (in pixels) of the image that
    you are passing as input to this ITK pipeline ?


2) Does this error happens the first time you run
    the pipeline, or does it happens when you have
    ran it once and then you change the input image
    for another  one  ?



Please let us know,


    Thanks


      Luis


------------------------
Rodrigo Trujillo wrote:

> Hi,
>  
> I'm testing itkShapeDetectionLevelSetImageFilter in Python. With 2D 
> images (1 slice) my program
> works very well, but when i tried to expand the code for 3D images, the 
> following error occurs:
>  
> "
> Traceback (most recent call last):
>   File "C:\Meus 
> Documentos\Rodrigo\Segmentacao\ShapeDetection_3D_MRI.py", line 84, in ?
>     shapeD.Update()
>   File "C:/ITK/Release\itkShapeDetectionLevelSetImageFilter.py", line 
> 692, in Update
>     def Update(*args): return 
> apply(_itkShapeDetectionLevelSetImageFilter.itkShapeDetectionLevelSetImageFilterF3F3_Pointer_Update,args)
> RuntimeError: 
> C:\Compilacao\InsightToolkit-1.6.0\Code\Common\itkDataObject.cxx:376:
> Requested region is (at least partially) outside the largest possible 
> region. 
> "
>  
> Does somebody know the reason of the error ? And how can i solve ?
> My code is attached.
>  
>  
> Thanks,
>  
> Rodrigo Trujillo
> 
> 
> ------------------------------------------------------------------------
> 
> import vtk
> import ConnectVTKITK as itk
> from vtk.tk.vtkTkRenderWidget import *
> 
> reader = vtk.vtkImageReader2()
> reader.SetDataByteOrderToLittleEndian()
> reader.SetFilePrefix('C://TEMP//InVesalius//Jair')
> reader.SetDataSpacing(0.32,0.32,0.6)
> reader.SetDataExtent(0,255,0,255,2,27)
> reader.Update()
> 
> imageCast = vtk.vtkImageCast()
> imageCast.SetOutputScalarTypeToFloat()
> imageCast.SetInput(reader.GetOutput())
> imageCast.Update()
> 
> vtkExporter = vtk.vtkImageExport()
> vtkExporter.SetInput(imageCast.GetOutput())
> itkImporter = itk.itkVTKImageImportF3_New()
> itk.ConnectVTKToITKF3(vtkExporter, itkImporter.GetPointer())
> itkImporter.Update()  # Sem este update, dah pau !!!
> 
> # Passando um filtro de SMOOTHING
> filter  = itk.itkCurvatureAnisotropicDiffusionImageFilterF3F3_New()
> filter.SetInput(itkImporter.GetOutput())
> 
> filter.SetNumberOfIterations(6)
> filter.SetTimeStep(0.0625)
> filter.SetConductanceParameter(3.0)
> filter.Update()
> 
> # Criando o gradiente das imagem
> gradiente = itk.itkGradientMagnitudeRecursiveGaussianImageFilterF3F3_New()
> gradiente.SetInput(filter.GetOutput())
> gradiente.SetSigma(0.7)
> gradiente.Update()
> 
> # Passando sigmoid. Cria a funcao SPEED
> sigmoid = itk.itkSigmoidImageFilterF3F3_New()
> sigmoid.SetInput(gradiente.GetOutput())
> sigmoid.SetOutputMinimum(0)
> sigmoid.SetOutputMaximum(4095)
> sigmoid.SetAlpha(-1.0) # Deve ser negativo para inverter as cores
> sigmoid.SetBeta(14.0)
> sigmoid.Update()
> 
> # Pipeline do FastMarching
> fastMarching = itk.itkFastMarchingImageFilterF3F3_New()
> 
> sementes = itk.itkNodeContainerF3_New()     # Guardara as sementes (LevelSetNodes !!)
> node = itk.itkLevelSetNodeF3()              # Pega a semente e transforma num LevelSetNode
> 
> idx = itk.itkIndex3()                       # Selecionando uma semente
> idx.SetElement(0,130)
> idx.SetElement(1,114)
> idx.SetElement(2,20)
> 
> #node.SetValue(0.0)
> node.SetValue(-5.0)
> node.SetIndex(idx)
> 
> sementes.Initialize()                       # Inicializa as sementes
> sementes.InsertElement(0,node)
> 
> fastMarching.SetTrialPoints(sementes.GetPointer())
> fastMarching.SetOutputSize( gradiente.GetOutput().GetBufferedRegion().GetSize())
> fastMarching.SetSpeedConstant( 1.0 )
> fastMarching.SetStoppingValue(5.0)
> fastMarching.Update()
> 
> 
> # Implementacao do SHAPE
> 
> shapeD = itk.itkShapeDetectionLevelSetImageFilterF3F3_New()
> 
> shapeD.SetInput( fastMarching.GetOutput() )
> shapeD.SetFeatureImage( sigmoid.GetOutput() )
> 
> shapeD.SetPropagationScaling( 10.0 )
> shapeD.SetCurvatureScaling( 5.0 )
> 
> shapeD.SetMaximumRMSError( 0.02 )
> shapeD.SetNumberOfIterations( 400 )
> shapeD.Update()
> 
> # Threshold binario na imagem do Shape
> bt = itk.itkBinaryThresholdImageFilterF3US3_New()  # Retorna UShort
> bt.SetInput(shapeD.GetOutput())
> bt.SetLowerThreshold(0)
> bt.SetUpperThreshold(204)
> bt.SetOutsideValue(0)
> bt.SetInsideValue(2000)
> bt.Update()
> 
> itkExporter = itk.itkVTKImageExportUS3_New()
> itkExporter.SetInput(bt.GetOutput())
> vtkImporter = vtk.vtkImageImport()
> itk.ConnectITKUS3ToVTK(itkExporter.GetPointer(), vtkImporter)
> vtkImporter.Update()
> 
> 
> # Visualizacao
> 
> # Definindo a opacidade: Funcao de transferencia de opacidade
> 
> tfun = vtk.vtkPiecewiseFunction()
> tfun.AddPoint(0.0, 0)
> tfun.AddPoint(1.0, 0)
> #tfun.AddPoint(199.0, 0)
> #tfun.AddPoint(200.0, 0.5)
> tfun.AddPoint(2.0, 1.0)
> tfun.AddPoint(4095.0, 1.0)
> 
> # Definindo a funcao de transferencia de cores
> 
> ctfun = vtk.vtkColorTransferFunction()
> ctfun.AddRGBPoint(0.0, 1.0, 0.1, 0.1)
> ctfun.AddRGBPoint(4095.0, 1.0, 1.0, 1.0)
> 
> 
> # Criando o volume
> compositeFunction = vtk.vtkVolumeRayCastCompositeFunction()
> 
> volumeMapper = vtk.vtkVolumeRayCastMapper()
> volumeMapper.SetInput(vtkImporter.GetOutput())
> volumeMapper.SetVolumeRayCastFunction(compositeFunction)
> 
> volumeProperty = vtk.vtkVolumeProperty()
> volumeProperty.SetColor(ctfun)
> volumeProperty.SetScalarOpacity(tfun)
> volumeProperty.SetInterpolationTypeToLinear()
> volumeProperty.ShadeOn()
> 
> Ator = vtk.vtkVolume()
> Ator.SetMapper(volumeMapper)
> Ator.SetProperty(volumeProperty)
> 
> ren = vtk.vtkRenderer()
> ren.AddVolume(Ator)
> 
> # Criando a grade em volta do volume
> out = vtk.vtkOutlineFilter()
> out.SetInput(vtkImporter.GetOutput())
> 
> isoMapper2 = vtk.vtkOpenGLPolyDataMapper()
> isoMapper2.SetInput(out.GetOutput())
> isoMapper2.ScalarVisibilityOff()
> 
> isoActor2 = vtk.vtkOpenGLActor()
> isoActor2.SetMapper(isoMapper2)
> 
> ren.AddActor(isoActor2)
> 
> janela = Tkinter.Tk()
> 
> imageframe = Frame(janela)
> imageframe.pack(side=LEFT,fill = BOTH)
> 
> pane = vtkTkRenderWidget(imageframe,width = 512, height = 512)
> pane.GetRenderWindow().AddRenderer(ren)
> pane.pack()
> 
> janela.mainloop()