vtkIterativeClosestPointTransform - number of points in source and target

Support
,
1

Hi,

I am using vtkIterativeClosestPointTransform in a pipeline that detects radiopaque stereotactic frame fiducials in a CT/MRI. The geometry of the stereotactic frame is know. For a single CT volume there could be ~700-800 markups points, which forms the stereotactic frame localizers. The image below shows all the points that ae found within a CT volume.

leksell_mean_clusters
leksell_mean_clusters1536Ă—1291 147 KB

When running the ICP registration I convert all the CT markups points to polydata and generate an equal number of points for the frame. The top and bottom of each localizer bar have known coordinates and the bars have a known length (points on the same line in the image above).

My question: Do the source and target point clouds need have the same number of points? Realistically, I could have many more points in the target point cloud as it is generate in real-time based on the frame geometry. I have tried setting both point clouds to the same size as well as setting the target to have 4x as many points. I notice a substantial reduction in RMS error when the target cloud had 4x more points.

Is super-sampling the target point cloud frowned upon? I know it can be an issue if the target cloud has fewer points than the source but I couldn’t find anything about the target having many more points.

Thank you!

2

ICP do not require matching number of fixed and moving points. In fact, it registers a point cloud with a surface (I think line cells work, too). So, it should not be necessary to represent the template with many points, you just need to make sure that you have lines in the template’s polygonal mesh, not just points.

However, you probably don’t need to implement Z-frame registration from scratch. Slicer has an extension for automated Z-frame detection, subpixel-accuracy intersection point position estimattion, and optimal fitting: GitHub - SlicerProstate/SlicerZFrameRegistration: This extension can be used for the registration between intra-procedural image frame of reference with the transperineal biopsy or cryoablation templates You may need to tune it a bit for your images and the geometry of your frame.

If this extension is not flexible enough to specify different Z-frame geometries then it could be updated to use the general linear registration method developed by @mholden8. It can register any pattern made up of points, lines, and planes and guarantees to find the optimal match with a non-iterative method. As far as I remember, it is already available in one of the IGT extensions in Slicer.

1 Like
3

This was the exact solution. I modified the code to use vtkLine instead and it performs much better now. I just passed the vtkPolyData object to the ICP transform instead.

tubePolyData = vtk.vtkPolyData()
tubePolyData.SetPoints(framePoints)
tubePolyData.SetLines(frameLines)

icpTransform = vtk.vtkIterativeClosestPointTransform()
icpTransform.SetSource(sourceModel.GetPolyData())
icpTransform.SetTarget(tubePolyData)
...
icpTransform.Modified()
icpTransform.Update()

transform_matrix = vtk.vtkMatrix4x4()
transform_matrix.DeepCopy(icpTransform.GetMatrix())

Thank you for the quick response. I tried playing around with the ZFrameRegistration module but the stereotactic frames I am using are quite different. I have my own module I will be contributing for human stereotactic frame detection very soon.

1 Like
4

I have a follow-up question after performing the registration for multiple patients. Following the registration, I am computing the registration error and storing the euclidean distance in a list object:

cellId = vtk.mutable(0)
subId = vtk.mutable(0)
dist2 = vtk.mutable(0.0)
cellLocator = vtk.vtkCellLocator()
cellLocator.SetDataSet(inputTargetModel.GetPolyData())
cellLocator.SetNumberOfCellsPerBucket(1)
cellLocator.BuildLocator()
        
totalDistance = 0.0
accumPoints = []

sourcePoints = inputSourceModel.GetPoints()
for sourcePointIndex in range(sourcePoints.GetNumberOfPoints()):
    sourcePointPos = [0]*3
    sourcePoints.GetPoint(sourcePointIndex, sourcePointPos)
    tSourcePointPos = [0]*4
    sourcePointPos = sourcePointPos+[1]
    icpTransform.MultiplyPoint(sourcePointPos, tSourcePointPos)
    surfacePoint = [0]*3
    transformedSourcePointPos = tSourcePointPos[:3]
    cellLocator.FindClosestPoint(tSourcePointPos, surfacePoint, cellId, subId, dist2)
    totalDistance = totalDistance + math.sqrt(dist2)
    accumPoints.append(np.array(transformedSourcePointPos)-np.array(surfacePoint))

When I then examine the euclidean distance error for all the fiducials I found that only the fiducials that produce an oblique artifact in the CT volume have error in the z-axis (fiducial 2,5,8). If you look at the previous picture I posted, the diagonal bars in the N would be the bars producing the oblique artifact. The vertical bars that run dorsal to ventral have an error of 0 in the z-axis:

Fiducial X coordinate Y coordinate Z coordinate Root mean square error
1 -0.094 (0.209) 0.144 (0.067) 0.000 (0.000) 0.253
2 0.034 (0.133) -0.027 (0.069) -0.027 (0.069) 0.203
3 -0.041 (0.083) 0.108 (0.065) 0.000 (0.000) 0.184
4 -0.006 (0.069) -0.085 (0.090) 0.000 (0.000) 0.181
5 -0.074 (0.070) -0.224 (0.066) 0.074 (0.070) 0.296
6 -0.023 (0.063) -0.133 (0.075) 0.000 (0.000) 0.192
7 0.139 (0.076) 0.096 (0.059) 0.000 (0.000) 0.215
8 0.029 (0.134) -0.044 (0.068) -0.044 (0.068) 0.193
9 0.030 (0.202) 0.167 (0.077) 0.000 (0.000) 0.267

axial_frame
axial_frame1920Ă—1976 80.5 KB

Could this have something to do with the partial volume effect of CT? or have I done something incorrect during ICP registration error calculation?

Thank you in advance,
Greydon

5

Initializing mutable arrays using * operator is a bad practice because it does not create independent copies (it does not cause a problem for simple values, but causes very hard to detect errors when applied to arrays). Also, [0]*3 creates an integer array, so computation results will have quantization errors (rounded to integer values). Where did you see such initialization of point position vectors?

A safer and simpler syntax is np.zeros(3). You could use [0.0, 0.0, 0.0], too, but it is more verbose and does not allow you to perform numpy operations on the array directly (you need to convert to numpy array first).

6

I apologize, I wrote [0]*3 in my previous post but meant [0.]*3. However, I replaced this method of initialization with np.zeros(3). Thank you for this suggestion.

for sourcePointIndex in range(sourcePoints.GetNumberOfPoints()):
    sourcePointPos = np.zeros(3)
    sourcePoints.GetPoint(sourcePointIndex, sourcePointPos)
    tSourcePointPos = np.append(np.zeros(3),1)
    icpTransform.MultiplyPoint(np.append(sourcePointPos,1), tSourcePointPos)
    surfacePoint = np.zeros(3)
    cellLocator.FindClosestPoint(tSourcePointPos[:3], surfacePoint, cellId, subId, dist2)
    totalDistance = totalDistance + math.sqrt(dist2)
    accumPoints.append(tSourcePointPos[:3]-surfacePoint)

After re-running, the error did not change. I obtained the same registrations errors as reported in the table above. I am trying to rationalize the registration error.

The stereotactic frame is placed into a holder that is attached to the CT table. Assuming the gantry tilt is zero, the vertical fiducial bars would run parallel to the acquisition plane and the position of the artifact they produce would remain relatively consistent across CT slices. However, the position of the artifact produced by the diagonal fiducial bars would change between adjacent CT slices. Am I correct in concluding that the z-axis registration error is the result of the change in position of the diagonal bar artifact?