Defining a new coordinate system using markups plane node

Support
, ,
2

The axes of a MarkupsPlane node are always orthogonal, even if the triangle of points used to define the plane is not a right triangle. If you have a MarkupsFiducial node with your three defining points, this function will produce a MarkupsPlane node

def makePlaneMarkupFromFiducial(self, FNode, planeName):
    """Create MarkupsPlane using first three control points of the input fiducial node.
    """
    if FNode.GetNumberOfControlPoints()<3:
      logging.warning('Not enough control points to create plane markup!')
      return
    planeNode = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLMarkupsPlaneNode', planeName)
    for cpIdx in range(3):
      pos = vtk.vtkVector3d()
      FNode.GetNthControlPointPositionWorld(cpIdx, pos)
      planeNode.AddControlPointWorld(pos)

    return planeNode

The resulting orthogonal axes, X, Y and Z are determined from the three input points as follows:

  • the positive X axis direction is the direction from point 1 to point 2
  • the positive Z axis is the cross product of the X-axis direction vector with the direction vector from point 1 to point 3
  • the positive Y axis is the cross product of the Z-axis and the X-axis

Conceptually, the Z-axis is perpendicular to the plane defined by the three points, and the positive Z-axis points such that, looking back at the three points from the positive Z side, point 1, point 2, and point 3 are counterclockwise (because the cross product follows a right hand rule). The positive Y-axis direction is the projection of the point1 to point2 vector onto the vector perpendicular to the X-Z plane; the cross product of the positive X-axis and positive Z-axis is another way to find this direction.

In terms of the display of the MarkupsPlane in the 3D view, it is shown as a rectangle in space. Point1 is the center of the rectangle, and the edges of the rectangle are parallel to the X and Y axes as defined above. The rectangle extends in the positive X direction far enough to touch Point2, which will always be at the midpoint of this edge because the X direction points directly from Point1 to Point2. The rectangle extends the same distance on the other side of Point1, keeping Point1 in the center of the rectangle. In the positive Y direction, the rectangle extends far enough so that the edge touches Point3. Note that Point3 will NOT necessarily be at the midpoint of this edge, because the angle Point2-Point1-Point3 may not be a right angle. The rectangle extends in the negative Y direction exactly as far as it did in the positive Y direction, keeping Point1 in the center of the rectangle.

I had to work these details out by looking at the code for MarkupsPlane nodes, because, as far as I can tell there isn’t really documentation for them yet.

Here is sample code which will work in the Python interactor using the point of interest you listed in the other thread:

fiducialNodeForBoneBlock1 = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLMarkupsFiducialNode', 'BoneBlock1')
fiducialNodeForBoneBlock1 .SetNthControlPointPositionWorld(0, -110.676, -105.632, 31.2008)

The 0 which is the first input to SetNthControlPointPositionWorld() is the point index, which tells which control point you want to set the position of, numbered starting at 0. So, the first control point has index 0, the second point has index 1, the third has index 2, etc. In my description above, Point1, the center point for the plane, would have index 0; Point2 would have index 1, and Point3 would have index 2.

In terms of a workflow, you could run the two lines above defining a fiducial node for the first bone block and placing the center control point, and then interactively place point 2 and point 3 on landmarks using the Markups module. Then you could define fiducialNodeForBoneBlock2 in the same way, just using the centroid for block 2 instead of the one for block 1, place the point 2 and point 3 landmarks interactively (note that these must be clicked in the same order as for the first block, or else the plane node derived from them will be inverted relative to block1). Lastly, you can use the makePlaneMarkupFromFiducial() function defined above to make plane nodes for each block:

planeNodeForBlock1 = makePlaneMarkupFromFiducial( fiducialNodeForBoneBlock1, 'Block1_Plane')
planeNodeForBlock2 = makePlaneMarkupFromFiducial( fiducialNodeForBoneBlock2, 'Block2_Plane')

Hope that helps!

5 Likes