Deformable Registration and .vtk file Help

1

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

2

Hi Shane,

To read the mesh, use itk::MeshFileReader. To register the meshes, use itk::PointSetToPointSetMetricv4.

HTH,
Matt

1 Like
3

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

4

Here is an example. You need to modify it to read a mesh, then get mesh’s vertices (instead of reading a list of points “manually”).

1 Like
5

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

6

ITK does not support VTK’s XMLPolyData file format. We only support the legacy file format (usual extension .vtk). I think vtkPolyDataWriter is the right class to use.

1 Like
7

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
8

I did find this: 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

9

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

1 Like
10

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;
}
11

You might want to look at ITKMeshToPolyData remote module. Possibly also of interest itkwidgets vtkPolyData.ipynb. Maybe also this VTK discussion. But ITKVtkGlue remote module might be of most interest.

12

Hi Dzenan,
This ITKMeshToPolyData remote module. 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

13

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.