Simple Muscle Architecture Analysis (SMA): An ImageJ macro tool to automate measurements in B-mode ultrasound scans

In vivo measurements of muscle architecture (i.e. the spatial arrangement of muscle fascicles) are routinely included in research and clinical settings to monitor muscle structure, function and plasticity. However, in most cases such measurements are performed manually, and more reliable and time-efficient automated methods are either lacking completely, or are inaccessible to those without expertise in image analysis. In this work, we propose an ImageJ script to automate the entire analysis process of muscle architecture in ultrasound images: Simple Muscle Architecture Analysis (SMA). Images are filtered in the spatial and frequency domains with built-in commands and external plugins to highlight aponeuroses and fascicles. Fascicle dominant orientation is then computed in regions of interest using the OrientationJ plugin. Bland-Altman plots of analyses performed manually or with SMA indicate that the automated analysis does not induce any systematic bias and that both methods agree equally through the range of measurements. Our test results illustrate the suitability of SMA to analyse images from superficial muscles acquired with a broad range of ultrasound settings.


Introduction
Musculoskeletal ultrasound imaging is used in a wide range of fields, including the study of muscle and tendon function [1,2], the effects of training on muscle architecture [3], and in the study of architectural parameters in different clinical populations [4]. The images generated by this method are complex and require a great deal of time and effort from practitioners to interpret and extract meaning from them. In recent years, efforts have been made to automate parts of this process [5][6][7][8][9][10]. However, these efforts have been fragmented, and suffer from a number of limitations: they often focus on analysing a single parameter of interest; most publications do not reveal specific details of how to implement the method; some methods rely on software that require expensive licence fees; the majority of methods are only semi-automated, requiring manual, subjective interpretation of at least some images; and often tracking methods involve complex mathematics and require computer programming experience. Collectively, these factors limit the widespread use of existing methods. In this study we present Simple Muscle Architecture Analysis (SMA), a fully automated method of analysing muscle architectural parameters from individual images or collections of images. The method consists of a single macro in Fiji [11], which is a distribution of ImageJ [12,13], and is open source software that is commonly used to process ultrasound images. We provide full instructions for using the method, and no previous programming experience is required. In the following sections we present an overview of the method, as well as examples of analysis output and a range of test metrics.
A preprint version of this paper was first published at the following address: https://arxiv. org/abs/1905.09490v2.

Overview
In addition to the connective tissue surrounding bundles of muscle fibres (perimysium) and whole muscles (aponeuroses), ultrasound scans of muscle architecture acquired in brightness (B) mode show several types of echogenic tissues (e.g. fat, other types of connective tissue, blood vessels). A challenge associated with automatic segmentation of fascicles and aponeuroses is therefore to discriminate between signals from all echogenic regions (and artefactual echoes) and those from objects of interest. One of the originalities of the present approach is the implementation of filters in the frequency domain to isolate objects of interest (i.e. aponeuroses or fascicles) and improve their detection and registration. The sample images used to test the present method were collected in previous projects, for which the subjects gave permission to use images by signing an informed consent and ethical approval was granted by the ethical committee of the Norwegian School of Sport Sciences. The analysis method is summarised in Fig 1. The analysis script is written in ImageJ 1.x macro language and runs in Fiji [11]. Although this scripting language is less powerful than others (e.g. Python), it is easy to learn and use, and was chosen with the goal of developing a tool customisable by most users. The script includes a header, a macro and auxiliary functions called within the macro. The header contains information about necessary plugins, the GNU General Public License, and the script parameters selectable by the user. The macro contains a main function SingleImageAnalysis and sub-scripts to display the results and run the macro on a single image or in batch mode (i.e. on series of images). The SingleImageAnalysis function is divided into four main sections, which perform the following operations: i) detect the ultrasound field of view (FoV), ii) segment and register superficial and deep aponeuroses, iii) measure dominant orientation of fascicles, and iv) calculate and compile variables of interest. The first three sections include a pre-processing step to isolate objects of interest and a function or a plugin (for fascicle orientation) to segment these objects. Several spatial filters are used in pre-processing. Although the purpose of their implementation is specified, their description would extend beyond the scope of this article. The reader is therefore referred to the relevant references for additional information.

Detection of the ultrasound field of view
Ultrasound images nearly always come with a frame containing various information about institution/patient name, scanning parameters and scaling. To proceed with the analysis of the ultrasound FoV, this section optionally (see the "Program description" section) identifies and crops the frame out. This step is achieved by filtering out most of the text from the frame by applying a convolution and a median filter, and binarizing the whole image with a median local threshold. The resulting segmentation enables the selection and deletion of the frame from a duplicate scan but does not alter the original image. The duplicate image of the FoV obtained at the end of this step is used in the three following sections.

Detection of the aponeuroses
The pre-processing of this section begins with a series of spatial filters using built-in commands or plugins (Enhance Local Contrast, CLAHE [14]; Tubeness [15]) in ImageJ/Fiji and a separately installed plugin (Non-Local Means Denoising [16,17]), to remove ultrasound noise and enhance muscle aponeuroses and fascicles. The Tubeness plugin is an implementation of multiscale vessel enhancement filtering first proposed by Frangi et al. [18]. The plugin scores each point on the basis of eigenvalues of the Hessian matrix, to determine an index of "tubeness". Although this approach was first developed to detect blood vessels, it can also effectively enhance muscle fascicles in ultrasound images [5]. A Fourier transform (FFT) is then applied to the filtered image. The resulting power spectrum is thresholded to retain the dominant signal from muscular aponeuroses and fascicles, and a mask is applied to the thresholded image to suppress features orientated in the expected direction of fascicles (see the "Program description" section). The inverse FFT of the filtered power spectrum is subsequently computed, and the edges of the aponeuroses are segmented (Canny Edge Detector plugin [19]) and registered using a customwritten function. This function uses a simple contrast threshold to search for parallel 'lines' in the horizontal direction, which typify the appearance of aponeuroses (see Fig 2B).

Measurement of the dominant fascicle orientation
At least two overlapping regions of interest (ROI) are defined along the width of non-filtered duplicates of the FoV, excluding aponeuroses (see the "Program description" section). Multiple overlapping ROIs were shown in pilot experiments to yield the dominant fascicle orientation more consistently than single or non-overlapping ROIs. Similarly to the aponeurosis detection, ROIs are first filtered spatially and in the frequency domain. However, to better preserve fascicular features, noise filtering is "lighter" than for aponeuroses (median filter and Non-Local Means Denoising plugin). Unless set manually (see the "Program description" section), a thresholding proportional to the brightness of the filtered image is applied to the power spectrum to mostly retain the frequencies corresponding to fascicles in the ROI. Following an inverse FFT, fascicles are enhanced with the Tubeness plugin. Their dominant orientation is measured within each ROI with the OrientationJ plugin [20,21]. OrientationJ is a plugin developed by Daniel Sage [20] to compute orientation and isotropy properties within a ROI, based on the evaluation of the gradient structure tensor in a local neighbourhood. The local

Fig 2. Examples of raw images of the ultrasound field of view, aponeuroses edge detection and overlaid ROIs (yellow) on a suitable scan (A-C), and on a scan where the lower aponeurosis is imaged with insufficient contrast (D-G).
In C and G, the paths of the superficial and deep aponeuroses are also indicated as green lines. When the aponeurosis is imaged with insufficient contrast it is not detected (E). In this case using a higher Tubeness sigma may help detection (F-G), although too high values may cause a slight distortion of the aponeurosis edge. window is characterized by a 2D Gaussian function of standard deviation σ, and here the gradient is computed with a (quasi-isotropic) cubic spline filter. The standard deviation of the Gaussian filter σ is defined by the user (see the "Program description" section), and should be proportional to the thickness of the fascicles (i.e. thicker fascicles are associated with a larger σ). The dominant angles from all ROIs are collected and either the greatest, the mean or the median angle-as selected by the user-is retained for each analysed image.

Calculation of architecture parameters and display
The angle of fascicle pennation is computed as the sum of the angles characterising the orientations of the deeper aponeurosis and fascicles, measured in previous steps. Fascicle length is computed as the length of a straight line running between aponeuroses, at an angle corresponding to the dominant fascicle orientation (either greatest, mean or median angle). Thus, the current version of the script estimates linear fascicle paths. A previous study comparing the length of gastrocnemius medialis fascicles when taking their curvature into account versus when assuming a straight path and parallel aponeuroses found a 6% difference during maximal contraction [22]. The difference seen in resting muscle was negligible. Because the current approach takes the orientation of aponeuroses into account, we expect an error negligible at rest and smaller than that reported by Muramatsu and colleagues during contraction. Future implementations of the script may estimate curved paths. Although the analysis was developed for scans displaying entire fascicles, the script enables the extrapolation of aponeuroses if the composite fascicle runs outside the FoV. Muscle thickness is calculated as the mean distance between aponeuroses. The SingleImageAnalysis function ends after this step; if batch mode was selected the macro loops back to the beginning of the function and repeats the analysis for all other images. Numerical results are displayed in a table and a composite fascicular line is overlaid on corresponding original images, arranged in a stack.

Program description
The script can be opened and launched in Fiji without installation, but end users are advised to install it in the Fiji.app/plugins folder. It can then be found in the Plugins menu of the Fiji software. When launched, the script opens a user interface divided into 5 sections. A brief description is displayed when the mouse pointer hovers over a parameter. A message at the top reminds the user that the analysis script was designed for ultrasound scans displaying the proximal side to the left of the image, where fascicles would typically be angled upward to the left and the lower aponeurosis would be orientated horizontally or downwards to the left. The check box below this message allows the user to flip analysed images horizontally if required. In addition to the orientation of the image, we recommend following some basic scanning guidelines. The transducer should be used in the highest frequency range and the focal zone set to the deeper aponeurosis. When possible, the ultrasound beam can be steered towards the perpendicular to the direction of the fascicles, as long as aponeuroses still appear clearly. Using compound imaging is a preferred alternative to changing the beam orientation, when this feature is available. Any other feature (e.g. tissue harmonic imaging, speckle reduction) capable of reducing noise is also recommended. The dynamic range (or compression) should be moderate or low. The following sections detail the choices offered in the user interface.
Type of analysis. The user chooses here to analyse single images, or, if multiple images are being analysed, selects target input and output folders. In the second case the extension corresponding to the file format must be specified. A check box enables the user to print all analysis parameters in the results table.
Image cropping. Images can be cropped manually. In this case the user is prompted to draw a selection around the field of view. In the case of folder analysis, the same dimensions of the cropping selection are used for all files, so it is advisable to only include images taken at the same scanning depth to avoid accidentally cropping the FoV. Manual cropping is recommended when aponeurosis detection fails, in particular when bright structures (e.g. other aponeuroses, bone edges) are visible under the muscle of interest.
Aponeurosis detection. The only parameter that can be adjusted at this step is the standard deviation-sigma-of the Gaussian used in the Tubeness plugin to convolve the image (see section "Detection of the aponeuroses"). The default value is set to 10 on the basis of our sample images, but different spatial resolutions may require another value. A lower sigma (e.g. 8) is recommended whenever possible, as larger values may result in inaccurate detection of aponeuroses edges. Aponeuroses imaged with an insufficient contrast may not be detected. In such cases, an error message is displayed and the user is invited to crop the image manually or to use a higher value of Tubeness sigma (Fig 2D-2G). Note that choosing too high values may cause a slight distortion of the aponeurosis edge. The alternative suggestion of cropping the image manually may help to exclude areas of the image (i.e. subcutaneous fat, blood vessels or connective tissue) that are too close to the aponeuroses and prevent their accurate detection.
Fascicle orientation. In this section the user chooses the number and the size of the ROIs used to measure dominant fascicle orientation. These parameters define the degree of overlap and the depth of the ROIs. The choice of ROI height enables the computation of fascicle orientation within the whole image or in regions closer to the deeper aponeurosis. The latter is often favoured by researchers because of the better alignment between the deeper aponeurosis and the direction of muscle force, but this theoretical argument should be weighed against the smaller number and consistency of visible fascicles within smaller ROIs.
As explained in section "Measurement of the dominant fascicle orientation", structures aligned in a different direction than the fascicles are filtered out by thresholding the power spectrum of ROIs (Fig 3). The threshold value is set automatically but the user can uncheck this option and set the threshold manually. In the next step, the Laplacian of Gaussian filter value (referred to as σ in "Measurement of the dominant fascicle orientation") can be set according to the spatial resolution of the image and expected fascicle thickness. The default value (4) was found to be suitable in most of our tests. However, the value of this parameter will influence the measurements and using a σ consistently is necessary to obtain reliable results. A "Test" function is available here, to display the output measurements with σ values computed between 0 and 7. The displayed table indicates notably the dominant orientation and coherency, the latter being an index (0-to-1) representing the anisotropy level. In addition to informing the decision of σ value, this test can be useful to identify scans of insufficient quality for this analysis: in case of an insufficient fascicle resolution the orientation will increase continuously with σ. For such scans, the dominant orientation will appear lower than that of the fascicles upon visual inspection of the final overlaid image. Finally, the user can select the method used to obtain the dominant fascicle orientation out of the values obtained in each ROI.
Pixel scaling. For single images or series of images acquired at the same depth, measurements can be scaled if this option is checked. The user should then select the scanning depth, and at the onset of the analysis when prompted, draw a straight line over the scale bar usually included in the frame of ultrasound scans (as is normally done when scaling images in ImageJ/Fiji).
System requirements. Hardware/software specifications are based on Fiji/ImageJ requirements. Currently, the following systems are supported: • Windows XP, Vista, 7, 8 and 10 • Mac OS X 10.8 "Mountain Lion" or later • Linux on amd64 and x86 architectures However, the capability to use Java 8 runtime is the only real requirement. Mode of availability. The SMA script and the necessary (Non Local Means Denoise, Canny Edge Detector and OrientationJ) or optional (FFMPEG) plugins are installed via the ImageJ Updater interface, by adding relevant update sites. This relatively simple procedure has the advantage of automatically updating all of the required components and libraries. Generic steps for adding an update site are explained on the ImageJ website (https://imagej. net/Following_an_update_site), but detailed instructions for installing SMA are provided in the S1 Supporting Information. The script can also be accessed directly: https://sites.imagej. net/SMA/plugins/.
Performance measure. To assess the performance of the automated analysis, scans of muscle architecture were analysed with SMA and manually, by the two investigators. For all variables, measurements obtained from SMA and manual analyses had equal variance (Bartlett criterion) and were normally distributed (Shapiro-Wilk test). The agreement between manual and automated analysis was assessed with Bland-Altman plots. Measures of reliability were not included since repeated SMA analyses of the same scans yield identical results. Likewise, we chose here not to include any comparison between SMA analyses performed on repeated scans. The repeatability of probe placement and generally, the skills of the operator of the ultrasound are critical factors affecting repeated measures of muscle architecture but this extends beyond the scope of this article as it does not relate to the performance of the algorithm.

Results
The analysis is currently optimised for superficial muscle architecture images taken at rest, with entire fascicles visible within the FoV. However, aponeuroses can be extrapolated linearly and the length of elongated fascicles (e.g. when the muscle is being stretched) can be estimated with satisfactory accuracy. The primary purpose of the macro is therefore to obtain reliable estimates of muscle architecture in repeated measurements from the same individuals. In the first example below, we compared the reliability and validity of this type of analysis against manual analysis (Figs 4 & 5). To illustrate the effect of variations in image parameters due to the use of different equipment and ultrasound settings, sample scans of gastrocnemius medialis muscle were acquired from 15 individuals (Sample A), at 9 MHz, with a 96-element transducer (60mm, LV7.5/60/96Z, LogicScan 128 EXT-1Z, Telemed, Lithuania), and 15 different sample scans (Sample B) were acquired at 12 MHz, with a 128-element transducer (50 mm, 5-12 MHz HD11XE, Phillips, Bothell, Washington, USA). Ultrasound parameters were adjusted to visualise fascicles and aponeuroses appropriately but differently for each sample (e.g. contrast, smoothing), to reflect a broader range of image quality than obtained with hardware alone. Therefore, the results of the analysis of these samples do not reflect brand-specific capabilities of the ultrasound systems used here. Both samples were analysed with the SMA script, and manually (and independently) by the two authors. Manual analysis consisted of digitising 3 fascicles per image and the distance between superficial and deep aponeuroses at 3 different sites (approximately 25, 50 and 75% of image width). An average of these measurements was used to represent pennation angle, fascicle length and muscle thickness.
Tracking of architecture parameters in image series is also possible with the SMA script, as long as fascicles are imaged consistently, as in during slow contractions or passive stretching. If scans were acquired as movies, the user should first install FFMPEG plugins [23] to be able to import them into Fiji, before saving them as image sequences. An example of analysis of images (Sample C) acquired at 12 MHz, with a 128-element transducer (60mm, HL9.0/ 60/128Z-2, LogicScan 128 EXT-1Z, Telemed, Lithuania) on the tibialis anterior muscle of an individual is illustrated in Fig 6. Scans were acquired at 15 frames per second, during an unconstrained dorsiflexion-plantarflexion movement. The analysis was performed with the following parameters: automatic cropping, Tubeness sigma (aponeurosis detection) = 7, 3 ROIs, ROI width = 60, ROI height = 90, automatic thresholding of the ROIs power spectrum, OrientationJ σ = 4 and maximal value of all dominant orientations detected in ROIs. The mean analysis time per frame was 8.4s on the computer used for this example.
The three datasets analysed during the current study are available as separate compressed files at the following address: http://dx.doi.org/10.17632/dpmf9bz8pt.2.

Discussion
The measurements of muscle architecture parameters obtained with SMA are within the range of expected values obtained manually. Bland-Altman plots show that the automated analysis does not induce any substantial bias, in particular when taking the inter-rater comparison into account. The bias in pennation angle is in all cases less than or equal to 1˚. Image quality (i.e. spatial resolution and contrast) probably influenced the estimation of fascicle length, with a bias of up to 5.6 mm in Sample A and less than 1 mm in Sample B. Importantly, the bias magnitude was not specific to manual vs. automated analysis. The bias in muscle thickness was negligible in all comparisons (<1 mm). None of the measured biases were found to be proportional to averaged values, indicating that manual and automatic analyses agree equally through the range of measurements. The validity of muscle architecture measurements performed manually has previously been shown [24]. The present comparisons demonstrate that the results obtained with SMA and manual measurements are comparable and, indirectly, equally valid. In addition, automated measurements improve the reliability of muscle architecture measurements by removing the variability induced by manual segmentation, particularly when multiple raters are involved [25]. As shown in examples of scans processed with the SMA script (Fig 7), the implemented filtering and user adjustments allow the processing of scans acquired with a broad range of ultrasound settings.
However, the analysis relies heavily on sufficient contrast and homogeneity of grey values when imaging fascicles and aponeuroses. For instance, large differences between grey values of the lower and upper aponeuroses may cause the detection to fail. Likewise, inconsistent grey values along the length of an aponeurosis may affect its segmentation. Users are advised Automated analysis of muscle architecture to test the suitability of new image settings prior to data collection. This step is important since different settings will influence the results of any analysis. For instance, different sigma values chosen to detect the aponeurosis (see sections "Detection of the aponeuroses" and "Aponeurosis detection") will yield slight differences in the computed orientation of these structures. Similarly, the parameters used to compute the dominant fascicle orientation (e.g. ROI size and number, σ, see section "Fascicle orientation") will influence the results. In addition to testing these parameters on sample data, the user is advised to keep the same settings when analysing subsequently collected data. Except in cases where parts of the image need to be excluded (e.g. due to excessive connective tissue) and cropping must be done manually, the above precautions will generally enable the analysis process to be fully automated if desired. Finally, the script was designed for superficial muscles. Deeper muscles may be targeted by manually cropping the FoV but the above requirements may limit their analysis. The automated analysis proposed here yields results comparable to those obtained manually. The SMA script, the first implementation of such an analysis to be proposed as a free and open source software, is therefore a suitable tool to reduce some of the variability introduced by manual analyses, and to improve the overall quality of muscle architecture measurements. Since it relies on sufficient contrast and homogeneity of grey values, it is mainly intended for use with images from resting muscles, but our test results also show that scans acquired during motion can be analysed, as long as qualitative image criteria are met.