[ITK-users] Why does GeodesicActiveContourShapePriorLevelSetImageFilter not take into account the PCA model?

Fri Jun 19 03:35:40 EDT 2015

I'm trying to find a vertebra in an image. I decided to use
GeodesicActiveContourShapePriorLevelSetImageFilter. The code is right below.

typedef
itk::GeodesicActiveContourShapePriorLevelSetImageFilter<InternalImageType,
InternalImageType > GeodesicActiveContourFilterType;
GeodesicActiveContourFilterType::Pointer geodesicActiveContour
=GeodesicActiveContourFilterType::New();

geodesicActiveContour->SetPropagationScaling( propagationScaling );
geodesicActiveContour->SetShapePriorScaling( shapePriorScaling );
geodesicActiveContour->SetCurvatureScaling( curvatureScaling);
geodesicActiveContour->SetAdvectionScaling( advectionScaling);

In order to check my model, I cleaned an image and painted the background
in grey so the vertebra limits were pretty clear. But, although the
vertebra should be easy to find, I can't do it unless I choose a very low
shapePriorScaling parameter.
For example:
propagationScaling=7
shapePriorScaling=0.01
curvatureScaling=5
advectionScaling=15

This is why I imagine that my PCA model is wrong.

To create it I read two sources:
https://github.com/midas-journal/midas-journal-812/blob/master/ActiveShapeModel/Create2DActiveShapeModel/2dasm.cxx
http://itk.org/Wiki/ITK/Examples/WishList/Segmentation/EstimatePCAModel

I realized that the way both author normalize PCs are different, and since
this is an important issue in
GeodesicActiveContourShapePriorLevelSetImageFilter (it prefers normalized
eigen vectors), I think that maybe I do it in a wrong way.

...Segmentation/EstimatePCAModel: here the author strangely divides every
output of the filter (even the mean image) by the first eigen value and
multiply it by the eigen value corresponding to the previous eigen vector.

  my_Estimatortype::VectorOfDoubleType v=filter->GetEigenValues();
  double sv_mean=sqrt(v[0]);

  for ( unsigned int k = 0; k < nb_modes; k++ )
    {
    double sv=sqrt(v[k]);
    double sv_n=sv/sv_mean;
    std::cout << "writing: " << outFileNames[k] << std::endl;

    std::cout << "svd[" << k << "]: " << sv << " norm: " << sv_n <<
std::endl;
    WriterType::Pointer writer = WriterType::New();
    writer->SetFileName( outFileNames[k].c_str() );
    scaler->SetInput(filter->GetOutput(k));
    scaler->SetConstant(sv_n);
    writer->SetInput( scaler->GetOutput() );
    writer->Update();
    }

...Create2DActiveShapeModel/2dasm.cxx: here the author divides every eigen
vector by the sqrt of its eigen value.
    for(unsigned int i = 0; i < NUMPC; i++)
    {
    //normalizing the images
    DivideFilterType::Pointer divider = DivideFilterType::New();
    divider->SetInput(model->GetOutput(i+1));
    divider->SetScale(1.0/sqrt(eigenValues(i)));
    WriterType::Pointer myWriter = WriterType::New();
    myWriter->SetFileName(outputImageNames[i].c_str());
    myWriter->SetInput(divider->GetOutput());
    myWriter->Update();
    }

My final code was:
    double sv_biggest=sqrt(autoVal(0));

    for(unsigned int i=1;i<=numPC;i++){

        ShiftScaleFilterType::Pointer escale=ShiftScaleFilterType::New();
        escale->SetInput(model->GetOutput(i));

        double sv=sqrt(autoVal(i-1));
        double sv_n=sv/sv_biggest;
        escale->SetScale(sv_n);

        resultPCA.push_back(escale->GetOutput());
        escale->Update();
    }//for i.

Please, can anyone tell me whether (and what) I'm doing wrong?

Thanks in advance,

Viki.

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////// FULL
FUNCTIONS
///////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

FixedImageType::Pointer
ActiveShapeModel::modifyThreshold(FixedImageType::Pointer
imagen){
    //Defining types
    typedef itk::BinaryThresholdImageFilter<FixedImageType,FixedImageType>
FiltroUmbral;
    typedef
itk::SignedDanielssonDistanceMapImageFilter<FixedImageType,FixedImageType>
FiltroSigno;

    //Modifying thresholding
    FiltroUmbral::Pointer umbral=FiltroUmbral::New();
    umbral->SetLowerThreshold(255); umbral->SetUpperThreshold(255);
umbral->SetInsideValue(255); umbral->SetOutsideValue(0);
    umbral->SetInput(imagen);
    umbral->Update();

    //Applying signed distance map.
    FiltroSigno::Pointer signo=FiltroSigno::New();
    signo->SetInput(umbral->GetOutput());
    signo->Update();

    return(signo->GetOutput());

}//modifyThreshold.

void ActiveShapeModel::PCAanalisys(int numPC,VectorImgs
imgsEntreno,VectorImgs &resultadoPCA, vnl_vector<double> &autoValores){
    //Defining types
    typedef itk::ImagePCAShapeModelEstimator<FixedImageType,FixedImageType>
EstimadorPCA;

    //Instanciate model estimator.
    EstimadorPCA::Pointer modelo=EstimadorPCA::New();
    modelo->DebugOn();
    modelo->SetNumberOfTrainingImages(imgsEntreno.size());
    modelo->SetNumberOfPrincipalComponentsRequired(numPC);

    try{
    int k=0;
    for (VectorImgs::iterator it = imgsEntreno.begin() ; it !=
imgsEntreno.end(); ++it){
        //Picking the image from vector
        FixedImageType::Pointer imagen=(FixedImageType::Pointer) *it;

        //Translating to signed distance map image.
        imagen=modifyThreshold(imagen);

        //Adding it to the estimator.
        modelo->SetInput(k,imagen);
        k++;
    }//for vector.
    }catch( itk::ExceptionObject & excep ){
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
    }

    try{
    //Invoking estimator.
    modelo->Update();
    modelo->Print(std::cout);
    autoValores=modelo->GetEigenValues();
    }catch( itk::ExceptionObject & excep ){
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
    }

    try{
    resultadoPCA.reserve(numPC+1);

    //Saving mean image.
    resultadoPCA.push_back(modelo->GetOutput(0));
    double sv_mayor=sqrt(autoValores(0));

    for(unsigned int i=1;i<=numPC;i++){
        //Normalizing PCs
        ShiftScaleFilterType::Pointer
filtroEscala=ShiftScaleFilterType::New();
        filtroEscala->SetInput(modelo->GetOutput(i));

         double sv=sqrt(autoValores(i-1));
         double sv_n=sv/sv_mayor;
         filtroEscala->SetScale(sv_n);

        resultadoPCA.push_back(filtroEscala->GetOutput());
        filtroEscala->Update();
    }//for i.
    }catch( itk::ExceptionObject & excep ){
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
    }

}//PCAanalisys.

VectorImgs  ActiveShapeModel::fitModel(
    FixedImageType::Pointer imgEntrada, FixedImageType::Pointer imgMedia,
VectorImgs imgsModo, const int numIter,
    const double conductSuavizado,const double initialDistance, const
double sigma,
    const double propagationScaling, const double shapePriorScaling,const
double curvatureScaling, const double advectionScaling,
    std::vector<int> x,std::vector<int>y, const double startX, const double
startY){

    VectorImgs imgsSalida;

typedef FixedImageType InternalImageType;
typedef itk::ChangeInformationImageFilter<InternalImageType >
CenterFilterType;
CenterFilterType::Pointer center = CenterFilterType::New();
center->CenterImageOn();
typedef
itk::CurvatureAnisotropicDiffusionImageFilter<InternalImageType,InternalImageType
> SmoothingFilterType;

SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
smoothing->SetTimeStep( 0.125 );
smoothing->SetNumberOfIterations( 5 );
smoothing->SetConductanceParameter( conductSuavizado ); //9.0 por defecto
para 2D

typedef
itk::GradientMagnitudeRecursiveGaussianImageFilter<InternalImageType,InternalImageType
> GradientFilterType;
GradientFilterType::Pointer gradientMagnitude = GradientFilterType::New();
gradientMagnitude->SetSigma(sigma);

typedef
itk::BoundedReciprocalImageFilter<InternalImageType,InternalImageType >
ReciprocalFilterType;
ReciprocalFilterType::Pointer reciprocal = ReciprocalFilterType::New();

typedef itk::FastMarchingImageFilter<InternalImageType,InternalImageType >
FastMarchingFilterType;
FastMarchingFilterType::Pointer fastMarching =
FastMarchingFilterType::New();
configurarConjuntoInicial<InternalImageType>(initialDistance,x,y,fastMarching);

fastMarching->SetSpeedConstant( 1.0 );

typedef
itk::GeodesicActiveContourShapePriorLevelSetImageFilter<InternalImageType,
InternalImageType > GeodesicActiveContourFilterType;
GeodesicActiveContourFilterType::Pointer geodesicActiveContour
=GeodesicActiveContourFilterType::New();

geodesicActiveContour->SetPropagationScaling( propagationScaling );
geodesicActiveContour->SetShapePriorScaling( shapePriorScaling );
geodesicActiveContour->SetCurvatureScaling( curvatureScaling); //1.0 );
geodesicActiveContour->SetAdvectionScaling( advectionScaling);//1.0 );

geodesicActiveContour->SetMaximumRMSError( 0.005 );
geodesicActiveContour->SetNumberOfIterations( numIter);

geodesicActiveContour->SetNumberOfLayers( 4 );

center->SetInput(((InternalImageType::Pointer)imgEntrada));
smoothing->SetInput( center->GetOutput() );
gradientMagnitude->SetInput( smoothing->GetOutput() );
reciprocal->SetInput( gradientMagnitude->GetOutput() );
center->Update();
fastMarching->SetOutputRegion(center->GetOutput()->GetBufferedRegion() );
fastMarching->SetOutputSpacing(center->GetOutput()->GetSpacing() );
fastMarching->SetOutputOrigin(center->GetOutput()->GetOrigin() );

geodesicActiveContour->SetInput( fastMarching->GetOutput() );
geodesicActiveContour->SetFeatureImage( reciprocal->GetOutput() );

const unsigned int numberOfPCAModes = imgsModo.size();
typedef itk::PCAShapeSignedDistanceFunction<double,Dimension,FixedImageType
> ShapeFunctionType;
ShapeFunctionType::Pointer shape = ShapeFunctionType::New();
shape->SetMeanImage(imgMedia);
shape->SetNumberOfPrincipalComponents( numberOfPCAModes );
shape->SetPrincipalComponentImages( imgsModo );

ShapeFunctionType::ParametersType pcaStandardDeviations( numberOfPCAModes );
pcaStandardDeviations.Fill( 1.0 );
shape->SetPrincipalComponentStandardDeviations( pcaStandardDeviations );

typedef itk::Euler2DTransform<double> TransformType;
TransformType::Pointer transform = TransformType::New();
shape->SetTransform( transform );

typedef itk::ShapePriorMAPCostFunction<InternalImageType,PixelType >
CostFunctionType;
CostFunctionType::Pointer costFunction = CostFunctionType::New();
CostFunctionType::WeightsType weights;
weights[0] = 1.0;
weights[1] = 20.0;
weights[2] = 1.0;
weights[3] = 1.0;
costFunction->SetWeights( weights );

CostFunctionType::ArrayType mean( shape->GetNumberOfShapeParameters() );
CostFunctionType::ArrayType stddev( shape->GetNumberOfShapeParameters() );
mean.Fill( 0.0 );
stddev.Fill( 1.0 );
costFunction->SetShapeParameterMeans( mean );
costFunction->SetShapeParameterStandardDeviations( stddev );

typedef itk::OnePlusOneEvolutionaryOptimizer OptimizerType;
OptimizerType::Pointer optimizer = OptimizerType::New();
typedef itk::Statistics::NormalVariateGenerator GeneratorType;
GeneratorType::Pointer generator = GeneratorType::New();
generator->Initialize( 20020702 );
optimizer->SetNormalVariateGenerator( generator );
OptimizerType::ScalesType scales( shape->GetNumberOfParameters() );
scales.Fill( 1.0 );
for( unsigned int k = 0; k < numberOfPCAModes; k++ ){
scales[k] = 20.0; // scales for the pca mode multiplier
}
scales[numberOfPCAModes] = 350.0; // scale for 2D rotation
optimizer->SetScales( scales );

double initRadius = 1.05;
double grow = 1.1;
double shrink = pow(grow, -0.25);
optimizer->Initialize(initRadius, grow, shrink);
optimizer->SetEpsilon(1.0e-6); // minimal search radius
optimizer->SetMaximumIteration(15);

ShapeFunctionType::ParametersType parameters(shape->GetNumberOfParameters()
);
parameters.Fill( 0.0 );
parameters[numberOfPCAModes + 1] = startX;
parameters[numberOfPCAModes + 2] = startY;

geodesicActiveContour->SetShapeFunction( shape );
geodesicActiveContour->SetCostFunction( costFunction );
geodesicActiveContour->SetOptimizer( optimizer );
geodesicActiveContour->SetInitialParameters( parameters );
typedef CIUSegmentacion<GeodesicActiveContourFilterType> CommandType;
CommandType::Pointer observer = CommandType::New();
geodesicActiveContour->AddObserver( itk::IterationEvent(), observer );

try{
    geodesicActiveContour->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: " <<
geodesicActiveContour->GetNumberOfIterations() << std::endl;
std::cout << "Max. RMS error: " <<
geodesicActiveContour->GetMaximumRMSError() << std::endl;
std::cout << std::endl;
std::cout << "No. elpased iterations: " <<
geodesicActiveContour->GetElapsedIterations() << std::endl;
std::cout << "RMS change: " << geodesicActiveContour->GetRMSChange() <<
std::endl;
std::cout << "Parameters: " <<
geodesicActiveContour->GetCurrentParameters() << std::endl;

imgsSalida.push_back(smoothing->GetOutput() );
imgsSalida.push_back(gradientMagnitude->GetOutput());
imgsSalida.push_back(reciprocal->GetOutput());
imgsSalida.push_back(fastMarching->GetOutput() );
imgsSalida.push_back(geodesicActiveContour->GetOutput());

typedef
itk::SpatialFunctionImageEvaluatorFilter<ShapeFunctionType,InternalImageType,InternalImageType
> EvaluatorFilterType;
EvaluatorFilterType::Pointer evaluatorIni = EvaluatorFilterType::New();
evaluatorIni->SetInput(geodesicActiveContour->GetOutput());
evaluatorIni->SetFunction(shape);

shape->SetParameters(geodesicActiveContour->GetInitialParameters());
evaluatorIni->Update();
imgsSalida.push_back(evaluatorIni->GetOutput());

EvaluatorFilterType::Pointer evaluatorFin = EvaluatorFilterType::New();
evaluatorFin->SetInput(geodesicActiveContour->GetOutput());
evaluatorFin->SetFunction(shape);

shape->SetParameters( geodesicActiveContour->GetCurrentParameters() );
evaluatorFin->Update();//evaluator->Modified();
imgsSalida.push_back(evaluatorFin->GetOutput());

return imgsSalida;
}//fitModel
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/insight-users/attachments/20150619/97e7f443/attachment.html>