ITK Geodesic Active Contour 3D Example (Java)

So like promised in the last posting the sample source code for a three dimensional Geodesic Active Contour. (Please note that you need the patches from the previous posting for this to work)

But first I have added links to some good books covering the Segmentation Algorithm used in this example and the whole field of either PDE Segementation, Level Sets, or Visualization -- after all Christmas is around the corner. I just got Shapiro's book in the mail (left) and I might post a thorough review later;-)



And here comes the code. Look at p. 537 of the ITK User Guide for more information on the parameters and the algorithms:


private void myITKCalls ( String inputFile, String OutputFile, float propagationScaling,
            double sigma, double alpha, double beta, int seedX, int seedY, int seedZ,
            float initialDistance ) {
    /* # We're going to build the following pipelines:
    # 1. reader -> smoothing -> gradientMagnitude -> sigmoid -> FI
    # 2. fastMarching -> geodesicActiveContour ( FI ) -> thresholder -> writer
    # The output of pipeline 1 is a feature image that is used by the
    # geodesicActiveContour object.  Also see figure 9.18 in the ITK
    # Software Guide.*/



        itkImageFileReaderUS3_Pointer reader = itkImageFileReaderUS3.itkImageFileReaderUS3_New ( ) ;
            itkImageFileWriterUC3_Pointer writer = itkImageFileWriterUC3.itkImageFileWriterUC3_New ( ) ;

            itkCastImageFilterUS3F3_Pointer inputCast = itkCastImageFilterUS3F3.itkCastImageFilterUS3F3_New ( ) ;

            itkCurvatureAnisotropicDiffusionImageFilterF3F3_Pointer smoothing = itkCurvatureAnisotropicDiffusionImageFilterF3F3.itkCurvatureAnisotropicDiffusionImageFilterF3F3_New ( ) ;

            itkRescaleIntensityImageFilterUS3UC3_Pointer outputCast = itkRescaleIntensityImageFilterUS3UC3.itkRescaleIntensityImageFilterUS3UC3_New ( ) ;

            itkGradientMagnitudeRecursiveGaussianImageFilterF3F3_Pointer gradientMagnitude = itkGradientMagnitudeRecursiveGaussianImageFilterF3F3.itkGradientMagnitudeRecursiveGaussianImageFilterF3F3_New ( ) ;

            itkSigmoidImageFilterF3F3_Pointer sigmoid = itkSigmoidImageFilterF3F3.itkSigmoidImageFilterF3F3_New ( ) ;
            
            itkFastMarchingImageFilterF3F3_Pointer fastMarching = itkFastMarchingImageFilterF3F3.itkFastMarchingImageFilterF3F3_New ( ) ;

            itkGeodesicActiveContourLevelSetImageFilterF3F3_Pointer geodesicActiveContour = itkGeodesicActiveContourLevelSetImageFilterF3F3.itkGeodesicActiveContourLevelSetImageFilterF3F3_New ( ) ;

            itkBinaryThresholdImageFilterF3US3_Pointer thresholder = itkBinaryThresholdImageFilterF3US3.itkBinaryThresholdImageFilterF3US3_New ( ) ;

            reader.SetFileName ( inputFile ) ;
            writer.SetFileName ( OutputFile ) ;

            inputCast.SetInput ( reader.GetOutput ( ) ) ;
            smoothing.SetInput ( inputCast.GetOutput ( ) ) ;

            smoothing.SetNumberOfIterations (   5 ) ;
            smoothing.SetTimeStep ( 0.0625 ) ;
            smoothing.SetConductanceParameter ( 3.0 ) ;
            smoothing.Update ( ) ;
            System.out.println ( "Smoothing done" ) ;

            gradientMagnitude.SetInput ( smoothing.GetOutput ( ) ) ;
            gradientMagnitude.SetSigma ( sigma ) ;
            gradientMagnitude.Update ( ) ;
            System.out.println ( "Gradient Magnitude done" ) ;

            sigmoid.SetInput ( gradientMagnitude.GetOutput ( ) ) ;
            sigmoid.SetOutputMinimum ( ( float ) 0.0 ) ;
            sigmoid.SetOutputMaximum ( ( float ) 1.0 ) ;
            sigmoid.SetAlpha ( alpha ) ;
            sigmoid.SetBeta ( beta ) ;
            sigmoid.Update ( ) ;
            System.out.println ( "Sigmoid done" ) ;
            
            //Fast Marching
            itkIndex3 seedPosition = new itkIndex3 ( ) ; //seed
            seedPosition.SetElement ( 0, seedX ) ;
            seedPosition.SetElement ( 1, seedY ) ;
            seedPosition.SetElement ( 2, seedZ ) ;

            float seedValue = -initialDistance;
            itkLevelSetNodeF3 node = new itkLevelSetNodeF3 ( ) ;
            node.SetValue ( seedValue ) ;
            node.SetIndex ( seedPosition ) ;

            itkNodeContainerF3_Pointer seeds = itkNodeContainerF3.itkNodeContainerF3_New ( ) ;
            seeds.Initialize ( ) ;
            seeds.InsertElement ( 0, node ) ;

            fastMarching.SetTrialPoints ( seeds.GetPointer ( ) ) ;
            fastMarching.SetSpeedConstant ( 1.0 ) ;
            fastMarching.SetOutputSize ( reader.GetOutput ( ) .GetBufferedRegion ( ) .GetSize ( ) ) ;
            fastMarching.Update ( ) ;
            System.out.println ( "Fast Marching done" ) ;

            geodesicActiveContour.SetInput ( fastMarching.GetOutput ( ) ) ;
            geodesicActiveContour.SetFeatureImage ( sigmoid.GetOutput ( ) ) ;
            
            geodesicActiveContour.SetPropagationScaling ( propagationScaling ) ;
            geodesicActiveContour.SetCurvatureScaling ( ( float ) 1.0 ) ;
            geodesicActiveContour.SetAdvectionScaling ( ( float ) 1.0 ) ;
            geodesicActiveContour.SetMaximumRMSError ( 0.02 ) ;
            geodesicActiveContour.SetNumberOfIterations ( 800 ) ;
            geodesicActiveContour.Update ( ) ;
            System.out.println ( "Geodesic done" ) ;

            thresholder.SetLowerThreshold ( ( float ) -1000.0 ) ;
            thresholder.SetUpperThreshold ( ( float ) 0.0 ) ;
            thresholder.SetOutsideValue ( 0 ) ;
            thresholder.SetInsideValue ( 65535 ) ;
            
            thresholder.SetInput ( geodesicActiveContour.GetOutput ( ) ) ;
            System.out.println ( "Threshold done" ) ;

            outputCast.SetInput ( thresholder.GetOutput ( ) ) ;
            writer.SetInput ( outputCast.GetOutput ( ) ) ;

            outputCast.SetOutputMinimum (    ( short ) 0   ) ;
            outputCast.SetOutputMaximum ( ( short ) 255 ) ;

            writer.Update ( ) ;

     }