*2.1. Virtual Model*

Virtual models of the proximal femur were used, obtained by downloading the CT scans of two male patients from the open-source virtual library *"The Cancer Imaging Archive*" with references TCGA-VP-A878 [11] and Pelvic-Ref-009 [12] corresponding to the geometric case 1 (GC1) and 2 (GC2). The medical images had a slice thickness of 2 and 3 mm in

the axial plane, and 0.909 and 1 mm in the coronal and sagittal planes for GC1 and GC2, respectively, each image being 512 × 512 pixels. Both CT scans were imported into 3D Slicer® 4.10.2 (https://www.slicer.org, accessed on 3 January 2022) to segment the right femur and its cortical part (Figure 1A) using the threshold, level tracing, paint, erase and smooth tools. Then, the trabecular part was obtained through the logical operation of subtraction between the femur and its cortical bone (Figure 1B). The 3D Slicer® allowed us to export the segmentation of the femur, cortical and trabecular bone as meshes in STL format, and these files were imported into Meshmixer® 3.5 (Autodesk Inc., Mill Valley, CA, USA) for inspection, repair and smoothing. Finally, the virtual models were processed as solids in NX® 10 (Siemens PLM Software Solutions, Plano, TX, USA), matching their coordinate system with the femur's coordinate system so the axial, coronal and sagittal planes were the XY, XZ and YZ planes, respectively (Figure 1C). mm in the axial plane, and 0.909 and 1 mm in the coronal and sagittal planes for GC1 and GC2, respectively, each image being 512 × 512 pixels. Both CT scans were imported into 3D Slicer® 4.10.2 (https://www.slicer.org, accessed on 04 January 2022) to segment the right femur and its cortical part (Figure 1A) using the threshold, level tracing, paint, erase and smooth tools. Then, the trabecular part was obtained through the logical operation of subtraction between the femur and its cortical bone (Figure 1B). The 3D Slicer® allowed us to export the segmentation of the femur, cortical and trabecular bone as meshes in STL format, and these files were imported into Meshmixer® 3.5 (Autodesk Inc., Mill Valley, CA, US) for inspection, repair and smoothing. Finally, the virtual models were processed as solids in NX® 10 (Siemens PLM Software Solutions, Plano, TX, US), matching their coordinate system with the femur's coordinate system so the axial, coronal and sagittal planes were the XY, XZ and YZ planes, respectively (Figure 1C).

Virtual models of the proximal femur were used, obtained by downloading the CT scans of two male patients from the open-source virtual library *"The Cancer Imaging Archive*" with references TCGA-VP-A878 [11] and Pelvic-Ref-009 [12] corresponding to the geometric case 1 (GC1) and 2 (GC2). The medical images had a slice thickness of 2 and 3

*Materials* **2022**, *15*, x FOR PEER REVIEW 3 of 32

**2. Materials and Methods** 

*2.1. Virtual Model* 

**Figure 1.** (**A**) Segmentation of the right femur using 3D Slicer®. (**B**) Process to obtain the trabecular bone. (**C**) Femoral coordinate system. **Figure 1.** (**A**) Segmentation of the right femur using 3D Slicer®. (**B**) Process to obtain the trabecular bone. (**C**) Femoral coordinate system.

### *2.2. Elliptical adjustment app 2.2. Elliptical Adjustment App*

Trapezoidal, oval, elliptical and circular cross-sections have been used for the design of femoral stems. Previous studies [13,14] determined that the elliptical section produces a good stress distribution along the stem, allows its primary stability and improves its adaptability to different bone sections with changes in shape and size. Therefore, using Streamlit®, an open-source Python® library, an elliptical adjustment app was created to obtain the ellipse that best fits the bone section (https://github.com/solor5/elliptical\_adjustment\_app/blob/main/app.py, accessed on 04 January 2022). To utilize the application, the user must sample the bone section to be adjusted by employing the NX® point tool, and the coordinates of each point are exported in a DAT file (Figure 2), which is inserted into the app file uploader. The mathematical fundamentals that enable the elliptical fitting are explained below. Trapezoidal, oval, elliptical and circular cross-sections have been used for the design of femoral stems. Previous studies [13,14] determined that the elliptical section produces a good stress distribution along the stem, allows its primary stability and improves its adaptability to different bone sections with changes in shape and size. Therefore, using Streamlit®, an open-source Python® library, an elliptical adjustment app was created to obtain the ellipse that best fits the bone section (https://github.com/solor5/elliptical\_ adjustment\_app/blob/main/app.py, accessed on 3 January 2022). To utilize the application,the user must sample the bone section to be adjusted by employing the NX® point tool, and the coordinates of each point are exported in a DAT file (Figure 2), which is inserted into the app file uploader. The mathematical fundamentals that enable the elliptical fitting are explained below. *Materials* **2022**, *15*, x FOR PEER REVIEW 4 of 32

**Figure 2.** (**A**) Bone section sampling and DAT file. (**B**) Representation of the fitted curve and the elliptical adjustment of its orthogonal projection. diaphyseal part, have an elliptical shape, the conic provided by the app will be an ellipse. **Figure 2.** (**A**) Bone section sampling and DAT file. (**B**) Representation of the fitted curve and the elliptical adjustment of its orthogonal projection.

As a consequence of the position of the femur in NX® (Figure 1C), the bone section is located in an oblique plane perpendicular to the XZ; therefore, the X and Y coordinates of

If = ሼଵ, ଶ,…,ሽ is the set of points obtained from sampling the orthogonal projection of the bone section, it only includes the X and Y coordinates of the DAT file loaded in the app. The vector of coefficients must be adjusted to ; for this purpose, the algebraic distance () was used. This distance is widely employed because it simplifies the calculations and needs less computational resources [15]. Mathematically, it is obtained by replacing the coordinates of a point = (, ) in the polynomial; hence, if belongs

The least squares technique optimizes the fit by minimizing the square of the algebraic distance between the points and the curve, and it can be expressed as the squared norm of the product between the design matrix , which contains information

ଶ ଵ ଵ ଵ

ଶ ଶ ଶ ଶ

ଶ

⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮

= (, )ଶ 

ୀଵ

To avoid the trivial solution of = 0ത, the vector of coefficients is bounded [15]. Paton [16] analyzed the chromosome shape using a conic fit with a constraint of ‖‖ଶ = 1, avoiding all coefficients being zero. Therefore, using this constraint, the solution can be an ellipse, hyperbola or parabola; however, because the bone sections, especially in the

() = ሼ ∈ ℝଶ | (, ) = 0ሽ (1)

<sup>ଶ</sup> + +

ଶ ଵ ଵ 1

ଶ ଶ ଶ 1

ଶ 1 ⎦

⎥ ⎥ ⎤

= <sup>ଶ</sup> + + <sup>ଶ</sup> + + + = 0 (2)

<sup>ଶ</sup> + + + (3)

= ‖‖ଶ (5)

(4)

⎣ ⎢ ⎢ ⎢ ⎢ ⎡ ⎦ ⎥ ⎥ ⎥ ⎥ ⎤

൫, ()൯ = (, ) =

=

൫, ()൯

ୀଵ

⎣ ⎢ ⎢ ⎡ ଵ

ଶ

ଶ

to the ellipse, its distance will be 0.

of , and the vector.

(, ) = [ଶ ଶ 1] .

ynomial (), defined by a vector of coefficients ( = [ ]்):

As a consequence of the position of the femur in NX® (Figure 1C), the bone section is located in an oblique plane perpendicular to the XZ; therefore, the X and Y coordinates of the points in the DAT file (*p*) allow the elliptical adjustment of the orthogonal projection of this section (Figure 2A). This conic (*W*) is represented by an implicit second-order polynomial (*Q*), defined by a vector of coefficients (*v* = [*A B C D E F*] *T* ):

$$\mathcal{W}(v) = \left\{ p \in \mathbb{R}^2 \mid \mathcal{Q}(p, v) = 0 \right\} \tag{1}$$

$$Q(p,v) = \begin{bmatrix} \mathbf{x}^2 & \mathbf{x}y & y^2 & \mathbf{x} & y & 1 \end{bmatrix} \cdot \begin{bmatrix} A \\ B \\ C \\ D \\ E \\ F \\ F \end{bmatrix} = A\mathbf{x}^2 + B\mathbf{x}y + Cy^2 + D\mathbf{x} + Ey + F = 0 \tag{2}$$

If *P* = {*p*1, *p*2, . . . , *pn*} is the set of points obtained from sampling the orthogonal projection of the bone section, it only includes the X and Y coordinates of the DAT file loaded in the app. The vector of coefficients must be adjusted to *P*; for this purpose, the algebraic distance (*DA*) was used. This distance is widely employed because it simplifies the calculations and needs less computational resources [15]. Mathematically, it is obtained by replacing the coordinates of a point *p<sup>i</sup>* = (*x<sup>i</sup>* , *yi*) in the *Q* polynomial; hence, if *p<sup>i</sup>* belongs to the ellipse, its distance will be 0.

$$D\_A(p\_{\rm i}, \mathcal{W}(v)) = Q(p\_{\rm i}, v) = Ax\_i^2 + Bx\_i y\_i + Cy\_i^2 + Dx\_i + Ey\_i + F \tag{3}$$

The least squares technique optimizes the fit by minimizing the square of the algebraic distance between the *P* points and the *W* curve, and it can be expressed as the squared norm of the product between the design matrix *DP*, which contains information of *P*, and the *v* vector.

$$D\_p = \begin{bmatrix} x\_1^2 & x\_1 & y\_1 & y\_1^2 & x\_1 & y\_1 & 1\\ x\_2^2 & x\_2 & y\_2 & y\_2^2 & x\_2 & y\_2 & 1\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ x\_n^2 & x\_n & y\_n & y\_n^2 & x\_n & y\_n & 1 \end{bmatrix} \tag{4}$$

$$\min \sum\_{i=1}^{n} D\_A(P, W(v))^2 = \min \sum\_{i=1}^{n} Q(P, v)^2 = \min \left\| D\_P v \right\|^2 \tag{5}$$

To avoid the trivial solution of *v* = 06, the vector of coefficients is bounded [15]. Paton [16] analyzed the chromosome shape using a conic fit with a constraint of k*v*k <sup>2</sup> = 1, avoiding all coefficients being zero. Therefore, using this constraint, the solution can be an ellipse, hyperbola or parabola; however, because the bone sections, especially in the diaphyseal part, have an elliptical shape, the conic provided by the app will be an ellipse. The Lagrange multipliers allowed us to minimize the distance considering the constraint. Consequently, *L* is the Lagrange function to be optimized.

$$L = \|D\_P v\|^2 - \lambda \left(\|v\|^2 - 1\right) = v^T D\_P^T D\_P v - \lambda \left(v^T v - 1\right) \tag{6}$$

Equating the gradient of *L* with respect to *v* to 0 for minimizing the function gives:

$$
\nabla\_{\upsilon} L = 0 \iff 2D\_P^T D\_P \upsilon - 2\lambda \upsilon = 0 \tag{7}
$$

$$D\_P^T D\_P v = \lambda v \tag{8}$$

The optimization leads to the eigenvector problem; then, *λ* and *v* must be an eigenvalue and an eigenvector of *D<sup>T</sup> <sup>P</sup>DP*. If *<sup>D</sup><sup>T</sup> <sup>P</sup>DPv* = *λv*, Equation (5) will be:

$$\min \left\| D\_P \upsilon \right\|^2 = \min \upsilon^T D\_P^T D\_P \upsilon = \min \lambda \left\| \upsilon \right\|^2 = \min \lambda \tag{9}$$

As a result, the coefficient vector (*v*) that minimizes the algebraic distance will be the eigenvector of *D<sup>T</sup> <sup>P</sup>D<sup>P</sup>* corresponding to the smallest eigenvalue (*λ*). Once *v* is found, *Q* is defined; however, although CAD programs allow users to enter functions to draw a curve, this operation is often tedious. Therefore, the *W* ellipse can be defined with five parameters: the coordinates of its center (*xc*, *yc*), the largest (*R*) and smallest (*r*) radius and the angle (*α*) that rotates the curve counterclockwise. The coefficients of *v* are renamed:

$$v = \begin{bmatrix} ABCDEF \end{bmatrix}^T = \begin{bmatrix} a' \ \mathbf{2}b' \ c' \ \mathbf{2}d' \ \mathbf{2}e' \ f' \end{bmatrix}^T \tag{10}$$

The five parameters can be calculated using the following equations [17]:

$$\propto\_{c} = \frac{c'd' - b'e'}{b'^2 - a'c'} \tag{11}$$

$$y\_c = \frac{a'e' - b'd'}{b'^2 - a'c'} \tag{12}$$

$$R = \sqrt{\frac{2\left(a'e'^2 + c'd'^2 + f'b'^2 - 2b'd'e' - a'c'f'\right)}{\left(b'^2 - a'c'\right)\left[\sqrt{\left(a' - c'\right)^2 + 4b'^2} - \left(a' + c'\right)\right]}}\tag{13}$$

$$\sigma = \sqrt{\frac{2\left(a'e'^2 + c'd'^2 + f'b'^2 - 2b'd'e' - a'c'f'\right)}{\left(b'^2 - a'c'\right)\left[-\sqrt{\left(a' - c'\right)^2 + 4b'^2} - \left(a' + c'\right)\right]}}\tag{14}$$

$$\alpha = \begin{cases} \begin{array}{c} 0; \text{if } b' = 0 \text{ and } a' < c'\\ \frac{\pi}{2}; \text{if } b' = 0 \text{ and } a' > c' \end{array} \\\begin{array}{c} \frac{\arctan\left(\frac{2b'}{a'-c'}\right)}{2}; \text{if } b' \neq 0 \text{ and } a' < c' \end{array} \\\frac{\pi}{2} + \frac{\arctan\left(\frac{2b'}{a'-c'}\right)}{2}; \text{if } b' \neq 0 \text{ and } a' > c' \end{array} \end{cases} \tag{15}$$

The fitted curve of the bone section is the intersection between an elliptical cylinder, with the *W* curve as its directrix and the Z-axis as its generatrix, and the section plane *z* = *Gx* + *H* (Figure 2B), where the constants (*G*, *H*) are fitted using linear regression from the X and Z coordinates of each point of the DAT file. There are two ways to export the fitted curve from the app to NX®. The first one allows the user to obtain the points of the curve in DAT format by clicking on *Download DAT file* (Figure 2B), and then import these points to NX® and, with the spline tool, obtain the fitted curve. Likewise, the application provides a graph of the *W* curve containing the parameters that define it (*xc*, *yc*, *R*, *r*,*α*); these are introduced in the NX® ellipse tool and *W* is projected to the section plane to obtain the fitted curve (Figure 2B). The elliptical adjustment app allows the user to download the graph of the *W* curve (Figure 2B) and provides a 3D view of the fitted curves because the app can make several adjustments at the same time.
