# VMTK Centerline Tree Model Implementation Problem

**URL:** <https://discourse.slicer.org/t/vmtk-centerline-tree-model-implementation-problem/41760>\
**Category:** Development\
**Tags:** python, vmtk, vmtk-centerline\
**Created:** [February 18, 2025, 4:39pm UTC](https://discourse.slicer.org/t/vmtk-centerline-tree-model-implementation-problem/41760 "2025-02-18T16:39:04Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![yoonjh119](https://sea2.discourse-cdn.com/flex002/user_avatar/discourse.slicer.org/yoonjh119/32/79677_2.png) [@yoonjh119](https://discourse.slicer.org/u/yoonjh119)\
**Post date:** [February 18, 2025, 4:39pm UTC](https://discourse.slicer.org/t/vmtk-centerline-tree-model-implementation-problem/41760/1 "2025-02-18T16:39:04Z")

</div>

Hello, I am implementing the centerline function of VMTK so that it can be executed with Python code. Recently, I have completed the implementation of the network model and am implementing the tree model. After loading the nifty file, I am creating polydata and obtaining the endpoint coordinates from the 3D slicer and then trying to get the centerline tree model that is the same as the slicer by inputting it.

 ![image](https://us1.discourse-cdn.com/flex002/uploads/slicer/original/3X/1/3/13ba3e780188da1fedbfb09977a1d0c545e795ba.png)

The code used is as follows.

```auto
import nibabel as nib
import numpy as np
import vtk
from vmtk import vtkvmtk
import vmtk.vtkvmtkMiscPython as vtkvmtkMisc
from vmtk import vmtkscripts
from vtk.util import numpy_support
import matplotlib.pyplot as plt
import pandas as pd
from vtk.util.numpy_support import vtk_to_numpy
from vtk.numpy_interface import dataset_adapter as dsa
import json
import logging
import vmtk.vtkvmtkComputationalGeometryPython as vtkvmtkComputationalGeometry

class endPointsMarkupsNode:
    def __init__ (self):
        self.controlPoints = []
        self.selectedControlPoints = []

    def GetNumberOfControlPoints(self):
        return len(self.controlPoints)

    def AddControlPoint(self, position):
        self.controlPoints.append(position)
        self.selectedControlPoints.append(False) # Not selected by default

    def GetNthControlPointPosition(self, index, pos):
        if index < 0 or index >= len(self.controlPoints):
            raise IndexError("Control point index out of range")
        pos[:] = self.controlPoints[index]

    def GetNthControlPointSelected(self, index):
        if index < 0 or index >= len(self.controlPoints):
            raise IndexError("Control point index out of range")
        return self.selectedControlPoints[index]

    def SetNthControlPointSelected(self, index, selected):
        if index < 0 or index >= len(self.controlPoints):
            raise IndexError("Control point index out of range")
        self.selectedControlPoints[index] = selected

    def RemoveNthControlPoint(self, index):
        if index < 0 or index >= len(self.controlPoints):
            raise IndexError("Control point index out of range")
        del self.controlPoints[index]
        del self.selectedControlPoints[index]

def extractCenterline(surfacePolyData, endPointsMarkupsNode, curveSamplingDistance=1.0):
    """Compute centerline.
    This is more robust and accurate but takes longer than the network extraction.
    :param surfacePolyData:
    :param endPointsMarkupsNode:
    :return:
    """
    capDisplacement = 0.0
    surfaceCapper = vtkvmtkComputationalGeometry.vtkvmtkCapPolyData()
    surfaceCapper.SetInputData(surfacePolyData)
    surfaceCapper.SetDisplacement(capDisplacement)
    surfaceCapper.SetInPlaneDisplacement(capDisplacement)
    surfaceCapper.Update()

    if not endPointsMarkupsNode or endPointsMarkupsNode.GetNumberOfControlPoints() < 2:
        raise ValueError("At least two endpoints are needed for centerline extraction")

    tubePolyData = surfaceCapper.GetOutput()
    pos = [0.0, 0.0, 0.0]

    numberOfControlPoints = endPointsMarkupsNode.GetNumberOfControlPoints()
    foundStartPoint = False
    for controlPointIndex in range(numberOfControlPoints):
        if not endPointsMarkupsNode.GetNthControlPointSelected(controlPointIndex):
            foundStartPoint = True
            break

    sourceIdList = vtk.vtkIdList()
    targetIdList = vtk.vtkIdList()

    pointLocator = vtk.vtkPointLocator()
    pointLocator.SetDataSet(tubePolyData)
    pointLocator.BuildLocator()

    for controlPointIndex in range(numberOfControlPoints):
        isTarget = endPointsMarkupsNode.GetNthControlPointSelected(controlPointIndex)
        if not foundStartPoint and controlPointIndex == 0:
            isTarget = False
        endPointsMarkupsNode.GetNthControlPointPosition(controlPointIndex, pos)
        pointId = pointLocator.FindClosestPoint(pos)
        if isTarget:
            targetIdList.InsertNextId(pointId)
        else:
            sourceIdList.InsertNextId(pointId)

    centerlineFilter = vtkvmtkComputationalGeometry.vtkvmtkPolyDataCenterlines()
    centerlineFilter.SetInputData(tubePolyData)
    centerlineFilter.SetSourceSeedIds(sourceIdList)
    centerlineFilter.SetTargetSeedIds(targetIdList)
    centerlineFilter.SetRadiusArrayName("Radius")
    centerlineFilter.SetCostFunction('1/R') # this makes path search prefer go through points with large radius
    centerlineFilter.SetFlipNormals(False)
    centerlineFilter.SetAppendEndPointsToCenterlines(0)
    centerlineFilter.SetSimplifyVoronoi(False)
    centerlineFilter.SetCenterlineResampling(False)
    centerlineFilter.SetResamplingStepLength(curveSamplingDistance)
    centerlineFilter.Update()

    if not centerlineFilter.GetOutput():
        raise ValueError("Failed to compute centerline (no output was generated)")
    centerlinePolyData = vtk.vtkPolyData()
    centerlinePolyData.DeepCopy(centerlineFilter.GetOutput())

    if not centerlineFilter.GetVoronoiDiagram():
        raise ValueError("Failed to compute centerline (no Voronoi diagram was generated)")
    voronoiDiagramPolyData = vtk.vtkPolyData()
    voronoiDiagramPolyData.DeepCopy(centerlineFilter.GetVoronoiDiagram())

    return centerlinePolyData, voronoiDiagramPolyData

# Path to NIfTI file
input_nifti_01 = 'combined_mask_02.nii.gz'
msk_name_01 = 'output_mask_poly.vtk'

nifti_img = nib.load(input_nifti_01)
data = nifti_img.get_fdata()

reader = vtk.vtkNIFTIImageReader()
reader.SetFileName(input_nifti_01)
reader.Update()

affine = nifti_img.affine

# Create numpy arrays for start and end points (LPS coordinates)
startPoint = [-1.96, -153.29446411132812, 379.3500061035156]
endPoint = [-10.934102058410645, -171.96652221679688, 319.4288635253906]

# Read NIfTI file into VTK
reader = vtk.vtkNIFTIImageReader()
reader.SetFileName(input_nifti_01)
reader.Update()

cast = vtk.vtkImageCast()
cast.SetInputData(reader.GetOutput())
cast.SetOutputScalarTypeToFloat()
cast.Update()

mc = vtk.vtkMarchingCubes()
mc.SetInputData(cast.GetOutput())
mc.ComputeNormalsOn()
mc.SetValue(0, 0.5)
mc.Update()

transform = vtk.vtkTransform()
transform.Scale(-1.0, -1.0, 1.0) # RAS to LPS
if reader.GetQFormMatrix():
    transform.Concatenate(reader.GetQFormMatrix())
elif reader.GetSFormMatrix():
    transform.Concatenate(reader.GetSFormMatrix())

tpd = vtk.vtkTransformPolyDataFilter()
tpd.SetInputData(mc.GetOutput())
tpd.SetTransform(transform)
tpd.Update()

# Clean and triangulate
cleaner = vtk.vtkCleanPolyData()
cleaner.SetInputData(tpd.GetOutput())

triangleFilter = vtk.vtkTriangleFilter()
triangleFilter.SetInputConnection(cleaner.GetOutputPort())
triangleFilter.Update()
simplifiedPolyData = triangleFilter.GetOutput()

writer = vtk.vtkPolyDataWriter()
writer.SetInputData(simplifiedPolyData)
writer.SetFileName(msk_name_01)
writer.Update()

# Add start and end points to the markup node
endPointsNode = endPointsMarkupsNode()
endPointsNode.AddControlPoint(startPoint)
endPointsNode.SetNthControlPointSelected(0, False)
endPointsNode.AddControlPoint(endPoint)
endPointsNode.SetNthControlPointSelected(1, True)

# Extract centerline
centerlinePolyData, voronoiDiagramPolyData = extractCenterline(simplifiedPolyData, endPointsNode)

# Save results (e.g., as VTK files)
writer = vtk.vtkPolyDataWriter()
writer.SetFileName("centerline.vtk")
writer.SetInputData(centerlinePolyData)
writer.Write()

writer.SetFileName("voronoi_diagram.vtk")
writer.SetInputData(voronoiDiagramPolyData)
writer.Write()

print("Centerline extraction completed and saved to 'centerline.vtk' and 'voronoi_diagram.vtk'")

```

However, the following error occurs and it does not proceed.

 ![image](https://us1.discourse-cdn.com/flex002/uploads/slicer/original/3X/b/4/b407213c7c87d544102d063f5cb6520b0f93ce11.jpeg)  
Is there a solution to this? The same problem occurs not only in the extractCenterline function but also in the previous version, computecenterline. Thank you.