<div dir="ltr"><div>I'm trying to find a vertebra in an image. I decided to use GeodesicActiveContourShapePriorLevelSetImageFilter. The code is right below.<br><br>typedef itk::GeodesicActiveContourShapePriorLevelSetImageFilter<InternalImageType, InternalImageType > GeodesicActiveContourFilterType;<br>GeodesicActiveContourFilterType::Pointer geodesicActiveContour =GeodesicActiveContourFilterType::New();<br><br>geodesicActiveContour->SetPropagationScaling( propagationScaling );<br>geodesicActiveContour->SetShapePriorScaling( shapePriorScaling );<br>geodesicActiveContour->SetCurvatureScaling( curvatureScaling); <br>geodesicActiveContour->SetAdvectionScaling( advectionScaling);<br><br><br><br>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.<br>For example:<br>propagationScaling=7<br>shapePriorScaling=0.01<br>curvatureScaling=5<br>advectionScaling=15<br> <br>This is why I imagine that my PCA model is wrong.<br><br>To create it I read two sources:<br><a href="https://github.com/midas-journal/midas-journal-812/blob/master/ActiveShapeModel/Create2DActiveShapeModel/2dasm.cxx" target="_blank">https://github.com/midas-journal/midas-journal-812/blob/master/ActiveShapeModel/Create2DActiveShapeModel/2dasm.cxx</a><br><a href="http://itk.org/Wiki/ITK/Examples/WishList/Segmentation/EstimatePCAModel" target="_blank">http://itk.org/Wiki/ITK/Examples/WishList/Segmentation/EstimatePCAModel</a><br><br>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.<br><br><span style="color:rgb(39,78,19)">...Segmentation/</span><span style="color:rgb(39,78,19)">EstimatePCAModel:</span>
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.<br><br> my_Estimatortype::VectorOfDoubleType v=filter->GetEigenValues();<br><span style="color:rgb(102,0,0)"> double sv_mean=sqrt(v[0]);</span><br> <br> for ( unsigned int <span style="color:rgb(102,0,0)">k = 0</span>; k < nb_modes; k++ )<br> {<br><span style="color:rgb(102,0,0)"> double sv=sqrt(v[k]);<br> double sv_n=sv/sv_mean;</span><br> std::cout << "writing: " << outFileNames[k] << std::endl;<br> <br> std::cout << "svd[" << k << "]: " << sv << " norm: " << sv_n << std::endl;<br> WriterType::Pointer writer = WriterType::New();<br> writer->SetFileName( outFileNames[k].c_str() );<br> scaler->SetInput(filter->GetOutput(k));<br> scaler->SetConstant(sv_n);<br> writer->SetInput( scaler->GetOutput() );<br> writer->Update();<br> }<br><br><span style="color:rgb(39,78,19)">...Create2DActiveShapeModel/</span><span style="color:rgb(39,78,19)">2dasm.cxx:</span> here the author divides every eigen vector by the sqrt of its eigen value.<br> for(unsigned int i = 0; i < NUMPC; i++)<br> {<br> //normalizing the images<br> DivideFilterType::Pointer divider = DivideFilterType::New();<br> divider->SetInput(<span style="color:rgb(102,0,0)">model-></span><span style="color:rgb(102,0,0)">GetOutput(i+1)</span>);<br> divider->SetScale(<span style="color:rgb(102,0,0)">1.0/sqrt(</span><span style="color:rgb(102,0,0)">eigenValues(i)</span>));<br> WriterType::Pointer myWriter = WriterType::New();<br> myWriter->SetFileName(outputImageNames[i].c_str());<br> myWriter->SetInput(divider->GetOutput());<br> myWriter->Update();<br> }<br><br><span style="color:rgb(39,78,19)">My final code was:</span><br> double sv_biggest=sqrt(autoVal(0));<br> <br> for(unsigned int i=1;i<=numPC;i++){<br> <br> ShiftScaleFilterType::Pointer escale=ShiftScaleFilterType::New();<br> escale->SetInput(model->GetOutput(i));<br> <br> double sv=sqrt(autoVal(i-1));<br> double sv_n=sv/sv_biggest;<br> escale->SetScale(sv_n);<br><br> resultPCA.push_back(escale->GetOutput());<br> escale->Update(); <br> }//for i.<br><br></div><div>Please, can anyone tell me whether (and what) I'm doing wrong? <br><br></div><div>Thanks in advance,<br><br></div><div>Viki.<br></div><div><br><br><br>////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////<br>////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////<br></div><div>////////////////////////////////////////////////////////////// FULL FUNCTIONS ///////////////////////////////////////////////////////////////////////////////////////<br></div><div>////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////<br>////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////<br><br><br><span style="color:rgb(39,78,19)">FixedImageType::Pointer ActiveShapeModel::modifyThreshold(FixedImageType:</span><span style="color:rgb(39,78,19)">:Pointer imagen){</span><br></div> //Defining types<br><div> typedef itk::BinaryThresholdImageFilter<FixedImageType,FixedImageType> FiltroUmbral;<br> typedef itk::SignedDanielssonDistanceMapImageFilter<FixedImageType,FixedImageType> FiltroSigno;<br> <br> <br> //Modifying thresholding<br> FiltroUmbral::Pointer umbral=FiltroUmbral::New(); <br> umbral->SetLowerThreshold(255); umbral->SetUpperThreshold(255); umbral->SetInsideValue(255); umbral->SetOutsideValue(0);<br> umbral->SetInput(imagen);<br> umbral->Update();<br><br> //Applying signed distance map.<br> FiltroSigno::Pointer signo=FiltroSigno::New();<br> signo->SetInput(umbral->GetOutput());<br> signo->Update();<br> <br> return(signo->GetOutput());<br> <br></div><div>}//modifyThreshold.<br></div><div><br><br><span style="color:rgb(39,78,19)">void ActiveShapeModel::PCAanalisys(</span><span style="color:rgb(39,78,19)">int numPC,VectorImgs imgsEntreno,VectorImgs &resultadoPCA, vnl_vector<double> &autoValores){</span><br> //Defining types<br> typedef itk::ImagePCAShapeModelEstimator<FixedImageType,FixedImageType> EstimadorPCA;<br> <br> //Instanciate model estimator.<br> EstimadorPCA::Pointer modelo=EstimadorPCA::New();<br> modelo->DebugOn();<br> modelo->SetNumberOfTrainingImages(imgsEntreno.size());<br> modelo->SetNumberOfPrincipalComponentsRequired(numPC);<br> <br> try{<br> int k=0; <br> for (VectorImgs::iterator it = imgsEntreno.begin() ; it != imgsEntreno.end(); ++it){ <br> //Picking the image from vector<br> FixedImageType::Pointer imagen=(FixedImageType::Pointer) *it;<br><br> //Translating to signed distance map image.<br> imagen=modifyThreshold(imagen);<br> <br> //Adding it to the estimator.<br> modelo->SetInput(k,imagen); <br> k++;<br> }//for vector. <br> }catch( itk::ExceptionObject & excep ){<br> std::cerr << "Exception caught !" << std::endl;<br> std::cerr << excep << std::endl;<br> }<br> <br> try{<br> //Invoking estimator.<br> modelo->Update();<br> modelo->Print(std::cout); <br> autoValores=modelo->GetEigenValues();<br> }catch( itk::ExceptionObject & excep ){<br> std::cerr << "Exception caught !" << std::endl;<br> std::cerr << excep << std::endl;<br> }<br> <br> try{<br> resultadoPCA.reserve(numPC+1);<br> <br> //Saving mean image.<br> resultadoPCA.push_back(modelo->GetOutput(0));<br> double sv_mayor=sqrt(autoValores(0));<br> <br> for(unsigned int i=1;i<=numPC;i++){<br></div><div> //Normalizing PCs<br></div><div> ShiftScaleFilterType::Pointer filtroEscala=ShiftScaleFilterType::New();<br> filtroEscala->SetInput(modelo->GetOutput(i));<br><br> double sv=sqrt(autoValores(i-1));<br> double sv_n=sv/sv_mayor;<br> filtroEscala->SetScale(sv_n);<br><br> resultadoPCA.push_back(filtroEscala->GetOutput());<br> filtroEscala->Update(); <br> }//for i.<br> }catch( itk::ExceptionObject & excep ){<br> std::cerr << "Exception caught !" << std::endl;<br> std::cerr << excep << std::endl;<br> }<br> <br>}//PCAanalisys.<br><br><br><span style="color:rgb(39,78,19)">VectorImgs ActiveShapeModel::fitModel(<br> FixedImageType::Pointer imgEntrada, FixedImageType::Pointer imgMedia, VectorImgs imgsModo, const int numIter,<br> const double conductSuavizado,const double initialDistance, const double sigma, <br>
const double propagationScaling, const double shapePriorScaling,const
double curvatureScaling, const double advectionScaling,<br> std::vector<int> x,std::vector<int>y, const double startX, const double startY){</span><br> <br> VectorImgs imgsSalida;<br><br><br>typedef FixedImageType InternalImageType;<br>typedef itk::ChangeInformationImageFilter<InternalImageType > CenterFilterType;<br>CenterFilterType::Pointer center = CenterFilterType::New();<br>center->CenterImageOn();<br>typedef itk::CurvatureAnisotropicDiffusionImageFilter<InternalImageType,InternalImageType > SmoothingFilterType;<br><br>SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();<br>smoothing->SetTimeStep( 0.125 );<br>smoothing->SetNumberOfIterations( 5 );<br>smoothing->SetConductanceParameter( conductSuavizado ); //9.0 por defecto para 2D<br><br><br>typedef itk::GradientMagnitudeRecursiveGaussianImageFilter<InternalImageType,InternalImageType > GradientFilterType;<br>GradientFilterType::Pointer gradientMagnitude = GradientFilterType::New();<br>gradientMagnitude->SetSigma(sigma);<br><br>typedef itk::BoundedReciprocalImageFilter<InternalImageType,InternalImageType > ReciprocalFilterType;<br>ReciprocalFilterType::Pointer reciprocal = ReciprocalFilterType::New();<br><br>typedef itk::FastMarchingImageFilter<InternalImageType,InternalImageType > FastMarchingFilterType;<br>FastMarchingFilterType::Pointer fastMarching = FastMarchingFilterType::New();<br>configurarConjuntoInicial<InternalImageType>(initialDistance,x,y,fastMarching); <br><br>fastMarching->SetSpeedConstant( 1.0 );<br><br>typedef itk::GeodesicActiveContourShapePriorLevelSetImageFilter<InternalImageType, InternalImageType > GeodesicActiveContourFilterType;<br>GeodesicActiveContourFilterType::Pointer geodesicActiveContour =GeodesicActiveContourFilterType::New();<br><br>geodesicActiveContour->SetPropagationScaling( propagationScaling );<br>geodesicActiveContour->SetShapePriorScaling( shapePriorScaling );<br>geodesicActiveContour->SetCurvatureScaling( curvatureScaling); //1.0 );<br>geodesicActiveContour->SetAdvectionScaling( advectionScaling);//1.0 );<br><br>geodesicActiveContour->SetMaximumRMSError( 0.005 );<br>geodesicActiveContour->SetNumberOfIterations( numIter);<br><br>geodesicActiveContour->SetNumberOfLayers( 4 );<br><br>center->SetInput(((InternalImageType::Pointer)imgEntrada));<br>smoothing->SetInput( center->GetOutput() );<br>gradientMagnitude->SetInput( smoothing->GetOutput() );<br>reciprocal->SetInput( gradientMagnitude->GetOutput() );<br>center->Update();<br>fastMarching->SetOutputRegion(center->GetOutput()->GetBufferedRegion() );<br>fastMarching->SetOutputSpacing(center->GetOutput()->GetSpacing() );<br>fastMarching->SetOutputOrigin(center->GetOutput()->GetOrigin() );<br><br>geodesicActiveContour->SetInput( fastMarching->GetOutput() );<br>geodesicActiveContour->SetFeatureImage( reciprocal->GetOutput() );<br><br>const unsigned int numberOfPCAModes = imgsModo.size();<br>typedef itk::PCAShapeSignedDistanceFunction<double,Dimension,FixedImageType > ShapeFunctionType;<br>ShapeFunctionType::Pointer shape = ShapeFunctionType::New();<br>shape->SetMeanImage(imgMedia);<br>shape->SetNumberOfPrincipalComponents( numberOfPCAModes );<br>shape->SetPrincipalComponentImages( imgsModo );<br><br>ShapeFunctionType::ParametersType pcaStandardDeviations( numberOfPCAModes );<br>pcaStandardDeviations.Fill( 1.0 );<br>shape->SetPrincipalComponentStandardDeviations( pcaStandardDeviations );<br><br>typedef itk::Euler2DTransform<double> TransformType;<br>TransformType::Pointer transform = TransformType::New();<br>shape->SetTransform( transform );<br><br>typedef itk::ShapePriorMAPCostFunction<InternalImageType,PixelType > CostFunctionType;<br>CostFunctionType::Pointer costFunction = CostFunctionType::New();<br>CostFunctionType::WeightsType weights;<br>weights[0] = 1.0; <br>weights[1] = 20.0; <br>weights[2] = 1.0;<br>weights[3] = 1.0;<br>costFunction->SetWeights( weights );<br><br>CostFunctionType::ArrayType mean( shape->GetNumberOfShapeParameters() );<br>CostFunctionType::ArrayType stddev( shape->GetNumberOfShapeParameters() );<br>mean.Fill( 0.0 );<br>stddev.Fill( 1.0 );<br>costFunction->SetShapeParameterMeans( mean );<br>costFunction->SetShapeParameterStandardDeviations( stddev );<br><br>typedef itk::OnePlusOneEvolutionaryOptimizer OptimizerType;<br>OptimizerType::Pointer optimizer = OptimizerType::New();<br>typedef itk::Statistics::NormalVariateGenerator GeneratorType;<br>GeneratorType::Pointer generator = GeneratorType::New();<br>generator->Initialize( 20020702 );<br>optimizer->SetNormalVariateGenerator( generator );<br>OptimizerType::ScalesType scales( shape->GetNumberOfParameters() );<br>scales.Fill( 1.0 );<br>for( unsigned int k = 0; k < numberOfPCAModes; k++ ){<br>scales[k] = 20.0; // scales for the pca mode multiplier<br>}<br>scales[numberOfPCAModes] = 350.0; // scale for 2D rotation<br>optimizer->SetScales( scales );<br><br>double initRadius = 1.05;<br>double grow = 1.1;<br>double shrink = pow(grow, -0.25);<br>optimizer->Initialize(initRadius, grow, shrink);<br>optimizer->SetEpsilon(1.0e-6); // minimal search radius<br>optimizer->SetMaximumIteration(15);<br><br>ShapeFunctionType::ParametersType parameters(shape->GetNumberOfParameters() );<br>parameters.Fill( 0.0 );<br>parameters[numberOfPCAModes + 1] = startX;<br>parameters[numberOfPCAModes + 2] = startY;<br><br>geodesicActiveContour->SetShapeFunction( shape );<br>geodesicActiveContour->SetCostFunction( costFunction );<br>geodesicActiveContour->SetOptimizer( optimizer );<br>geodesicActiveContour->SetInitialParameters( parameters );<br>typedef CIUSegmentacion<GeodesicActiveContourFilterType> CommandType;<br>CommandType::Pointer observer = CommandType::New();<br>geodesicActiveContour->AddObserver( itk::IterationEvent(), observer );<br><br><br>try{<br> geodesicActiveContour->Update();<br>}catch( itk::ExceptionObject & excep ){<br> std::cerr << "Exception caught !" << std::endl;<br> std::cerr << excep << std::endl;<br>}<br><br><br>std::cout << std::endl;<br>std::cout << "Max. no. iterations: " << geodesicActiveContour->GetNumberOfIterations() << std::endl;<br>std::cout << "Max. RMS error: " << geodesicActiveContour->GetMaximumRMSError() << std::endl;<br>std::cout << std::endl;<br>std::cout << "No. elpased iterations: " << geodesicActiveContour->GetElapsedIterations() << std::endl;<br>std::cout << "RMS change: " << geodesicActiveContour->GetRMSChange() << std::endl;<br>std::cout << "Parameters: " << geodesicActiveContour->GetCurrentParameters() << std::endl;<br><br>imgsSalida.push_back(smoothing->GetOutput() );<br>imgsSalida.push_back(gradientMagnitude->GetOutput());<br>imgsSalida.push_back(reciprocal->GetOutput());<br>imgsSalida.push_back(fastMarching->GetOutput() );<br>imgsSalida.push_back(geodesicActiveContour->GetOutput()); <br><br>typedef itk::SpatialFunctionImageEvaluatorFilter<ShapeFunctionType,InternalImageType,InternalImageType > EvaluatorFilterType;<br>EvaluatorFilterType::Pointer evaluatorIni = EvaluatorFilterType::New();<br>evaluatorIni->SetInput(geodesicActiveContour->GetOutput());<br>evaluatorIni->SetFunction(shape); <br><br>shape->SetParameters(geodesicActiveContour->GetInitialParameters());<br>evaluatorIni->Update();<br>imgsSalida.push_back(evaluatorIni->GetOutput()); <br><br>EvaluatorFilterType::Pointer evaluatorFin = EvaluatorFilterType::New();<br>evaluatorFin->SetInput(geodesicActiveContour->GetOutput());<br>evaluatorFin->SetFunction(shape); <br><br>shape->SetParameters( geodesicActiveContour->GetCurrentParameters() );<br>evaluatorFin->Update();//evaluator->Modified();<br>imgsSalida.push_back(evaluatorFin->GetOutput()); <br><br>return imgsSalida;<br></div>}//fitModel</div>