Defining a new coordinate system using markups plane node

Continuing the discussion from Creating a new coordinate system:

Hello all, I would like to continue this discussion on defining a coordinate system using the markups plane node, as suggested by Andras in the forum linked above.
I had a couple of questions about this process:

  1. How do I make sure that my axes that I set are orthogonal to each other?
  2. How do I use the markupsNode.SetNthControlPointPositionWorld(i, r, a, s) method to place a fiducial at a specific coordinate point in the global system?
    a. It would be great if I could see a screenshot of the Python interactor, where the code has run successfully.
    b. What is a pointIndex?

Thank you very much for your help!

Michele

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

Hello,

Sorry to revive this topic, however i have been trying to create a plane from 3 fiducial points using the code you’ve shown here. Right now i’m creating a 3 point plane and editing the control points manually in the properties.

I have an extremely limited knowledge of python, please bear with me.
I have a fiducial node with 3 points. What are the steps i should take from this point?

When i run the code i get the following error:

File “”, line 1
return planeNode
IndentationError: unexpected indent

If someone could help me it would be greatly appreciated

If you cut and paste into the python interactor, it doesn’t like blank lines in the middle of function definitions, so that’s what gave you the indentation error above. Also, the default plane type changed since that post was written (it used to be 3points and is now pointNormal), so the function needed to be updated anyway before it would work now.

def makePlaneMarkupFromFiducial(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)
    planeNode.SetPlaneType(planeNode.PlaneType3Points)
    for cpIdx in range(3):
      pos = vtk.vtkVector3d()
      FNode.GetNthControlPointPositionWorld(cpIdx, pos)
      planeNode.AddControlPointWorld(pos)
    return planeNode

This should now work just fine from the python interactor, just cut and paste and hit enter twice. If your fiducial point list is named ‘myPoints’, you could then use the first three to define a markups plane node by typing

makePlaneMarkupFromFiducial(getNode('myPoints'), 'myPlane')

into the interactor and hitting enter.

1 Like

I follow the above process, but the results are as follows. Can you give me some help
mmexport1667962279153
mmexport1667962284060

I success ,thank you

1 Like

I think the issue here was again line breaks. The python console won’t allow blank lines within a function definition (because it interprets it as ending the function, and then thinks the next indented line is an unexplained and incorrect indentation), and it also wants at least one full blank line at the end of a function definition before you try to use the function (I don’t know why this makes sense, but it seems to be how it works). In this case, you didn’t leave a blank line between the return statement and the call to the function, so it gave an error. It sounds like you got things working, I just figured I’d leave an explanation to help the next person who might try this.

1 Like