A review of biomechanically informed breast image registration

John H Hipwell; Vasileios Vavourakis; Lianghao Han; Thomy Mertzanidou; Björn Eiben; David J Hawkes

doi:10.1088/0031-9155/61/2/R1

1. Introduction

Breast cancer is the most commonly diagnosed cancer and the leading cause of cancer death in women worldwide. With more than 1.3 million women diagnosed each year, it is the second most common cancer, after lung cancer, in men and women combined (Jemal et al 2010). Incidence and mortality rates vary internationally however (Ferlay et al 2013), with incidence rising whilst mortality exhibits a stable or negative trend in western countries, due to widespread screening and improved treatment respectively, but mortality displaying an upward trajectory in Eastern Europe, Asia, Latin America and Africa (Jemal et al 2010). Whilst early detection remains a key factor in maximising probability of survival, the risk of overdiagnosis is also a concern (Marmot et al 2012). Prognostic grading of breast cancer, the importance of which has been known for more than 50 years (Bloom and Richardson 1957), has been thrown into sharp relief in recent years with the discovery of specific molecular sub-groups of breast cancer (Curtis et al 2012) and associated stroma (Farmer et al 2009) which exhibit a range of prognostic characteristics. These developments, together with advances in detection, diagnosis, staging, and treatment monitoring via radiological imaging, suggest that personalised stratification of breast cancer, with respect to treatment response and overall prognosis, may be possible if all the available imaging and histopathological information can be combined (de Abreu et al 2013).

Relating the appearance of the breast across the range of imaging modalities used in breast radiology; from routine imaging via x-ray mammography, MRI (with and without contrast) and ultrasound (both 2D and 3D); to those becoming available more recently such as digital breast tomosynthesis (DBT), and dedicated breast imaging systems for PET and tomographic ultrasound; is a challenging task due to the highly deformable nature of the breast. In each of these modalities the subject's pose typically changes between one of either upright, prone or supine positions. In screening or diagnostic breast imaging, such as x-ray mammography and DBT, or invasive interventions, such as biopsy, the breast is immobilised between compression plates (Highnam and Brady 1999) and may be subject to interventional forces. During breast conserving surgery the patient is usually lying supine, possibility with the bed at an angle, with her arm extended. In an MRI acquisition the patient is usually lying in the prone position with the breast pendulous, however during an ultrasound acquisition or surgery, the patient lies in the supine position with the breast compressed due to gravity against the chest wall. During ultrasound the breast may be further subjected to compression forces due to pressure of the ultrasound probe. Tissue sampling by biopsy or fine needle aspiration would benefit from image guidance yet, for images acquired before the procedure, the breast is invariably in a very different position to that during the intervention. The problem of breast deformation is even more acute when pre-operative imaging is used to plan and guide lumpectomy in breast conserving surgery, where it is critical that clear surgical margins are achieved, free from tumour spread.

Due to these large scale deformations there has been considerable interest in deformation modelling to capture the transformation of the breast between different imaging sessions and between imaging and biopsy or surgical resection. Many of these methods have focused on biomechanical modelling and how this computational tool can be utilised to create non-rigid image registration methods to solve the inter- and intra-modality correspondence problem as well as correspondence between images and the physical space of an intervention.

In this review we investigate current research in the field of breast biomechanical modelling and identify publications where the resulting models have been incorporated into breast image registration algorithms. We begin with an introduction to relevant concepts in the fields of non-rigid image registration and deformation modelling, placing particular emphasis on biomechanical modelling and finite element methods which have been widely applied in this field. Following this we review the state of the art in breast biomechanical modelling and soft tissue material characterisation, and perform a comprehensive survey of breast image registration methods, broken down by the specific breast deformation combinations addressed.

2. An overview of non-rigid breast image registration

The goal of image registration is to calculate the correspondence between points in an image, I, and those in a fixed reference image or the physical space of an intervention, R. Registration algorithms typically consist of three components: (i) a transformation model, ${{T}_{\theta}}(x)$ , to constrain the deformation between the images, or image and physical space, over the spatial domain, x, specified via a set of parameters, θ; (ii) a cost function, $\Psi$ , to capture the similarity (or distance), D, between the images, and optionally regularise the result (λ);

$\begin{eqnarray}&&\Psi=D\left[I\left({{T}_{\theta}}(x)\right),R(x)\right]+\lambda \left({{T}_{\theta}}(x)\right)\,,\end{eqnarray}$

and (iii) an optimisation strategy to minimise the cost function: $\arg {{\min}_{\theta}}\,\left(\Psi\right)$ (Hill and Batchelor 2001, Modersitzki 2004). Whilst this classification is a gross simplification of the range of algorithms applied to the problem of medical image registration, it provides a convenient context for the current discussion of breast image registration.

2.1. Breast image registration topics

Breast image registration has been utilised for mammographic breast density measurement, e.g. Pereira et al (2010), and as an aid to computer aided detection of mammographic lesions, e.g. van Engeland et al (2003). It has been applied to MRI when estimating dynamic contrast enhancement (DCE), e.g. Rueckert et al (1999), and diffusion weighted imaging, e.g. Arlinghaus et al (2011), to eliminate the effects of patient motion of the breast. A confounding factor in DCE-MRI, from the point of view of registration, is the change in intensity over time generated by the contrast enhancement. This issue has been tackled by a number of authors who propose a variety of methods to separate this enhancement from the effects of patient motion (Ebrahimi and Martel 2009, Melbourne et al 2011a). The application to radiotherapy planning was proposed by Alderliesten et al (2013), by computing the large deformation between prone pre-operative and supine planning scans.

Ultrasound imaging is a particularly challenging modality to align with itself or other modalities due to the directional nature of the received signal. For this reason a number of groups have investigated acquiring co-registered images by building dedicated multi-modality scanners combining, for instance, mammography and ultrasound for the purposes of biopsy guidance e.g. Irwin et al (2008), Kapur et al (2004) and Kapur et al (2002). Despite the technical challenges, attempts have been made to compensate for patient motion during and between ultrasound scans, e.g. Boehler and Peitgen (2008), and compound image volumes using elastic registration, e.g. Krucker et al (2002).

Multi-modality breast imaging has been achieved by combining medical imaging techniques such as positron emission mammography and x-ray mammography, e.g. Bergman et al (1998). When these modalites are not co-registered at the point of acquisition, image registration has been used to fuse them. Examples include fusion of the functional information present in MRI with the x-ray attenuation representation of mammography, Behrenbruch et al (2003) and Krueger et al (2013), registration of MRI and ultrasound to guide biopsy, Causer et al (2008), and registration of MRI and PET-CT to aid breast tumour characterisation, Dmitriev et al (2013).

Image guided biopsy or fine needle aspiration procedures currently in routine clinical practice either rely on real-time imaging, such as ultrasound, in which the lesion is located directly in the frame of reference of the tissue sampling needle, or the breast is immobilised between imaging and tissue sampling, in for example MR guided biopsy, Orel et al (1994), or stereotactic x-ray biopsy, Parker et al (1991). Accurate co-localisation of other imaging data enables much more effective use of the information they contain.

Image guided breast surgery is a particularly challenging task that has received limited attention, however preliminary research suggests that a biomechanical approach could be of benefit (Carter et al 2005, 2008).

2.2. Transformation models

Transformation models vary greatly in complexity and the range of deformations they attempt to capture. At the simplest level are rigid or affine transformation matrices, composed of combinations of translation, rotation, scaling and shear. With a maximum of twelve parameters (for a 3D affine transformation), they are commonly used to initialise a more complex non-rigid registration. However, despite their simplicity, affine transformations have found application in breast image registration, with reported accuracies for mammography registration of around 6 mm (Pereira et al 2010) or 8 mm (van Engeland et al 2003) and for MRI to x-ray registration of 13 mm (Mertzanidou et al 2012b). Arlinghaus et al (2011) applied affine registration to the alignment of DWI-MRI and reported a reduction in the variance of ADC measurements as a consequence. Nonetheless, feature based approaches have broad applicability to mammogram registration via breast specific coordinate systems that utilise landmarks such as the pectoral muscle boundary, nipple location and parabolic breast shape (Brandt et al 2011).

One of the earliest non-rigid breast image registration algorithms, which has also found wide application to other medical image registration problems, was that developed by Rueckert et al (1999). They used a free-form deformation model based on B-splines in which the coordinates of a set of control points were optimised to produce a smoothly varying deformation. This method has also been adapted to register MRI and PET-CT (Dmitriev et al 2013). Application of this transformation model to mammogram registration has generated conflicting results (Diez et al 2010, Pereira et al 2010). Closely related thin-plate splines (Bookstein 1989) have also been used to model the compression deformation of the breast during x-ray mammography (Behrenbruch et al 2003).

The class of non-parametric, variational registration methods typically involves computing a dense displacement field, at points arranged in a regular grid superimposed on the images. This grid will generally be derived from the source or target image's voxel matrix (or some sub-sampling thereof) and results in a large number of unknowns: the number of grid points multiplied by the dimensionality of the images. Given the large number of unknowns, and to ensure that the optimisation is well-posed, these algorithms are frequently characterised by the regularisation employed i.e. diffusion, elastic (Fischer et al 1999), fluid (Crum et al 2005), curvature (Modersitzki 2004) etc. Linear elastic registration has been used to register both DCE-MRI (Bruckner et al 2000) and follow-up MRI with images from a previous visit (Boehler et al 2010). A combination of diffusion and curvature regularisation has also been used to register DCE-MRI (Melbourne et al 2011b). Optical flow is a related method that has been applied to MRI motion correction (Froh et al 2006, Botella et al 2010). Constraining a registration transformation to be smoothly varying, to avoid folding and so-called non-diffeomorphisms, is an issue addressed by the div-curl transformation model adopted by Chu et al (2009) to register DCE-MRI. A recent approach applied to prone-supine registration of breast CT images uses a 'dual-dynamic' grid transformation model, specified over an irregular grid and incorporating an elastic regulariser, to tackle the problem of large deformations and disappearing structures (Alderliesten et al 2013).

The recent review of Boehler et al (2012) describes the range of breast image registration techniques in more detail. The current review focuses on biomechanical deformation models as a means of generating physically realistic and well-constrained deformations to address the fundamental challenge inherent in breast image registration problems: gross, non-rigid, soft-tissue deformation of the breast.

3. Breast biomechanical modelling

Biomechanical modelling and simulation of biological tissues has been the subject of intense research within the past three decades by various research groups. It has been investigated for various medical applications, such as surgical procedure training, pre-operative planning, diagnosis and clinical biopsy, image guided surgery and image registration, material parameter estimation etc.

For example, in surgical planning, biomechanical modelling could provide realistic and accurate evaluation of virtual biological tissue deformations, enabling potential outcomes of various surgical strategies, or alteratives that fall within the range of available surgical options, to be compared. The goal of surgical training simulators is to provide physicians with an integrated interface that can support force feedback through reliable biomechanical simulation, in order to practice and plan surgical procedures of various medical tasks. For diagnosis and clinical biopsy, reliable and efficient simulation tools for tissue deformation and motion estimation, are of particular interest. For example, a biomechanical model could be used to define the three-dimensional location of a lesion or model its temporal evolution within the human body. To that end, modelling can be used standalone or in conjunction with imaging procedures for treatment purposes.

Likewise, biomechanical models can be used for image guided surgery and image registration to realistically warp pre-operative image data to match the surgical situation. This overcomes the limitations inherent in acquiring pre-operative images in an alternative pose to that used in surgery and accounts for the fact that human organs shift and deform during surgery by updating images concurrent with surgery using physically realistic models. Moreover, modelling tools in conjuction with non-invasive imaging technologies (e.g. shear-wave elastography, dimensional time resolved 3D ultrasound) and in vivo experiments (e.g. identation measurements, aspiration techniques) can be also used for tissue health monitoring and biomechanical material response characterization.

In summary, soft tissue simulators aim to model living biological tissues to an appropriate level of accuracy with respect to physical, mechanical, and anatomical properties, to address the following tasks effectively:

modelling of patient-specific data;
representation of the internal structure of specific biological tissues and/or organs of interest;
prediction of the biomechanical behaviour and deformation of tissues and organs according to the relevant clinical conditions (e.g. multi-modal diagnosis, surgical planning and image guided biopsy);
computational efficiency to provide users with results in a timeframe that is compatible with the relevant clinical workflow.

Various physics based computational methodologies have been proposed for the purposes of biomechanical modelling and soft tissue deformation simulation. The most common numerical procedures are based on linear or nonlinear biomechanical models, which include the mass-spring method (MSM) (Roose et al 2005, Chang et al 2010, Patete et al 2013), the mass-tensor method (MTM) (Cotin et al 2000, Picinbono et al 2003, Roose et al 2006), the point-associated finite-field approach (De et al 2006), the boundary element method (Greminger and Nelson 2003, James and Pai 2005), and the conventional finite element (FE) method (see cited papers in sections 3.1 and 3.2).

3.1. Biomechanical constitutive models

The tissue composition and the macro-/micro-structure of the breast directly contributes to its mechanical behaviour, therefore accurate and realistic biomechanical modeling of the breast requires understanding of its anatomical characteristics. The adult female breast is predominately composed of glandular lobules, milk ducts, adipose, connective tissues, and skin (figure 1). The glandular lobules and milk ducts—usually referred to as fibroglandular tissue—are surrounded by dense connective tissues that maintain the breast integrity. Each breast is attached to the chest by the pectoral fascia over the pectoralis major muscles, while its shape is established and maintained by the skin. The breast changes under hormonal influence and due to ageing. For example, the ductal and glandular structure increases in size during the monthly premenstrual phase—contrary to the postmenstrual phase—while adipose is being almost completely replaced by ductal and glandular tissue during pregnancy. After menopause, fibroglandular tissues shrink and fat becomes predominant. Evidently, due to these changes, breast tissue biomechanics are expected to change. For example, in vivo magnetic resonance elastography (MRE) measurements have revealed that the distensibility of breast tissues in healthy volunteers depends on the menstrual cycle (Lorenzen et al 2003).

**Figure 1.** Adult female breast anatomy illustration.
Download figure:
Standard image High-resolution image

In general, experimental studies measuring mechanical properties of breast tissues have been mainly motivated by diagnostic requirements, i.e. identifying abnormal tissues by assessing tissue stiffness/elasticity (Sarvazyan 1993, Sarvazyan et al 1995, Krouskop et al 1998, van Houten et al 2003, Sinkus et al 2005, Samani et al 2007, O'Hagan and Samani 2009, Parker et al 2011, Doyley 2012). A summary of the most commonly used constitutive models and the associated tissue properties proposed for the biomechanical description of female breast tissues is outlined in table 1.

Table 1. Exemplary biomechanical models and tissue properties of the female breast.

Reference	Model	Procedure	Material properties (kPa)
Reference	Model	Procedure	Adipose	Fibroglandular	DCIS
Krouskop et al (1998)	Bi-linear Elastic (E₁/E₂)	Ex vivo Dynamic Compression	$18.\pm 7.$ ^a/ $20.\pm 8.$ ^b	$28.\pm 14.$ ^a/ $48.\pm 15.$ ^b	$22.\pm 8.$ ^a/ $291.\pm 67.$ ^b
Wellman et al (1999)	Linear Elastic (E)	Ex vivo indentation tests	$4.8\pm 2.5$ ^c/ $6.6\pm 7$ ^d	$17.5\pm 8.6$ ^c/ $33.0\pm 12.$ ^d	71.2^c/188.7^d
			$10.4\pm 7.9$ ^e/ $17.4\pm 8.4$ ^f	$88.1\pm 66.7$ ^e/ $271.8\pm 167.7$ ^f	638.7^e/2162.1^f
	Exponential ( ${{b}^{\star}}$ ; ${{m}^{\star}}$ )		$4460.\pm 2345.6$ ; $7.4\pm 4.$	$15\,174.5\pm 6750.7$ ; $12.3\pm 7.4$	$55\,776.$ ; 24.4
Azar et al (2001)	Exponential (b_n; m_n)	Adapted from Wellman et al (1999)	4.46; 7.4	15.1; 10.	N/A
Samani and Plewes (2004)	Mooney-Rivlin (C₁₁) (C₁₀; C₀₁) (C₂₀; C₀₂)	Ex vivo indentation tests and inverse analysis	$2.25\pm 0.3$ $0.31\pm 0.03$ ; $0.3\pm 0.02$ $3.80\pm 0.6$ ; $4.72\pm 0.7$	$4.49\pm 0.8$ $0.33\pm 0.04$ ; $0.28\pm 0.03$ $7.72\pm 1.1$ ; $9.45\pm 1.3$	N/A

O'Hagan and Samani (2009)	Polynomial Hyperelastic (C₁₁) (C₁₀; C₀₁; C₂₀; C₀₂) Ogden (μ; α) Yeoh (C₁₀; C₂₀; C₃₀) Arruda–Boyce (μ; λ) Veronda-Westmann (C₁; C₂; α)	Ex vivo indentation tests	N/A N/A N/A N/A N/A	N/A N/A N/A N/A N/A	52.6 0.975; 0.975; 0.65; 0.65 4.84; 51.228 2.09; 33.2; 0.002 5.; 4.774 2.034; −0.506; 1.51
		Ex vivo indentation tests	N/A N/A N/A N/A N/A	N/A N/A N/A N/A N/A
Han et al (2012)	Neohookean Iso. (μ)	MRI and material Optimisation	$3.398\pm 0.047$	$56.859\pm 69.491$	247.5
	Neohookean Aniso. (μ; η)		$3.371\pm 0.046$ ; $30.8\pm 29.$	$51.8\pm 94.2$ ; $1289.5\pm 2635.3$	328.9; $13\,239.4$

^a5% pre-compression, 0.1 Hz loading frequency. ^b20% pre-compression, 0.1 Hz loading frequency. ^cat 1% strain. ^dat 5% strain. ^eat 10% strain. ^fat 15% strain.

3.1.1. Ex vivo biomechanical properties.

Mechanical properties of ex vivo breast tissues have previously been measured through tension and compression experiments. Uniaxial compression experiments of 150 ex vivo specimens of healthy, cancerous and fibroadenomatous tissues showed that fibroadenomas are four times stiffer than normal tissues, and cancerous tissues are as much as seven times stiffer than normal ones (Sarvazyan et al 1995). Recognizing the potential effect of strain level, Krouskop et al (1998) measured the elastic moduli of 142 ex vivo small samples of healthy and pathological breast tissues at two different strain levels: low-strain level (5%–15%) and high-strain level (20%–30%). With low frequency sinusoidal loads imposed by a circular indenter, experiments were performed on breast tissue slices, which were assumed to be in a semi-infinite elastic solid domain, so as to evaluate Young's modulus of elasticity. Krouskop and his colleagues observed that the elastic modulus of fat tissue is constant over the strain range, while elastic moduli of the other tissue groups (glandular tissue, fibrous tissue, ductal carcinoma, invasive and infiltrating ductal carcinoma) are strain dependent. Adipose gave the lowest modulus measurements, while fibrous tissues were one to two orders of magnitude stiffer than fat. Glandular tissue had an elastic modulus similar to that of fat at low-strain levels but one order of magnitude stiffer at a higher-strain level. Carcinomas were stiffer than the other tissues at the higher-strain level; intra-ductal in situ carcinomas were similar to fat at the low-strain level and much stiffer than glandular tissues at the high-strain level. Furthermore, infiltrating ductal carcinoma was much stiffer than any of the other breast tissues, giving twice the stiffness of fibrous tissue at the high-strain level. Thus, the dependence of modulus of elasticity on the two strain levels reveals the nonlinear elastic behaviour of breast tissues.

Wellman et al (1999) experimentally measured Young's modulus of ex vivo breast tissues. They found a significant difference in the stiffness distribution, and an increased rate of change of stiffness with strain, between cancerous and benign breast tissues. Infiltrating ductal cancer was more than ten times as stiff as normal fat tissue and two and a half times as stiff as glandular tissue at 1% strain, and more than seventy times as stiff as normal fat tissue and five times as stiff as glandular tissue at 15% strain. However, fat tissue only showed a three-fold increase in stiffness from 1% to 15% strain regimes.

Since ex vivo breast tissue specimens tend to be small and no analytical expressions are currently available to directly extract the elastic modulus, E, from force-displacement measurements, finite element (FE) based inverse analysis methods have been proposed to estimate E through indentation testing on small tissue specimens (Samani et al 2003, Samani and Plewes 2007). The developed methods have been used to extract elasticity measurements of 169 ex vivo breast tissues and tumours. Under small deformation conditions, Young's modulus of normal breast fat and fibroglandular tissues was found nearly equivalent (around 3 kPa), while fibroadenomas are approximately twice as stiff as in fibrocystic disease. Malignant tumours exhibited a three to six fold increased stiffness with high-grade, invasive, ductal carcinoma exhibiting up to a thirteen fold increase in stiffness compared to fibroglandular tissue. Along these lines, Samani and Plewes (2004) reported an inverse problem, FE solution strategy to extract the hyperelastic parameters of ex vivo breast tissues through indentation testings. Based on the data of ex vivo compression experiments published in the literature, Young's modulus of adipose varied from 0.5 kPa to 25 kPa and in glandular tissue from 2 kPa to 66 kPa (Gefen and Dilmoney 2007).

3.1.2. In vivo biomechanical properties.

Mechanical properties of in vivo breast tissues have primarily been measured with elasticity imaging techniques, such as, sono-elastography, magnetic resonance elastography, shear wave elasticity imaging and mechanical imaging (Parker et al 2011, Sarvazyan et al 2011, Doyley 2012). Elasticity parameters of in vivo breast tissues were extracted from displacement/strain fields based on inverse problem solution techniques (Doyley 2012), or from shear-wave velocity measurements. These experimental methods were typically performed with breast tissues subjected to small deformations, whereas in Han et al (2003) the possibility of using clinical ultrasound probes as indenters for large deformation measurements was demonstrated.

In an early study involving six healthy volunteers and six patients with malignancies (McKnight et al 2002), experimental evidence of the mean shear modulus measured from MRE showed that breast carcinoma was 4.2 times as stiff as surrounding breast tissues. Similarly, Lorenzen et al (2002) used MRE procedures to assess the viscoelastic properties of breast tumours and surrounding tissues in a clinical study involving fifteen patients with malignant breast tumours, five patients with benign breast tumours and fifteen healthy volunteers. They found that in general, malignant invasive tumours (median 15.9 kPa) were significantly stiffer than benign lesions (7 kPa). However, an overlap in the elasticity ranges of soft malignant tumours and stiff benign lesions was observed, hence obscuring tissue differentiation through MRE mapping. However, the results of fifteen patients reported in Sinkus et al (2005) showed that there was a good separation between breast cancer and benign fibroadenoma when measuring the shear modulus through MRE. More recently, Athanasiou et al (2010) employed a supersonic shear imaging technique to quantitatively assess the correlation between the stiffness of breast lesions and pathologic results in a study of 46 women with 48 breast lesions (28 benign, 20 malignant). The stiffness of malignant lesions (146.6 kPa) was found to be higher than benign ones (45.3 kPa).

Besides pathological breast tissues, elastic parameters of normal breast tissues were also investigated with elasticity imaging techniques (McKnight et al 2002, van Houten et al 2003, Sinkus et al 2005, Tanter et al 2008). The results revealed that gland/fibroglandular tissues were stiffer than fat, and elastic moduli of in vivo breast tissues measured from elasticity imaging techniques have relatively good correlations with those measured through ex vivo mechanical testing procedures (Krouskop et al 1998, Wellman et al 1999). Recently, Gamage et al (2011) proposed a non-invasive in vivo material parameter identification technique, where multiple gravity loading deformation modes were used in combination with nonlinear optimization techniques and finite element modelling. Their approach was validated by performing tests on a two-layer heterogeneous silicone gel phantom geometry. Although they demonstrated that three gravity loading orientations were sufficient to identify the heterogeneous constitutive parameters and maximize model predictability, their modelling was based on a rather simple constitutive material description.

The methods for characterizing the mechanical properties of healthy and pathological breast tissues have been reviewed above. Based on these methods, both ex vivo and in vivo estimations of breast tissue elastic properties have been reported in the literature. Nonetheless, this existing data is insufficient for high fidelity biomechanical simulations, i.e. very accurate numerical predictions. This is because in vivo measurements were typically limited to small-strain experiments, and mechanical properties of ex vivo samples are likely to differ from in vivo tissues, e.g. due to fixing of tissue samples, stress relaxation effects, etc. Moreover, it is generally accepted that the mechanical properties of living tissues vary widely across the population and over time for an individual. However, these data, and their relative values, are still useful for quantitative biomechanical analyses and could be used as an initial estimate for material parameter optimization in patient-specific biomechanical modelling (Rajagopal et al 2008, Han et al 2012), and in in vivo estimation of breast tissue elastic constants (Lago et al 2015). In contrast to breast tissues, experimental studies of in vivo and in vitro human skin have been more extensive. Mechanical measurement methods and mechanical properties of human skin can be found in the cited papers (Edward and Marks 1995, Elsner et al 2001, Pierre and Philippe 2004, Hendriks 2005, Geerligs 2006, Groves et al 2013).

3.1.3. Constitutive models of breast tissues.

In general, biological soft tissues, and more precisely breast tissues, are observed to exhibit nonlinear, anisotropic and time-dependent mechanical response (Han et al 2003). However, in Krouskop et al (1998) and Wellman et al (1999), it was found that viscous effects were negligible when the mechanical loads are applied within short time scales. To model the mechanical response of breast tissues to applied loadings, various material constitutive models have been used, i.e. linear or nonlinear elastic relationships, depending on deformation magnitude and the purpose of specific applications. Amongst these, hyperelastic (or Green-elastic) material descriptions are most commonly used in breast biomechanical modelling (Bakic 2000, Tanner et al 2006), as shown also in table 1. For hyperelastic materials, the constitutive relationship is derived from a strain-energy function (Malvern 1977). Particular forms of the strain-energy potential include the polynomial form and its variations (the reduced polynomial form, the Neohookean form, the Mooney–Rivlin and the Yeoh form), the Ogden, the Arruda–Boyce, and the Van der Waals model (Ogden 1984, Holzapfel 2000). Anisotropic linear elastic models and anistropic hyperelastic models have also been used to describe the anisotropic elastic behaviour due to the existence of Cooper's ligaments (Han et al 2012).

The published experimental data of ex vivo breast tissues have been used to determine model parameters of the material constitutive models. For example, Azar et al (2001) fitted the experimental curves of Wellman et al (1999) in order to determine the modulus of elasticity of glandular tissues, in which it is described as an exponential function of strain. Yin et al (2004) used corresponding moduli of fibroglandular and adipose tissues measured by Krouskop et al (1998) to define two parameter Mooney–Rivlin type hyperelastic models for adipose tissue and for fibroglandular tissue. Furthermore, Pathmanathan et al (2008) used a second-order polynomial hyperelastic model to describe fibroglandular and adipose tissues, while the model parameters were experimentally determined by Samani and Plewes (2004), and an incompressible exponential hyperelastic material law for skin was employed (Veronda and Westmann 1970).

3.2. 'State-of-the-art' in breast biomechanics

3.2.1. Breast deformation modelling.

To date, various numerical procedures have been proposed to model the biomechanics of adult female breast and simulate breast deformation during mammography screening, MR guided biopsy techniques and image registration. Amongst the pioneering works in patient-specific finite element (FE) modelling of breast tissue deformations was that of Samani et al (2001). They proposed two hexahedral mesh generation methods: one that was voxel based where FE meshing is accomplished by down-sampling MR images, and a surface-based meshing procedure by transfinite interpolation (Knupp and Steinberg 1993). In a numerical example, they presented hexahedral meshes of a breast from MR images with the two meshing techniques and demonstrated breast deformation under plate compression explicitly modelled as a contact problem. Elastic moduli of fat and fibroglandular tissues were experimentally fitted to polynomial functions of strain, while skin was considered to be elastic and modelled using membrane elements. Although their numerical scheme performed well, the authors mentioned that more reliable material properties are required to validate their results quantitatively and increase the accuracy of their model.

Motivated by the requirement to estimate the location and extent of a tumour during mammographic compression associated with needle biopsy, Azar et al (2001) and Azar et al (2002) proposed an FE-based methodology to model breast deformations and predict the tumour position during plate compression. Three-dimensional hexahedral FE meshes were generated from multiple two-dimensional slice contours. The breast was deformed using virtual compression plates by applying displacements to the surface nodes. To reduce the total computational time, the large deformation compression was divided into a number of small increments using small-strain considerations. Non-constant Young's modulus, as a function of the strain tensor, was used for modelling the nonlinear behaviour of fat and glandular tissues while the presence of skin was neglected. Compression simulations of three patient-specific breasts produced reasonable lesion estimation results.

In the works of Ruiter et al (2002) and Ruiter (2003), a method for automatic registration of x-ray mammograms and MR volumes of the female breast was developed using a biomechanical model of the breast. This study used the commercial finite element platform ANSYS for the mammographic compression simulation of MR reconstructed volumes of the undeformed breast. In Ruiter et al (2002), different tissue models (i.e. inhomogenous nonlinear, homogenous nonlinear, and homogenous linear models) for the breast tissues and boundary conditions (BCs) were considered in a parametric modelling study. They observed that the inhomogenous nonlinear material model assumption is the least accurate, while the rest perform quite similarly, although their simulations were confined to less than 25% breast compression. In their analyses the authors assumed identical material parameters for the lesion and the healthy breast tissues, although it is known that some carcinomas are an order of magnitude stiffer.

Rajagopal et al (2004) reported a biomechanical modelling framework of patient-specific geometries of the breast anatomy. They formulated a semi-automatic procedure to generate patient-specific finite element grids from magnetic resonance images, while they utilized a Total Langrangian formulation in order to numerically solve the equilibrium equations. The FE solution procedure was validated experimentally on silicon gel samples subjected to natural gravitational loading. Additionally, they attempted to investigate the material assumptions for soft tissues (i.e. isotropy, homogeneity and incompressibility) made in the literature at that time.

In del Palomar et al (2008), an FE-based approach to construct models that could be employed to predict patient-specific breast shape after lumpectomy interventions was presented. To achieve this, computed tomography images were used to reconstruct the shape of the breast. In their model, they didn't differentiate fatty from fibroglandular tissues; thus, breast tissues were assumed to be a region of homogeneous material by assigning average material properties for the breast tissues. Two patients were chosen in order to adjust the material properties and to further validate their FE model. A woman with a high proportion of fat tissue was considered and the material constants were fitted using several landmarks located on the patient's breast. Then, these material constants were used for the breast deformation analysis of the second patient with a different proportion of fat/fibroglandular tissue, and the final shape of the breast was compared with 3D scanner images of the patient in standing position. Their numerical results demonstrated the efficiency of their model to predict breast deformations at various positions of the patient with reasonable accuracy.

Chung et al (2008a) developed a biomechanical model of the breast to simulate compression during mammographic screening. MR images of the uncompressed breast were segmented into skin and pectoral muscle, in order to reconstruct the Hermite element grid of a 36-year old volunteer's breast geometry using a nonlinear geometric fitting process; while special numerical treatment was applied to the FE model to account for contact. Their model was reported to reliably reproduce the compressive deformation of the breast during mammography (up to 32% compression). The accuracy was assessed by comparing the predicted surface data and the locations of three internal features within the compressed breast, which has been further tested and validated through experimental evidence on soft gel phantoms. In this regard, the FE predictions of the surface deformation yielded a root-mean-square error of 1.5 mm, while the Euclidean error in evaluating the locations of three internal landmarks ranged between 4 and 6.5 mm. Furthermore, the FE numerical framework of the above mentioned authors was later refined to account for frictional contact mechanics when modelling the interaction between the breast skin and compression plates (Chung et al 2008b). In this work, the authors investigated two frictional contact formulations based on Coulomb's frictional law, while the simulation results were compared with experiments on a silicon gel phantom. Both numerical approaches produced similar accuracy results of the breast surface deformation (root-mean-square errors were less than 3 mm approximately for both cranio-caudal and mediolateral compression).

In the same year, Pathmanathan et al (2008) proposed a FE numerical procedure to model breast deformation during cranio-caudal and mediolateral-oblique mammographic screening. Their breast model was constructed from MR images, where fibroglandular tissues were segmented from adipose tissues into their model, while the skin was also accounted in their analysis. Nonlinear constitutive relations were used for the biomechanical description of breast tissues and the skin (Samani and Plewes 2004). The numerical procedure involved two consecutive steps: a backward analysis in order to evaluate the zero-gravity configuration of the breast, followed by a forward analysis of breast deformation due to the compression pads. In the latter analysis step, the interaction between breast skin and compression pads assumed zero friction. Their approach, nonetheless, differs from the previous works in the fact that it recovers the undeformed configuration of the breast. However, as the authors argue, their numerical scheme was not validated using clinical data and they raised a concern about the material parameters utilized in their analysis.

In the work of Tanner et al (2011), individual-specific breast compression computer simulations of eight female volunteers were investigated with particular emphasis on breast tissue material parameters specification. The realism of mammography simulations was evaluated by comparing the breast shapes of simulated and real mammograms. In the numerical simulations, the performance of isotropic and transverse isotropic material models to predict the displacement of internal landmarks was compared. The mean displacement error of several landmarks was observed to reduce for isotropic materials after optimizing the material parameters with respect to breast surface alignment and image similarity. However, when the authors specified transverse isotropic breast tissues (with an increased stiffness in the anterior-posterior direction), they obtained more realistic mammogram simulations from MR images. They also observed that the amount of anisotropy had a significant impact on the numerical predictions, whereas homogeneous material models under-performed significantly.

Subsequently, Han et al (2012) proposed an integrated biomechanical modelling framework for surgical simulations, supported by a fast graphics processing unit based explicit nonlinear FE solver. Patient-specific FE model generation was accomplished using a semi-automatic segmentation method for tissue classification and a fully automated FE grid generation. The main contribution of this work was the development of an optimization algorithm for the numerical determination of breast tissues' material parameters through FE simulations. The optimized deformation prediction was obtained by iteratively updating the material parameters, to maximize the geometry similarity (i.e. via the location of various landmarks) between the FE predicted breast shape and the experimentally acquired MRI scans. Their numerical methodology was validated and tested through breast compression simulations for five patient-specific cases, and their results revealed that both the heterogeneity and the anisotropy of soft tissues are essential in predicting accurate large breast deformations under plate compression, despite the fact that the skin wasn't explicitly modelled in their analyses.

Costa (2012) proposed an interesting alternative numerical approach for large deformation analysis of soft tissues. Breast tissue bulk was modelled as an incompressible fluid medium, while a dense network of distensible elastic fibers connecting the surface (boundary) vertices was also considered. The deformed state of the mesh was computed by taking the quasi-static equilibrium of the internal forces, due to fluid pressure and fiber tension, with external forces acting in an area associated with each superficial vertex node. In this approach, as opposed to the MSM, the mass is distributed in the entire object by the fluid density. The author argues that the proposed methodology can easily incorporate anisotropic and nonlinear behaviour of soft tissues, and performs independently of the mesh resolution. However, this approach is limited in modelling objects filled with fluids and no dynamic simulations could be performed. Costa also attempted to simulate breast deformation during mammography, while his methodology was validated through an ex vivo compression experimentation on bovine liver.

Recently, Kuhlmann et al (2013) presented a breast biomechanical model to predict soft tissue deformation during imaging. Their numerical method was founded on a coupled Eulerian-Lagrangian FE approach that, as argued by the authors, would allow for a more adequate representation of very large deformations experienced by soft living tissues. The internal tissues of the breast were represented as an Eulerian material and enclosed by a Lagrange-membrane representation for the skin. Nonlinear constitutive relations were adopted from the literature for mechanically describing the breast tissues and the skin. MR images were utilized for generating patient-specific breast geometries. The authors validated their model by comparing the numerical results with experimental laser-scan measurements of a subject through a microwave breast imaging system. Furthermore, they simulated cranio-caudal mammographic compression using two rigid plates. However, as reported in their paper, although this numerical approach provides significant improvements for stable, large deformation simulations even for highly incompressible materials, the computational time cost is remarkably high due to the mesh resolution requirements, when compared to pure biomechanical models.

Patete et al (2013) developed and validated a three-dimensional female breast deformation model that was based on MSM for computer assisted breast surgery. MR scans of a patient—for breast biopsy preparatory analysis—were acquired using two different datasets: the uncompressed breast, and the compressed breast due to lateral single breast compression. The acquired images were segmented using a semi-automatic procedure, and then a tetrahedral based grid representing the skin, fat and glandular tissues was generated. Tissue modelling was performed through a two-step algorithm. First, an iterative approximation algorithm was implemented in order to estimate the springs' rest length and evaluate their stiffness due to the inherent gravitational loading on the uncompressed breast grid. Subsequently, the resulting spring parameters were used to carry out the final deformation simulation and hence recover the original deformed breast grid. Numerical results revealed that the applied gravity force was compensated by the internal elastic forces of the springs' network, thus leading to maximum relative distance error between the two configurations of less than half a millimeter. However, it was reported by the authors that further investigation on a larger dataset is required to assess the accuracy and robustness of the model to predict breast deformations.

3.2.2. Predicting the undeformed breast shape.

As discussed in section 1, conventional image acquisition and analysis techniques are unable to fully and accurately predict large breast deformations at various positions of the patient. A prerequisite of this approach, however, is to assess the loading state of the individual-specific breast geometry in order to ensure the validity and correctness of the mathematical modelling procedures. This is especially true for prone MRI, which delivers good soft tissue contrast and is hence a very good basis for a three-dimensional biomechanical model, but shows the breast in a freely pendulous position with gravity pulling the tissues in the anterior direction. Thus, it is vital to evaluate breast tissue pre-stressing due to the in vivo conditions (i.e. gravity) and accurately predict the unloaded (i.e. gravity-free) configuration of an individual-specific breast geometry. This section attempts to review the relevant research into this topic.

Rajagopal et al (2007) reported a pure Lagrangian iterative scheme that numerically evaluates the unloaded configuration of a solid body. The proposed formulation is founded on a standard forward finite deformation approach, where the only unknown for the problem is the undeformed (or reference) configuration (as opposed to the deformed state in conventional forward analysis). In this approach, an estimated undeformed state is perturbed and the residual of the unbalanced forces is evaluated repeatedly until convergence is achieved. The authors validated their FE methodology with an analytical solution and verified their predictions with experimental results on phantom geometries.

Similarly, Carter et al (2008) proposed an iterative FE methodology to evaluate the unloaded shape of individual-specific breast geometries. Initially, an estimate of the reference state is assumed having zero internal stresses and applying gravity in the anterior direction. Then, in the posterior direction, gravity is re-applied to the reference state model, and the nodes' spatial location is compared between the deformed reference model and the supine model. This process is repeated until the difference between the two configurations falls below a predefined tolerance. This method was later refined by Eiben et al (2013), in which an initial guess of the zero-gravity state is established by inverting the direction of gravity and relaxing the stresses on each breast model. The corresponding configurations are then re-loaded, and an iterative procedure is followed where the estimated zero-gravity configuration is updated. The update is performed by pulling the difference vector between the original loaded model and the re-loaded model into the unloaded configuration and adding a scaled version of this vector to the node position in the unloaded configuration. The iterative process ends when convergence is met between the original breast (prone or supine) model and the estimated gravity-free model. Eder et al (2014) proposed a similar iterative approach to predict the unloaded configuration from a prone breast MRI but instead of using the displacement vectors directly to update the zero gravity state estimate, they applied the reaction forces that arise when the re-loaded configuration is displaced to the real loaded MRI configuration. These forces are then used in the unloaded configuration to calculate the corresponding node update.

Eiben et al (2014) assessed the performance and accuracy of three algorithms for recovering the unloaded state of patient-specific breast geometries. These algorithms were the simple inversion of the direction of gravity, an inverse finite deformation approach (Govindjee and Mihalic 1996) and a fixed-point type iterative approach (Eiben et al 2013). They showed that the simple inversion of the direction of gravity does not provide high accuracy estimates for the unloaded breast shape especially for larger breast volumes. The iterative approach on the other hand generated results comparable to those produced with the inverse finite deformation approach. Furthermore a sensitivity analysis study using a phantom geometry, having as known the ground-truth unloaded breast configuration and tissue properties, was carried out. They reported that the simple inversion of gravity approach is more sensitive than the 'pure' inverse analysis methodologies with respect to the material parameters assumed in biomechanical modelling.

More recently, Vavourakis et al (2015) proposed a generalized mixed FE formulation for the evaluation of the unloaded configuration of biological soft tissues and organs, which exhibit incompressible behaviour under finite deformations (figure 2). In their study, in vivo patient-specific geometries—acquired through medical imaging—were utilised in simulating the recovery of the pressure-free configuration of human aortas and the gravity-free shape of the female breast; where as far as breast unloaded shape is concerned, they considered in their analysis the spatial-varying structure of the organ and the biomechanics of the skin.

4. Biomechanically informed breast image registration for multi-modality imaging and image guided interventions

In this section we review the state-of-the-art in biomechanically informed breast image registration. These methods are categorised according to the specific transformations addressed, which in turn correspond to one or more clinical applications. Section 4.1 introduces prone to prone image registration for the purposes of aligning (a) the temporal frames of a dynamic contrast enhanced MR sequence and (b) MRI and PET images, both of which are acquired in the prone position. For completeness we include a discussion of the validation of prone-prone breast image registration algorithms (section 4.2) in which the methods discussed, whilst not utilising biomechanical transformations directly in their methods, were however validated using biomechanical simulations. In section 4.3 prone to supine breast registration is considered. This covers registration of prone MRI to supine US, CT or MRI. CT is routinely used in radiotherapy planning and relating this image to the pre-operative MRI (if available) would help direct the radiotherapy dose to the most critical regions within the breast. Whilst not currently standard clinical practice, an additional supine MRI acquisition is proposed, in the context of surgical planning, as a means of providing a pre-operative estimate of the initial breast position at the start of surgery. Finally section 4.4 reviews methods developed to address applications in which the breast was imaged with and without compression. This equates to registration of MRI or Ultrasound Computed Tomography (USCT), in which the breast is imaged prone and freely pendulous, to x-ray mammography, in which the breast is immobilised via compression between parallel plates.

4.1. Prone to prone breast image registration

The fact that the magnitude of the deformation expected in the prone-to-prone registration problem is relatively small might explain why few biomechanical approaches exist to inform this type of registration problem. Roose et al (2008) for example investigate how boundary conditions can be applied to a biomechanical mass-tensor model in order to obtain the alignment of prone MR images. Displacement boundary conditions applied on the skin surface are an obvious choice to push the model into the position of potential alignment. The interior of the breast in this case would be deformed, completely determined by the potential energy minimisation or relaxation of the breast tissue. Pure displacement boundary conditions are in general, however, problematic, as they require a point to point correspondence to be established in advance. This correspondence is difficult to establish between two surfaces as it is mostly unstructured. In Roose et al (2008) three different types of boundary conditions are proposed: i) sliding surface contact, ii) image intensity guided surface matching, and iii) combined surface attraction and image intensity guided matching. The results look best for the last approach which allows sliding of the surface but is driven by image intensities. A registration accuracy is given only in terms of an improvement in image similarities. However they clearly show that the biomechanically guided alignment procedure results in significantly less volume changes, compared to free form deformation registration, and therefore offers a more plausible solution.

Image registration between MRI and PET is inherently difficult since both modalities capture fundamentally different biochemical and physical tissue properties. In particular, registration algorithms have to account for the sparse, locally confined information of PET images. One approach to align highly different image information is to use a known commonality such as skin fiducial markers which are visible in both modalities. Based on this additional information Coman et al (2004), Unlu et al (2005) and Krol et al (2006) proposed a method whereby first corresponding markers are detected and from this sparsely distributed displacement vectors are calculated. To generate the deformation vector field which aligns the MRI with the PET image, the breast is segmented from the MR image and a finite element mesh is generated to discretise the domain. Three displacement components are regarded as three separate local temperatures which are imposed on the surface of the model. To interpolate the displacements across the whole domain, the Steady State Heat Transfer problem is solved. The authors highlight that no biomechanical model is required for this kind of registration. Subsequent work from this group (Unlu et al 2010) proposed an optional extension by incorporating the CT images from combined PET/CT acquisitions, if these are available. The CT image is in this case used to perform additional iterative surface matching between the MRI and the CT.

PET/CT acquisitions provide a desirable combination of imaging modalities as CT images—being acquired in the same patient position and imaging session as the PET images—provide essential information for the PET image reconstruction in terms of attenuation correction. Thus it is logical for image registration techniques to use the anatomically more informative CT images, in place of the PET images, to solve the PET-MR registration task.

Dmitriev et al (2013) overcome the highly different nature of DCE-MRI and PET images by using CT images as an intermediate modality. They fuse DCE-MRI and PET/CT images by multi-scale registration of the T1 weighted pre-contrast image of the DCE-MR image and the anatomical CT image. During both acquisitions the patient position is kept as similar as possible, i.e. prone. The PET image is assumed to be registered to the CT image as both are captured at the same time. They report that 94.2% of the 140 breast images from 70 patients were aligned with an error of less than 4 mm. The target registration error was evaluated by using five corresponding, manually picked landmarks within each breast. The authors employ a conventional B-spline transformation and list the lack of a biomechanical model as a limitation of their approach.

4.2. Validation of prone to prone breast image registration

Quantitative evaluation of image registration algorithms requires a known ground truth deformation which an algorithm aims to recover. This ground truth deformation can be provided via a set of spatially distributed landmarks or as a dense deformation vector field. Tanner et al (2007) chose biomechanical simulations to generate known deformations in order to quantify the performance of the well known B-spline registration algorithm proposed by Rueckert et al (1999) for dynamic contrast enhanced MR images. Intensity changes in the dynamic contrast series caused by the contrast agent are likely to be removed by intensity based image registration algorithms and thus constraints such as volume preservation are required. Hence, Tanner and her colleagues selected ten pre- and post-contrast image pairs from different patients which did not show any displacement of the tissue between the acquisition time points. The images were manually segmented into fat, fibroglandular and cancerous tissue and the breast region was discretised using second order tetrahedral elements. These elements were assigned with isotropic, linear elastic material properties according to the segmentation. Displacement boundary conditions were developed and imposed to simulate different types of motion with two amplitudes of 5 mm and 10 mm.

Hill et al (2009) followed a similar approach to examine the effect of motion on the evaluation of DCE-MRI enhancement characteristics. To this end, a patient-specific heterogeneous model was built. They observed a relationship between the maximum enhancement value from the DCE series, with the different tissue types (fat, fibroglandular and lesion) and used a linear mapping to convert maximum enhancement into Young's modulus (Hill et al 2008). Thus an explicit segmentation step is only required for the gross shape of the breast. They derived displacement boundary conditions from a dataset with visible movement at the breast-pectoralis interface and imposed this motion onto a different dataset, and corresponding model, where no motion was apparent. They conclude that small motion artefacts can influence the enhancement curve of a lesion significantly and thus image alignment of the contrast series should be performed in order to avoid misleading results.

4.3. Prone to supine breast image registration

This section reviews the methods developed to address prone to supine breast image registration; in particular, prone MRI to supine US, CT or MRI.

Prone-to-supine breast image registration is a particularly challenging task due to the large deformations which occur due to the change in the direction of gravity. Pure intensity based image registration algorithms are likely to fail as the initial overlap is too small to drive the registration in the correct direction. Thus biomechanical models have been used in registration strategies to incorporate the direction of gravity reversal as prior knowledge.

Early work by Carter et al (2008) aimed to make pre-surgical information from prone dynamic contrast enhanced Magnetic-Resonance-Imaging (DCE-MRI) available for image guided surgery. In a surgical setting, the facility to perform imaging is limited. Three dimensional surfaces from stereoscopic reconstruction techniques can be captured relatively easily as they are based on optical systems which only require an unobstructed view of the patient's chest. The main limitation of this modality however is that it provides surface information only. Carter et al (2008) thus follow a two step approach of (1) aligning prone and supine breast MR images of the same patient, where the latter image acts as an intermediate which is an addition to current standard clinical practice and (2) by aligning the supine MRI with the reconstructed chest surface.

As discussed in section 3.2.2 a number of authors have addressed the issue of computing the unloaded shape of the breast when developing prone to supine registration methods because pre-stresses in either pose of the breast are generally unknown (Rajagopal et al 2007, Carter et al 2008, Eiben et al 2013).

All biomechanical models and corresponding simulations will—to some extent—be inaccurate due to limited knowledge about the exact in vivo tissue response to loading, among other factors. The alignment accuracy achieved by simulation alone can be improved by intensity based image registration. Carter et al (2008) used a fluid registration approach whereas Lee et al (2010) incorporated a B-Spline registration scheme.

Eiben et al (2013) also used a B-Spline registration and showed that the registration can also be performed in the unloaded zero-gravity frame, claiming that this method is advantageous, especially for larger breasts, as folding is reduced in the virtual unloaded state. In this approach it is necessary to build a biomechanical model from the prone and the supine position.

Similarly, prone images are deformed into the supine position by Han et al (2014) by using FE simulations and a subsequent image registration step. While the pre-stresses from gravity are not taken into account, motion of the breast tissue along the chest wall is permitted by using frictionless boundary conditions. The approach further refines previous approaches since the material properties of the biomechanical model as well as the direction of gravity is optimised at the same time, where the optimisation criterion is the similarity measure of the final registration step.

In this work, the computational cost was relatively high (in the range of 3–4 h), despite the use of a high performance GPU-based FE solver. This is clearly incompatible with real-time/fast clinical applications.

The target information in all publications referenced in this section so far were spatially resolved images. In the current clinical work-flow however, supine images, to which a prone MRI could be deformed, are not usually acquired. In the context of surgical guidance, advances in affordable surface scanning techniques could potentially provide alternative target information to deform prone images into the supine position. Preliminary work by Lago et al (2012) imposed boundary conditions from simulated range data mimicking the depth field from a time-of-flight camera. Unfortunately no quantitative evaluation of this technique is included and results are only presented on the level of visual quality.

In addition to significant effort and much progress in the field of prone-to-supine image registration of the breast, and related biomechanical modelling, recent research has also considered options to avoid the large scale deformation registration problem by acquiring a pre-surgical contrast enhanced MRI in the supine position (Siegler et al 2012, 2011). However, it remains to be seen if the supine breast acquisition is accepted in clinical diagnostic imaging.

4.4. Image registration of the breast under compression

In this section we review breast image registration developments in which biomechanical modelling has been used to capture large deformations caused by compression of the breast, such as that observed in mammography (figure 3). For this purpose, patient-specific breast models have been developed and plate compression simulations performed to predict the deformed breast shape.

Methods that use a simplified generic breast model have also been proposed. For the CC to MLO compression (Kita et al 2001) proposed the use of quadratic equations to model the deformation of each curve within the breast from the compressed to the uncompressed state, while (van Schie et al 2011) used a semi-spherical model. A population Statistical Deformation Model (SDM) (Tanner et al 2009) was used for multimodal registration and finally a volume-preserving affine (Mertzanidou et al 2012a) and an ellipsoidal SDM (Mertzanidou et al 2011) were employed for MRI to x-ray registration.

We first review patient-specific breast compression modelling approaches and then refer in more detail to the utility of FE modelling within an MRI to x-ray mammography registration framework.

4.4.1. Patient-specific breast modelling for plate compression simulations.

A number of approaches to FE modelling of the breast for plate compression simulation have been proposed for a range of applications:

prediction of 3D deformations (Samani et al 2001, Pathmanathan et al 2004, Chung et al 2008b, Tanner et al 2011, Han et al 2012),
breast MRI to x-ray registration (Ruiter et al 2006, Mertzanidou et al 2012a, Hopp et al 2013a, Lee et al 2013),
the validation of x-ray mammography registration algorithms (Hipwell et al 2007),
biomechanical software phantom compressions for evaluation of tomosynthesis image processing and reconstruction algorithms (Ruiter et al 2008, Bakic et al 2011, Pokrajac et al 2012),
ultrasound CT to x-ray mammography registration (Hopp et al 2012b), and
validation of mammographic density methods (Alonzo-Proulx et al 2010).

In this section, we review the various patient-specific modelling techniques, in the context of image registration, according to the biomechanical modelling components presented in section 1. Sensitivity analysis of the various modelling parameters with respect to compression simulations has been performed using a three-dimensional breast phantom (Hsu et al 2011) and real patient data (Mertzanidou et al 2014).

Geometry and meshing.

As discussed previously, the geometry of a biomechanical breast model is frequently defined using MRI. The first step when building a biomechanical model therefore consists of segmenting the breast volume from the image in order to extract the surface and the volumetric mesh. This can be a challenging segmentation task, due to the poorly-defined boundary between the breast tissue and the pectoral muscle. Therefore, authors often used either manual or interactive approaches (Samani et al 2001, Tanner et al 2011, Han et al 2012, Mertzanidou et al 2014). However, automated approaches have also been proposed (Ruiter et al 2006, Hopp et al 2013a, Lee et al 2013, Solves-Llorens et al 2014).

The choice of the optimal number of nodes and elements in the volumetric mesh is an open research topic. Fine meshes, with a large number of nodes and elements, have the advantage of representing the breast anatomy and shape more accurately, but have the disadvantage of requiring longer computational times by the FE solver. Also, large deformations, such as mammographic compression, can cause element topology problems when applied to very fine meshes, especially in cases where folding is present. Folding can occur, particularly for large breasts, due to contact with the MRI breast coil or in cases where the subject is clothed during scanning.

In the literature authors have proposed a wide range of values for mesh node density, from very fine to very coarse meshes, but no experimental study has determined the optimal value for this application. When modelling breast plate compressions, the number of elements used varied from a few hundred (Ruiter et al 2006, Chung et al 2008b, Lee et al 2013), to several thousand (Mertzanidou et al 2014, Pathmanathan et al 2004) and finally tens of thousands (Samani et al 2001, Tanner et al 2011, Han et al 2012, Hopp et al 2013a, 2013b, Solves-Llorens et al 2014). The majority of FE analyses have been performed on tetrahedral meshes, however hexahedral meshes have been proposed by Chung et al (2008b), Lee et al (2013), Pathmanathan et al (2004) and Ruiter et al (2006).

Material properties.

When modelling large breast compressions authors have chiefly proposed non-linear material models, as linear elastic models are considered less accurate for large deformations. More specifically, a hyperelastic, (nearly) incompressible material was proposed by Chung et al (2008b), Han et al (2012), Hopp et al (2013a), (2013b), Lee et al (2013), Mertzanidou et al (2014), Ruiter et al (2006), Samani et al (2001) and Solves-Llorens et al (2014). Linear elastic models where used by Tanner et al (2011) and Alonzo-Proulx et al (2010), however in the experiments described by Tanner et al (2011), the displacements of all the surface nodes were known and constrained. Therefore for this particular application the effect of the material model is expected to have less influence on the results.

In Alonzo-Proulx et al (2010), Han et al (2012), Ruiter et al (2006), Samani et al (2001), Solves-Llorens et al (2014) and Tanner et al (2011) three different tissue classes were used: fibro-glandular, fat and skin, whereas Chung et al (2008b), Hopp et al (2012b, 2013a), Mertzanidou et al (2014) adopted a single homogeneous tissue type and Hopp et al (2013b) vary Young's modulus across the breast according to the speed of sound measured using ultrasound CT. The experimental work of Ruiter (2003) showed no significant effect on the results when different tissue models were used. The work of Ruiter et al (2006) and Tanner et al (2011) reported an anisotropic behaviour of the breast when applying mammographic compressions, with a reduced elongation in the anterior-posterior direction and an increased expansion in the medial-lateral direction, for a CC-view compression. This was incorporated into the modelling performed by Han et al (2012), Mertzanidou et al (2014) and Tanner et al (2011).

In most approaches the material parameters of the breast tissue were taken from the literature on studies of ex vivo tissue samples, however a patient-specific in vivo parameter estimation method was proposed in Han et al (2012) and Mertzanidou et al (2014).

Compression simulation and boundary conditions.

The plate compression that occurs in mammography has been modelled in the literature, primarily using two different techniques. The first involves applying displacements to the surface nodes of the breast mesh in the compression direction (Ruiter et al 2006, Hipwell et al 2007). The compression plates and their interaction with the breast are not explicitly modelled in this case, but it is assumed that the displacements applied on the surface nodes have the same effect. It was reported by Ruiter (2003) that this technique leads to artifacts on the breast surface (breast swelling) around the nodes that are adjacent to the ones that the displacements are applied to. An alternative method is to explicitly model the interaction between the contact plates and the breast tissue using either a friction (Alonzo-Proulx et al 2010, Hopp et al 2013a, Lee et al 2013) or a frictionless model (Pathmanathan et al 2004, Chung et al 2008b, Han et al 2012, Hopp et al 2012b, Mertzanidou et al 2014).

The behaviour of the breast tissue under compression is influenced by the presence of the pectoral muscle. This was modelled in the literature either by constraining the nodes close to the pectoral muscle to be fixed (Samani et al 2001, Chung et al 2008b, Hopp et al 2012b, 2013a), or by allowing them to slide along the chest wall (Alonzo-Proulx et al 2010, Tanner et al 2011, Han et al 2012, Lee et al 2013, Mertzanidou et al 2014, Solves-Llorens et al 2014). Although the surface of the pectoral muscle can be extracted from the MRI, for simplicity this has sometimes been neglected (Hopp et al 2013a) or has been approximated as a plane (Mertzanidou et al 2014, Solves-Llorens et al 2014).

4.4.2. MRI to x-ray mammography registration.

A patient-specific FE modelling approach was first used as part of an MRI to x-ray registration framework by Ruiter et al (2006) and Ruiter (2003). This implementation used the breast outline for alignment rather than, for instance, the intensities in the two images. The registration was performed in two stages: in the first step a plate compression was applied and in the second a breast outline alignment was achieved by applying additional displacements on the surface nodes of the breast model. Hopp et al (2013a) extended this approach by introducing an additional step in which the rotation of the breast about the anterior-posterior axis was estimated using an intensity-based optimisation. Similarly, another FEM-based approach with a contact model was proposed (Lee et al 2013) which also employed an iterative intensity-based registration framework. However, this was limited to a 2D rigid-body registration of the simulated projection to the x-ray mammogram. Similarly, Solves-Llorens et al (2014) followed an approach that also used a single breast compression of the MRI, and then employed a 2D registration between the simulated projection image and the real mammogram. Finally, Mertzanidou et al (2014) proposed an intensity-based registration, where both pose and FEM parameters were optimised in the same framework, and multiple FEM-based simulations were tested at each iteration of the optimisation.

Regarding validation, quantitative results on clinical cases showed initially a mean error of 4.3 mm on six cases (Ruiter et al 2006). A more recent semi-automated implementation of the same approach had a mean error of 11.8 mm, on CC-view mammograms only, of eleven patients (Hopp et al 2012a). Finally, in the latest automated implementation (Hopp et al 2013a), where the results were improved by optimising the roll angle using intensity-based registration, the mean error was 13.2 mm on CC-view only mammograms of a substantial 78 patients. The FEM-based technique followed by a 2D rigid intensity-based registration of Lee et al (2013) was tested on both CC- and MLO-view mammograms of six patients and gave a mean error of 19 mm. The method followed by Solves-Llorens et al (2014) reported errors of 4.5 mm in fourteen datasets. It is worth noting that for this method a two-dimensional, non-rigid transformation was employed to warp the simulated projection image to match the real mammogram, which may introduce physically non-realistic deformations. Finally, a simultaneous optimisation of pose and FEM parameters (Mertzanidou et al 2014) produced a mean error of 11.3 mm when tested using CC- and MLO-view mammograms of ten patients. Of note is that this value is only marginally lower than the value of 12 mm obtained using the same data registered via a volume-preserving affine transformation. The latter was previously validated on 113 CC and MLO mammograms and a median error of 13 mm obtained (Mertzanidou et al 2012a).

Direct comparison between these methodologies is problematic, due to the different datasets used and variations in the features used for validation. Of note is the additional challenge associated with registering MLO- versus CC-view mammograms. Registration of the former is often more challenging, due to greater uncertainties regarding the positioning of the breast before compression and also the effect of the pectoral muscle on the compression simulation.

Given the 3D-to-2D nature of this registration task, an accuracy approaching that of the MR voxel resolution of one to two millimetres may be unrealistic, however the above results suggest that the combined benefits of a physically realistic biomechanical transformation, robust optimisation framework, appropriate cost function and accurate MRI to x-ray simulation have yet to be fully realised for this application.

4.4.3. Ultrasound computed tomography to x-ray mammography registration.

Hopp et al (2012a) and Hopp et al (2013a) applied a similar approach to that used for MRI/x-ray mammography alignment to Ultrasound Computed Tomography (USCT) to mammography registration (Hopp et al 2010, 2012b, Hopp et al 2013b). Initially, using a homogeneous model and manual estimation of the rotation in 2D of the projected volume, a mean Euclidean distance between corresponding lesion centers in 2D of 7.3 mm was reported using four data sets (Hopp et al 2010). In a subsequent study, estimation of this 2D rotation was automated and a number of similarity measures compared for this purpose (Hopp et al 2012b). Gradient correlation was found to reduce the error in a data set of thirteen subjects from 15.8 to 10.4 mm but this configuration was not applied to an independent test set. In their most recent work (Hopp et al 2013b), Young's modulus was computed from the USCT speed of sound distribution across the breast and each element of the biomechanical model assigned the value for its corresponding region using four squared or exponential relationships between the speed of sound and Young's modulus. Compared to the mean error of 12.2 mm obtained for five data sets using a homogeneous model, the error using each of the heterogeneous models varied between 11.8 and 13.1 mm.

5. Conclusion and outlook

Biomechanically informed breast image registration methods have been developed for a number of applications, namely:

registration of prone MRI (Roose et al 2008),
registration of prone MRI to PET (Coman et al 2004, Unlu et al 2005, Krol et al 2006) or PET/CT (Unlu et al 2010, Dmitriev et al 2013),
registration of prone MRI to supine MRI (Carter et al 2008, Lee et al 2010, Eiben et al 2013, Han et al 2014),
registration of prone MRI to supine surgical surface/depth data (Carter et al 2008, Lago et al 2012),
MRI to x-ray registration (Ruiter 2003, Ruiter et al 2006, Mertzanidou et al 2011, Hopp et al 2013a, Lee et al 2013, Mertzanidou et al 2014, Solves-Llorens et al 2014), and
registration of ultrasound CT and x-ray mammography (Hopp et al 2010, 2012b, Hopp et al 2013b).

Biomechanical modelling has also been used to validate more conventional breast image registration methods (see for example the papers of Tanner et al (2007) and Hill et al (2009)). Despite these developments there remain a number of issues that limit clinical applications of biomechanical modelling, e.g. the accuracy of biomechanical modelling, failure to meet clinically acceptable levels of computational cost and time, and the complexity of applying patient-specific biomechanical modelling methods to routine clinical practice. These still need to be addressed.

In many cases the perceived benefit of ensuring 'physically plausible' deformations offered by biomechanical models is not bourne out by the reported registration accuracy. This suggests that the state-of-the-art in mathematical modelling techniques adopted by current studies lack the sophistication required to capture the biomechanical complexities of highly deformable breast tissue. A case in point is MRI to x-ray registration (section 4.4) in which a number of groups have developed a range of methods, none of which has been able to demonstrate an accuracy of less than 10 mm on a sizeable dataset. The quality of the biomechanical modelling is undoubtedly at least partially responsible for this disappointing performance, however the extent to which other factors, such as the optimisation strategy, cost function and boundary constraints etc are also complicit, is as yet unclear. The required accuracy of these methods is obviously application dependent but a suitable goal would be significantly greater than the size of the lesion of interest (or margin for surgical applications), and ideally approaching the resolution of the specific imaging modality involved.

Obtaining accurate material properties is an on-going research topic in tissue engineering and constitutive modelling. More exact properties could be provided, however, by in vivo measurements offered by emerging techniques such as magnetic resonance elastography (section 3.1). Additionally, boundary conditions on the chest wall are typically specified from a limited set of extreme, and hence crude, constraints, i.e. fixed boundary conditions or free-sliding interfaces (e.g. section 4.4.1). Interaction of the breast and skin surface with mammographic compression plates is an outstanding modelling issue and the introduction of tilting paddles has implications for patient comfort and benefit which should be investigated.

The computational cost of obtaining numerical solutions is high. Therefore extremely scalable algorithmic implementations are required to ensure completion of the computations in useable timeframes. Whilst progress has been made in developing fast solvers (Taylor et al 2008, Johnsen et al 2015), such explicit dynamic methods are not universally appropriate for all applications. This becomes increasingly critical when high-fidelity breast model discretisation is required, for instance in real-time utilisation of surgical simulators and in iterative registration methods.

Finally, the construction of individual-specific biomechanical models is challenging to automate. Intervention from an experienced technical user is typically required to ensure, for instance, the correctness and integrity of a discretised breast geometry. This issue needs to be addressed to enable computational methods to be introduced into routine clinical practice (section 3.2).

**Figure 2.** Unloaded (or gravity-free) configuration, numerical prediction of a subject-specific breast geometry, using an inverse problem analysis approach (Vavourakis *et al* 2015). (a) Arrow vectors illustrating the displacement magnitude distribution predictions of the inverse analysis. (b) Superposition of the prone (transparent orange) and the gravity-free (blue surface grid) breast geometry.
Download figure:
Standard image High-resolution image

**Figure 3.** FE models for breast plate compression deformation analysis (reproduced with permission from Han *et al* (2012)). (a) Before compression. (b) After compression.
Download figure:
Standard image High-resolution image

Acknowledgments

This work was funded by the European 7th Framework Program grants VPH-PRISM (FP7-ICT-2011-9, 601040), VPH-PICTURE (FP7-ICT-2011-9, 600948), a Marie-Curie Fellowship (FP7-PEOPLE-2013-IEF, 627025), and the Engineering and Physical Sciences Research Council grant MIMIC (EP/K020439/1).

A review of biomechanically informed breast image registration

Article metrics

Author e-mails

Author affiliations

Dates

Abstract

1. Introduction