How to fill the big holes for further processing of having a Riemannian surface?

Algorithms
1

Hi,

I am facing an issue for creating the SimplexMesh from the TriangleMesh. It seems the surface does not have Riemannian geometry.
I filled the holes by the following commands

typedef itk::BinaryFillholeImageFilter<ImageType> FillHoleImageFilterType;
 FillHoleImageFilterType::Pointer fillHoleFilter=FillHoleImageFilterType::New();
 fillHoleFilter->SetInput(image);
 fillHoleFilter->SetFullyConnected(true);
 fillHoleFilter->SetForegroundValue(itk::NumericTraits<PixelType>::max());
 fillHoleFilter->Update();

However, it seems there are some holes as I have shown with circles in the images:
snapshot0001 snapshot0002
According to some published work, they are using Conditional Random field (CRFs) to fill the holes.

Is there any suggestion for parameter setting of FillHoleImageFilterType or any other algorithm that can be used?

Is there any itk library for CRFs that can be used?

Thanks

2

I would be interested to hear more about the CRF algorithm, not sure if it is implemented in ITK or in any external module.

ITK has this VotingBinaryIterativeHoleFillingImageFilter that works fine, test it for your purposes.
Example from ITK docs

Edit: I would also recommend to use Slicer or an interactive ITK jupyter notebook for finding the right parameters during segmentation. Slicer uses ITK algorithms, but you get a 3D visualization, and faster feedback loop when changing parameters than doing it in the command line.
Once you find parameters that work, you can set them in your script to automatize the process.

1 Like
3

We have itk::MRFImageFilter with some background in the ITK Software Guide, Book 2, Section 5.4.5.

1 Like
4

Thanks for your suggestion,

DeepMedic has provided Dense3DCRF tool (however with python wrapper). Is it possible to include the libraries into my C++ code for itk?

5

Thank you,

Could I ask what are the mean values that is taking as input arguments in this line:

std::cerr << " mean1 mean2 … meanN " << std::endl;

and it seems that the means are being used in this line

centroid[0] = std::stod(argv[i + numberOfArgumentsBeforeMeans]);

Since I did not know, I could not proceed after the line mentioned above. And Also it is not clearly explained in the book.

Thanks

6

As noted in the book:

Figure 5.8 illustrates the effect of this filter with three classes. In this example the filter was run with
a smoothing factor of 3. The labeled image was produced by ScalarImageKmeansClassifier.cxx and
the means were estimated by ScalarImageKmeansModelEstimator.cxx.

1 Like
7

Thanks again,

Since I have used my algorithm to have the initial segmentation, how can I assign means for MRF (as they have used KMeans for the getting the centroid of each cluster (i.e., class here)?

8

One option is to use `itk::LabelStatisticsImageFilter.

1 Like
9

@matt.mccormick I wrote a function that takes two inputs (float 3D ImageType is the intensity input image and a unsigned char 3D image (the prediction output of my algorithm). Then, I calculated the means of predicted labels using itk::LabelStatisticsImageFilter. I set the parameters numberOfIterations=1000, smothingFactor=2.5, and SetErrorTolerance(1e-7), and the weights for 3D images I set according to this link. However, the output image is a blank image, which I do not know what is the reason causing this issue.

CharImageType::Pointer ApplyMRF(ImageType::Pointer image,CharImageType::Pointer prediction)
{
    
/*
     * This function applies MRF for the correction of
     * Since the Markov Random Field algorithm is defined in general for images whose pixels have multiple
       components, that is, images of vector type, we must adapt our scalar image in order to satisfy the
       interface expected by the MRFImageFilter. We do this by using the itk::ComposeImageFilter.
       With this filter we will present our scalar image as a vector image whose vector pixels contain a
       single component.
     */

    LabelStatisticsImageFilterType::Pointer labelStatisticsImageFilter=LabelStatisticsImageFilterType::New();
    labelStatisticsImageFilter->SetLabelInput(prediction);
    labelStatisticsImageFilter->SetInput(image);
    labelStatisticsImageFilter->Update();
    std::cout << "Number of labels: " << labelStatisticsImageFilter->GetNumberOfLabels() << std::endl;

    //
    // getting the label Statistics
    // Link: https://itk.org/Doxygen/html/classitk_1_1LabelStatisticsImageFilter.html
    //
    //using LabelPixelType=LabelStatisticsImageFilterType::LabelPixelType ;  //-
    //LabelStatisticsImageFilterType::LabelPixelType LabelPixelType;  //-

    float labelMeanValue [labelStatisticsImageFilter->GetNumberOfLabels()];  //an array for saving the means of labels

    for (auto vIt = labelStatisticsImageFilter->GetValidLabelValues().begin();
         vIt != labelStatisticsImageFilter->GetValidLabelValues().end();
         ++vIt) {
        if (labelStatisticsImageFilter->HasLabel(*vIt)) {
            //LabelPixelType labelValue = *vIt;  //-
            labelMeanValue[(int)*vIt]=labelStatisticsImageFilter->GetMean(*vIt);  //+
            std::cout << "Label: " << (int)*vIt << std::endl;
            //std::cout << "\tmean: " <<labelStatisticsImageFilter->GetMean(labelValue)<< std::endl;   //-
            std::cout << "\tmean: " <<labelMeanValue[(int)*vIt]<< std::endl;   //+

        }
    }

    /*  printing the sorted means of labels
     * int len=sizeof(labelMeanValue)/sizeof(labelMeanValue[0]);

    for(int i=0;i<len;++i){                 //0         1           2           3           4
        std::cout<<labelMeanValue[i]<<"\t";  //0.278904	0.115523	0.496039	0.121493	0.43653
    }*/

    ScalarToArrayFilterType::Pointer scalerToArrayFilter=ScalarToArrayFilterType::New();
    scalerToArrayFilter->SetInput(image);
    scalerToArrayFilter->Update();

    /*
     * With the input image type ImageType and labeled image type CharImageType we instantiate the
       type of the itk::MRFImageFilter that will apply the Markov Random Field algorithm in order to
       refine the pixel classification
     */

    MRFFilterType::Pointer mrfFilter=MRFFilterType::New();
    mrfFilter->SetInput(scalerToArrayFilter->GetOutput());

    /*
     * We set now some of the parameters for the MRF filter. In particular, the number of classes to be
       used during the classification, the maximum number of iterations to be run in this filter and the error
       tolerance that will be used as a criterion for convergence.
     */
    int numberOfIterations=1000;

    mrfFilter->SetNumberOfClasses(numberOfClasses);
    mrfFilter->SetMaximumNumberOfIterations(numberOfIterations);
    mrfFilter->SetErrorTolerance(1e-7);
    // The smoothing factor represents the tradeoff between fidelity to the observed image and the smoothness
    // of the segmented image. Typical smoothing factors have values between 1-5
    float smoothingFactor=2.5;
    mrfFilter->SetSmoothingFactor(smoothingFactor);

    SupervisedClassifierType::Pointer classifier=SupervisedClassifierType::New();

    /*
     *  The classifier need a decision rule to be set by the user. Note that we must use GetPointer() in the
        call of the SetDecisionRule() method because we are passing a SmartPointer, and smart pointer
        cannot perform polymorphism, we must then extract the raw pointer that is associated to the smart
        pointer. This extraction is done with the GetPointer() method.
     */
    typedef itk::Statistics::MinimumDecisionRule DecisionRuleType;
    DecisionRuleType::Pointer classificationDecisionRule=DecisionRuleType::New();
    classifier->SetDecisionRule(classificationDecisionRule);
    //classifier->SetClassifiedImage(prediction); //+



    // Instantiating membership function
    typedef itk::Statistics::DistanceToCentroidMembershipFunction<ArrayPixelType> MembershipFunctionType;
    typedef MembershipFunctionType::Pointer MembershipFunctionPointer;

    double meanDistance=0;

    MembershipFunctionType::CentroidType centroid(1);

    for (unsigned int i=0;i<numberOfClasses;i++) {
        MembershipFunctionPointer membershipFunction = MembershipFunctionType::New();

        centroid[0] = labelMeanValue[i];
        membershipFunction->SetCentroid(centroid);
        classifier->AddMembershipFunction(membershipFunction);

        meanDistance += static_cast<double>(centroid[0]);
    }
    if (numberOfClasses>0)
    {
        meanDistance/=numberOfClasses;
    }
    else
    {
        std::cerr <<"Error: numberOfClasses is 0"<<std::endl;

    }

    //
    // we set the neighborhoodradius that will define the size of the clique
    //
    mrfFilter->SetNeighborhoodRadius(1);


    //Start with below slice, setting it as the example mentions
    /*
     * As far as I understand it, it goes along the rows from the top or bottom slice upwards/downwards…
     * The weights have to be the same size as the neighbourhood. For a 3d image with a radius of 1 (3x3x3) It needs 9x3 weights
     * Link: https://discourse.itk.org/t/solved-copying-headers-after-performing-operation/1111/3
     */
    std::vector<double>weights;    //The array is arranged and then passed to the filter by using the method SetMRFNeighborhoodWeight().
    weights.push_back(1.3);
    weights.push_back(1.3);
    weights.push_back(1.3);
    weights.push_back(1.3);
    weights.push_back(1.5); //center
    weights.push_back(1.3);
    weights.push_back(1.3);
    weights.push_back(1.3);
    weights.push_back(1.3);

    weights.push_back(1.7);
    weights.push_back(1.7);
    weights.push_back(1.7);
    weights.push_back(1.7);
    weights.push_back(1.7);
    weights.push_back(1.7);
    weights.push_back(1.7);
    weights.push_back(1.7);
    weights.push_back(1.7);

    weights.push_back(1.3);
    weights.push_back(1.3);
    weights.push_back(1.3);
    weights.push_back(1.3);
    weights.push_back(1.5);
    weights.push_back(1.3);
    weights.push_back(1.3);
    weights.push_back(1.3);
    weights.push_back(1.3);

    /*
     * We now scale weights so that the smoothing function and the image fidelity functions have comparable value.
     * This is necessary since the label image and the input image can have different dynamic ranges.
     * The fidelity function is usually computed using a distance function, such as the
    itk::DistanceToCentroidMembershipFunction or one of the other membership functions. They
        tend to have values in the order of the means specified
     */

    double totalWeight=0;
    for(std::vector<double>::const_iterator wcIt=weights.begin();
        wcIt!=weights.end(); ++wcIt)
    {
        totalWeight+=*wcIt;
    }
    for(std::vector<double>::iterator wIt=weights.begin();
        wIt!=weights.end();++wIt)
    {
        *wIt=static_cast<double>((*wIt)*meanDistance/(2*totalWeight));
    }
    mrfFilter->SetMRFNeighborhoodWeight(weights);
    mrfFilter->SetClassifier(classifier);
    mrfFilter->Update();

    //
    // The output image produced by the \doxygen{MRFImageFilter} has the same pixel
    // type as the labeled input image. In the following lines we use the
    // \code{OutputImageType} in order to instantiate the type of a
    // \doxygen{ImageFileWriter}. Then create one, and connect it to the output of
    // the classification filter after passing it through an intensity rescaler
    // to rescale it to an 8 bit dynamic range
    //
    using OutputImageType=MRFFilterType::OutputImageType ;
    using RescaledOutputImageType=itk::Image<unsigned char, Dimension>;   //Similar to the CharImageType
    typedef itk::RescaleIntensityImageFilter<RescaledOutputImageType ,OutputImageType > RescalerType;

    RescalerType::Pointer intensityRescaler=RescalerType::New();
    intensityRescaler->SetOutputMinimum(0);
    intensityRescaler->SetOutputMaximum(255);
    intensityRescaler->SetInput(mrfFilter->GetOutput());

    // Software Guide : BeginCodeSnippet

    using WriterType = itk::ImageFileWriter<OutputImageType>;

    WriterType::Pointer writer = WriterType::New();

    writer->SetInput(intensityRescaler->GetOutput());

    writer->SetFileName("MRFOutput.nii.gz");
    // Software Guide : EndCodeSnippet


    // Software Guide : BeginLatex
    //
    // We are now ready for triggering the execution of the pipeline. This is done
    // by simply invoking the \code{Update()} method in the writer. This call will
    // propagate the update request to the reader and then to the MRF filter.
    //
    // Software Guide : EndLatex


    writer->Update();


    std::cout << "Number of Iterations : ";
    std::cout << mrfFilter->GetNumberOfIterations() << std::endl;
    std::cout << "Stop condition: " << std::endl;
    std::cout << "  (1) Maximum number of iterations " << std::endl;
    std::cout << "  (2) Error tolerance:  " << std::endl;
    std::cout << mrfFilter->GetStopCondition() << std::endl;

    return intensityRescaler->GetOutput();



}// end of ApplyMRF function
10

Hi @Sara_Caffe,

Inspecting the code, a few items stick out:

  • This is commented: //classifier->SetClassifiedImage(prediction);
  • It may help to reduce the number of iterations and smoothing factor to understand what is happening.
1 Like
11

@matt.mccormick Hi and thank you for your reply.

Once I un-commented the
classifier->SetClassifiedImage(prediction); //+
classifier->Update(); //+
It is showing the following error:

terminate called after throwing an instance of 'itk::ExceptionObject'
  what():  /usr/local/include/ITK-5.0/itkClassifierBase.hxx:67:
itk::ERROR: ImageClassifierBase(0x17adc50): Zero class

Then I added the classifier->SetNumberOfClasses(numberOfClasses); and the output is still blank image. Once I add the classifier->Update(); for updating the classifier, it still shows an error in the output:

WARNING: In /usr/local/include/ITK-5.0/itkImageSink.hxx, line 88
LabelStatisticsImageFilter (0x20fed30): Unable to convert input "LabelInput" to type N3itk5ImageIfLj3EEE

Number of labels: 5
Label: 4
	mean: 0.43653
Label: 0
	mean: 0.278904
Label: 1
	mean: 0.115523
Label: 2
	mean: 0.496039
Label: 3
	mean: 0.121493

Process finished with exit code 139 (interrupted by signal 11: SIGSEGV)

I printed out the MRF and the output is:

MRFImageFilter (0xfe4020)
  RTTI typeinfo:   itk::MRFImageFilter<itk::Image<itk::FixedArray<unsigned char, 1u>, 3u>, itk::Image<unsigned char, 3u> >
  Reference Count: 1
  Modified Time: 544
  Debug: Off
  Object Name: 
  Observers: 
    none
  Inputs: 
    Primary: (0xfeaf80) *
  Indexed Inputs: 
    0: Primary (0xfeaf80)
  Required Input Names: Primary
  NumberOfRequiredInputs: 1
  Outputs: 
    Primary: (0xfea1e0)
  Indexed Outputs: 
    0: Primary (0xfea1e0)
  NumberOfRequiredOutputs: 1
  Number Of Work Units: 32
  ReleaseDataFlag: Off
  ReleaseDataBeforeUpdateFlag: Off
  AbortGenerateData: Off
  Progress: 0
  Multithreader: 
    RTTI typeinfo:   itk::PoolMultiThreader
    Reference Count: 1
    Modified Time: 528
    Debug: Off
    Object Name: 
    Observers: 
      none
    Number of Work Units: 32
    Number of Threads: 8
    Global Maximum Number Of Threads: 128
    Global Default Number Of Threads: 8
    Global Default Threader Type: PoolMultiThreader
    SingleMethod: 0
    SingleData: 0
  DynamicMultiThreading: On
  CoordinateTolerance: 1e-06
  DirectionTolerance: 1e-06
   MRF Image filter object 
   Number of classes: 5
   Maximum number of iterations: 100
   Error tolerance for convergence: 1e-07
   Size of the MRF neighborhood radius:[1, 1, 1]
   Number of elements in MRF neighborhood :27
   Neighborhood weight : (0.00385276, 0.00385276, 0.00385276, 0.00385276, 0.00444549, 0.00385276, 0.00385276, 0.00385276, 0.00385276, 0.00503822, 0.00503822, 0.00503822, 0.00503822, 0.00503822, 0.00503822, 0.00503822, 0.00503822, 0.00503822, 0.00385276, 0.00385276, 0.00385276, 0.00385276, 0.00444549, 0.00385276, 0.00385276, 0.00385276, 0.00385276)
   Smoothing factor for the MRF neighborhood:0.8

Then, I checked the minimum and maximum value of intensity rescaled mrfFilter output image.

    ImageMinMaxCalculatorFilterType::Pointer imageMinMaxCalculator=ImageMinMaxCalculatorFilterType::New();
    imageMinMaxCalculator->SetImage(intensityRescaler->GetOutput());
    imageMinMaxCalculator->Compute();
    std::cout<<"Min value of image:  "<<(int)imageMinMaxCalculator->GetMinimum()<<" and Max: "<<(int)imageMinMaxCalculator->GetMaximum()<<std::endl;

Min value of image: 255 and Max: 0, in which the setting for intensityRescaler is the opposite. I am wondering what is happening here.

shows:

StopCondition: 1
   Number of iterations: 0
Min value of image:  255 and Max: 0
Number of Iterations : 1
Stop condition: 
  (1) Maximum number of iterations 
  (2) Error tolerance:  
2
12

You might be using ITK as shared libraries (DLLs), and some version is mismatched, or some other linking issue is manifesting itself.

2 Likes
13

:+1: or the LabelStatisticImageFilter template arguments are not correct.

1 Like
14

@dzenanz thank you for mentioning the point. I do not have access to my device at the moment due to the pandemic. So I chose a lazy way and changes the unsigned char type of label images to the float. Could I ask do you mean that two different versions of ITK are available on my device? How can I resolve this?
Sorry for basic questions; I am not a professional programmer.

15

@matt.mccormick What I did was to convert the unsigned char CharType image to float. The output of MRF is like this for all 5 labels (including 4 objects background as 0 label):

Then I removed the other labels, calculate the mean for two classes, and the input to the mrfFilter with binary values (0: background, 1: object of interest that the holes need to be filled) gives this:

  • Is this the correct output of MRF or I need to change some parameters?

  • Should I feed the probability map (obtained from my algorithm) to the MRF?

  • Why MRF has not filled the holes in my object of interest, which is the label: 1? should I give the input with all labels of only the image with object that needs to be corrected by MRF?

Thanks

16

Looks good, but adjusting parameters may help.

Yes, isolated that label that you want to filter may help.

Also of interest is itk::ClosingByReconstructionImageFilter.

1 Like