<div dir="ltr"><div style="font-size:12.8000001907349px">Hi everyone,</div><div style="font-size:12.8000001907349px"><br></div><div style="font-size:12.8000001907349px">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.</div><div style="font-size:12.8000001907349px"><br></div><div style="font-size:12.8000001907349px">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...</div><div style="font-size:12.8000001907349px"><br></div><div style="font-size:12.8000001907349px">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 (<a href="https://www.slicer.org/publications/item/view/2037" target="_blank">https://www.slicer.org/publications/item/view/2037</a>).</div><div style="font-size:12.8000001907349px"><br></div><div style="font-size:12.8000001907349px">Any insight about what is possibly causing these results would be appreciated. Thanks in advance.</div><div style="font-size:12.8000001907349px"><br></div><div style="font-size:12.8000001907349px"><br></div><div style="font-size:12.8000001907349px">#include <iostream><br></div><div style="font-size:12.8000001907349px"><div>#include <vector></div><div>#include <cstdlib></div><div>#include <cstring></div><div>#include <cmath></div><div>#include <ctime></div><div><br></div><div>#include <itkImage.h></div><div>#include <itkImageFileReader.h></div><div>#include <itkImageFileWriter.h></div><div>#include <itkMesh.h></div><div>#include <itkMeshFileReader.h></div><div>#include <itkMeshFileWriter.h></div><div>#include <itkBSplineTransform.h></div><div>#include <itkResampleImageFilter.h></div><div>#include <itkTransformMeshFilter.h></div><div>#include <itkLinearInterpolateImageFunction.h></div><div><br></div><div>int main(int argc, char const *argv[])</div><div>{</div><div> if (argc < 3)</div><div> {</div><div> std::cerr << "Usage: " << argv[0]</div><div> << " gridNodesInOneDimension inputImage [inputMesh1, ...]"</div><div> << std::endl;</div><div><br></div><div> return EXIT_FAILURE;</div><div> }</div><div><br></div><div> srand(time(NULL));</div><div><br></div><div> const int splineOrder = 3;</div><div> const int dimensions = 3;</div><div><br></div><div> typedef itk::Image<double, dimensions> ImageType;</div><div> typedef itk::ImageFileReader<ImageType> ImageReaderType;</div><div> typedef itk::ImageFileWriter<ImageType> ImageWriterType;</div><div><br></div><div> typedef itk::Mesh<double, dimensions> MeshType;</div><div> typedef itk::MeshFileReader<MeshType> MeshReaderType;</div><div> typedef itk::MeshFileWriter<MeshType> MeshWriterType;</div><div><br></div><div> typedef itk::BSplineTransform<double, dimensions, splineOrder></div><div> BSplineTransformType;</div><div> typedef itk::ResampleImageFilter<ImageType, ImageType, double></div><div> ResampleImageFilterType;</div><div> typedef itk::TransformMeshFilter<MeshType, MeshType, BSplineTransformType></div><div> TransformMeshFilterType;</div><div> typedef itk::LinearInterpolateImageFunction<ImageType, double></div><div> LinearInterpolatorType;</div><div><br></div><div> // ============================================================================================</div><div> // Reading the input image</div><div> // ============================================================================================</div><div><br></div><div> ImageReaderType::Pointer imageReader = ImageReaderType::New();</div><div> imageReader->SetFileName(argv[2]);</div><div> imageReader->Update();</div><div><br></div><div> ImageType::Pointer inputImage = imageReader->GetOutput();</div><div><br></div><div> // ============================================================================================</div><div> // Reading the input meshes</div><div> // ============================================================================================</div><div><br></div><div> std::vector<MeshType::Pointer> inputMeshVector;</div><div><br></div><div> for (int i = 3; i < argc; ++i)</div><div> {</div><div> MeshReaderType::Pointer meshReader = MeshReaderType::New();</div><div> meshReader->SetFileName(argv[i]);</div><div> meshReader->Update();</div><div><br></div><div> MeshType::Pointer inputMesh = meshReader->GetOutput();</div><div> inputMeshVector.push_back(inputMesh);</div><div> }</div><div><br></div><div> // ============================================================================================</div><div> // Generating the parameters of the B-Spline transform</div><div> // ============================================================================================</div><div><br></div><div> const int gridNodesInOneDimension = atoi(argv[1]);</div><div><br></div><div> BSplineTransformType::Pointer bsplineTransform = BSplineTransformType::New();</div><div> BSplineTransformType::PhysicalDimensionsType physicalDimensions;</div><div> BSplineTransformType::MeshSizeType meshSize;</div><div><br></div><div> for (int i = 0; i < dimensions; ++i)</div><div> {</div><div> physicalDimensions[i] = inputImage->GetSpacing()[i] *</div><div> static_cast<double>(inputImage->GetLargestPossibleRegion().GetSize()[i] - 1);</div><div> }</div><div> meshSize.Fill(gridNodesInOneDimension - splineOrder);</div><div> </div><div> bsplineTransform->SetTransformDomainOrigin(inputImage->GetOrigin());</div><div> bsplineTransform->SetTransformDomainDirection(inputImage->GetDirection());</div><div> bsplineTransform->SetTransformDomainPhysicalDimensions(physicalDimensions);</div><div> bsplineTransform->SetTransformDomainMeshSize(meshSize);</div><div><br></div><div> BSplineTransformType::ParametersType parameters(bsplineTransform->GetNumberOfParameters());</div><div><br></div><div> double minValue = -20.0;</div><div> double maxValue = 20.0;</div><div> for (int i = 0; i < bsplineTransform->GetNumberOfParameters(); ++i)</div><div> {</div><div> parameters[i] = (double) rand() / RAND_MAX * (maxValue - minValue) + minValue;</div><div> std::cout << parameters[i] << " ";</div><div> }</div><div> std::cout << std::endl;</div><div><br></div><div> bsplineTransform->SetParameters(parameters);</div><div><br></div><div> // ============================================================================================</div><div> // Applying the B-Spline transform to the input image</div><div> // ============================================================================================</div><div><br></div><div> LinearInterpolatorType::Pointer linearInterpolator = LinearInterpolatorType::New();</div><div><br></div><div> ResampleImageFilterType::Pointer resampleFilter = ResampleImageFilterType::New();</div><div> resampleFilter->SetTransform(bsplineTransform);</div><div> resampleFilter->SetInterpolator(linearInterpolator);</div><div> resampleFilter->SetInput(inputImage);</div><div> resampleFilter->SetSize(inputImage->GetLargestPossibleRegion().GetSize());</div><div> resampleFilter->SetOutputStartIndex(inputImage->GetLargestPossibleRegion().GetIndex());</div><div> resampleFilter->SetOutputOrigin(inputImage->GetOrigin());</div><div> resampleFilter->SetOutputSpacing(inputImage->GetSpacing());</div><div> resampleFilter->SetOutputDirection(inputImage->GetDirection());</div><div> resampleFilter->Update();</div><div><br></div><div> // ============================================================================================</div><div> // Applying the B-Spline transform to the meshes</div><div> // ============================================================================================</div><div><br></div><div> std::vector<MeshType::Pointer> outputMeshVector;</div><div><br></div><div> for (int i = 0; i < inputMeshVector.size(); ++i)</div><div> {</div><div> TransformMeshFilterType::Pointer transformFilter = TransformMeshFilterType::New();</div><div> transformFilter->SetInput(inputMeshVector[i]);</div><div> transformFilter->SetTransform(bsplineTransform);</div><div> transformFilter->Update();</div><div><br></div><div> outputMeshVector.push_back(transformFilter->GetOutput());</div><div> }</div><div><br></div><div> // ============================================================================================</div><div> // Writing the transformed image</div><div> // ============================================================================================</div><div><br></div><div> char outputImageFileName[256];</div><div> int length1 = strlen(argv[2]) - 4;</div><div><br></div><div> strncpy(outputImageFileName, argv[2], length1);</div><div> outputImageFileName[length1] = '\0';</div><div> strcat(outputImageFileName, "_transformed.nii");</div><div><br></div><div> ImageWriterType::Pointer imageWriter = ImageWriterType::New();</div><div> imageWriter->SetFileName(outputImageFileName);</div><div> imageWriter->SetInput(resampleFilter->GetOutput());</div><div> imageWriter->Update();</div><div><br></div><div> // ============================================================================================</div><div> // Writing the transformed meshes</div><div> // ============================================================================================</div><div><br></div><div> for (int i = 0; i < outputMeshVector.size(); ++i)</div><div> {</div><div> char outputMeshFilename[256];</div><div> int length2 = strlen(argv[i + 3]) - 4;</div><div><br></div><div> strncpy(outputMeshFilename, argv[i + 3], length2);</div><div> outputMeshFilename[length2] = '\0';</div><div> strcat(outputMeshFilename, "_transformed.vtk");</div><div><br></div><div> MeshType::Pointer outputMesh = outputMeshVector[i];</div><div><br></div><div> MeshWriterType::Pointer meshWriter = MeshWriterType::New();</div><div> meshWriter->SetFileName(outputMeshFilename);</div><div> meshWriter->SetInput(outputMesh);</div><div> meshWriter->Update();</div><div> }</div><div><br></div><div> return EXIT_SUCCESS;</div><div>}</div></div><div style="font-size:12.8000001907349px"><br></div><div style="font-size:12.8000001907349px"><br></div><div style="font-size:12.8000001907349px">[]s</div><div><br></div>-- <br><div class="gmail_signature"><div dir="ltr">Carlos Henrique Villa Pinto<div>Graduate Student in Computer Science</div><div>Federal University of São Carlos - Brazil</div><div>XCS</div></div></div>
</div>