The program seems suspended (or halted) when it reached to the line Updating the Simplex filter

Hi everyone,

I followed all the comments here, and I am trying to apply the deformation as a post-processing step for the segmentation results. I have applied this algorithm to my data. I have the following files:

  1. a mesh created from the segmentation result of an algorithm (mesh created by Cuberille filter)
  2. an intensity image with normalized values between [0,1].

I am wondering how long the processing takes time, because the algorithm, after almost 30 minutes, still has not translated the mesh to the origin [0,0,0]. I followed the same code, however, I do not know why when it reaches to the following line to convert the triangle mesh to simplex mesh, it halts the program without any output, then I have to stop running since nothing is shown and keeps showing as running.

TriangleToSimplexFilterType::Pointer cast=TriangleToSimplexFilterType::New();
cast->SetInput(triMesh);;
cast->Update();

This is the function that I have defined my main file:

void ApplyDeformableModel(ImageType::Pointer  image,TriangleMeshType::Pointer triMesh){
    /*
     * This function takes:
     *              1) the intensity image and triangle mesh as input,
     *              2) it converts the triangle mesh into Simplex mesh,
     *              3) applies 3D Simplex Deformable model
     * Link:
     *      https://discourse.itk.org/t/simplex-mesh-deformable-model-doesnt-change/386
     */

    float dspl_x, dspl_y, dspl_z;

    dspl_x=image->GetOrigin()[0];
    dspl_y=image->GetOrigin()[1];
    dspl_z=image->GetOrigin()[2];
    std::cout<<"Displacement x: "<<dspl_x<<"  Displacement y: "<<dspl_y<<"  Displacement z: "<<dspl_z<<std::endl;

    MeshTransformType::Pointer translation=MeshTransformType::New();
    MeshTransformType::OutputVectorType displacement;

    displacement[0]=-dspl_x;
    displacement[1]=-dspl_y;
    displacement[2]=-dspl_z;

    translation->Translate(displacement);

    MeshFilterType::Pointer filter=MeshFilterType::New();
    filter->SetInput(triMesh);
    filter->SetTransform(translation);

    triMesh->DisconnectPipeline();  //+
    try{
        filter->Update();
    }
    catch(itk::ExceptionObject& exp){
        std::cerr<<"Exception thrown during Update() "<<std::endl;
        std::cerr<<exp<<std::endl;
    }
    TriangleMeshType::Pointer translatedMesh=filter->GetOutput();
    translatedMesh->DisconnectPipeline();

    MeshWriterType::Pointer meshWriter=MeshWriterType::New();
    meshWriter->SetFileName("translatedMesh.vtk");
    meshWriter->SetInput(translatedMesh);//conversionResult);
    meshWriter->Update();



    TriangleToSimplexFilterType::Pointer cast=TriangleToSimplexFilterType::New();
    cast->SetInput(triMesh);;
    cast->Update();

    SimplexMeshType::Pointer simplexMesh=cast->GetOutput();
    simplexMesh->DisconnectPipeline();
    std::cout<<" Simplex Mesh was Created !"<<std::endl;

    // Correct Image Origin
    std::cout<<"Changing image origin "<<std::endl;

    ImageType::PointType newOrigin;
    newOrigin[0]=0;
    newOrigin[1]=0;
    newOrigin[2]=0;

    typedef itk::ChangeInformationImageFilter<ImageType> ChangeImageInfoType;
    ChangeImageInfoType::Pointer changeInfo=ChangeImageInfoType::New();
    changeInfo->SetInput(image);
    changeInfo->SetOutputOrigin(newOrigin);
    changeInfo->ChangeOriginOn();
    try{
        changeInfo->Update();
    }
    catch(itk::ExceptionObject&exp){
        std::cerr <<"exception thrown during Update() "<<std::endl;
        std::cerr <<exp<<std::endl;
    }

    image=changeInfo->GetOutput();

    GradientAnisotropicImageType::Pointer gradientAnisotropicFilter=GradientAnisotropicImageType::New();
    gradientAnisotropicFilter->SetInput(image);
    gradientAnisotropicFilter->SetNumberOfIterations(5);
    gradientAnisotropicFilter->SetTimeStep(0.0625);
    gradientAnisotropicFilter->SetConductanceParameter(3);
    gradientAnisotropicFilter->Update();
    std::cout<< "GradientAnisotropicDiffusion is DONE! "<<std::endl;



    GradientFilterType::Pointer gradientFilter=GradientFilterType::New();
    gradientFilter->SetInput(gradientAnisotropicFilter->GetOutput() );
    //gradientFilter->SetInput(sigmoidImageFilter->GetOutput());
    gradientFilter->SetSigma(1.0);
    gradientFilter->Update();
    std::cout<<"Gradient Filter DDNE!"<< std::endl;



    GradientType::Pointer gradientImage=gradientFilter->GetOutput();

    DeformFilterType::Pointer deformFilter=DeformFilterType::New();

    const unsigned int numberOfCycles=100;
    for(unsigned int i=0;i<numberOfCycles;i++){
        // Must disconnect the pipline

        deformFilter->SetInput(simplexMesh);
        deformFilter->SetGradient(gradientImage);
        deformFilter->SetAlpha(0.1);
        deformFilter->SetBeta(0.9);
        deformFilter->SetIterations(500);
        deformFilter->SetRigidity(1);
        deformFilter->Update();
    }
    SimplexMeshType::Pointer deformResult=deformFilter->GetOutput();

    SimplexToTriangleFilterType::Pointer simplexToTriangleFilter=SimplexToTriangleFilterType::New();
    simplexToTriangleFilter->SetInput(deformResult);
    simplexToTriangleFilter->Update();
    TriangleMeshType::Pointer conversionResult=simplexToTriangleFilter->GetOutput();

    conversionResult->DisconnectPipeline();

    // return the coordinates back
    MeshTransformType::Pointer returnTranslation=MeshTransformType::New();
    MeshTransformType::OutputVectorType returnDisplacement;

    returnDisplacement[0]=dspl_x;
    returnDisplacement[1]=dspl_y;
    returnDisplacement[2]=dspl_z;

    returnTranslation->Translate(returnDisplacement);

    MeshFilterType::Pointer filter2=MeshFilterType::New();
    filter2->SetInput(conversionResult);
    filter2->SetTransform(returnTranslation);
    try{
        filter2->Update();
    }
    catch(itk::ExceptionObject & exp){
        std::cerr<<"Exception thrown during Update() " <<std::endl;
        std::cerr<<exp<<std::endl;
    }

    TriangleMeshType::Pointer outputMesh=TriangleMeshType::New();
    outputMesh=filter2->GetOutput();
    outputMesh->DisconnectPipeline();


    meshWriter->SetFileName("deformedMesh.vtk");
    meshWriter->SetInput(outputMesh);//conversionResult);
    meshWriter->Update();

}

I have attached the inputs to the function and the translated mesh.

Your expernt opinion is really appreciated

Your problem is probably that your mesh is not a 2-manifold (it is not a closed mesh, like representing surface of a 3D object), and therefore it cannot be converted into a simplex mesh. Trying cleaning up the mesh and closing it.

@dzenanz thanks for responding, how a mesh can be cleaned up and closed in ITK?

Pad the image with 0 or -1024, to make sure your surface is closed. As for the noise inside the object, you can try fill holes. Both need to be applied before converting to mesh.

1 Like

@dzenanz Thanks for your help. I wrote a function to clean the data as follows:
the input to the function is the output of largest connected component filter:

using LargestCCFilterType=itk::LargestConnectedComponentFilter<ImageType ,ImageType >; 
LargestCCFilterType::Pointer largestCCFilter=LargestCCFilterType::New();
        largestCCFilter->SetInputBackgroundValue(0);
        largestCCFilter->SetOutputForegroundValue(1);
        largestCCFilter->SetInput(pred);
       pred=largestCCFilter->GetOutput();

pred is the ImageType defined with float PixelType and dimention of 3, and this is the function:

ImageType::Pointer PrepareData(ImageType::Pointer pred, int padSize=3){
    /*
     * This function prepares the prediction file for surface extraction.
     * It cleans up the data for mesh generation
     *  1) padding
     *  Link: https://itk.org/Doxygen/html/classitk_1_1ConstantPadImageFilter.html
     *
     *  2) filling the wholes
     *  Link: https://itk.org/Doxygen/html/classitk_1_1BinaryFillholeImageFilter.html
     */
    PadFilterType::Pointer padFilter=PadFilterType::New();
    FillholeFilterType::Pointer fillHoleFilter=FillholeFilterType::New();

    ImageType::SizeType lowerExtendRegion;
    lowerExtendRegion.Fill(padSize);    // pad from up-left side    lowerBound[0]: pad from left side, lowerbound[1]: pad from up

    ImageType::SizeType upperExtendRegion;
    upperExtendRegion.Fill(padSize);      // pad from the bottom side, upperBound[0]: pad from right side, lowerbound[1]: pad from bottom
    padFilter->SetInput(pred);
    padFilter->SetPadLowerBound(lowerExtendRegion);
    padFilter->SetPadUpperBound(upperExtendRegion);
    padFilter->SetConstant(0);   // fill with zero values, default=0

    fillHoleFilter->SetInput(padFilter->GetOutput());
    fillHoleFilter->SetForegroundValue(itk::NumericTraits<PixelType>::max());


    ImageType::Pointer outputImage=fillHoleFilter->GetOutput();

    //==== calculate the minimum and maximum value of the image
    using ImageCalculatorFilterType = itk::MinimumMaximumImageCalculator<ImageType>;

    ImageCalculatorFilterType::Pointer imageCalculatorFilter = ImageCalculatorFilterType::New();
    imageCalculatorFilter->SetImage(outputImage);
    imageCalculatorFilter->Compute();
    std::cout<<"Min value of image:  "<<imageCalculatorFilter->GetMinimum()<<" and Max: "<<imageCalculatorFilter->GetMaximum()<<std::endl;
 ImageType::RegionType region = pred->GetLargestPossibleRegion();   //ImageType::RegionType inputRegion = image->GetBufferedRegion();
    ImageType::SizeType size = region.GetSize();
    ImageType::IndexType start = region.GetIndex();
    cout << size << endl;
    cout << start << endl;
    cout << "===========================================" << endl;


    return outputImage;


} // End PrepareData

However, there is an error:

Min value of image: 3.40282e+38 and Max: -3.40282e+38
[0, 0, 0]
[0, 0, 0]
===========================================
terminate called after throwing an instance of 'itk::ExceptionObject’
** what(): /usr/local/include/ITK-5.0/itkConstNeighborhoodIterator.h:304:**
In method IsAtEnd, CenterPointer = 0x34 is greater than End = 0
** Neighborhood:**
** Radius:[1, 1, 1]**
** Size:[3, 3, 3]**
** DataBuffer:NeighborhoodAllocator { this = 0x7ffe17e772e8, begin = 0x3d97800, size=27 }**

I am struggling for whole day, unfortuantely could not solve the issue and did not find out what could be the reason.
The min and maximum in the input file is 0 and 1, however, it is showing different value, and the size of image is shown as [0, 0, 0].

I changed the information pred input image as the information of original image while reading it into C++, however, I am still getting the same error, should I change the information of image after padding? and what is the reason the values of min and max is different.

This is the pred file after taking the largestCC and it is input to the funtion after changing information as original input:
pred_input_origin.nii.gz (68.0 KB)

I don’t think you need LargestConnectedComponentFilter, because you don’t have islands. It looks like you have holes in your main component, so you need fill hole filter.

You are missing fillHoleFilter->Update();, that’s why you get Min value of image: 3.40282e+38 and Max: -3.40282e+38.