As musculoskeletal illnesses continue to increase, practical computerised muscle modelling is crucial. This paper addresses this concern by proposing a mathematical model for a dynamic 3D geometrical surface representation of muscles using a Radial Basis Function (RBF) approximation technique. The objective is to obtain a smoother surface while minimising data use, contrasting it from classical polygonal (e.g. triangular) surface mesh models or volumetric (e.g. tetrahedral) mesh models. The paper uses RBF implicit surface description to describe static surface generation and dynamic surface deformations based on its spatial curvature preservation during the deformation. The novel method is tested on multiple data sets, and the experiments show promising results according to the introduced metrics.

Computerised muscle modelling garners increasing attention with the rising prevalence of musculoskeletal illnesses (Cieza

Muscle modelling has a rich history, with Hill (

While these methods often operate in one dimension, moving to higher dimensions, especially in 2D space, evolved crucial because the 1D model cannot accurately define the original shape, producing errors up to 75% (Valente

The open, ultimate problem is developing an accurate, simple, fast, and smooth musculoskeletal model using the least possible parameters. All of the state-of-the-art methods lack one or the other or give some compromises.

The paper primarily emphasises a detailed description of two-dimensional techniques, allowing faster calculations while keeping the ability to infer internal muscle composition. Two main objectives include achieving more immediate model deformation with fewer updated parameters and data reduction. The goal is achieved using a mathematically described RBF implicit surface geometric model, providing infinitely differentiable smooth surfaces and using fewer parameters than the “traditional” triangular mesh approaches.

The contribution of this paper lies not just in the overall RBF model but also in a novel metric to determine the approximation accuracy presented. This metric can also be used in different research fields where two surfaces must be compared. The RBF solver is another contribution, and its usage can also be extended to applications where some smooth objects (without sharp edges) need to be approximated, such as fluid dynamics, aerodynamics, computer graphics and more.

Section

Having established the theoretical background, we now advance to a detailed investigation of some techniques relative to muscle modelling.

The origin of a musculoskeletal model is deeply rooted in empirical data. The primary raw materials for these models are medical scans, such as Magnetic Resonance Imaging and Computed Tomography screenings, originating from living subjects and cadavers. Recently, the processing has evolved, transitioning from semi-automatic to fully automatic techniques. This progression involves the segmentation of images, a critical step where the relevant structures, such as bones, joints and muscles, are isolated and extracted from the rest of the image.

Following segmentation, surface models of bones are created. This process is not just about assembling a skeletal structure; it involves the intricate specification of joints, setting constraints on the degrees of freedom within these joints. Here, e.g. the STAPLE algorithm by Modenese and Renault (

Triangular meshes are the most common way to approximate bone surfaces by connecting triangles to define the object’s overall shape. The methodology for obtaining a triangular muscle model from a person entails a comprehensive process that begins with acquiring high-resolution, segmented medical images, commonly from MRI or CT scans. These images provide a detailed view of the bones (and surrounding tissues). The next step involves the extraction of a triangular mesh from these segmented images, a task typically accomplished using algorithms such as Marching Cubes (Lorensen and Cline,

Following the initial mesh extraction, the mesh undergoes a series of refinement procedures to enhance its quality and anatomical accuracy. These procedures include the removal of non-manifold vertices and edges to eliminate anomalies that do not adhere to the criteria of a well-defined surface. Concurrently, any gaps or holes in the mesh are meticulously filled to ensure continuity of the bone surface, a critical step for maintaining anatomical fidelity. The mesh is further refined by optimising the shapes of the triangles to more accurately align with the bone’s contours and by reducing the mesh to lower its complexity while keeping the model’s overall shape and intricate details.

Additionally, Laplacian smoothing is applied to the mesh. This process adjusts the position of each vertex based on the average of its neighbouring vertices, effectively smoothing out irregularities and noise, resulting in a uniform and more realistic representation of the bone. The refined mesh is then rigorously reviewed and compared against the original segmented images by medical professionals to validate its accuracy, ensuring the model precisely mirrors the anatomical structure.

While the forces used on bones during the movement can induce deformations, these deformations are generally so minute that, for practical purposes, they can be disregarded. This assumption yields the bones as rigid bodies in motion. This assumption is in contrast to the muscle-tendon units.

Joint data about the muscles and tendons (together form muscle-tendon units – MTU) is acquired because it is difficult to distinguish the muscle and the tendon apart from the imaging; otherwise, the data can be obtained in the same fashion as in the case of bones. However, the problem is that these tissues are less visible in the imaging. Fortunately, the segmentation can also be automated, but more complex strategies are required, with machine learning methods emerging as promising approaches (Goyanes

If the deformation is discussed, muscles behave differently than bone structures; they exhibit elastic deformations, presenting a significant challenge in modelling their behaviour accurately to the known movement of bones. Various models attempt to address this, often employing many oversimplifications.

Also, in the real-world scenario, muscle deformation results from the contraction or relaxation of muscle fibres, driven by the sliding of actin and myosin filaments within muscle fibres. This contraction leads to changes in muscle shape, causing deformation. In inverse kinematics, muscle deformation is pretended to be caused by adjacent bone movement, posing the challenge of figuring out the shape of unknown parts of the muscle model. If we do not consider measurements of the physiological signals (e.g. electromyography), it is required to perform the inverse kinematics because the bone displacements are known; however, the shape of the muscles during the movement is not (due to the problematic data acquisition).

The first models created were one-dimensional and were formed by a straight line, polyline, or curve and their multiples.

The simplest models might represent a muscle as a straight line connecting two points on different bones, blatantly ignoring any potential intersection with the bone itself. This approach relies heavily on the assumption that the attachment points are accurately chosen (Kohout and Cervenka,

It’s also imperative to represent muscles not just as mere lines or curves but as substantial higher-dimensional models (Kedadria

Geometrical 3D models have been explored, utilising the fact that a large set of muscle fibres creates the muscle. Modelling approaches using mass-spring systems (Janák and Kohout,

Despite muscles lying in the 3D space, modelling them in lower dimensions is viable. Straight lines or polylines can be employed in one dimension but inaccurately. In two dimensions, already obtained triangular surfaces are limited in the original object description but can be accurate enough. Each “well-behaved” 3D surface model constructed by a single component manifold can be parameterised in two dimensions or described utilising a 2D representation.

In the case of model dynamics, numerous approaches have been developed, primarily working directly on the triangular meshes that have already been obtained. Here, we provide a brief overview of the most significant ones.

An extended version, XPBD (Macklin

In 2019, Angles

ARAP minimises shape deformation by discouraging non-rigid transformation via a cost function. Mathematically, this is characterised as the search for a solution to a nonhomogeneous system of linear equations. The matrix corresponds to the discrete Laplace operator of the mesh, while the right-hand side vector encompasses the second differences of each vertex concerning its neighbours.

The main drawback of ARAP and similar approaches lies in recalculating each mesh parameter in every iteration, proving time-consuming for fine triangular meshes. Also, the triangular mesh has the smoothness issue mentioned before. Therefore, a continuous 2D model better approximates smoothness and is better suited for the deformation approaches.

The Non-Uniform Rational Basis Spline (NURBS) is an interpolation and approximation method (Nie

D-NURBS (Terzopoulos and Qin,

Ni

While these approaches offer valuable insights, a common challenge is recalculating every parameter in each iteration, which can be time-consuming. Also, the neighbourhood of vertices needs to be defined for all the techniques mentioned.

Even the previously described NURBS techniques require some relations between the vertices (neighbourhood), but Radial Basis Function (RBF) techniques do not require that knowledge. The task of estimating and fitting scattered data is widespread across various fields of engineering and scientific research. To list only a limited subset of the RBF usage possibilities, Zhang

The RBF method was presented by Hardy (

The radial basis function approach uses “centre points” instead of surface vertices, serving as descriptors for spatial function locations. Unlike NURBS, modifying a single RBF centre location can affect the entire surface because of its “blending” behaviour with other RBFs nearby, propagating recursively throughout the volume. To address the RBF placement possibilities, Majdisova and Skala (

The RBF approach offers advantages over triangular mesh representations of muscle surfaces, including:

There is no need to define the connectivity between control points.

RBFs can generate

RBF may use fewer parameters compared to triangular meshes in general.

Deforming a muscle may involve changing only a portion of the parametric space while maintaining smoothness.

For each application, the correct RBF must be chosen first to reach desirable properties and satisfactory results. Two global RBFs are commonly considered: the Gaussian RBF or thin-plate spline (TPS). The Gaussian RBF is defined as

A single RBF is a mathematical function defined solely by the distance from a designated point, denoted as the centre. Mathematically, it can be expressed as

The expression of RBF approximation, defining the function value at any specific input points, is as follows:

The idea can be extended using polynomial conditions, which can improve accuracy. The extended system of equations, considering these conditions, is expressed as:

Suppose we disregard the polynomial conditions for a while. In that case, RBF approximation involves solving the equation (

The Ordinary Least Squares (OLS) method can choose weights

Due to the often intricate shape of muscles/bones (e.g. with many triangles, quadrilaterals, etc.), approximating the muscle/bone rather than interpolating, is advantageous for subsequent calculations. When approximating, spatial placement of individual RBF is critical. A naive approach involves uniformly sampling the input shape (scalar distance field of the volumetric data or triangular mesh in the case of muscle modelling), but this may not adequately capture the muscle/bone underlying properties.

Even though utilising the polynomial conditions improves the outcomes and enhances the matrix conditionality in some cases, the presumption is that the input data behaves like the selected polynomial. However, this is not the case in muscle modelling, where the data’s overall shape resembles somewhat an ellipsoidal shape rather than a polynomial form, further diminishes the relevance of polynomial constraints, which will not be discussed further.

The number of parameters of the triangular mesh is often large. To address the issue of reducing parameters during deformation, a new mathematical model is proposed, encompassing both static muscle shape and its dynamic configuration during bone movement.

Even though the number of parameters will be reduced, many will still be present. Therefore, the extensive emerging RBF linear equation system would be ill-conditioned. The static model is constructed using a greedy approach, adding one RBF after another to solve the issue. In each iteration, the algorithm identifies the position and shape of a single RBF to minimise the approximation error effectively.

The objective is to minimise the squared error of the approximation to the original triangle mesh across all vertices throughout the bounded space

Here,

Creating the geometrical static model initiates the generation of a scalar distance field (SDF) on the triangular mesh. This field is later sampled using a (quasi-)random sampling method. Throughout the research, various distributions, including uniform, random uniform, Gaussian, and Halton, were tested. The findings indicated that the Halton distribution emerged as the most effective option, which is following other research (see, e.g. Majdisova and Skala,

The geometrical static model depicts a gluteus maximus muscle. The red represents the original triangular mesh, while the blue represents its RBF approximation.

In this study, we followed the placement of each RBF by Carr

Greedy placement is a sub-optimal approach because two independently and subsequently placed centres independently select a shape parameter for each RBF compared to two “cooperating” centres, which are not required to be optimal at the time of their selection. However, this approach avoids the stability issue by using an extensive system of linear equations as a solution.

The optimal position is for each centre selected, and the RBF function is subtracted from the volume bounded by, e.g. the AABB box and the process is iterated until either the allotted number of RBFs is utilised or the predetermined level of accuracy is achieved. This centre point distribution approach was designed to avoid the OLS stability issue altogether.

Choosing the radius of action, also called a shape parameter, is pivotal in achieving precise interpolation or approximation. While one option is to independently select a shape parameter for each RBF, this can lead to many stored parameters. On the other hand, determining an optimal global shape parameter for all RBFs leads to inaccuracies and using it as a compromise remains an open question. Different approaches have been suggested, such as those by Wang

The task at hand for evaluating the placement of individual RBFs is to assess how accurately the placement represents the original shape. One possible metric for this assessment is the Jaccard index; the second is the mean square error.

Strictly speaking, it is a disparity between the area inside and outside both objects and the areas inside one object but outside the other.

between two objects: the original mesh and the surface formed by RBF. It is mathematically defined as follows:In this context,

The computation of the MSE objective function entails integrating the squared L2 norm across the subset Ω of the original space

Here, the function

It’s important to note that these objective functions may encounter challenges when dealing with the translational motion of the shape. This paper assumes that the muscle movement occurs only in part of the muscle while the rest remains static (for instance, the muscle is attached to two bones, but only one moves). This assumption holds as long as the deformation of a muscle is done in the local scope, with one static and one moving bone (meaning both bones cannot translate together during deformation).

Delving further into the challenge of selecting the appropriate objective function, opting for the MSE tends to yield a smoother surface. At the same time, the Jaccard index aligns more closely with the original model at the expense of introducing a less smooth surface. Hence, a combination of both proves to be a favourable approach.

The metric used to compare the original surface with a new one is a weighted sum between the Jaccard index and MSE. This metric is proposed because the smooth surface provided by MSE and the surface location accuracy provided by the Jaccard index are required. The weight

The equation calculates the final metric

Regularisation optimises how the RBFs are placed, increasing the likelihood of placing RBF centres in a predefined position or formation. When forming the RBF surface, it may be advantageous to position RBFs within the muscle volume and less outside, resulting in a more accurate curvature field, including fewer local extrema near the surface. The central concept involves penalising the RBFs placed significantly outside the desired region more heavily than the others.

When utilising the Signed Distance Field (SDF), the sign indicates whether the target vertex is inside or outside. Let’s denote the signed distance as

Additionally, to enable the use of wider RBFs (avoiding the potential replacement by less advantageous numerous RBFs with minimal overall influence), a similar penalisation for the RBF shape parameters can be implemented:

To summarise, the parameters for the following need to be determined. The regularisation factor

In our studied scenario, the bones undergo a rigid transformation in inverse kinematics, typically following the desired real-world movement. Subsequently, in an inverse manner, the shape of all the muscles related to the bone needs to be reconstructed. The portions of the muscle attached to the bones are mandated to move in conjunction with these bones, while the remaining parts of all muscles require reconstruction. In the following text, the muscle is meant to be the whole muscle-tendon unit to ensure that the “muscle” is correctly attached to a set of bones. Also, for the simplicity of the following text, just a singular muscle involved in deformation is considered.

The main idea of the novel dynamic model is to find a new shape according to the curvature of the original one, preserving its initial shape as much as possible. The mathematical description to maintain the initial shape as much as possible can be described in many ways. In the case of muscle modelling, the muscle’s initial curvature throughout the whole volume is a suitable shape descriptor. Mathematically speaking, we define the cost function as the total difference of the original (

Then we use the definition of the 3D RBF approximation (and for the sake of simplicity, use

To calculate the curvature, we need to find the eigenvalues of the Hessian matrix of

For this paper, only global RBFs were considered at first due to their influence over the entire space, allowing the change of the whole model by altering fewer parameters. Moreover, for that purpose, the global Gaussian RBF is chosen for this paper due to its simplicity if differentiation is considered, which is beneficial when creating this dynamical model. If the Gaussian RBF is considered, then the mean of all of the eigenvalues can be calculated as:

Combining (

The calculation, including all the computational steps, can be found in Appendix

The new mathematical model for the RBF shape deformation defines where all centre points should be translated (in the direction specified by its gradient). At this point, we can, according to the broader shape parameters and larger weights, decide which RBFs are more significant and calculate the deformation more often. In contrast, the less critical parameters can be recalculated only a few times, which may reduce the computational complexity significantly.

It’s crucial to re-emphasise the centrality of our mathematical model in this study. Even though this paper does not adequately address any collision detection and response (CD/R) issues, similar techniques to Cani-Gascuel and Desbrun (

The proposed approach has been tested to ensure its usefulness for musculoskeletal modelling. Here, we detail the critical experiments conducted to test the relevance and efficacy of the radial basis function approach in muscle modelling, concentrating on processes and their implications for our theoretical description.

The mathematical model of the curvatures can be evaluated using actual data from the gluteus maximus, gluteus medius, iliacus and adductor brevis muscles. The initial experiment will be conducted without regularisation, solely relying on the greedy placement approach.

Before we dive deeper into the muscle dynamics experiments, the first experiment shows the results of the static model. These experiments prove the correct centre placements.

The original gluteus maximus muscle in red and its RBF counterpart in blue.

Using a greedy approach, we employ a naive RBF placement strategy in the initial experiment. The results of such placement can be seen in Figs.

The original gluteus medius muscle in red and its RBF counterpart in blue.

The original iliacus muscle in red and its RBF counterpart in blue.

The original adductor brevis muscle in red and its RBF counterpart in blue.

By including the regularisation term, higher curvatures will be concentrated at the borders of the original triangular mesh. The regularisation creates more promising conditions for the gradient descent method to achieve a superior approximation of the muscle in motion, mitigating the problem of numerous local minima in the field. Three examples are presented here: one with a regularisation factor of 0.7, another with a factor of 0.3, and the last with a regularisation factor of 0.05. The outcomes can be observed in the Appendix

The weighted metrics (between the MSE and JI) have been tested with two weights (0.15 and 0.85) and on two different sample resolutions (10k and 100k vertices). The evaluation for the first hundred RBF centres is shown in Fig.

The metrics evaluation on gluteus maximus muscle up to 100 RBF centres, with different sampling density and smoothness weight.

The metrics evaluation on gluteus maximus muscle, with different sampling density and smoothness weight.

One may tell from the results that both metrics converge similarly, and in the end, they tend to go to a fixed value. However, if you continue further, it will be evident that this is not the case because the penalisation of the Jaccard index holds near a fixed value. Still, after a while, it goes faster towards zero (on a non-logarithmic scale). The results for more (7000) centres are visualised in Fig.

A well-known problem with the static surface RBF approximation is its ineffectiveness in accurately approximating sharp edges. This drawback is less critical in muscle modelling, as most soft tissues are relatively smooth and do not possess sharp edges, often resembling spherical shapes. However, this issue can be demonstrated using a simple 3D volumetric shape like a tetrahedron. A tetrahedron has six sharp edges (each with a dihedral angle exceeding 70 degrees) and four corners, which present significant challenges for RBF approximation.

The fundamental issue with RBFs is that each RBF inherently forms a spherical isosurface at a given value. To create a sharp corner would theoretically require an infinite number of infinitesimal RBFs. An example of attempting to develop an RBF static surface for a tetrahedron is shown in Fig.

An attempt to create an RBF surface for a tetrahedron demonstrates the unsuitability of RBF for approximating sharp edges. RBF is ineffective for shapes like a tetrahedron due to its inability to model the sharp edges accurately.

Both options involving and neglecting the regularisation factor were tested. For testing, the centre points were shifted randomly, with a magnitude of 5% of the AABB box diagonal, and the expectation is that the centre points return to/near their initial positions. The successful experiment with regularisation is shown in Figs.

The initialisation, where the randomly deformed green muscle should return to its initial shape in blue.

The restoration copies the original blue shape with the green one, despite minor differences on the right side.

This paper shows a novel approach for modelling a surface with the radial basis function approach. The proposed method offers a way to model a static and dynamic muscle surface, altogether avoiding recalculating the geometrical model’s parameters to deform it. Even though there are some examples where the approach is not applicable (e.g. the tetrahedron), the suitability for the muscle reconstruction is apparent. Regularisation helps reduce the number of local extrema. However, it also reduces the resulting precision; luckily, a compromise can be reached using a suitable regularisation term selection.

A deformation methodology respecting muscle properties (muscle volume, shape preservation, and bone avoidance) is the most obvious candidate for future work. A good starting point seems to be the work from Cani-Gascuel and Desbrun (

It is not just the novelty but also the oversight of the approach that is apparent; the metric discussed in the article extends beyond its use in muscle modelling or general modelling applications. It is versatile enough to be applied in any situation that requires the comparison of two surfaces. Similarly, the concept of regularisation, with its mathematical underpinnings suitable for RBF surface modelling, can be adapted for broader applications. Its fundamental principle of constraining centres to remain inside applies to more general contexts, such as clustering methodologies.

The curvature of the RBF surface of the gluteus maximus muscle is depicted in the image. The prominent centres of the RBFs are represented by green spheres in the space. The curvature values are displayed in a logarithmic scale, as the differences in curvature are not linear but rather exponential. This result stems from the same experiment shown in Fig.

On the left, we have the curvature field, while on the right, we see the approximation result with a regularisation factor of 0.7. Although there is a preference for curvature fluctuation within the muscle, there is also substantial fluctuation outside. Decreasing the local minima outside the muscle volume results in a less precise approximation of the original red muscle to the blue one.

On the left, we have the curvature field, while on the right, we see the approximation result with a regularisation factor of 0.3. The fluctuation in the curvature field has been further reduced, but as a consequence, the resulting geometrical surface model has become more rough in texture.

The curvature field (on the left) and approximation result (on the right) with the regularisation factor of 0.05. While the fluctuation in the curvature field has been nearly eliminated, the roughness of the muscle surface has become quite pronounced.

The initialisation of the simulation without the regularisation factor, where the randomly deformed green muscle should return to its initial shape in blue.

The progress of the shape restoration. The green surface seems to diverge from the blue one due to omitting the regularisation factor.

The final part of the shape restoration. The green shape found other local minima than the muscle’s blue (original) shape.

The initialisation of the simulation with the regularisation factor, where the randomly deformed green muscle should return to its initial shape in blue.

The green shape seems to converge to the blue shape. The main difference is the restored top part of the green shape.

The final restoration. Excluding the minor differences in the right part of the muscle, the green shape quite accurately approximates the original blue one.

Let us assume the function

Let’s assume that Gaussian Radial Basis Functions (RBFs) are utilised exclusively. Additionally, we can define

The shape of the geometrical model can be well described using its curvature. To compute the curvature, we need to determine the Hessian matrix, consisting of the second partial derivatives with respect to the function

By employing the chain rule, where the expression

We will now construct the Hessian matrix, which will be used to calculate the curvature. As previously mentioned, its computation relies on knowing the second partial derivatives. In an

To compute the second partial derivatives, we rely on the gradient of a function

Applying the product and chain rules, the second partial derivative with respect to

The calculation of the mixed second partial derivatives is performed similarly to before. Upon examining the outcome in (

The Hessian matrix can be populated with the second partial derivatives. Both sets of second partial derivatives involve multiplication by the distances from the centre point in their respective directions. This process can be encapsulated using an outer product. The value of

For curvature calculation, one option is to employ the mean curvature, which is defined as the average of all the eigenvalues

There are zero eigenvalues with the multiplicity of

The subsequent step involves specifying the cost function, which is essentially the squared L2 norm between the new and original curvatures across the subspace

One may employ the gradient descent method to discover the optimal values of

The resulting gradient takes the following form:

Now, we need to compute the gradient of the cost function

Given the definition of the partial derivatives of

The authors thank their colleagues at the University of West Bohemia and the University of Maribor. Thanks to Vaclav Skala and Ivana Kolingerova for their valuable insight and discussion.