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, 2022). A slight improvement to this model is to replace the straight line with a polyline or a curve that either passes through predetermined points relative to a bone or traces the surface of a parametric body, aiming to minimise the curve length. However, fine-tuning these models to mirror reality is arduous (Hájková and Kohout, 2014).

It’s also imperative to represent muscles not just as mere lines or curves but as substantial higher-dimensional models (Kedadria et al., 2023) to enhance the realism of musculoskeletal models.

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, 2014) and the finite element (FE) method are currently gaining prominence (Delp and Blemker, 2005). The FE method, utilising a template of internal structure projected onto the surface shape, proves superior accuracy. The main issue of those methods is their large number of parameters, which need to be figured out, and their computational complexity. Therefore, this paper concentrates more on the two-dimensional (surface) models instead, which is a reasonable compromise between a simple but inaccurate 1D modelling and 3D accurate but too complex 3D modelling (Macklin et al., 2016).

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.

Position-Based Dynamics.  Position-based dynamics, introduced by Müller et al. (2007), is a rapid technique widely utilised in the animation industry for simulating deformations of elastic objects. Initially created for generic shapes, the method inputs a smooth manifold surface mesh and generates its deformed counterpart.

An extended version, XPBD (Macklin et al., 2016), introduces the concept of elastic potential energy, removing the need to determine time steps and iteration counts. XPBD has been applied to muscle modelling challenges by Romeo et al. (2018), who addressed limitations regarding convergence speed, setup simplicity, intuitive controls, and artistic control. They constructed an internal structure inside the surface mesh, accounting for muscle anisotropy. However, the approach lacked attention to collision detection and resolution, later addressed by Cervenka et al. (2023).

In 2019, Angles et al. (2019) introduced an XPBD-based method for muscle modelling, virtually decomposing the muscle into flexible “rods” approximating muscle fibres. These rods adjust their diameter to maintain volume, enabling real-time simulation. The main issue is the same: the problem of muscle-bone penetration is not adequately addressed or apparent in the paper or its supplementary material.

As-Rigid-as-Possible Surface Modelling.  The as-rigid-as-possible approach (ARAP), developed by Sorkine and Alexa (2007), is a technique for determining minimal non-rigid transformations in a surface mesh. Unlike the method by Kellnhofer and Kohout (2012), which involves volume constraints and an internal skeleton, “it distinguishes itself by not requiring an internal structure”. Fasser et al. (2021) applied ARAP in a medical context to morph the template of a pelvic bone onto subject-specific landmarks, though overlooking certain critical features like volume preservation. Wang et al. (2021) investigated ARAP but chose not to use it due to non-smooth shapes and spikes in results. ARAP has also found application in time-varying meshes by Dvořák et al. (2021).

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 et al., 2022) describing a surface via a set of B-spline curves. NURBS is useful for interactive applications, offering higher smoothness than triangular meshes, allowing intuitive deformation of specific parts of the geometrical model (see, e.g. Clapés et al., 2008 or Sánchez-Reyes and Chacón, 2020). For muscle modelling, the NURBS is not required to form a structure as triangular meshes when creating the data and then approximate those data, meaning that after applying the Marching Cubes or Triangles step, the vertices will not be connected by triangles but rather by individual NURBS patches and the smoothness is accomplished automatically. Also, if a part of the segmentation is missing, the NURBS techniques can restore those parts more precisely than in the triangular mesh case.

D-NURBS (Terzopoulos and Qin, 1994) (Dynamic-NURBS) can be used for dynamics. It extends the NURBS approach by incorporating physical phenomena for more intuitive interactive shape deformation relevant to various purposes such as high-resolution digital terrain (Ye et al., 2020), virtual reality (Lavoie et al., 2006), or CAD applications (Zhang and Qin, 2001).

Ni et al. (2023) used D-NURBS for muscle shape parametrisation and deformation. They achieved satisfactory results by defining muscle shape with a set of neighbouring 2D NURBS plates, reaching a mean square difference of approximately 1.5 mm and a maximum volume error of 0.75%.

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.

In this study, we followed the placement of each RBF by Carr et al. (2001), placing centres at the three level-sets of the isosurface. However, in the case of muscle modelling, only a subset of the proposed set is required. This subset is determined through a straightforward trial-and-error method. It involves evaluating an objective function and dynamically minimising it. Therefore, this placement approach can be described as greedy. The process involves checking the objective function of each RBF for every vertex across all samples within the input space (the bounding box of a mesh), as well as for each predefined shape parameter (refer to Section 4.1.2 for more details).

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 et al. (2021), Afiatdoust and Esmaeilbeigi (2015), and Sarra and Sturgill (2009). Still, they may not always identify the optimal shape parameter for each RBF. During testing, we discovered that the single global shape parameter is not exhaustive. Search from a limited list of possible shape parameters is suitable because the number of parameters is low. The condensed set of parameters reduces the storage capacity and computational complexity while searching for the most suitable one.

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.

Jaccard Index.  The Jaccard index quantifies the disparity in volume1¹

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

In this context, i represents the intersection volume, u symbolises the union of both meshes and d represents the difference between both meshes. Calculating the exact volume of intersections can be complex so that a rough approximation might be a more practical choice. Rather than using volumes, the number of samples (using the already described Halton distribution) in each group can be used. A higher number of samples leads to more significant computational accuracy. It is worth noting that while the Jaccard index helps specify the difference between two shapes, it does not provide an understanding of how much is “inside” or “outside” something, which is a feature to consider.

Mean Square Error.  The Mean Square Error (MSE) provides a more detailed insight into both muscles’ interior and exterior rather than a binary assessment. Calculating MSE requires knowledge of interior information. Interestingly, no additional information is necessary for the RBF representation, as it can be evaluated at any point within the volume. For the original surface mesh, interior information can be acquired using a scalar distance field (SDF) already obtained for its creation, allowing us to determine the distance of each vertex from the surface mesh.

The computation of the MSE objective function entails integrating the squared L2 norm across the subset Ω of the original space R d and can be expressed as: (5) M s e = ∫ Ω ‖ s ( x ) − h ( x i ) ‖ 2 2 d x .

Here, the function s represents the Scalar Distance Field (SDF) constructed over the mesh, while h represents the sum of all currently utilised RBFs. When computed, a more efficient approach is subtracting each RBF from the function s and then calculating the updated s norm.

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.

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 d v i . To obtain relative values rather than absolute distances, we can divide d v by the maximum possible distance d max . This yields r d = d v i d max , which falls within the interval ( − ∞ , 1 ).

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: r s = α i α max . Since ∀ i , α i > 0, the r s value will fall within the interval of ( 0 , 1 ). The wider RBFs are then useful in forcing fewer parameters to move to deform the overall shape.

To summarise, the parameters for the following need to be determined. The regularisation factor γ represents the ratio of how much of the original cost function is employed. (7) C r = γ C m + ( 1 − γ ) d v i d max α i α max .

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. 2, 3, 4 and 5. While this approach is straightforward to implement and relatively fast in execution, as depicted in Fig. 11, it exhibits issues with numerous small (but still different in size) RBFs distributed across the entire space. The approach forms a complex curvature field with numerous potential local extrema. Additionally, the curvature outside the muscle boundaries (though barely visible in the image) is cluttered with unneeded curvature differences, further exacerbating the proliferation of local extrema. The problem is solved with regularisation, and the next section describes the experiment.

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 A in Figs. 12, 13, and 14, respectively.

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. 6, which describes the initial part of the metrics convergence.

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. 7. The other takeaway is that the finer the sampling (green and black dotted curves), the longer it takes to lower the penalisation of the Jaccard index (red and blue). Also, the lower the smoothing factor (red and green), the higher the MSE (dashed lines), confirming that both metrics describe the shape difference significantly differently.

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. 8. Let’s also note that in the case of the tetrahedron and similar shapes, often no data reduction is achieved, but rather the opposite.