doi:10.1155/2010/389356 Research Article Soft Tissue Structure Modelling for Use in Orthopaedic

We present our methodology for the three-dimensional anatomical and geometrical description of soft tissues, relevant for orthopaedic surgical applications and musculoskeletal biomechanics. The technique involves the segmentation and geometrical description of muscles and neurovascular structures from high-resolution computer tomography scanning for the reconstruction of generic anatomical models. These models can be used for quantitative interpretation of anatomical and biomechanical aspects of different soft tissue structures. This approach should allow the use of these data in other application fields, such as musculoskeletal modelling, simulations for radiation therapy, and databases for use in minimally invasive, navigated and robotic surgery.


Introduction
In the last decade, technology revolutionized medical imaging, biomechanical modelling and surgical techniques in the field of orthopaedics. These advances have demonstrated the necessity and feasibility of supportive technologies in clinical practice, including image processing technologies, computer-assisted preoperative planning, image-guided and robotic assisted surgery.
Knowledge of the anatomical-geometrical manipulation of bone, muscles, and neighbouring nervous and vascular structures is essential for safe computed preoperative planning, during navigated and robotic-assisted surgical applications, and for the correct interpretation of postoperative outcomes. All this requires the development of anatomical models that provide digitized data that can be used for geometrical visualization, reconstruction and biomechanical analysis of both preoperative and postoperative anatomy, including the origin, insertion, location and lapses of the muscle fascicules and neurovascular structures surrounding the different joints, as well as detailed three-dimensional data of all osseous anatomical structures involved. The resulting data should obviously not only be anatomically correct but also generalizing, simple in format and easy to communicate, handle and manipulate.
Currently, musculoskeletal imaging techniques such as magnetic resonance imaging (MRI) and conventional computer tomography (CT) can potentially provide a source for such complete anatomical models. Among the various imaging modalities for visualization of bony anatomy, CT is by far the ultimate medium, with the highest resolution. Soft tissue structures such as muscles, vascular tissues and nerves are significantly better visualized on MRI, and image-fusion techniques have been developed to combine the advantages in reconstruction of different imaging modalities into a single specimen [2,3]. However, computational restraints and constraints of clinical reality, for example, time, cost, radiation exposure, and last but not least image distortion when large metallic implants are used, currently preclude the creation of full complex customized models of a region of interest. This gap can be bridged by creating idealized generic musculoskeletal models that can be scaled, morphed and fitted into a patient-specific model, using limited imaging and morphometric data of the patient and a database created from cadaver studies, where bone and functionally 2 EURASIP Journal on Advances in Signal Processing relevant soft tissues such as muscles and neurovascular specifies are collected in detail [4][5][6]. To create generic models that include detailed bony as well as soft tissue anatomy, CT scanning of cadaver models was performed. Soft tissues were contrasted to allow semi-automated data retrieval in a format that would permit data processing for geometrical visualization in different applications and for biomechanical simulations or treatment planning. Nerve and vascular tissue were mathematically described as tubular structures. Muscle tissue was represented as single fibres, surfaces or solid volumes. This allows these data to be used for the construction of subject-specific data on generic anatomically based models. Alternatively, subject-specific anatomically based models can be constructed by customization/transformation of generic finite element geometric models [7,8].

Preparation of the Anatomical Model
Muscle fascicle positions and orientations were predefined on a cadaver model by suturing 0.7-mm flexible copper wires to the muscles, from origin to insertion, following the fibre paths. The brachial plexus, as an example of anatomical complex neural tissue, was carefully dissected and injected with an iodine contrast mixture (Visipaque). A CT multislice scan (Siemens Somatom Volume Zoom) with a slice thickness of 1 mm was performed. The scan started above the occiput and continued down below the hip joint, resulting in about 750 slices. Each slice has 512 × 512 pixels, and each pixel had a grey value in the Hounsfield scale of 4096 grey scale values, meaning that it is represented with a 12-bit value. Voxel size was 0.88 mm 3 . A total data set from a single scan was therefore 512 × 512 × 750 × 12 = 2.36 GBit or about 300 MB. (Figure 1)

Segmentation and 3D Reconstruction of Muscle Fibre Paths and Surface Anatomy
The commercially available Mimics software package (Materialise NV, Heverlee, Belgium) was used for a density based segmentation and reconstruction of the different metallic markers and neurovascular structures. The software package was chosen because it allows for semi-automated segmentation using thresholding, dynamic region growing, multislice editing, Boolean operations and hole filling. Contrary the neurovascular structures, the segmentation could be performed more easily and fully automated due to the uniform and high density of the metallic markers.
Postprocessing of the segmented volumes was then performed for the mathematical definition of the orientation and position of the metallic markers used to define the muscle fibre paths. The algorithm to describe fibre positions for use in biomechanical simulations was based on the original description by Van der Helm et al. [20] A point cloud, generally consisting of more than 2000 points representing each marked muscle fibre, was sorted in x, y or z according to the dominant anatomical direction of the fibre path.  A centreline defining the muscle fibre path was created based on a cluster method algorithm. The clustering algorithm was designed to provide a piecewise segmentation of the structure at interest from beginning to ending, orthogonal to its dominant anatomical direction. Subsequently, a variable t (0 ≤ t ≤ 1) was generated in such way that d(l)/d(t) =l, where l is the total length of the centreline polygon defining the muscle fibre path. This allows to define equidistant parts within the muscle fibre path, a feature that will be used later in the error estimation of the fitting procedure and further surface modelling of the muscle. Furthermore, the insertion (t = 1) and origin (t = 0) are easily deductible. Polynomials were fitted by a least-squares criterion using vectors t x , t y , t z , which are polynomial functions of the variable t (0 ≤ t ≤ 1), to represent all muscle fibre paths: x = a 0 + a 1 t + · · · + a n t n EURASIP Journal on Advances in Signal Processing The mean resultant error in x i , y i and z i was expressed as where N is the number of points defining the centreline, x i,m , y i,m and z i,m are the predicted points on the tpolynomial, and x c,m , y c,m , z c,m the corresponding points on the centreline. Finally to visualize the position and geometry of muscles, a surface patch was generated using the fitted t-polynomials representing the segmented muscle fibres and a cubic spline interpolation of n equidistant points on the t-polynomials. (Figure 2) This method allows to represent each muscle as a limited number of fibres, described by a small n by 3 matrix, from which fibre origin, insertion, 3D positions and muscle surfaces can be easily obtained for use in scalable generic anatomical models. Figure 3 demonstrates the above as the pectoralis major muscle was represented both as a surface patch and a limited number of separated muscle fibres defining the muscles origin and insertional area.

Parametric Reconstruction of Neurovascular Tissues
Numerous vessel extraction techniques and algorithms have been described in the literature. Some of those are applicable to tubular objects that show similar characteristics to vessels, in our case nervous tissue. An interesting overview and classification of most available methods was provided by Kirbas and Quek [9]. The choice of a specific algorithm mainly depends on the intended use of its output parameters.
As we intend to generate orientation and general morphology, for use in generic models, we need an extraction technique that provides a simple description of curvature and vessel diameter. We therefore opted for a generalized cylinder model. The generalized cylinder is a volume created by cross-section swept along a path, the vessel centreline. The centreline is represented by a 3D cubic B-spline, in our case simplified by a polynomial approximation. Similar to the muscle fibre paths, the Mimics software package (Materialise NV, Heverlee, Belgium) was used for initial segmentation of the neurovascular tissue. As the nervous tissue was injected with a contrasting mixture and vessels are air-filled structures in a cadaveric model, semiautomated segmentation of the tissues could be performed using thresholding, dynamic region growing, multislice editing, Boolean operations and hole filling.
Postprocessing of the segmented data of the different nervous and vascular structures was performed to define the orientation, morphology and position of these structures. A centreline of the different nervous and vascular components was generated as previously described for the muscle fibre paths. Nerves and vessels were approximated as tubular structures, of which the local radius was defined using the previously described clustering algorithm. Each cluster represented a piecewise elliptical section of the tube corresponding with the shape of the vessel or nerve. The norm of the projection orthogonal to the centreline of the vector defined by the maximum distance (r max ) between the centreline (v) and the elliptical section, represented the approximated local radius of a cylinder sectioned by that specific cluster: where l is the total length of the centreline polygon defining the anatomical structure. A polynomial function of t (0 ≤ t ≤ 1) was fitted to the centreline and its corresponding radii by a least-squares criterion. The output of the algorithm was a set of directed, 4-dimensional points indicating the (x, y, z) spatial position of each sequential vessel or nerve skeleton point with an associated radius at each point. The total generated equation consisted of x ,y, z representing the centreline, and r representing the radius: x = a 0 + a 1 t + · · · + a n t n The surface of the vessels and nerves were then modelled as curved tubes with variable radii in three dimensional space, based on estimating Frenet-Serret frames along the t-polynomial as originally described by Zerroug and Nevatia [10]. Although easy to implement, the Frenet-Serret formulation model and tube model are known to suffer from serious drawbacks of discontinuities and non-intuitive twisting behaviour at infliction points along the curve [9].  One way around this is to calculate the frame not using the second derivative of the curve (which becomes parallel to the tangent at these infliction points) but using an arbitrary vector that is never parallel to the tangent. (Figure 2) We realize that the above description to obtain the vessel or nerve centreline and radius is approximative and by no means comparable to the advanced algorithms for vessel extraction used in today's automated radiological diagnostic systems. The method obviously is not intended for the diagnosis or description of local morphological abnormalities, but for the generation of generic anatomical models. The algorithm is easily reproducible and computationally fast. Moreover, it delivers anatomy as a generalized model in a format that is easy to communicate and manipulate. The entire subclavian artery for example is represented by a 4 by 4 matrix. Because of these features, our method fulfils the requirements of a generic model aimed at applications in orthopaedics and musculoskeletal biomechanics

Anatomical Model Reconstructions, Challenging Cases
To demonstrate the robustness of our method the following muscles were chosen for their complex anatomical structure: Latissimus dorsi, Trapezius and Deltoid muscle. In general 3rd or 4th-order t-polynomials were sufficient to approximate the different muscle fibre paths. Muscles of less complex morphology will probably be adequately represented by even lower degree functions. The use of t-polynomials EURASIP Journal on Advances in Signal Processing 5 Figure 4: Surface reconstruction of the brachial plexus and intersecting subclavian artery.
facilitates communication on the data and allows for easy deduction of lengths, moment arms,... for use in anatomical reconstructions and biomechanical applications. The errors generated by the fitting algorithm applied to the metallic markers describing the different muscle fibre paths, were on average 1.7 mm (± 0.9 mm) and were mainly due to smoothing of local bends and kinks in the markers from manipulation during dissections and fixation. The use of higher-order t-polynomials fitting did not significantly decrease the error of estimate nor did it affect the overall shape of the reconstructed fibres.
Its complex geometry and close relation to the subclavian artery, make the geometrical reconstruction of the brachial plexus particularly challenging. Because of its tortuous anatomy, higher-degree polynomial fits were necessary to represent the different branches of the brachial plexus ( Figure 3). Instead of these high order polynomial fittings the use of B-splines might be considered to describe such complex neural structures, although these would obscure the anatomical model description which was originally meaned to remain uncomplicated. An anatomical model including the provided examples is shown in Figures 4 and 5.

Transition to Volumetric Muscle Models
Muscles are geometrically positioned between bone and other neighbouring muscles. Once the outer surfaces of all possible delimiting structures are derived, the remaining volume described by each muscle can therefore be predicted on a layer by layer basis. (Figure 6) However, the transition to volumetric data requires though a muscle-specific approach. Muscles close to the bony skeleton such as the rotator cuff muscles, for example, the supraspinatus, infraspinatus and subscapularis muscle, fill discrete cavities onto the scapular body. The muscle's volume is therefore defined by the bony surface delimiting the muscle's inner surface, the scapular body, and the muscle's outer surface as previously obtained. To describe the path of muscle fibres on the inner surface of these muscles,  the problem can be approached as a minimizing geodesic path on the delimiting bony surfaces [1].
Superficial muscle layers on the other hand tend to cover the deeper layers as an elastic membrane. For these, the previously obtained outer surface can be defined as a conformation of this elastic membrane of high potential energy, which is then progressively released to a position of lesser potential energy without penetration of any delimiting structure, for example, bony surfaces or underlying muscles. This can be achieved by defining a finite number of spring elements connecting the surface nodes on the previously obtained muscle's outer surface.
This surface is composed of N-by-N nodes, each describing (N − 1)-by-(N − 1) tetragons. For each tetragon the respective centroid is defined. These are then used to define 6 EURASIP Journal on Advances in Signal Processing  the basic element in the further optimization procedure. The movement of each original node is iteratively calculated by using the locations of its adjacent previously defined centroids only. Each edge connecting the central node with its neighbouring centroids can be seen as a linear spring with an initial length of zero.
Let v i be the vector from the central node to the i the neighbouring centroid: The sum of the spring forces acting on the central node is: where K is the spring constant, and k is the number of neighbouring centroids. Considering that all the springs have initial lengths of zero, we can compute the potential energy of the system as: The cost function to be minimized is the sum of the squared lengths of the edges shared by the same centroid: We can obtain position (x, y, z) that minimizes the cost function by simply finding the geometric centre of the neighbouring centroids: In case the algorithm results in a node displacement that would cause penetration of any local obstacle, the surface point on the obstacle closest to the local minimum is withheld. The process is repeated until a steady state is obtained, defining the inner surface of muscle. Upon closure of both the inner and outer surfaces, the volume describing the muscle at interest remains. (Figure 7)

Application in Orthopaedic Simulations
An implementation of the technique and its applications can be illustrated in reversed shoulder arthroplasty or nonanatomical shoulder replacement. Compared to normal shoulder anatomy, this prosthesis is a reverse ball-andsocket design, which results in important anatomical and biomechanical changes around the shoulder joint. Not only does the design produce significant changes in muscle moment arms, several anatomical structures are translated or stretched [11][12][13].
A cadaver model was prepared as outlined above for the description of normal anatomy. Next, a plastic model of a reversed shoulder prosthesis obtained by rapid prototyping was surgically implanted. The model was then scanned a second time in order to allow model reconstruction following nonanatomical shoulder replacement. (Figure 8) The preoperative and postoperative specimens were studied with use of a helical CT scan (Siemens/volume zoom). Scanning parameters were similar to those described earlier. The shoulder of the specimen was positioned in adduction-internal rotation and the elbow in approximately 90 degrees of flexion. The cervical spine was placed in a neutral position and both wrists were placed and strapped down on the lower abdomen. (Figure 9) The CT images of the first and the second session were uploaded separately into the Mimics software package for further segmentation and 3D reconstruction (Materialise N.V., Heverlee, Belgium). Thresholding for bone, air, the iodine-contrasted tissues and the metallic markers was the first action performed to create a segmentation mask. Then, each region of interest was further selected.
Our method allows for a detailed analysis of changes in geometrical and biomechanical parameters, caused by the surgical procedure. For example measurement of excursion, elongation and displacement of the brachial plexus nerves after reversed prosthesis surgery of the shoulder joint following model preparation as outlined in the present paper has previously been described in detail [12].

Conclusions
Technological assistance in orthopaedic surgery has progressed mainly from the combination of pre-and intraoperative imaging modalities on rigid osseous structures, sometimes combined with the use of tracking systems [14]. Few studies have focused on functional and surgically relevant deformable soft tissues other than skin and fat, such as muscle, vascular and nervous structures, and the postoperative outcome related to their surgical manipulation in orthopaedic applications [12]. Mainly due to the different nature of the treated pathologies, other surgical disciplines have evolved in a completely different, and compared to orthopaedics, mainly soft tissue focused direction, for example, augmented reality, soft tissue navigation and haptic technology applications in abdominal and pelvic surgery and soft tissue models for use in virtual reality training for minimally invasive and laparoscopic surgery [15,16]. Parallel to the technological advances in computer and robotic aided orthopaedic surgery, musculoskeletal modelling has evolved greatly and appears to be close on the verge of integrating biomechanical simulations of soft tissue structures with intraoperative image guidance on bony anatomy [17][18][19]. Such near future improvements require general data that can bridge the gap between the two technological modalities and can be used for the development and validation of the fused endproduct.
Although a variety of anatomical data sets and reports on anthropometric scaling of bone and muscle attachment sites, usually within the context of biomechanical modelling, are available, few studies in the literature have focused on the digitization and parametric description of complete soft tissue anatomical models for direct application in modelling or computer-assisted surgical techniques in orthopaedics [5,6,12,[20][21][22].
Anatomical datasets for biomechanical analysis of an orthopaedic intervention usually describe the anatomy of a normal model, from which surgical manipulation of muscles and the corresponding biomechanical effect can be calculated [20,23]. The classical example and gold standard for comparison remains the work published by Van der Helm et al. [20] who used a 3D-palpator and digitizer for the polynomial description of position and geometry of muscle and muscle attachment sites in the upper and lower limb. A variety of wrapping algorithms are then available for the full description of muscle fibre positions for biomechanical simulation of a specific surgical procedure [21,[24][25][26]. However, currently no data or method exists that can be used for the validation of such mathematically reconstructed anatomy, nor are there data on the choice and position of wrapping objects used to create the resulting simulation environment. Finally, functional and surgically relevant structures such as vessels and nerves have never been included in these simulations, although surgery affects their position and length and can compromise the postoperative outcome [12]. For nervous structures for example, it has been shown that a nerve strain of 5-10% already impairs axonal transport and nerve conduction [27].
The present method allows for such analysis. Following anatomical preparation, estimation of muscle and joint parameters necessary for biomechanical analysis and 3D imaging of the cadaver model, a specific surgical procedure can be performed and the resulting deflection and positions of relevant soft tissues can be visualized and geometrically analyzed. Wrapping objects for use in biomechanical simulations can be estimated and validated from both the reconstructed bony surfaces and the resulting deflections of the reconstructed muscle fascicles and surface models. In the field of navigated and computer-assisted surgery, the geometrical description of soft tissues is particularly useful to create a virtual environment that is reassuring and familiar to the operating surgeon, and to alert him that caution must be exerted during surgery in the proximity of fragile structures such as vessels and nerves.
The technique and its applications obviously have limitations. In the transition from generic to customized models, care should be taken in extrapolating the results; not because the validity of the model is being questioned, but because of the subject specificity, ethnic and racial variations, and the not uncommon occurrence of anatomical variants [6].
Despite this important limitation we believe that, for the current status of technology in imaging, surgical navigation and biomechanical analysis, the described technique offers a number of advantages that might aid in the further development of biomechanical models and computer-assisted surgical applications as well as for the closer integration of both.