# Deformable Registration and .vtk file Help

**URL:** <https://discourse.itk.org/t/deformable-registration-and-vtk-file-help/2699>\
**Category:** Uncategorized\
**Created:** [February 10, 2020, 10:58pm UTC](https://discourse.itk.org/t/deformable-registration-and-vtk-file-help/2699 "2020-02-10T22:58:21Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![ShaneKirkbride](https://discourse.itk.org/user_avatar/discourse.itk.org/shanekirkbride/32/1225_2.png) [@ShaneKirkbride](https://discourse.itk.org/u/ShaneKirkbride)\
**Post date:** [February 10, 2020, 10:58pm UTC](https://discourse.itk.org/t/deformable-registration-and-vtk-file-help/2699/1 "2020-02-10T22:58:21Z")

</div>

I have a set of similar DICOM files. I need to convert them to a 3D mesh and then register them. Then I need to transform the meshes so they are both on the same coordinate system. After this I need to calculate the area of each triangle in the mesh and do things with these values.

I’ve used VTK to create the 3D isoSurface similiar to marching cubes. Then I write to a .vtk ascii file the isoSurfaces.

I’d like to read these into the DeformableRegistration7 example but I keep getting an error:  
VTKImageIO(0000025A4187AD90): Not structured points, can’t read

Here are the main changes I’ve made

```
    FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
    MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
    ImageIOType::Pointer vtkIO = ImageIOType::New();
    vtkIO->SetFileTypeToASCII();

    fixedImageReader->SetFileName(targetFilename);
    fixedImageReader->SetImageIO(vtkIO);

    movingImageReader->SetFileName(sourceFilename);
    movingImageReader->SetImageIO(vtkIO);

    FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();

    registration->SetFixedImage(fixedImage);
    registration->SetMovingImage(movingImageReader->GetOutput());

```

What type of 3D Images are good for this example? Is there a better way to do this? I like using VTK but if there is a way to do this in ITK i’d go for it.

I’m working in C++.

Thanks,  
~Shane

---

<div class="post-metadata">

**Author:** ![matt.mccormick](https://discourse.itk.org/user_avatar/discourse.itk.org/matt.mccormick/32/7_2.png) [@matt.mccormick](https://discourse.itk.org/u/matt.mccormick)\
**Post date:** [February 10, 2020, 11:32pm UTC](https://discourse.itk.org/t/deformable-registration-and-vtk-file-help/2699/2 "2020-02-10T23:32:26Z")

</div>

Hi Shane,

To read the mesh, use [itk::MeshFileReader](https://itk.org/ITKExamples/src/IO/Mesh/ReadMesh/Documentation.html). To register the meshes, use [itk::PointSetToPointSetMetricv4](https://itk.org/Doxygen/html/classitk_1_1PointSetToPointSetMetricv4.html).

HTH,  
Matt

---

<div class="post-metadata">

**Author:** ![ShaneKirkbride](https://discourse.itk.org/user_avatar/discourse.itk.org/shanekirkbride/32/1225_2.png) [@ShaneKirkbride](https://discourse.itk.org/u/ShaneKirkbride)\
**Post date:** [February 11, 2020, 12:38am UTC](https://discourse.itk.org/t/deformable-registration-and-vtk-file-help/2699/3 "2020-02-11T00:38:00Z")

</div>

Hi Matt,  
Thanks for the pointers!  
MeshFileReader makes sense.

Any examples for the PointToPoint Classes? I can probably figure it out but guidance is always helpful.  
-Shane

---

<div class="post-metadata">

**Author:** ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)\
**Post date:** [February 11, 2020, 3:20pm UTC](https://discourse.itk.org/t/deformable-registration-and-vtk-file-help/2699/4 "2020-02-11T15:20:42Z")

</div>

Here is an [example](https://github.com/InsightSoftwareConsortium/ITK/blob/master/Examples/RegistrationITKv4/IterativeClosestPoint3.cxx). You need to modify it to read a mesh, then get mesh’s vertices (instead of reading a list of points “manually”).

---

<div class="post-metadata">

**Author:** ![ShaneKirkbride](https://discourse.itk.org/user_avatar/discourse.itk.org/shanekirkbride/32/1225_2.png) [@ShaneKirkbride](https://discourse.itk.org/u/ShaneKirkbride)\
**Post date:** [February 11, 2020, 4:49pm UTC](https://discourse.itk.org/t/deformable-registration-and-vtk-file-help/2699/5 "2020-02-11T16:49:35Z")

</div>

Hi Dzenan,  
Thank you for the reference. I tried reading in .vtp and .vtk files like this:

```
    constexpr unsigned int Dimension = 3;
    using CoordinateType = float;
    using MeshType = itk::Mesh<CoordinateType, Dimension>;
    typedef itk::MeshFileReader<MeshType> MeshReaderType;
    MeshReaderType::Pointer Meshreader = MeshReaderType::New();
    Meshreader->SetFileName("C:\\Users\\shane\\source\\repos\\Registration-1\\Registration-1\\targetMesh.vtp");
    Meshreader->Update();

```

And I keep getting an error saying that the file suffix is wrong or there is something wrong with the file.

I’m writing the file in VTK like this, which I believe is correct:

```
// Write the file
vtkSmartPointer<vtkXMLPolyDataWriter> writer = vtkSmartPointer<vtkXMLPolyDataWriter>::New();
writer->SetFileName(filename.c_str());
writer->SetInputData(mesh);

// Optional - set the mode. The default is binary.
//writer->SetDataModeToBinary();
//writer->SetDataModeToAscii();

writer->Write();
return 0;

```

Here is the first few lines of the header from the file generated:

```
<?xml version="1.0"?>
<VTKFile type="PolyData" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
  <PolyData>
    <Piece NumberOfPoints="28245" NumberOfVerts="0" NumberOfLines="0" NumberOfStrips="0" NumberOfPolys="56486" >
      <PointData Scalars="DICOMImage" Normals="Normals">
        <DataArray type="Int16" Name="DICOMImage" format="appended" RangeMin="150" RangeMax="150" offset="0" />
        <DataArray type="Float32" Name="Normals" NumberOfComponents="3" format="appended" RangeMin="0.99999984582" RangeMax="1.0000001505" offset="168" />
      </PointData>

```

What am I doing wrong here?  
Thanks,  
~Shane

---

<div class="post-metadata">

**Author:** ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)\
**Post date:** [February 11, 2020, 4:53pm UTC](https://discourse.itk.org/t/deformable-registration-and-vtk-file-help/2699/6 "2020-02-11T16:53:40Z")

</div>

ITK does not support VTK’s `XMLPolyData` file format. We only support the legacy file format (usual extension `.vtk`). I think [vtkPolyDataWriter](https://vtk.org/doc/nightly/html/classvtkPolyDataWriter.html) is the right class to use.

---

<div class="post-metadata">

**Author:** ![ShaneKirkbride](https://discourse.itk.org/user_avatar/discourse.itk.org/shanekirkbride/32/1225_2.png) [@ShaneKirkbride](https://discourse.itk.org/u/ShaneKirkbride)\
**Post date:** [February 11, 2020, 5:10pm UTC](https://discourse.itk.org/t/deformable-registration-and-vtk-file-help/2699/7 "2020-02-11T17:10:06Z")

</div>

Hi Dzenan,  
Tried the polyDataWriter and still have the same issue (both Ascii and binary formats)  
Could it be a .vtk version issue? I read that ITK only supports version 3.0 but not sure if this is still accurate.  
Here is the .vtk header with polyDataWriter:

```
# vtk DataFile Version 5.0
vtk output
ASCII
DATASET POLYDATA
POINTS 28245 float
```

---

<div class="post-metadata">

**Author:** ![ShaneKirkbride](https://discourse.itk.org/user_avatar/discourse.itk.org/shanekirkbride/32/1225_2.png) [@ShaneKirkbride](https://discourse.itk.org/u/ShaneKirkbride)\
**Post date:** [February 11, 2020, 5:18pm UTC](https://discourse.itk.org/t/deformable-registration-and-vtk-file-help/2699/8 "2020-02-11T17:18:46Z")

</div>

I did find this: [https://github.com/InsightSoftwareConsortium/ITKApps/blob/master/Auxiliary/vtk/vtkPolyDataToITKMesh.cxx](https://github.com/InsightSoftwareConsortium/ITKApps/blob/master/Auxiliary/vtk/vtkPolyDataToITKMesh.cxx)

However it seems somewhat outdated. I basically need to get my polyData into an ITK Mesh for the registration.

Thanks,  
~Shane

---

<div class="post-metadata">

**Author:** ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)\
**Post date:** [February 11, 2020, 6:53pm UTC](https://discourse.itk.org/t/deformable-registration-and-vtk-file-help/2699/9 "2020-02-11T18:53:57Z")

</div>

The extension also needs to be `.vtk`, I think. If it is `.vtp`, ITK thinks it is an unsupported format.

---

<div class="post-metadata">

**Author:** ![ShaneKirkbride](https://discourse.itk.org/user_avatar/discourse.itk.org/shanekirkbride/32/1225_2.png) [@ShaneKirkbride](https://discourse.itk.org/u/ShaneKirkbride)\
**Post date:** [February 11, 2020, 8:03pm UTC](https://discourse.itk.org/t/deformable-registration-and-vtk-file-help/2699/10 "2020-02-11T20:03:58Z")

</div>

Hi Dzenanz,  
Both .vtk and .vtp fail. However I was able to do this; and it seems to work:

```
SmartPointer<MeshType> convertToITKData(string filename)
{
    vtkPolyDataReader* reader = vtkPolyDataReader::New();
    reader->SetFileName(filename.c_str());
    //reader->SetFileName("C:\\Users\\shane\\source\\repos\\Registration-1\\Registration-1\\targetMesh.vtk");
    reader->Update();

    vtkPolyData* polyData = reader->GetOutput();

    MeshType::Pointer mesh = MeshType::New();

    //
    // Transfer the points from the vtkPolyData into the itk::Mesh
    //
    const unsigned int numberOfPoints = polyData->GetNumberOfPoints();

    vtkPoints* vtkpoints = polyData->GetPoints();

    mesh->GetPoints()->Reserve(numberOfPoints);

    for (unsigned int p = 0; p < numberOfPoints; p++)
    {

        vtkFloatingPointType* apoint = (float*)vtkpoints->GetPoint(p);

        mesh->SetPoint(p, MeshType::PointType(apoint));

    }

    //
    // Transfer the cells from the vtkPolyData into the itk::Mesh
    //
    vtkSmartPointer<vtkCellArrayIterator> triangleStripsIter = vtkSmartPointer<vtkCellArrayIterator>::New();
    vtkCellArray* triangleStrips = polyData->GetStrips();
    int numberOfCells = triangleStrips->GetNumberOfCells();

    vtkIdType* cellPoints = new vtkIdType;
    vtkIdType numberOfCellPoints = numberOfCells*2;

    //
    // First count the total number of triangles from all the triangle strips.
    //
    unsigned int numberOfTriangles = 0;
    triangleStrips->InitTraversal();
    triangleStripsIter = vtk::TakeSmartPointer(triangleStrips->NewIterator());
    
    for (triangleStripsIter->GoToFirstCell(); !triangleStripsIter->IsDoneWithTraversal(); triangleStripsIter->GoToNextCell())
    {
        numberOfTriangles += numberOfCellPoints - 2;
    }
    
    vtkSmartPointer<vtkCellArrayIterator> polygonsIter = vtkSmartPointer<vtkCellArrayIterator>::New();
    vtkCellArray* polygons = polyData->GetPolys();
    polygons->InitTraversal();
    polygonsIter = vtk::TakeSmartPointer(polygons->NewIterator());
    for (polygonsIter->GoToFirstCell(); !polygonsIter->IsDoneWithTraversal(); polygonsIter->GoToNextCell())
    {

        if (numberOfCellPoints == 3)
        {
            numberOfTriangles++;
        }
    }
          
    //
    // Reserve memory in the itk::Mesh for all those triangles
    //
    mesh->GetCells()->Reserve(numberOfTriangles);
          
    // 
    // Copy the triangles from vtkPolyData into the itk::Mesh
    //
    //

    typedef MeshType::CellType CellType;

    typedef itk::TriangleCell< CellType > TriangleCellType;

    int cellId = 0;

    // first copy the triangle strips
    triangleStrips->InitTraversal();
    triangleStripsIter = vtk::TakeSmartPointer(triangleStrips->NewIterator());
    for (triangleStripsIter->GoToFirstCell(); !triangleStripsIter->IsDoneWithTraversal(); triangleStripsIter->GoToNextCell())
    {

        unsigned int numberOfTrianglesInStrip = numberOfCellPoints - 2;

        unsigned long pointIds[3];
        pointIds[0] = cellPoints[0];
        pointIds[1] = cellPoints[1];
        pointIds[2] = cellPoints[2];

        for (unsigned int t = 0; t < numberOfTrianglesInStrip; t++)
        {
            MeshType::CellAutoPointer c;
            TriangleCellType* tcell = new TriangleCellType;
            TriangleCellType::PointIdentifier itkPts[3];
            for (int ii = 0; ii < 3; ++ii)
            {
                itkPts[ii] = static_cast<TriangleCellType::PointIdentifier>(pointIds[ii]);
            }
            tcell->SetPointIds(itkPts);
            c.TakeOwnership(tcell);
            mesh->SetCell(cellId, c);
            cellId++;
            pointIds[0] = pointIds[1];
            pointIds[1] = pointIds[2];
            pointIds[2] = cellPoints[t + 3];
        }
    }

    // then copy the normal triangles
    polygons->InitTraversal();
    polygonsIter = vtk::TakeSmartPointer(polygons->NewIterator());
    for (polygonsIter->GoToFirstCell(); !polygonsIter->IsDoneWithTraversal(); polygonsIter->GoToNextCell())
    {
        if (numberOfCellPoints != 3) // skip any non-triangle.
        {
            continue;
        }
        MeshType::CellAutoPointer c;
        TriangleCellType* t = new TriangleCellType;
        TriangleCellType::PointIdentifier itkPts[3];
        for (int ii = 0; ii < numberOfCellPoints; ++ii)
        {
            itkPts[ii] = static_cast<TriangleCellType::PointIdentifier>(cellPoints[ii]);
        }
        t->SetPointIds(itkPts);
        c.TakeOwnership(t);
        mesh->SetCell(cellId, c);
        cellId++;
    }

    std::cout << "Mesh " << std::endl;
    std::cout << "Number of Points = " << mesh->GetNumberOfPoints() << std::endl;
    std::cout << "Number of Cells = " << mesh->GetNumberOfCells() << std::endl;
       
    // Release vtk objects
    //
    reader->Delete();
         
    return mesh;
}
```

---

<div class="post-metadata">

**Author:** ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)\
**Post date:** [February 11, 2020, 8:36pm UTC](https://discourse.itk.org/t/deformable-registration-and-vtk-file-help/2699/11 "2020-02-11T20:36:40Z")

</div>

You might want to look at [ITKMeshToPolyData remote module](https://github.com/InsightSoftwareConsortium/ITKMeshToPolyData). Possibly also of interest [itkwidgets vtkPolyData.ipynb](https://github.com/InsightSoftwareConsortium/itkwidgets/blob/master/examples/vtkPolyData.ipynb). Maybe also [this VTK discussion](https://discourse.vtk.org/t/vtk-js-api-to-load-vtk-file-for-volume-rendering/2069/8?u=dzenanz). But [ITKVtkGlue remote module](https://github.com/InsightSoftwareConsortium/ITKVtkGlue) might be of most interest.

---

<div class="post-metadata">

**Author:** ![ShaneKirkbride](https://discourse.itk.org/user_avatar/discourse.itk.org/shanekirkbride/32/1225_2.png) [@ShaneKirkbride](https://discourse.itk.org/u/ShaneKirkbride)\
**Post date:** [February 11, 2020, 9:16pm UTC](https://discourse.itk.org/t/deformable-registration-and-vtk-file-help/2699/12 "2020-02-11T21:16:00Z")

</div>

Hi Dzenan,  
This [ITKMeshToPolyData remote module](https://github.com/InsightSoftwareConsortium/ITKMeshToPolyData). did work but I’m doubtful of a good conversion. I’m going to try the ITK method of reading in Dicom files.  
Thanks for the help so far!  
~Shane

---

<div class="post-metadata">

**Author:** ![matt.mccormick](https://discourse.itk.org/user_avatar/discourse.itk.org/matt.mccormick/32/7_2.png) [@matt.mccormick](https://discourse.itk.org/u/matt.mccormick)\
**Post date:** [February 19, 2020, 1:54pm UTC](https://discourse.itk.org/t/deformable-registration-and-vtk-file-help/2699/13 "2020-02-19T13:54:18Z")

</div>

> Both .vtk and .vtp fail.

The `.vtk` file format data type should work. `itk::MeshFileReader` also support other formats, e.g. `OFF`, `BYU`, `OBJ` – these could also be used.