Get best fit sphere of humeral head

I am tyring to makes best fit sphere of humeral head.
So i makes 3 land marks and then import above code to python interactor.
but there is no sphere and some red text is shown…
plase some help…
thank you

slicer.util.arrayfromMarkupsControlPoints(markups) has an undefined “markups” object. This variable must be a reference to your point list.

markups=slicer.util.getNode("F") is example code of getting a point list with a current name of “F” which is a common default name for it.

1 Like

Thank you…
That`s works.
I have some more questions.

  1. Can i get center point of that sphere?
  2. Is there any method that one plane move to crossing the point which made by question 1. (Center ponints of sphere)
  3. Is there any method that calculate the distance between above 2 planes.

Yes, sure. You have set it yourself in the script: sphere.SetCenter(centerX, centerY, centerZ)

You can create a plane at any position. Add a new plane node, set its position to the sphere centerpoint, and set a normal.

You can use basic computational geometry methods for computing distances, angles between points, lines, planes, project points, lines, determine intersections, etc. There are a couple of examples in the script repository to see how to get inputs and show outputs in Slicer, and you can learn the rest by typing your question into Google search (for example python distance between two planes) and reading through the search results. You will get both detailed mathematical explanation and complete Python implementations (typically few lines of numpy code). If you really get stuck then you can ask here, but it could be useful for you to spend some with reading and experimenting.

thank you…
I have some more question…
I can`t find by google search…

  1. Can i get center point of that sphere? → I`m trying to find center point of sphere aleady made by humeral head fiduciel not set a center point before making sphere.

  2. I want make plane which are exactly parallel to one plane (Aleady made by GT foot print) and crossing center points of sphere.

  3. can i subtract plane by overlapped with humerus bone defect?
    I want measure the depth of defect. image below

The script computes the centerpoint position and stores it in centerX, centerY, centerZ.

The plane orientation is defined by the plane normal vector. You can create parallel planes by setting the same plane normal vector.

You can use Dynamic Modeler to cut a model with planes.

For making measurements you probably don’t need to actually cut a model, but it is enough to view a cross-section in a slice view. You can enable display of slice intersections and then rotate the slice intersections; or you can enable reformat widget to translate/rotate a slice plane in 3D.

I’m experiencing the same problem:

# Get markup node from scene
>>> pointListNode = slicer.util.getNode("F")
>>> from scipy.optimize import least_squares
>>> import numpy
>>> def fit_sphere_least_squares(x_values, y_values, z_values, initial_parameters, bounds=((-numpy.inf, -numpy.inf, -numpy.inf, -numpy.inf),(numpy.inf, numpy.inf, numpy.inf, numpy.inf))):
...   """
...   Source:
...   Uses scipy's least squares optimisor to fit a sphere to a set
...   of 3D Points
...   :return: x: an array containing the four fitted parameters
...   :return: ier: int An integer flag. If it is equal to 1, 2, 3 or 4, the
...           solution was found.
...   :param: (x,y,z) three arrays of equal length containing the x, y, and z
...           coordinates.
...   :param: an array containing four initial values (centre, and radius)
...   """
...   return least_squares(_calculate_residual_sphere, initial_parameters, bounds=bounds, method="trf", jac="3-point", args=(x_values, y_values, z_values))
>>> def _calculate_residual_sphere(parameters, x_values, y_values, z_values):
...   """
...   Source:
...   Calculates the residual error for an x,y,z coordinates, fitted
...   to a sphere with centre and radius defined by the parameters tuple
...   :return: The residual error
...   :param: A tuple of the parameters to be optimised, should contain [x_centre, y_centre, z_centre, radius]
...   :param: arrays containing the x,y, and z coordinates.
...   """
...   #extract the parameters
...   x_centre, y_centre, z_centre, radius = parameters
...   #use numpy's sqrt function here, which works by element on arrays
...   distance_from_centre = numpy.sqrt((x_values - x_centre)**2 + (y_values - y_centre)**2 + (z_values - z_centre)**2)
...   return distance_from_centre - radius
>>> # Fit a sphere to the markups fidicual points
>>> markupsPositions = slicer.util.arrayFromMarkupsControlPoints("F")
Traceback (most recent call last):
  File "<console>", line 1, in <module>
  File "D:\gsirc\AppData\Local\NA-MIC\Slicer 4.11.20210226\bin\Python\slicer\", line 1808, in arrayFromMarkupsControlPoints
    numberOfControlPoints = markupsNode.GetNumberOfControlPoints()
AttributeError: 'str' object has no attribute 'GetNumberOfControlPoints'
>>> import numpy as np
>>> # initial guess
>>> center0 = np.mean(markupsPositions, 0)
Traceback (most recent call last):
  File "<console>", line 1, in <module>
NameError: name 'markupsPositions' is not defined
>>> radius0 = np.linalg.norm(np.amin(markupsPositions,0)-np.amax(markupsPositions,0))/2.0
Traceback (most recent call last):
  File "<console>", line 1, in <module>
NameError: name 'markupsPositions' is not defined
>>> fittingResult = fit_sphere_least_squares(markupsPositions[:,0], markupsPositions[:,1], markupsPositions[:,2], [center0[0], center0[1], center0[2], radius0])
Traceback (most recent call last):
  File "<console>", line 1, in <module>
NameError: name 'markupsPositions' is not defined
>>> [centerX, centerY, centerZ, radius] = fittingResult["x"]
Traceback (most recent call last):
  File "<console>", line 1, in <module>
NameError: name 'fittingResult' is not defined
>>> # Create a sphere using the fitted parameters
>>> sphere = vtk.vtkSphereSource()
>>> sphere.SetPhiResolution(30)
>>> sphere.SetThetaResolution(30)
>>> sphere.SetCenter(centerX, centerY, centerZ)
Traceback (most recent call last):
  File "<console>", line 1, in <module>
NameError: name 'centerX' is not defined
>>> sphere.SetRadius(radius)
Traceback (most recent call last):
  File "<console>", line 1, in <module>
NameError: name 'radius' is not defined
>>> sphere.Update()
>>> # Add the sphere to the scene
>>> modelsLogic = slicer.modules.models.logic()
>>> model = modelsLogic.AddModel(sphere.GetOutput())
>>> model.GetDisplayNode().SetSliceIntersectionVisibility(True)
>>> model.GetDisplayNode().SetSliceIntersectionThickness(3)
>>> model.GetDisplayNode().SetColor(1,1,0)

Can anybody provide some help?

F is the a name of a node. slicer.util.arrayFromMarkupsControlPoints requires a node object as input. You can get the node object from a node name by calling getNode():

markupsPositions = slicer.util.arrayFromMarkupsControlPoints(slicer.util.getNode("F"))

1 Like

Thank you for the help. It worked smoothly.

1 Like

Hi, I succesfully extratced the coordinates and radius and I assume, they are in mm, right?

Yes, length unit is millimeters by default. It is always recommended to validate your entire workflow by testing on some object of known size.