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:
- a mesh created from the segmentation result of an algorithm (mesh created by Cuberille filter)
- 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