# How I can find the center of the humeroulnar joint using 3d slicer ?

Hello everyone,
I am currently working on the morphology of the ulna and I need to calculate the biomechanical length of some ulnas by placing a landmark at the center of the trochlear notch (the center of the joint). I thought I could do this by creating a sphere that passes through different landmarks at the center of the notch, I use this code :

``````# Get markup node from scene
markups = 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: https://github.com/thompson318/scikit-surgery-sphere-fitting/blob/master/sksurgeryspherefitting/algorithms/sphere_fitting.py
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: https://github.com/thompson318/scikit-surgery-sphere-fitting/blob/master/sksurgeryspherefitting/algorithms/sphere_fitting.py
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)

# Fit a sphere to the markups fidicual points
markupsPositions = slicer.util.arrayFromMarkupsControlPoints(markups)
import numpy as np
# initial guess
center0 = np.mean(markupsPositions, 0)
fittingResult = fit_sphere_least_squares(markupsPositions[:,0], markupsPositions[:,1], markupsPositions[:,2], [center0[0], center0[1], center0[2], radius0])
[centerX, centerY, centerZ, radius] = fittingResult["x"]

# Create a sphere using the fitted parameters
sphere = vtk.vtkSphereSource()
sphere.SetPhiResolution(30)
sphere.SetThetaResolution(30)
sphere.SetCenter(centerX, centerY, centerZ)
sphere.Update()

# Add the sphere to the scene
modelsLogic = slicer.modules.models.logic()
model.GetDisplayNode().SetSliceIntersectionVisibility(True)
model.GetDisplayNode().SetSliceIntersectionThickness(3)
model.GetDisplayNode().SetColor(1,1,0)
``````

and when I tried I got this :

My ultimate goal is to calculate the length from the center of the joint to the end of the ulnar head. Could you please help me achieve this ?

Would this fit your requirements ? You can draw an arbitrary sphere if you install ExtraMarkups with the Extensions manager. Please note that you can also transform an arbitrary sphere using IGT extension.

1 Like

Thanks a lot for your answer It can resolve my probleme but I have to try it first and I really can find it on the extensions manager, neither with the name ExtraMarkups or SlicerExtraMarkups, can you helps me ?
And IGT I find it but donât undertand how it can help me âŠ

Youâre probably using Slicer stable 5.2. ExtraMarkups is available in Slicer preview that you can download similarly.

As for IGT, you can use the âCreate modelsâ module in this extension. Then in âData moduleâ, you can create a transform on this object, show the transform in GUI and manipulate it.

Tinker with one at a time to be less confused.

1 Like

Thank you very much !

I try it :

As you can see it work pretty well !
But in fact Iâm not sure itâs really the best fit sphere and if Iâm really in the center joint but honestly I think itâs enough for my purpose.
Is it possible to define de sphere with more than one point or no ?

Yes, but it implies complicating the code, meaning risks of errors and risks of crash, with little benefit. So no.

If you have the source images, you may reslice on the humeral trochlea and get a reliable sphere.

Ok I understand

I donât have the source image, I just have the 3D surface model.

Well it works for now so thank you very much but if anyone else (or even you) have another idea to better fit my requirement, Iâm open.

You can get the optimal best fit sphere the way you tried it at the very beginning. The only mistake was that you only placed 8 points, mostly in a single plane. If you place dozens of points, well dispersed on the entire joint surface then the automatically fit sphere should provide the correct solution.

I tried again and you have right, itâs work, thank you very much so I have two methods :

And this one :

I donât know if itâs really a probleme that the sphere donât fit only the keeled part of the joint. I will try to get the center of that sphere and compare both mesures

Thank you very much again, if anyone know how to create a landmarks at the center of that sphere iâm interested in !

PS: I already looked here : Get best fit sphere of humeral head - #6 by lassoan, But I donât understand.

Thank you for your help !

I suppose you are referring to the markups Shape node, in which case, you can use this function to get the center :

``````def getShapeCenter(markupsNode):
center = None
if markupsNode.GetClassName() != "vtkMRMLMarkupsShapeNode":
return center
if markupsNode.GetShapeName( ) != slicer.vtkMRMLMarkupsShapeNode.Sphere:
return center
p1 = [ 0.0 ] * 3
p2 = [ 0.0 ] * 3
markupsNode.GetNthControlPointPositionWorld(0, p1)
markupsNode.GetNthControlPointPositionWorld(1, p2)
return p1
else:
center = [ 0.0 ] * 3
for i in range(3):
center[i] = (p1[i] + p2[i]) / 2.0

return center
``````
1 Like

@chir.set it would be nice to add a âbest fitâ mode for the sphere markup that would compute least square fit sphere for all the control points. It could provide the center and radius as measurements.

1 Like

Yes it could be awsome to add this mode !

I try this code, but anything happen when I run it :

What am I doing wrong ?

Iâll get back for more info in a few days, not possible right now.

Itâs a function that you must call, passing it a Sphere markups object :

``````# After pasting the function in the Python console :
sh = slicer.util.getNode("SH")
center = getShapeCenter(sh)
print(center)
``````

where âSHâ is the name of a markups Shape node of type Sphere.

This page would be quite helpful for this. The solution is still in Python and it should be implemented in C++ to fit in ExtraMarkups. I would need to call np.linalg.lstsq from C++. Is it doable and how ?

You can call Python from C++ using the Python manager class. However, it would a bit nicer to use the optimization methods that are directly available in C++. ITK has several optimizer, both linear and non-linear. For example lbfgs-b is a good non-linear one. Iâve found a simple Powell optimizer example here: PlusLib/vtkPlusProbeCalibrationOptimizerAlgo.cxx at master Â· PlusToolkit/PlusLib Â· GitHub

This task is out of reach to me.

I donât see how to input point coordinates to the ITK optimizer, where it processes the input and where/how to get output, despite the example you provide. After many hours, I conclude itâs time to stop.

Python has a great advantage of higher level functions here. Perhaps it could be a slicer.util.bestFitSphereFromPoints function, in the same line as dataFrameFrom*, arrayFrom* functions in slicer.util. Using the Python manager via a temporary file is obviously not elegant from C++.