At all levels of biological complexity, muscle structure and function are linked. Intermediate-scale structural properties, commonly known as muscle architecture, are defined as the geometric arrangement of a muscle’s fibers with respect to its mechanical line of action.1–4 Seminal work in the early 1990s introduced brightness-mode (B-mode) ultrasound to measure the structural properties of muscle fascicles.5,6 Ultrasound and other noninvasive imaging methods continue to contribute significantly to our understanding of muscle structure–function relationships. Like other noninvasive imaging methods, ultrasound cannot resolve individual fibers; rather, fascicle-scale measurements are taken as proxies for the fiber-scale properties of interest.

Several muscle architectural properties of interest include fascicle length, orientation, and curvature. Fascicle length and orientation characterize the phenomenon of pennation, a force-enhancing strategy in which shortened-length fibers are oriented obliquely to the line of action, placing more sarcomeres in parallel per unit volume. Pennation angle is commonly measured in ultrasound images by manually specifying tangent lines to the aponeurosis and a fascicle at its point of insertion into the aponeurosis and then calculating the angle between the 2 lines. Fascicle length can be measured by manually marking a fascicle of interest, potentially with additional mathematical modeling,7–9 if it is contained wholly within a panoramic imaging field of view or acquired in dual-transducer mode or the length has been estimated by measuring the muscle thickness and pennation angle and applying standard trigonometric relationships under a straight-line approximation.10 Fiber curvature, a requirement for muscle mechanical stability and a proposed determinant of intramuscular fluid pressure gradients during contraction,11–13 was first estimated from ultrasound images by assuming the fascicle to be an arc of uniform radius.14

The predominant method for estimating these properties from ultrasound images has been the manual definition of the structures of interest. However, manual analyses are intrinsically subjective, making them prone to random human error, the limitations of human image perception, and, when image analysts are aware of hypotheses and subject identity, the possibility of experimenter bias. In addition, there is an effective loss of resolution to the measurements of interest, as images obtained with spatial resolutions on the order of micrometers or tens of micrometers per pixel are analyzed with millimeter- or centimeter-scale methods.15

To overcome these limitations, several automated or semiautomated approaches to estimating architectural parameters from B-mode ultrasound images have been described. Published approaches for estimating fascicle orientation include vesselness filtering with orientation modeling,16–20 analysis of spatial gradients in image textural properties,21 feature extraction methods for elongated structures with orientation modeling,22–24 shape modeling based on prior information,25,26 integral transforms such as revoting Hough27 and Radon transforms,27–31 deep32 or machine33 learning approaches, and the implementation and comparison of multiple methods.34 Several of these methods have been released as freely available software.19,20,35 However, we are not aware of any automated or semiautomated method having the ability to estimate fascicle length, pennation angle, and curvature.

The objective of this project was to develop an in vivo ultrasound analysis method for estimating fascicles’ lengths, pennation angles, and curvatures across an entire imaging plane of a skeletal muscle. Our approach was based, in part, on diffusion tensor magnetic resonance imaging (DTMRI) tractography, a structural connectivity mapping method for estimating the trajectories taken by bundles of highly elongated cells between their beginning and ending points. The structural anisotropy of such cells causes the preferential diffusion of water parallel to the long axis of the cell of water diffusion. In DTMRI, a diffusion tensor is used to model this diffusion anisotropy, allowing the local tissue orientation to be represented as a series of unit-length vectors. In tractography, these vectors are integrated from starting points (“seed points”) until a logical stopping point is reached, producing streamlines that represent the underlying tissue structure. While DTMRI tractography has been most commonly applied to white matter tracts of the brain, it has also been applied in skeletal muscle (for a review, see Damon et al36). Here, we adapt a previously described fascicle orientation modeling algorithm16 and use DTMRI tractography methods for skeletal muscle37 to allow B-mode ultrasound-based skeletal muscle tractography. The resulting MATLAB toolbox, the MuscleUS_Toolbox, allows ultrasound-based analyses of fascicle length, pennation angle, and curvature across a region of interest within an ultrasound image. Building on other open science efforts for related methods,19,20,35,37 the package is freely available.

General Description

The MuscleUS_Toolbox was developed in MATLAB R2022a; custom-written functions are housed at GitHub.38 The toolbox includes functions for image preprocessing and fiber tract propagation and quantification, some of which use MATLAB’s Image Processing Toolbox. In addition, there are scripts that illustrate a typical workflow (preprocess_us.m and process_tracts_us.m). The preprocessing and tract propagation/quantification steps are described here and illustrated in Figure 1.

View Full Size Figure 1

—Flow diagram illustrating processing steps. Black-outlined boxes illustrate steps that function automatically after the function is called from the script or command line; a gray-outlined box is used for the step that requires manual intervention. (A) A DICOM-formatted ultrasound image, here showing the tibialis anterior muscle, is opened using the function read_dicom_us( ). (B) By calling the function define_muscleroi_us( ), the user (1) defines the portion of the muscle to be analyzed (dashed light blue lines) and (2) defines the aponeurosis of fascicle insertion (yellow line), which polynomial smoothed and interpolated to define seed points. (C) The grayscale image from (A) and the mask defined in (B) are input into the function bmode2angle_us( ), which models the fascicle orientations by (1) identifying fascicles using a vesselness filter (top panel; upper color bar shows the value of VMax), (2) convolving each pixel of the vesselness image with an anisotropic wavelet (second panel from the top) across a user-defined range of angles according to the angular convention illustrated, (3) taking the angle corresponding to the maximal convolution as the local orientation (third panel; lower color bar shows the angle in degrees), (4) masking out nonfascicular structures using a fascicle mask (fourth panel), producing the pixel-wise fascicle orientation image shown in the fifth panel, and (5) averaging over user-defined grid sections (bottom panel), here shown at 45 × 45-pixel spacing. (D) The gridded fascicle orientation image from (C) and the mask and aponeurosis definition from (B) are input into the function fiber_track_us( ), which integrates the local orientation vectors to form a set of unsmoothed fiber tracts. (E) The function fiber_smoother_us( ) is used to smooth the tracts according to a polynomial function of user-defined order. (F) The function fiber_quantifier_us( ) is used to estimate the fascicles’ length, pennation angle, and curvature. B-Mode images, masks, aponeurosis definitions, and tracts are displayed using the fiber_visualizer_us( ) function. Images have been cropped and zoomed in to highlight details. B-mode indicates brightness-mode; VMax, maximum value across all convolution steps.

Citation: Journal of Applied Biomechanics 39, 6; 10.1123/jab.2022-0270

Figure 1

—Flow diagram illustrating processing steps. Black-outlined boxes illustrate steps that function automatically after the function is called from the script or command line; a gray-outlined box is used for the step that requires manual intervention. (A) A DICOM-formatted ultrasound image, here showing the tibialis anterior muscle, is opened using the function read_dicom_us( ). (B) By calling the function define_muscleroi_us( ), the user (1) defines the portion of the muscle to be analyzed (dashed light blue lines) and (2) defines the aponeurosis of fascicle insertion (yellow line), which polynomial smoothed and interpolated to define seed points. (C) The grayscale image from (A) and the mask defined in (B) are input into the function bmode2angle_us( ), which models the fascicle orientations by (1) identifying fascicles using a vesselness filter (top panel; upper color bar shows the value of VMax), (2) convolving each pixel of the vesselness image with an anisotropic wavelet (second panel from the top) across a user-defined range of angles according to the angular convention illustrated, (3) taking the angle corresponding to the maximal convolution as the local orientation (third panel; lower color bar shows the angle in degrees), (4) masking out nonfascicular structures using a fascicle mask (fourth panel), producing the pixel-wise fascicle orientation image shown in the fifth panel, and (5) averaging over user-defined grid sections (bottom panel), here shown at 45 × 45-pixel spacing. (D) The gridded fascicle orientation image from (C) and the mask and aponeurosis definition from (B) are input into the function fiber_track_us( ), which integrates the local orientation vectors to form a set of unsmoothed fiber tracts. (E) The function fiber_smoother_us( ) is used to smooth the tracts according to a polynomial function of user-defined order. (F) The function fiber_quantifier_us( ) is used to estimate the fascicles’ length, pennation angle, and curvature. B-Mode images, masks, aponeurosis definitions, and tracts are displayed using the fiber_visualizer_us( ) function. Images have been cropped and zoomed in to highlight details. B-mode indicates brightness-mode; VMax, maximum value across all convolution steps.

Citation: Journal of Applied Biomechanics 39, 6; 10.1123/jab.2022-0270

The function read_dicom_us() opens a Digital Imaging and Communications in Medicine (DICOM)-formatted image, reads its metadata, and then prepares output structures suitable for use in subsequent functions. An example of a B-mode image is shown in Figure 1A. The second preprocessing step is muscle boundary and aponeurosis definition, performed using the function define_muscleroi_us() and illustrated in Figure 1B. Muscle boundary and aponeurosis definition are manual steps that are used to limit analysis to a region of interest and define the seed points for tractography. The X and Y coordinates of the user-selected points on the aponeurosis are fitted to third-order polynomials as functions of distance and interpolated at a user-defined spacing. Then, fascicle orientations are estimated without additional manual intervention, using procedures modified from Rana et al16 as implemented in the function bmode2angle_us().

The first step involves the use of a vesselness response (V[s]) function to detect elongated structures (ie, fascicles). Briefly, this involves convolving the B-mode image with Gaussian functions of increasing SD from an initial value, SD1. The Hessian matrix and V(s) are found at each convolution step: