Simplex Mesh Deformable Model doesn't change

Dear itk-developers,

I’m trying to use the itkDeformableSimplexMesh3D filter to make a local adjustment of a simplex mesh. The mesh that I want to deform using this filter is already initialized close to the target structure in the image. However, after playing with all parameter values of the filter, it seems that there is not significant differences between the input and output meshes.

I’m reading the mesh from a file (stl or vtk extension) and I want to deform it using a brain MRI (nifiti format) image as a reference. After the filter application, the result seems to not approximate the mesh to the local structure in the image (only occurs a small smoothness of the model). I’ve tried to change the parameters (alpha, beta, gamma, rigidity, number of iterations) values to check if some change would occur, but the result always seems the same. I would like that deformation to be able to move the mesh to closer to the structure in the image. My code is below. Is there something wrong with my code? Why it seems that the external force has no influence in the mesh deform result?

The image and the mesh I’m using in this test is available in the link below:

https://drive.google.com/drive/folders/0B3FtwZ6urVW3X01HbG5kQTgwdVE?usp=sharing

Here is my code:

=====================================================================
enum PARAMS {
PROGRAM_NAME = 0,
IMAGE_FILENAME,
INPUT_MESH_FILENAME,
OUTPUT_MESH_FILENAME,
NUM_ARGS
};

using namespace std;

int main( int argc, char *argv[] )
{
if( argc < NUM_ARGS )
{
std::cerr << “\n\n”;
std::cerr << “Usage: " << argv[PROGRAM_NAME] << " imageFileName inputMeshFileName ouputMeshFileName”;
std::cerr << “\n\n”;
return 1;
}

itk::STLMeshIOFactory::RegisterOneFactory();
const int dimension = 3;

typedef float PixelType;
typedef itk::Image<PixelType, dimension> ImageType;

typedef itk::DefaultDynamicMeshTraits<PixelType, dimension> TriangleMeshTraits;
typedef itk::DefaultDynamicMeshTraits<PixelType, dimension> SimplexMeshTraits;
typedef itk::Mesh<PixelType, dimension, TriangleMeshTraits> TriangleMeshType;
typedef itk::SimplexMesh<PixelType, dimension, SimplexMeshTraits> SimplexMeshType;

typedef itk::TriangleMeshToSimplexMeshFilter<TriangleMeshType, SimplexMeshType> TriangleToSimplexFilterType;
typedef itk::SimplexMeshToTriangleMeshFilter<SimplexMeshType, TriangleMeshType> SimplexToTriangleFilterType;

typedef itk::DeformableSimplexMesh3DFilter<SimplexMeshType, SimplexMeshType> DeformFilterType;
typedef DeformFilterType::GradientImageType GradientType;

typedef itk::GradientAnisotropicDiffusionImageFilter <ImageType, ImageType> DiffusionFilterType;

typedef itk::GradientAnisotropicDiffusionImageFilter < ImageType, ImageType > GradientAnisotropicImageType;
typedef itk::GradientMagnitudeRecursiveGaussianImageFilter < ImageType, ImageType > GradientMagnitudeType;
typedef itk::SigmoidImageFilter< ImageType, ImageType > SigmoidImageType;
typedef itk::GradientRecursiveGaussianImageFilter<ImageType,GradientType> GradientFilterType;

typedef itk::ImageFileReader ImageReaderType;
typedef itk::ImageFileWriter ImageWriterType;
typedef itk::MeshFileReader MeshReaderType;
typedef itk::MeshFileWriter MeshWriterType;

ImageReaderType::Pointer imageReader = ImageReaderType::New();
imageReader->SetFileName(argv[IMAGE_FILENAME]);
imageReader->Update();

MeshReaderType::Pointer meshReader = MeshReaderType::New();
meshReader->SetFileName(argv[INPUT_MESH_FILENAME]);
meshReader->Update();

std::cout << " starting to Filter Image" << std::endl;
GradientAnisotropicImageType::Pointer gradientanisotropicfilter = GradientAnisotropicImageType::New();
gradientanisotropicfilter->SetInput(imageReader->GetOutput());
gradientanisotropicfilter->SetNumberOfIterations(5);
gradientanisotropicfilter->SetTimeStep(0.0625);
gradientanisotropicfilter->SetConductanceParameter(3);
gradientanisotropicfilter->Update();
std::cout << “GradientAnisotropicDiffusion is DONE!” << std::endl;

GradientMagnitudeType::Pointer gradientmagnitudefilter = GradientMagnitudeType::New();
gradientmagnitudefilter->SetInput( gradientanisotropicfilter->GetOutput() );
gradientmagnitudefilter->SetSigma(1.0);
gradientmagnitudefilter->Update();
std::cout << “GradientMagnitude is DONE!” << std::endl;

SigmoidImageType::Pointer sigmoidimagefilter = SigmoidImageType::New();
sigmoidimagefilter->SetInput( gradientmagnitudefilter->GetOutput());
sigmoidimagefilter->SetOutputMinimum(0);
sigmoidimagefilter->SetOutputMaximum(1);
sigmoidimagefilter->SetAlpha(10);
sigmoidimagefilter->SetBeta(100);
sigmoidimagefilter->Update();
std::cout << “Sigmoid is DONE!” << std::endl;

GradientFilterType::Pointer gradientFilter = GradientFilterType::New();
gradientFilter->SetInput( sigmoidimagefilter->GetOutput() );
gradientFilter->SetSigma(1.0);
gradientFilter->Update();
std::cout << “GradientMagnitude is DONE!” << std::endl;

std::cout << “Gradient filter ok\n”;

TriangleToSimplexFilterType::Pointer triangleToSimplexFilter = TriangleToSimplexFilterType::New();
triangleToSimplexFilter->SetInput(meshReader->GetOutput());
triangleToSimplexFilter->Update();

GradientType::Pointer gradientImage = gradientFilter->GetOutput();
SimplexMeshType::Pointer simplexMesh = triangleToSimplexFilter->GetOutput();
DeformFilterType::Pointer deformFilter = DeformFilterType::New();

const unsigned int numberOfCycles = 100;

for (unsigned int i = 0; i < numberOfCycles; i++)
{
// must disconnect the pipeline
simplexMesh->DisconnectPipeline();
deformFilter->SetInput( simplexMesh );
deformFilter->SetGradient( gradientImage );
deformFilter->SetAlpha(0.1);
deformFilter->SetBeta(-0.1);
deformFilter->SetIterations(5);
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();

MeshWriterType::Pointer meshWriter = MeshWriterType::New();
meshWriter->SetFileName(argv[OUTPUT_MESH_FILENAME]);
meshWriter->SetInput(conversionResult);
meshWriter->Update();

return EXIT_SUCCESS;
}

=====================================================================

I really appreciate any help you can provide me.

Thank you very much,
Breno

Which version of ITK are you using Breno? The program you supplied (after minor modifications to make it compile, given below) crashes with the provided inputs, either during reading the mesh (debug mode) or writing it (release mode).

#include <iostream>
#include <itkDeformableSimplexMesh3DFilter.h>
#include <itkTriangleMeshToSimplexMeshFilter.h>
#include <itkSimplexMeshToTriangleMeshFilter.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkMeshFileReader.h>
#include <itkMeshFileWriter.h>
#include <itkSTLMeshIOFactory.h>
#include <itkDefaultDynamicMeshTraits.h>
#include <itkGradientAnisotropicDiffusionImageFilter.h>
#include <itkGradientMagnitudeRecursiveGaussianImageFilter.h>
#include <itkSigmoidImageFilter.h>
#include <itkGradientRecursiveGaussianImageFilter.h>

enum PARAMS {
    PROGRAM_NAME = 0,
    IMAGE_FILENAME,
    INPUT_MESH_FILENAME,
    OUTPUT_MESH_FILENAME,
    NUM_ARGS
};

using namespace std;

int main(int argc, char *argv[])
{
    if (argc < NUM_ARGS)
    {
        std::cerr << "\n\n";
        std::cerr << "Usage: " << argv[PROGRAM_NAME] << " imageFileName inputMeshFileName ouputMeshFileName";
        std::cerr << "\n\n";
        return 1;
    }

    itk::STLMeshIOFactory::RegisterOneFactory();
    const int dimension = 3;

    typedef float PixelType;
    typedef itk::Image<PixelType, dimension> ImageType;

    typedef itk::DefaultDynamicMeshTraits<PixelType, dimension> TriangleMeshTraits;
    typedef itk::DefaultDynamicMeshTraits<PixelType, dimension> SimplexMeshTraits;
    typedef itk::Mesh<PixelType, dimension, TriangleMeshTraits> TriangleMeshType;
    typedef itk::SimplexMesh<PixelType, dimension, SimplexMeshTraits> SimplexMeshType;

    typedef itk::TriangleMeshToSimplexMeshFilter<TriangleMeshType, SimplexMeshType> TriangleToSimplexFilterType;
    typedef itk::SimplexMeshToTriangleMeshFilter<SimplexMeshType, TriangleMeshType> SimplexToTriangleFilterType;

    typedef itk::DeformableSimplexMesh3DFilter<SimplexMeshType, SimplexMeshType> DeformFilterType;
    typedef DeformFilterType::GradientImageType GradientType;

    typedef itk::GradientAnisotropicDiffusionImageFilter <ImageType, ImageType> DiffusionFilterType;

    typedef itk::GradientAnisotropicDiffusionImageFilter < ImageType, ImageType > GradientAnisotropicImageType;
    typedef itk::GradientMagnitudeRecursiveGaussianImageFilter < ImageType, ImageType > GradientMagnitudeType;
    typedef itk::SigmoidImageFilter< ImageType, ImageType > SigmoidImageType;
    typedef itk::GradientRecursiveGaussianImageFilter<ImageType, GradientType> GradientFilterType;

    typedef itk::ImageFileReader<ImageType> ImageReaderType;
    typedef itk::ImageFileWriter<ImageType> ImageWriterType;
    typedef itk::MeshFileReader<SimplexMeshType> MeshReaderType;
    typedef itk::MeshFileWriter<TriangleMeshType> MeshWriterType;

    ImageReaderType::Pointer imageReader = ImageReaderType::New();
    imageReader->SetFileName(argv[IMAGE_FILENAME]);
    imageReader->Update();

    MeshReaderType::Pointer meshReader = MeshReaderType::New();
    meshReader->SetFileName(argv[INPUT_MESH_FILENAME]);
    meshReader->Update();

    std::cout << " starting to Filter Image" << std::endl;
    GradientAnisotropicImageType::Pointer gradientanisotropicfilter = GradientAnisotropicImageType::New();
    gradientanisotropicfilter->SetInput(imageReader->GetOutput());
    gradientanisotropicfilter->SetNumberOfIterations(5);
    gradientanisotropicfilter->SetTimeStep(0.0625);
    gradientanisotropicfilter->SetConductanceParameter(3);
    gradientanisotropicfilter->Update();
    std::cout << "GradientAnisotropicDiffusion is DONE!" << std::endl;

    GradientMagnitudeType::Pointer gradientmagnitudefilter = GradientMagnitudeType::New();
    gradientmagnitudefilter->SetInput(gradientanisotropicfilter->GetOutput());
    gradientmagnitudefilter->SetSigma(1.0);
    gradientmagnitudefilter->Update();
    std::cout << "GradientMagnitude is DONE!" << std::endl;

    SigmoidImageType::Pointer sigmoidimagefilter = SigmoidImageType::New();
    sigmoidimagefilter->SetInput(gradientmagnitudefilter->GetOutput());
    sigmoidimagefilter->SetOutputMinimum(0);
    sigmoidimagefilter->SetOutputMaximum(1);
    sigmoidimagefilter->SetAlpha(10);
    sigmoidimagefilter->SetBeta(100);
    sigmoidimagefilter->Update();
    std::cout << "Sigmoid is DONE!" << std::endl;

    GradientFilterType::Pointer gradientFilter = GradientFilterType::New();
    gradientFilter->SetInput(sigmoidimagefilter->GetOutput());
    gradientFilter->SetSigma(1.0);
    gradientFilter->Update();
    std::cout << "GradientMagnitude is DONE!" << std::endl;

    std::cout << "Gradient filter ok\n";

    TriangleToSimplexFilterType::Pointer triangleToSimplexFilter = TriangleToSimplexFilterType::New();
    triangleToSimplexFilter->SetInput(meshReader->GetOutput());
    triangleToSimplexFilter->Update();

    GradientType::Pointer gradientImage = gradientFilter->GetOutput();
    SimplexMeshType::Pointer simplexMesh = triangleToSimplexFilter->GetOutput();
    DeformFilterType::Pointer deformFilter = DeformFilterType::New();

    const unsigned int numberOfCycles = 100;

    for (unsigned int i = 0; i < numberOfCycles; i++)
    {
        // must disconnect the pipeline
        simplexMesh->DisconnectPipeline();
        deformFilter->SetInput(simplexMesh);
        deformFilter->SetGradient(gradientImage);
        deformFilter->SetAlpha(0.1);
        deformFilter->SetBeta(-0.1);
        deformFilter->SetIterations(5);
        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();

    MeshWriterType::Pointer meshWriter = MeshWriterType::New();
    meshWriter->SetFileName(argv[OUTPUT_MESH_FILENAME]);
    meshWriter->SetInput(conversionResult);
    meshWriter->Update();

    return EXIT_SUCCESS;
}

Sorry. I did not specify the “includes” in the source code. I guess that the test is failing because of the missing line #include “itkSTLMeshIO.h”. Here is my entire code. Moreover, I am able to read a stl mesh, but I am only able to save in the vtk format.
I usually give that input specification:
./main input.nii input_mesh.stl output_mesh.vtk

And I am using itk 4.12.2, the latest stable release.


#include <iostream>


#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkSTLMeshIOFactory.h"
#include "itkSTLMeshIO.h"

#include "itkMeshFileReader.h"
#include "itkMeshFileWriter.h"
#include "itkDeformableSimplexMesh3DFilter.h"
#include "itkDefaultDynamicMeshTraits.h"
#include "itkMesh.h"
#include "itkSimplexMesh.h"
#include "itkTriangleMeshToSimplexMeshFilter.h"
#include "itkSimplexMeshToTriangleMeshFilter.h"

#include "itkGradientAnisotropicDiffusionImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkSigmoidImageFilter.h"
#include "itkGradientRecursiveGaussianImageFilter.h"


enum PARAMS {
    PROGRAM_NAME = 0,
    IMAGE_FILENAME,
    INPUT_MESH_FILENAME,
    OUTPUT_MESH_FILENAME,
    NUM_ARGS
};


using namespace std;

int main( int argc, char *argv[] )
{

	if( argc < NUM_ARGS )
	{
		std::cerr << "\n\n";
		std::cerr << "Usage: " << argv[PROGRAM_NAME] << " imageFileName inputMeshFileName ouputMeshFileName";
		std::cerr << "\n\n";
		return 1;
	}
	itk::STLMeshIOFactory::RegisterOneFactory();

  const int dimension = 3;

  typedef float PixelType;
  typedef itk::Image<PixelType, dimension> ImageType;

  typedef itk::DefaultDynamicMeshTraits<PixelType, dimension>       TriangleMeshTraits;
  typedef itk::DefaultDynamicMeshTraits<PixelType, dimension>       SimplexMeshTraits;
  typedef itk::Mesh<PixelType, dimension, TriangleMeshTraits>       TriangleMeshType;
  typedef itk::SimplexMesh<PixelType, dimension, SimplexMeshTraits> SimplexMeshType;

  typedef itk::TriangleMeshToSimplexMeshFilter<TriangleMeshType, SimplexMeshType> TriangleToSimplexFilterType;
  typedef itk::SimplexMeshToTriangleMeshFilter<SimplexMeshType, TriangleMeshType> SimplexToTriangleFilterType;

  typedef itk::DeformableSimplexMesh3DFilter<SimplexMeshType, SimplexMeshType> DeformFilterType;
  typedef DeformFilterType::GradientImageType                                  GradientType;

  typedef itk::GradientAnisotropicDiffusionImageFilter <ImageType, ImageType> DiffusionFilterType;

  typedef itk::GradientAnisotropicDiffusionImageFilter < ImageType, ImageType >  GradientAnisotropicImageType;
  typedef itk::GradientMagnitudeRecursiveGaussianImageFilter < ImageType, ImageType >  GradientMagnitudeType;
  typedef itk::SigmoidImageFilter< ImageType, ImageType > SigmoidImageType;
  typedef itk::GradientRecursiveGaussianImageFilter<ImageType,GradientType> GradientFilterType;

  typedef itk::ImageFileReader<ImageType>       ImageReaderType;
  typedef itk::ImageFileWriter<ImageType>       ImageWriterType;
  typedef itk::MeshFileReader<TriangleMeshType> MeshReaderType;
  typedef itk::MeshFileWriter<TriangleMeshType>  MeshWriterType;


  ImageReaderType::Pointer imageReader = ImageReaderType::New();
  imageReader->SetFileName(argv[IMAGE_FILENAME]);
  imageReader->Update();

  MeshReaderType::Pointer meshReader = MeshReaderType::New();
  meshReader->SetFileName(argv[INPUT_MESH_FILENAME]);
  meshReader->Update();

  std::cout << " starting to Filter Image" << std::endl;
  GradientAnisotropicImageType::Pointer gradientanisotropicfilter = GradientAnisotropicImageType::New();
  gradientanisotropicfilter->SetInput(imageReader->GetOutput());
  gradientanisotropicfilter->SetNumberOfIterations(5);
  gradientanisotropicfilter->SetTimeStep(0.0625);
  gradientanisotropicfilter->SetConductanceParameter(3);
  gradientanisotropicfilter->Update();
  std::cout << "GradientAnisotropicDiffusion is DONE!" << std::endl;

  GradientMagnitudeType::Pointer gradientmagnitudefilter = GradientMagnitudeType::New();
  gradientmagnitudefilter->SetInput( gradientanisotropicfilter->GetOutput() );
  gradientmagnitudefilter->SetSigma(1.0);
  gradientmagnitudefilter->Update();
  std::cout << "GradientMagnitude is DONE!" << std::endl;

  SigmoidImageType::Pointer sigmoidimagefilter = SigmoidImageType::New();
  sigmoidimagefilter->SetInput( gradientmagnitudefilter->GetOutput());
  sigmoidimagefilter->SetOutputMinimum(0);
  sigmoidimagefilter->SetOutputMaximum(1);
  sigmoidimagefilter->SetAlpha(10);
  sigmoidimagefilter->SetBeta(100);
  sigmoidimagefilter->Update();
  std::cout << "Sigmoid is DONE!" << std::endl;

  GradientFilterType::Pointer gradientFilter = GradientFilterType::New();
  gradientFilter->SetInput( sigmoidimagefilter->GetOutput() );
  gradientFilter->SetSigma(1.0);
  gradientFilter->Update();
  std::cout << "GradientMagnitude is DONE!" << std::endl;

  std::cout << "Gradient filter ok\n";

  TriangleToSimplexFilterType::Pointer triangleToSimplexFilter = TriangleToSimplexFilterType::New();
  triangleToSimplexFilter->SetInput(meshReader->GetOutput());
  triangleToSimplexFilter->Update();

  GradientType::Pointer gradientImage    = gradientFilter->GetOutput();
  SimplexMeshType::Pointer simplexMesh   = triangleToSimplexFilter->GetOutput();
  DeformFilterType::Pointer deformFilter = DeformFilterType::New();

  const unsigned int numberOfCycles = 100;

  for (unsigned int i = 0; i < numberOfCycles; i++)
    {
    // must disconnect the pipeline
    simplexMesh->DisconnectPipeline();
    deformFilter->SetInput( simplexMesh );
    deformFilter->SetGradient( gradientImage );
    deformFilter->SetAlpha(0.1);
    deformFilter->SetBeta(-0.1);
    deformFilter->SetIterations(5);
    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();

  MeshWriterType::Pointer meshWriter = MeshWriterType::New();
  meshWriter->SetFileName(argv[OUTPUT_MESH_FILENAME]);
  meshWriter->SetInput(conversionResult);
  meshWriter->Update();

  return EXIT_SUCCESS;
}


The error I was getting on writing was “Found Non-Triangular Cell.”, which was a cellType == VERTEX_CELL. Simplex mesh contains vertices, edges and triangles (if I remember correctly). I guess STL format only supports triangles, that’s why this can be only written into .vtk format.

I found the reason why the mesh doesn’t get deformed. In itkDeformableSimplexMesh3DFilter.hxx, line 433 there is this condition:

  if ( ( coord[0] >= 0 ) && ( coord[1] >= 0 ) && ( coord[2] >= 0 )
       && ( coord2[0] < m_ImageWidth ) && ( coord2[1] < m_ImageHeight ) && ( coord2[2] < m_ImageDepth ) )

It means this filter assumes 0 origin (whereas your image had -128 or something similar). With both image and input mesh translated so the origin is 0, the filter deforms the mesh (at least with SetBeta(0.95) and a few other modifications I made during debugging).

The filter should be modified to be more general. As a stop-gap measure, I submitted a patch to warn users of these severe limitations.

Thank you for your help. I have a doubt. How can I change the mesh origin to 0 and still ensure that it will be aligned to the same structure in the translated image? Is there a way to set the mesh origin accordingly to the image origin?

Once the mesh is loaded, you could iterate through the points, and add an offset to their coordinates.

For example, if the image’s origin is -10,-100,+50, you should add +10,+100,-50 to the image’s origin (to make it 0,0,0) and the same offset to all the points in the mesh to keep the mesh aligned to the same structure in the image.

After the processing is done, you could translate your mesh by -10,-100,+50 to get it back into original coordinate system.

1 Like

I’m having a little bit of difficult to change this coordinates on the image and on the mesh the same way. You still has the code you used? Could you please post here, it would help me a lot to understand the problem and your solution.

I didn’t write the code to do it, I manually edited .nrrd file and transformed the mesh using Slicer.

Hello again…

Thank you for the help. I’m doing the conversions inside the code, changing the image origin to (0,0,0) and aligning the mesh to the image through a translation (the offset was chosen observing the image and the mesh at paraview…these values made the mesh align with the image in the (0,0,0) origin when seen in paraview). After this, I translate the mesh back to the original position (Am I doing something wrong in the conversions?).

When I compare the meshes, still there is no much difference. Moreover, the surface texture seems to be disconnected to the lines and points in the result (is this normal?). Is the simplex mesh deform filter able to move the mesh, or the filter deformation occurs around the center of the mesh (only expanding/contracting the model)?

Here is my code:


#include <iostream>


#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkSTLMeshIOFactory.h"
#include "itkSTLMeshIO.h"

#include "itkMeshFileReader.h"
#include "itkMeshFileWriter.h"
#include "itkDeformableSimplexMesh3DFilter.h"
#include "itkDefaultDynamicMeshTraits.h"
#include "itkMesh.h"
#include "itkSimplexMesh.h"
#include "itkTriangleMeshToSimplexMeshFilter.h"
#include "itkSimplexMeshToTriangleMeshFilter.h"

#include "itkGradientAnisotropicDiffusionImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkSigmoidImageFilter.h"
#include "itkGradientRecursiveGaussianImageFilter.h"

#include "itkTranslationTransform.h"
#include "itkTransformMeshFilter.h"

#include "itkChangeInformationImageFilter.h"


enum PARAMS {
    PROGRAM_NAME = 0,
    IMAGE_FILENAME,
    INPUT_MESH_FILENAME,
    OUTPUT_MESH_FILENAME,
    NUM_ARGS
};


using namespace std;

int main( int argc, char *argv[] )
{

	if( argc < NUM_ARGS )
	{
		std::cerr << "\n\n";
		std::cerr << "Usage: " << argv[PROGRAM_NAME] << " imageFileName inputMeshFileName ouputMeshFileName";
		std::cerr << "\n\n";
		return 1;
	}
	itk::STLMeshIOFactory::RegisterOneFactory();

  const int dimension = 3;

  typedef float PixelType;
  typedef itk::Image<PixelType, dimension> ImageType;

  typedef itk::DefaultDynamicMeshTraits<PixelType, dimension>       TriangleMeshTraits;
  typedef itk::DefaultDynamicMeshTraits<PixelType, dimension>       SimplexMeshTraits;
  typedef itk::Mesh<PixelType, dimension, TriangleMeshTraits>       TriangleMeshType;
  typedef itk::SimplexMesh<PixelType, dimension, SimplexMeshTraits> SimplexMeshType;

  typedef itk::TriangleMeshToSimplexMeshFilter<TriangleMeshType, SimplexMeshType> TriangleToSimplexMeshType;
  typedef itk::SimplexMeshToTriangleMeshFilter<SimplexMeshType, TriangleMeshType> SimplexToTriangleFilterType;

  typedef itk::DeformableSimplexMesh3DFilter<SimplexMeshType, SimplexMeshType> DeformFilterType;
  typedef DeformFilterType::GradientImageType                                  GradientType;

  typedef itk::GradientAnisotropicDiffusionImageFilter <ImageType, ImageType> DiffusionFilterType;

  typedef itk::GradientAnisotropicDiffusionImageFilter < ImageType, ImageType >  GradientAnisotropicImageType;
  typedef itk::GradientMagnitudeRecursiveGaussianImageFilter < ImageType, ImageType >  GradientMagnitudeType;
  typedef itk::SigmoidImageFilter< ImageType, ImageType > SigmoidImageType;
  typedef itk::GradientRecursiveGaussianImageFilter<ImageType,GradientType> GradientFilterType;

  typedef itk::ImageFileReader<ImageType>       ImageReaderType;
  typedef itk::ImageFileWriter<ImageType>       ImageWriterType;
  typedef itk::MeshFileReader<TriangleMeshType> MeshReaderType;
  typedef itk::MeshFileWriter<TriangleMeshType>  MeshWriterType;

  typedef itk::TranslationTransform< TriangleMeshType::PointType::CoordRepType, 3 >   MeshTransformType;
  typedef itk::TransformMeshFilter< TriangleMeshType, TriangleMeshType, MeshTransformType >   MeshFilterType;


  float dspl_x, dspl_y, dspl_z;


  std::cout << "Reading image - " << std::flush;

  ImageReaderType::Pointer reader = ImageReaderType::New();
  reader->SetFileName( argv[IMAGE_FILENAME] );
  reader->Update();

  ImageType::Pointer img = reader->GetOutput();

  dspl_x = img->GetOrigin()[0];
  dspl_y = img->GetOrigin()[1];
  dspl_z = img->GetOrigin()[2];

  std::cout << "finished\n";

    // Read Mesh File
  SimplexMeshType::Pointer simplexMesh = ITK_NULLPTR;
  {
    std::cout << "Reading mesh - " << std::flush;

    MeshReaderType::Pointer meshReader = MeshReaderType::New();
    meshReader->SetFileName( argv[INPUT_MESH_FILENAME] );

    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( meshReader->GetOutput() );
    filter->SetTransform( translation );
    try
    {
      filter->Update();
    }
    catch( itk::ExceptionObject & exp )
    {
      std::cerr << "Exception thrown during Update() " << std::endl;
      std::cerr << exp << std::endl;
      return EXIT_FAILURE;
    }

    TriangleMeshType::Pointer translatedMesh = filter->GetOutput();
    translatedMesh->DisconnectPipeline();

    MeshWriterType::Pointer mesh_writer = MeshWriterType::New();
 
      // Set the writer's input and file name, then update.
    mesh_writer->SetInput( translatedMesh );
    mesh_writer->SetFileName( "translatedMesh.vtk" );
    mesh_writer->Update();

    TriangleToSimplexMeshType::Pointer cast = TriangleToSimplexMeshType::New();
    cast->SetInput( translatedMesh );
    cast->Update();

    simplexMesh = cast->GetOutput();
    simplexMesh->DisconnectPipeline();  

    std::cout << "finished\n";
    }  

    //
    // Correct image origin
    //
    {
      std::cout << "Changing image origin - " << std::flush;

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

      typedef itk::ChangeInformationImageFilter< ImageType > ChangeImgInfoType;
      ChangeImgInfoType::Pointer changeInfo = ChangeImgInfoType::New();
      changeInfo->SetInput( img );
      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;
        return EXIT_FAILURE;
      }

      img = changeInfo->GetOutput();

      ImageWriterType::Pointer writerimg = ImageWriterType::New();
      writerimg->SetInput( img );
      writerimg->SetFileName( "img_origin000.nii" );
      writerimg->Update();

      std::cout << "finished\n";
    }


  std::cout << " starting to Filter Image" << std::endl;
  GradientAnisotropicImageType::Pointer gradientanisotropicfilter = GradientAnisotropicImageType::New();
  gradientanisotropicfilter->SetInput(img);
  gradientanisotropicfilter->SetNumberOfIterations(5);
  gradientanisotropicfilter->SetTimeStep(0.0625);
  gradientanisotropicfilter->SetConductanceParameter(3);
  gradientanisotropicfilter->Update();
  std::cout << "GradientAnisotropicDiffusion is DONE!" << std::endl;

  GradientMagnitudeType::Pointer gradientmagnitudefilter = GradientMagnitudeType::New();
  gradientmagnitudefilter->SetInput( gradientanisotropicfilter->GetOutput() );
  gradientmagnitudefilter->SetSigma(1.0);
  gradientmagnitudefilter->Update();
  std::cout << "GradientMagnitude is DONE!" << std::endl;

  SigmoidImageType::Pointer sigmoidimagefilter = SigmoidImageType::New();
  sigmoidimagefilter->SetInput( gradientmagnitudefilter->GetOutput());
  sigmoidimagefilter->SetOutputMinimum(0);
  sigmoidimagefilter->SetOutputMaximum(1);
  sigmoidimagefilter->SetAlpha(10);
  sigmoidimagefilter->SetBeta(100);
  sigmoidimagefilter->Update();
  std::cout << "Sigmoid is DONE!" << std::endl;

  GradientFilterType::Pointer gradientFilter = GradientFilterType::New();
  gradientFilter->SetInput( sigmoidimagefilter->GetOutput() );
  gradientFilter->SetSigma(1.0);
  gradientFilter->Update();
  std::cout << "GradientMagnitude is DONE!" << std::endl;

  std::cout << "Gradient filter ok\n";

  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 pipeline
    simplexMesh->DisconnectPipeline();
    deformFilter->SetInput( simplexMesh );
    deformFilter->SetGradient( gradientImage );
    deformFilter->SetAlpha(0.1);
    deformFilter->SetBeta(0.1);
    deformFilter->SetIterations(5);
    deformFilter->SetRigidity(1);
    deformFilter->Update();
    }
  SimplexMeshType::Pointer deformResult =  deformFilter->GetOutput();


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

	  conversionResult->DisconnectPipeline();

	  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;
	    return EXIT_FAILURE;
	  }

	  outputMesh = filter2->GetOutput();
	  outputMesh->DisconnectPipeline();
  }


  MeshWriterType::Pointer meshWriter = MeshWriterType::New();
  meshWriter->SetFileName(argv[OUTPUT_MESH_FILENAME]);
  meshWriter->SetInput(outputMesh);
  meshWriter->Update();

  return EXIT_SUCCESS;
}


With minor modifications, the mesh is deformed. The displacement caclulation was wrong. Gradient calculation was overly complex, and probably did not yield the result you wanted. You don’t need to run the filter a 100 times, rather increase the number of iterations 100-fold. That avoids overhead of a 100 initializations. Here is the updated code. Best visualized using a diff tool against your original:

#include <iostream>


#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkSTLMeshIOFactory.h"
#include "itkSTLMeshIO.h"

#include "itkMeshFileReader.h"
#include "itkMeshFileWriter.h"
#include "itkDeformableSimplexMesh3DFilter.h"
#include "itkDefaultDynamicMeshTraits.h"
#include "itkMesh.h"
#include "itkSimplexMesh.h"
#include "itkTriangleMeshToSimplexMeshFilter.h"
#include "itkSimplexMeshToTriangleMeshFilter.h"

#include "itkGradientAnisotropicDiffusionImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkSigmoidImageFilter.h"
#include "itkGradientRecursiveGaussianImageFilter.h"

#include "itkTranslationTransform.h"
#include "itkTransformMeshFilter.h"

#include "itkChangeInformationImageFilter.h"


enum PARAMS {
    PROGRAM_NAME = 0,
    IMAGE_FILENAME,
    INPUT_MESH_FILENAME,
    OUTPUT_MESH_FILENAME,
    NUM_ARGS
};


using namespace std;

int main( int argc, char *argv[] )
{

	if( argc < NUM_ARGS )
	{
		std::cerr << "\n\n";
		std::cerr << "Usage: " << argv[PROGRAM_NAME] << " imageFileName inputMeshFileName ouputMeshFileName";
		std::cerr << "\n\n";
		return 1;
	}
	itk::STLMeshIOFactory::RegisterOneFactory();

  const int dimension = 3;

  typedef float PixelType;
  typedef itk::Image<PixelType, dimension> ImageType;

  typedef itk::DefaultDynamicMeshTraits<PixelType, dimension>       TriangleMeshTraits;
  typedef itk::DefaultDynamicMeshTraits<PixelType, dimension>       SimplexMeshTraits;
  typedef itk::Mesh<PixelType, dimension, TriangleMeshTraits>       TriangleMeshType;
  typedef itk::SimplexMesh<PixelType, dimension, SimplexMeshTraits> SimplexMeshType;

  typedef itk::TriangleMeshToSimplexMeshFilter<TriangleMeshType, SimplexMeshType> TriangleToSimplexMeshType;
  typedef itk::SimplexMeshToTriangleMeshFilter<SimplexMeshType, TriangleMeshType> SimplexToTriangleFilterType;

  typedef itk::DeformableSimplexMesh3DFilter<SimplexMeshType, SimplexMeshType> DeformFilterType;
  typedef DeformFilterType::GradientImageType                                  GradientType;

  typedef itk::GradientAnisotropicDiffusionImageFilter <ImageType, ImageType> DiffusionFilterType;

  typedef itk::GradientAnisotropicDiffusionImageFilter < ImageType, ImageType >  GradientAnisotropicImageType;
  typedef itk::GradientMagnitudeRecursiveGaussianImageFilter < ImageType, ImageType >  GradientMagnitudeType;
  typedef itk::SigmoidImageFilter< ImageType, ImageType > SigmoidImageType;
  typedef itk::GradientRecursiveGaussianImageFilter<ImageType,GradientType> GradientFilterType;

  typedef itk::ImageFileReader<ImageType>       ImageReaderType;
  typedef itk::ImageFileWriter<ImageType>       ImageWriterType;
  typedef itk::MeshFileReader<TriangleMeshType> MeshReaderType;
  typedef itk::MeshFileWriter<TriangleMeshType>  MeshWriterType;

  typedef itk::TranslationTransform< TriangleMeshType::PointType::CoordRepType, 3 >   MeshTransformType;
  typedef itk::TransformMeshFilter< TriangleMeshType, TriangleMeshType, MeshTransformType >   MeshFilterType;


  float dspl_x, dspl_y, dspl_z;


  std::cout << "Reading image - " << std::flush;

  ImageReaderType::Pointer reader = ImageReaderType::New();
  reader->SetFileName( argv[IMAGE_FILENAME] );
  reader->Update();

  ImageType::Pointer img = reader->GetOutput();

  dspl_x = img->GetOrigin()[0];
  dspl_y = img->GetOrigin()[1];
  dspl_z = img->GetOrigin()[2];

  std::cout << "finished\n";

    // Read Mesh File
  SimplexMeshType::Pointer simplexMesh = ITK_NULLPTR;
  {
    std::cout << "Reading mesh - " << std::flush;

    MeshReaderType::Pointer meshReader = MeshReaderType::New();
    meshReader->SetFileName( argv[INPUT_MESH_FILENAME] );

    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( meshReader->GetOutput() );
    filter->SetTransform( translation );
    try
    {
      filter->Update();
    }
    catch( itk::ExceptionObject & exp )
    {
      std::cerr << "Exception thrown during Update() " << std::endl;
      std::cerr << exp << std::endl;
      return EXIT_FAILURE;
    }

    TriangleMeshType::Pointer translatedMesh = filter->GetOutput();
    translatedMesh->DisconnectPipeline();

    MeshWriterType::Pointer mesh_writer = MeshWriterType::New();
 
      // Set the writer's input and file name, then update.
    mesh_writer->SetInput( translatedMesh );
    mesh_writer->SetFileName( "translatedMesh.vtk" );
    mesh_writer->Update();

    TriangleToSimplexMeshType::Pointer cast = TriangleToSimplexMeshType::New();
    cast->SetInput( translatedMesh );
    cast->Update();

    simplexMesh = cast->GetOutput();
    simplexMesh->DisconnectPipeline();  

    std::cout << "finished\n";
    }  

    //
    // Correct image origin
    //
    {
      std::cout << "Changing image origin - " << std::flush;

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

      typedef itk::ChangeInformationImageFilter< ImageType > ChangeImgInfoType;
      ChangeImgInfoType::Pointer changeInfo = ChangeImgInfoType::New();
      changeInfo->SetInput( img );
      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;
        return EXIT_FAILURE;
      }

      img = changeInfo->GetOutput();

      ImageWriterType::Pointer writerimg = ImageWriterType::New();
      writerimg->SetInput( img );
      writerimg->SetFileName( "img_origin000.nii" );
      writerimg->Update();

      std::cout << "finished\n";
    }


  std::cout << " starting to Filter Image" << std::endl;
  GradientAnisotropicImageType::Pointer gradientanisotropicfilter = GradientAnisotropicImageType::New();
  gradientanisotropicfilter->SetInput(img);
  gradientanisotropicfilter->SetNumberOfIterations(5);
  gradientanisotropicfilter->SetTimeStep(0.0625);
  gradientanisotropicfilter->SetConductanceParameter(3);
  gradientanisotropicfilter->Update();
  std::cout << "GradientAnisotropicDiffusion is DONE!" << std::endl;

  //GradientMagnitudeType::Pointer gradientmagnitudefilter = GradientMagnitudeType::New();
  //gradientmagnitudefilter->SetInput( gradientanisotropicfilter->GetOutput() );
  //gradientmagnitudefilter->SetSigma(1.0);
  //gradientmagnitudefilter->Update();
  //std::cout << "GradientMagnitude is DONE!" << std::endl;

  //SigmoidImageType::Pointer sigmoidimagefilter = SigmoidImageType::New();
  //sigmoidimagefilter->SetInput( gradientmagnitudefilter->GetOutput());
  //sigmoidimagefilter->SetOutputMinimum(0);
  //sigmoidimagefilter->SetOutputMaximum(1);
  //sigmoidimagefilter->SetAlpha(10);
  //sigmoidimagefilter->SetBeta(100);
  //sigmoidimagefilter->Update();
  //std::cout << "Sigmoid is DONE!" << std::endl;

  GradientFilterType::Pointer gradientFilter = GradientFilterType::New();
  gradientFilter->SetInput(gradientanisotropicfilter->GetOutput() );
  gradientFilter->SetSigma(1.0);
  gradientFilter->Update();
  std::cout << "GradientMagnitude is DONE!" << std::endl;

  std::cout << "Gradient filter ok\n";

  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 pipeline
    simplexMesh->DisconnectPipeline();
    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();


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

	  conversionResult->DisconnectPipeline();

	  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;
	    return EXIT_FAILURE;
	  }

	  outputMesh = filter2->GetOutput();
	  outputMesh->DisconnectPipeline();
  }


  MeshWriterType::Pointer meshWriter = MeshWriterType::New();
  meshWriter->SetFileName(argv[OUTPUT_MESH_FILENAME]);
  meshWriter->SetInput(outputMesh);
  meshWriter->Update();

  return EXIT_SUCCESS;
}
1 Like

@dzenanz Thanks a lot for your helpful comments. 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. I do not know if I am doing some mistakes r it is because of my device that is slowing down the process?!

Thanks in advance