New algorithm for curved extrusion to design the custom titanium plate

[mauigna06]

Hi. I would like you to see the results of this new algorithm dedicated to making the mesh for the base of the custom titanium plate compatible with dental implants. Here is an example scene: https://gofile.io/d/7alZoU

You can load your sample mandible CT to see how well it fits the reconstruction. This is a prototype version and it can be improved. For example, the extrusion is not always normal to the mandible reconstruction (because it only depends on the GetCurvePointToWorldTransformAtPointIndex of the curve). In the example it works well but to guarantee that the algorithm works every time another user input would be needed.

I would like to know your opinions about how good it looks from the clinical point of view (@cmfsx) and the engineering point of view (@lassoan). Andras, I tried to smooth it using the surfaceToolbox but I was unsuccessful. Does the mesh need to be subdivided first? I tried that and it worked kind of better but the mesh smoothed heterogeneously, no good-enough results.

The algorithm needs as input a polygon to be extruded defined in code (it's name is startingPolygonPoints) and it needs a a curve called "MarkupsCurve" to extrude along created by the user

Here is the code to execute on the python console if you want to try it:

import numpy as np

originalCurve = getNode('MarkupsCurve')

myCurve = slicer.mrmlScene.CreateNodeByClass("vtkMRMLMarkupsCurveNode")
myCurve.SetName("C2")
slicer.mrmlScene.AddNode(myCurve)
slicer.modules.markups.logic().AddNewDisplayNodeForMarkupsNode(myCurve)

points = vtk.vtkPoints()
curvePointsArray = slicer.util.arrayFromMarkupsControlPoints(originalCurve)
vtkPointsData = vtk.util.numpy_support.numpy_to_vtk(curvePointsArray, deep=1)
points.SetNumberOfPoints(len(curvePointsArray))
points.SetData(vtkPointsData)
myCurve.SetControlPointPositionsWorld(points)

myCurve.ResampleCurveWorld(10)

displayNode = myCurve.GetDisplayNode()
displayNode.UseGlyphScaleOff()
displayNode.SetGlyphSize(1.5)
displayNode.SetCurveLineSizeMode(1)
displayNode.SetLineDiameter(0.5)

xLength = 2.5
yLength = 7

startingPolygonPoints = np.array(
            [
                [xLength/2, yLength/2, 0],
                [-xLength/2, yLength/2, 0],
                [-xLength/2, -yLength/2, 0],
                [xLength/2, -yLength/2, 0]
            ]
        )

cell_array = vtk.vtkCellArray()
points = vtk.vtkPoints()
point_id = 0

startIndex = 0
curveMatrix = vtk.vtkMatrix4x4()
myCurve.GetCurvePointToWorldTransformAtPointIndex(startIndex,curveMatrix)
myCurveX = np.array([curveMatrix.GetElement(0,0),curveMatrix.GetElement(1,0),curveMatrix.GetElement(2,0)])

curveToWorldTransform = vtk.vtkTransform()
curveToWorldTransform.PostMultiply()
curveToWorldTransform.Concatenate(curveMatrix)
curveToWorldTransform.Translate(-myCurveX*xLength/2)

firstTransformedPolygonPoints = []
for i in range(len(startingPolygonPoints)):
  transformedPolygonPoint = np.zeros(3)
  curveToWorldTransform.TransformPoint(startingPolygonPoints[i], transformedPolygonPoint)
  firstTransformedPolygonPoints.append(transformedPolygonPoint)

polygon = vtk.vtkPolygon()
polygon.GetPointIds().SetNumberOfIds(4)
for i in range(4):
    points.InsertNextPoint(firstTransformedPolygonPoints[i])
    polygon.GetPointIds().SetId(i, point_id)
    point_id += 1

cell_array.InsertNextCell(polygon)

curvePoints = slicer.util.arrayFromMarkupsCurvePoints(myCurve)

for j in range(1,len(curvePoints)):
  curvePoint = curvePoints[j]
  #
  secondTransformedPolygonPoints = []
  #
  closestCurvePoint = [0,0,0]
  closestCurvePointIndex = myCurve.GetClosestPointPositionAlongCurveWorld(curvePoint,closestCurvePoint)
  #
  curveMatrix = vtk.vtkMatrix4x4()
  myCurve.GetCurvePointToWorldTransformAtPointIndex(closestCurvePointIndex,curveMatrix)
  myCurveX = np.array([curveMatrix.GetElement(0,0),curveMatrix.GetElement(1,0),curveMatrix.GetElement(2,0)])
  #
  curveToWorldTransform = vtk.vtkTransform()
  curveToWorldTransform.PostMultiply()
  curveToWorldTransform.Concatenate(curveMatrix)
  curveToWorldTransform.Translate(-myCurveX*xLength/2)
  #
  for i in range(len(startingPolygonPoints)):
    transformedPolygonPoint = np.zeros(3)
    curveToWorldTransform.TransformPoint(startingPolygonPoints[i], transformedPolygonPoint)
    secondTransformedPolygonPoints.append(transformedPolygonPoint)

  for i in range(4):
    points.InsertNextPoint(secondTransformedPolygonPoints[i])
    point_id += 1

  for k in range(len(firstTransformedPolygonPoints)):
    polygon = vtk.vtkPolygon()
    polygon.GetPointIds().SetNumberOfIds(3)
    polygon.GetPointIds().SetId(0, k + point_id - len(secondTransformedPolygonPoints))
    polygon.GetPointIds().SetId(2, k + point_id - 2*len(secondTransformedPolygonPoints))
    if k!=3:
      polygon.GetPointIds().SetId(1, k + 1 + point_id - len(secondTransformedPolygonPoints))
    else:
      polygon.GetPointIds().SetId(1, point_id - len(secondTransformedPolygonPoints))

    cell_array.InsertNextCell(polygon)
    #
    polygon = vtk.vtkPolygon()
    polygon.GetPointIds().SetNumberOfIds(3)
    polygon.GetPointIds().SetId(0, k + point_id - 2*len(secondTransformedPolygonPoints))
    if k!=3:
      polygon.GetPointIds().SetId(1, k + 1 + point_id - len(secondTransformedPolygonPoints))
      polygon.GetPointIds().SetId(2, k + 1 + point_id - 2*len(secondTransformedPolygonPoints))
    else:
      polygon.GetPointIds().SetId(1, point_id - len(secondTransformedPolygonPoints))
      polygon.GetPointIds().SetId(2, point_id - 2*len(secondTransformedPolygonPoints))

    cell_array.InsertNextCell(polygon)

  firstTransformedPolygonPoints = secondTransformedPolygonPoints.copy()
  #
  if j == len(curvePoints)-1:
    polygon = vtk.vtkPolygon()
    polygon.GetPointIds().SetNumberOfIds(4)
    for i in range(4):
        polygon.GetPointIds().SetId(i, point_id-4)
        point_id += 1

    cell_array.InsertNextCell(polygon)

polydata = vtk.vtkPolyData()
polydata.SetPoints(points)
polydata.SetPolys(cell_array)

triangleFilter = vtk.vtkTriangleFilter()
triangleFilter.SetInputData(polydata)
triangleFilter.Update()

extrusionModel = slicer.mrmlScene.CreateNodeByClass("vtkMRMLModelNode")
slicer.mrmlScene.AddNode(extrusionModel)
extrusionModel.CreateDefaultDisplayNodes()
extrusionModel.SetAndObservePolyData(triangleFilter.GetOutput())

[cmfsx]

It looks good as an initial step to start. It needs to have a bevel (smooth edges) so the soft tissues are not irritated. And of course, planning how to put the holes. When it comes to placing screw holes ideally it should not be just straight holes but rather in this shape, this is because in the inside surface contact with the bone surface should be minimal.

[mauigna06]

It looks good as an initial step to start. It needs to have a bevel (smooth edges) so the soft tissues are not irritated.

Yeah. It needs to be smoothed but it would be better to smooth it after the screw holes are created. I had problems making smoothing filters work right till now.

And of course, planning how to put the holes. When it comes to placing screw holes ideally it should not be just straight holes but rather in this shape, this is because in the inside surface contact with the bone surface should be minimal.

I didn't understand what inside surface you are referring to. However, the screw holes shouldn't be difficult to create making boolean operations I think a boolean difference of a cylinder and a cone (or another cylinder) could get the shape you showed on the picture @cmfsx for the holes.

It's great that you approved the shape of the algorithm-created fixing plate that was the most difficult thing to do. The smoothing and the holes should be solved soon. Do you think will need guides to help the surgeon make the shape of the reconstruction with the fibula pieces and put the custom fixing plate?

[cmfsx]

Do you think will need guides to help the surgeon make the shape of the reconstruction with the fibula pieces and put the custom fixing plate? I am not sure I understand the question exactly. If I understand it properly, surgeons will not need a another guide top put the reconstruction plate (custom made).

[mauigna06]

Consider use of this example and this filter

[mauigna06]

Half-toroids with rounded corners (positioned over the neomandible) could be interpolated using a spline that joins them together with rounded rectangles

[mauigna06]

Predictive holes (that's the custom plate fixation holes) should be added to the mandible surgical guide

[mauigna06]

Half-toroids with rounded corners (positioned over the neomandible) could be interpolated using a spline that joins them together with rounded rectangles

About doing the interpolation see LERP and smoothMax function used on game development

[mauigna06]

Consider using vtkSmoothPolyDataFilter with setSource for neomandible as input

[mauigna06]

related #117