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)
return distance_from_centre - radius
# Fit a sphere to the markups fidicual points
markupsPositions = slicer.util.arrayFromMarkupsControlPoints(markups)
import numpy as np
# initial guess
center0 = np.mean(markupsPositions, 0)
radius0 = np.linalg.norm(np.amin(markupsPositions,0)-np.amax(markupsPositions,0))/2.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.SetRadius(radius)
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)
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 ?
Thank you for your answers !
