[ITK-users] [ITK] Applying a B-Spline transform to a 3D mesh
Carlos Henrique Villa Pinto
chvillap at gmail.com
Thu Sep 24 09:43:09 EDT 2015
Hi, Andras. Thank you very much for your answer.
I'm not very experienced with B-Spline transforms and 3D meshes, so maybe
this is a silly question, but... why exactly do I need inverse transforms
(forward for the mesh, inverse for the image), considering that I want to
deform both the same way? Could you elaborate, or perhaps give me some
links with explanations, please?
By the way, I generated a deformation field from the transform I created,
and visualized/played with it in 3D Slicer. As expected, applying such
transform to the original MR image in 3D Slicer results in the same
deformed image that I get using my ITK code. However, when I apply the
same transform to the original mesh (again in 3D Slicer), the result is a
deformed mesh that is neither the same result that I obtain with my ITK
code (using the itk::TransformMeshFilter class) nor the expected result (a
mesh aligned to the brain structure in the deformed image). And applying
the inverse transform (also calculated in 3D Slicer) doesn't give me better
results at all. So now I'm even more confused.
Any other insights would be greatly appreciated.
Thanks.
[]s
2015-09-23 10:56 GMT-03:00 Andras Lasso <lasso at queensu.ca>:
> You need forward transform to deform a mesh and inverse (resampling)
> transform to transform the corresponding image. Make sure the transform
> that you generate is invertible:
>
> · Displace control points just so much that there is no “folding”
> in the field (two different areas are transformed to the same area)
>
> · Make the grid larger than the image so that the displacements are
> all zero at the boundary (as ITK transforms cannot extrapolate to have
> smooth boundary conditions but they just cut the displacement to 0; if
> there is discontinuity in the field then inverse computation will fail)
>
>
>
> You can play with non-linear transforms, import/export from ITK, visualize
> them, invert them, combine them, apply to volumes, surfaces, points, etc.
> in 3D Slicer, see:
>
>
> http://www.slicer.org/slicerWiki/index.php/Documentation/Nightly/Modules/Transforms
>
>
>
> Andras
>
>
>
> *From:* Community [mailto:community-bounces at itk.org] *On Behalf Of *Carlos
> Henrique Villa Pinto
> *Sent:* September 23, 2015 9:38 AM
> *To:* insight-users at itk.org
> *Subject:* [ITK] [ITK-users] Applying a B-Spline transform to a 3D mesh
>
>
>
> Hi everyone,
>
>
>
> I have a 3D MR brain image and a set of 3D meshes of different brain
> structures aligned to it. I'm trying to write a very simple program in
> which I create a B-Spline transform (with random parameters) and apply such
> transform to the image and to the meshes in order to deform them all
> together.
>
>
>
> Am I wrong to assume that, if I apply exactly the same transform to an
> image and to a mesh which is perfectly aligned to a specific brain
> structure in such image, the resulting deformed mesh should be aligned to
> the same brain structure in the deformed image? Because that's not what
> it's happening in practice. The transform deforms the meshes to some
> extent, but the deformed meshes do not match the brain structures in the
> deformed image, even though the applied transform is the same for all. So I
> guess I must be doing something wrong...
>
>
>
> The source code is below. I'm using ITK 4.8, and the MR image (NIfTI) and
> meshes (VTK) used in my tests come from the NAC Brain Atlas (
> https://www.slicer.org/publications/item/view/2037).
>
>
>
> Any insight about what is possibly causing these results would be
> appreciated. Thanks in advance.
>
>
>
>
>
> #include <iostream>
>
> #include <vector>
>
> #include <cstdlib>
>
> #include <cstring>
>
> #include <cmath>
>
> #include <ctime>
>
>
>
> #include <itkImage.h>
>
> #include <itkImageFileReader.h>
>
> #include <itkImageFileWriter.h>
>
> #include <itkMesh.h>
>
> #include <itkMeshFileReader.h>
>
> #include <itkMeshFileWriter.h>
>
> #include <itkBSplineTransform.h>
>
> #include <itkResampleImageFilter.h>
>
> #include <itkTransformMeshFilter.h>
>
> #include <itkLinearInterpolateImageFunction.h>
>
>
>
> int main(int argc, char const *argv[])
>
> {
>
> if (argc < 3)
>
> {
>
> std::cerr << "Usage: " << argv[0]
>
> << " gridNodesInOneDimension inputImage [inputMesh1,
> ...]"
>
> << std::endl;
>
>
>
> return EXIT_FAILURE;
>
> }
>
>
>
> srand(time(NULL));
>
>
>
> const int splineOrder = 3;
>
> const int dimensions = 3;
>
>
>
> typedef itk::Image<double, dimensions> ImageType;
>
> typedef itk::ImageFileReader<ImageType> ImageReaderType;
>
> typedef itk::ImageFileWriter<ImageType> ImageWriterType;
>
>
>
> typedef itk::Mesh<double, dimensions> MeshType;
>
> typedef itk::MeshFileReader<MeshType> MeshReaderType;
>
> typedef itk::MeshFileWriter<MeshType> MeshWriterType;
>
>
>
> typedef itk::BSplineTransform<double, dimensions, splineOrder>
>
> BSplineTransformType;
>
> typedef itk::ResampleImageFilter<ImageType, ImageType, double>
>
> ResampleImageFilterType;
>
> typedef itk::TransformMeshFilter<MeshType, MeshType,
> BSplineTransformType>
>
> TransformMeshFilterType;
>
> typedef itk::LinearInterpolateImageFunction<ImageType, double>
>
> LinearInterpolatorType;
>
>
>
> //
> ============================================================================================
>
> // Reading the input image
>
> //
> ============================================================================================
>
>
>
> ImageReaderType::Pointer imageReader = ImageReaderType::New();
>
> imageReader->SetFileName(argv[2]);
>
> imageReader->Update();
>
>
>
> ImageType::Pointer inputImage = imageReader->GetOutput();
>
>
>
> //
> ============================================================================================
>
> // Reading the input meshes
>
> //
> ============================================================================================
>
>
>
> std::vector<MeshType::Pointer> inputMeshVector;
>
>
>
> for (int i = 3; i < argc; ++i)
>
> {
>
> MeshReaderType::Pointer meshReader = MeshReaderType::New();
>
> meshReader->SetFileName(argv[i]);
>
> meshReader->Update();
>
>
>
> MeshType::Pointer inputMesh = meshReader->GetOutput();
>
> inputMeshVector.push_back(inputMesh);
>
> }
>
>
>
> //
> ============================================================================================
>
> // Generating the parameters of the B-Spline transform
>
> //
> ============================================================================================
>
>
>
> const int gridNodesInOneDimension = atoi(argv[1]);
>
>
>
> BSplineTransformType::Pointer bsplineTransform =
> BSplineTransformType::New();
>
> BSplineTransformType::PhysicalDimensionsType physicalDimensions;
>
> BSplineTransformType::MeshSizeType meshSize;
>
>
>
> for (int i = 0; i < dimensions; ++i)
>
> {
>
> physicalDimensions[i] = inputImage->GetSpacing()[i] *
>
>
> static_cast<double>(inputImage->GetLargestPossibleRegion().GetSize()[i] -
> 1);
>
> }
>
> meshSize.Fill(gridNodesInOneDimension - splineOrder);
>
>
>
> bsplineTransform->SetTransformDomainOrigin(inputImage->GetOrigin());
>
>
> bsplineTransform->SetTransformDomainDirection(inputImage->GetDirection());
>
>
> bsplineTransform->SetTransformDomainPhysicalDimensions(physicalDimensions);
>
> bsplineTransform->SetTransformDomainMeshSize(meshSize);
>
>
>
> BSplineTransformType::ParametersType
> parameters(bsplineTransform->GetNumberOfParameters());
>
>
>
> double minValue = -20.0;
>
> double maxValue = 20.0;
>
> for (int i = 0; i < bsplineTransform->GetNumberOfParameters(); ++i)
>
> {
>
> parameters[i] = (double) rand() / RAND_MAX * (maxValue - minValue)
> + minValue;
>
> std::cout << parameters[i] << " ";
>
> }
>
> std::cout << std::endl;
>
>
>
> bsplineTransform->SetParameters(parameters);
>
>
>
> //
> ============================================================================================
>
> // Applying the B-Spline transform to the input image
>
> //
> ============================================================================================
>
>
>
> LinearInterpolatorType::Pointer linearInterpolator =
> LinearInterpolatorType::New();
>
>
>
> ResampleImageFilterType::Pointer resampleFilter =
> ResampleImageFilterType::New();
>
> resampleFilter->SetTransform(bsplineTransform);
>
> resampleFilter->SetInterpolator(linearInterpolator);
>
> resampleFilter->SetInput(inputImage);
>
>
> resampleFilter->SetSize(inputImage->GetLargestPossibleRegion().GetSize());
>
>
> resampleFilter->SetOutputStartIndex(inputImage->GetLargestPossibleRegion().GetIndex());
>
> resampleFilter->SetOutputOrigin(inputImage->GetOrigin());
>
> resampleFilter->SetOutputSpacing(inputImage->GetSpacing());
>
> resampleFilter->SetOutputDirection(inputImage->GetDirection());
>
> resampleFilter->Update();
>
>
>
> //
> ============================================================================================
>
> // Applying the B-Spline transform to the meshes
>
> //
> ============================================================================================
>
>
>
> std::vector<MeshType::Pointer> outputMeshVector;
>
>
>
> for (int i = 0; i < inputMeshVector.size(); ++i)
>
> {
>
> TransformMeshFilterType::Pointer transformFilter =
> TransformMeshFilterType::New();
>
> transformFilter->SetInput(inputMeshVector[i]);
>
> transformFilter->SetTransform(bsplineTransform);
>
> transformFilter->Update();
>
>
>
> outputMeshVector.push_back(transformFilter->GetOutput());
>
> }
>
>
>
> //
> ============================================================================================
>
> // Writing the transformed image
>
> //
> ============================================================================================
>
>
>
> char outputImageFileName[256];
>
> int length1 = strlen(argv[2]) - 4;
>
>
>
> strncpy(outputImageFileName, argv[2], length1);
>
> outputImageFileName[length1] = '\0';
>
> strcat(outputImageFileName, "_transformed.nii");
>
>
>
> ImageWriterType::Pointer imageWriter = ImageWriterType::New();
>
> imageWriter->SetFileName(outputImageFileName);
>
> imageWriter->SetInput(resampleFilter->GetOutput());
>
> imageWriter->Update();
>
>
>
> //
> ============================================================================================
>
> // Writing the transformed meshes
>
> //
> ============================================================================================
>
>
>
> for (int i = 0; i < outputMeshVector.size(); ++i)
>
> {
>
> char outputMeshFilename[256];
>
> int length2 = strlen(argv[i + 3]) - 4;
>
>
>
> strncpy(outputMeshFilename, argv[i + 3], length2);
>
> outputMeshFilename[length2] = '\0';
>
> strcat(outputMeshFilename, "_transformed.vtk");
>
>
>
> MeshType::Pointer outputMesh = outputMeshVector[i];
>
>
>
> MeshWriterType::Pointer meshWriter = MeshWriterType::New();
>
> meshWriter->SetFileName(outputMeshFilename);
>
> meshWriter->SetInput(outputMesh);
>
> meshWriter->Update();
>
> }
>
>
>
> return EXIT_SUCCESS;
>
> }
>
>
>
>
>
> []s
>
>
>
> --
>
> Carlos Henrique Villa Pinto
>
> Graduate Student in Computer Science
>
> Federal University of São Carlos - Brazil
>
> XCS
>
--
Carlos Henrique Villa Pinto
Graduate Student in Computer Science
Federal University of São Carlos - Brazil
XCS
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/insight-users/attachments/20150924/8a0cc2f6/attachment.html>
More information about the Insight-users
mailing list