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