[Insight-users] fast marching 3D problem

[Gokhan Mustafa Uzunbas]

Hi Luca,
Thanks for your suggestions. I checked Task manager. It shows % 50 CPU  
  usage with 135 Mb virtual memory. I think cpu is doing some  
prcessing. I am sending my code below. I modified from   
"GeodesicActiveContourImageFilter.cxx". I replaced GeodesicActive with  
  ChanVese. I am using my own level set filter instead of  
itk::GeodesicActiveContourImageFilter. Here is my code:
     Thank you in advance.
     Gokhan
     int main( int argc, char *argv[] )
     {
        if( argc < 11 )
          {
          std::cerr << "Missing Parameters " << std::endl;
          std::cerr << "Usage: " << argv[0];
          std::cerr << " inputImage  outputImage";
          std::cerr << " seedX seedY seedZ InitialDistance";
          std::cerr << " Sigma SigmoidAlpha SigmoidBeta";
          std::cerr << " PropagationScaling"  << std::endl;
          return 1;
          }
        ThresholdingFilterType::Pointer thresholder1 =
     ThresholdingFilterType::New();

        thresholder1->SetLowerThreshold( -1000.0 );
        thresholder1->SetUpperThreshold(     0.0 );

        thresholder1->SetOutsideValue( 0 );
        thresholder1->SetInsideValue(  255 );


        //GT 3D Reader Begin
             VolumeReaderType::Pointer apReader = VolumeReaderType::New();
             std::string strFileName = "TestMRIData.dcm";
             apReader->SetFileName(strFileName.c_str());

             typedef itk::GDCMImageIO ImageIOType;
             ImageIOType::Pointer gdcmImageIO = ImageIOType::New();
             apReader->SetImageIO(gdcmImageIO);
             try {
                     apReader->Update();
             } catch(itk::ExceptionObject & e) {
                     std::cout << "Cannot read file " <<strFileName <<  
std::endl;
                     exit(-1);
             }

             itk::CastImageFilter<UnsignedShortVolumeType,
     InternalImageType>::Pointer pCastFilter =
                     itk::CastImageFilter<UnsignedShortVolumeType,  
InternalImageType>::New();
             pCastFilter->SetInput(apReader->GetOutput());
             pCastFilter->Update();

             InternalImageType::Pointer pReadImage  = pCastFilter->GetOutput();


             pReadImage->Update();

        typedef itk::RescaleIntensityImageFilter<
                                     InternalImageType,
                                     OutputImageType >   CastFilterType;

        typedef  itk::FastMarchingImageFilter<
                                    InternalImageType,
                                    InternalImageType >     
FastMarchingFilterType;

        FastMarchingFilterType::Pointer  fastMarching =
     FastMarchingFilterType::New();

        typedef  itk::ChanVeseLevelSetImageFilter< InternalImageType,
                      InternalImageType >    ChanVeseFilterType;
        ChanVeseFilterType::Pointer ChanVese =
                                           ChanVeseFilterType::New();

        const double propagationScaling = atof( argv[10] );
        ChanVese->SetPropagationScaling( propagationScaling );
        ChanVese->SetCurvatureScaling( 1.0 );
        ChanVese->SetAdvectionScaling( 1 );

        ChanVese->SetMaximumRMSError( 0.02 );
        ChanVese->SetNumberOfIterations(10);
        ChanVese->SetInput(  fastMarching->GetOutput() );
        ChanVese->SetFeatureImage( sigmoid->GetOutput() );
        ChanVese->SetFeatureImage( pReadImage );

        thresholder1->SetInput( ChanVese->GetOutput() );

        typedef FastMarchingFilterType::NodeContainer  NodeContainer;
        typedef FastMarchingFilterType::NodeType       NodeType;
        NodeContainer::Pointer seeds = NodeContainer::New();
        InternalImageType::IndexType  seedPosition;
        NodeType node;
        const double initialDistance = atof( argv[6] );
        int initialPosX = atoi( argv[3] );
        int initialPosY = atoi( argv[4] );
        int initialPosZ = atoi( argv[5] );
        //const double seedValue = - initialDistance;
        double seedValue = 0;
        int Id = 0;
        seeds->Initialize();
        for(int nI = -initialDistance; nI <= initialDistance; nI++){
               seedPosition[0] = initialPosX + nI;
               double diffX = (seedPosition[0] - initialPosX);
               double disX2 = diffX * diffX;
               for(int nY = -initialDistance; nY <= initialDistance; nY++){
                       seedPosition[1] = initialPosY + nY;
                       double diffY = (seedPosition[1] - initialPosY);
                       double disY2 = diffY * diffY;
                       for(int nZ = -initialDistance; nZ <=  
initialDistance; nZ++){
                             seedPosition[2] = initialPosZ + nZ;
                             double diffZ = (seedPosition[2] - initialPosZ);
                             double disZ2 = diffZ * diffZ;
                             seedValue = sqrt(disX2 + disY2 + disZ2)-  
initialDistance;
                             node.SetValue( seedValue );
                             node.SetIndex( seedPosition );
                             seeds->InsertElement( Id, node );
                             Id++;
                       }
               }
         }

        itk::Object::GlobalWarningDisplayOn();
        fastMarching->SetTrialPoints(  seeds  );
        fastMarching->SetSpeedConstant( 1.0 );
        fastMarching->SetOutputSize(
                 pReadImage->GetBufferedRegion().GetSize() );
        fastMarching->SetCollectPoints(true);
         try
          {
        thresholder1->Update();
          }
        catch( itk::ExceptionObject & excep )
          {
          std::cerr << "Exception caught !" << std::endl;
          std::cerr << excep << std::endl;
          }
        std::cout << std::endl;
        std::cout << "Max. no. iterations: " <<
     ChanVese->GetNumberOfIterations() << std::endl;
        std::cout << "Max. RMS error: " << ChanVese->GetMaximumRMSError()
     << std::endl;
        std::cout << std::endl;
        std::cout << "No. elpased iterations: " <<
     ChanVese->GetElapsedIterations() << std::endl;
        std::cout << "RMS change: " << ChanVese->GetRMSChange() << std::endl;

        return 0;
     }