International Journal of Biomedical Imaging

Volume 2017 (2017), Article ID 7232751, 12 pages

https://doi.org/10.1155/2017/7232751

## Multitemporal Volume Registration for the Analysis of Rheumatoid Arthritis Evolution in the Wrist

DITEN, Università degli Studi di Genova, Via Opera Pia 11a, 16145 Genova, Italy

Correspondence should be addressed to Silvana G. Dellepiane

Received 21 November 2016; Revised 9 May 2017; Accepted 12 June 2017; Published 19 October 2017

Academic Editor: Jyh-Cheng Chen

Copyright © 2017 Roberta Ferretti and Silvana G. Dellepiane. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

This paper describes a method based on an automatic segmentation process to coregister carpal bones of the same patient imaged at different time points. A rigid registration was chosen to avoid artificial bone deformations and to allow finding eventual differences in the bone shape due to erosion, disease regression, or other eventual pathological signs. The actual registration step is performed on the basis of principal inertial axes of each carpal bone volume, as estimated from the inertia matrix. In contrast to already published approaches, the proposed method suggests splitting the 3D rotation into successive rotations about one axis at a time (the so-called basic or elemental rotations). In such a way, singularity and ambiguity drawbacks affecting other classical methods, for instance, the Euler angles method, are addressed. The proposed method was quantitatively evaluated using a set of real magnetic resonance imaging (MRI) sequences acquired at two different times from healthy wrists and by choosing a direct volumetric comparison as a cost function. Both the segmentation and registration steps are not based on a priori models, and they are therefore able to obtain good results even in pathological cases, as proven by the visual evaluation of actual pathological cases.

#### 1. Introduction

Diseases related to the wrist, such as rheumatoid arthritis (RA), have a major negative impact on quality of life because they are the most common causes of severe long-term pain and physical disability [1]. An early diagnosis of the pathology in such cases provides the patient with targeted care, a high probability of improvement, and a reduction in social costs [2]. The focus of this work was the erosive arthritis of wrist bones, which is the most disabling rheumatic disease [3].

Even though conventional radiography is considered the traditional gold standard for evaluating bone erosions, poor performances (i.e., overall sensitivity, specificity, and accuracy of 25.5%, 98.3%, and 70.1%, resp.) have been reported for detecting wrist erosions [4]. As clearly stated in [5], computed tomography (CT) has been shown to be more sensitive than radiography but, limited by the radiation dose, it is rarely used in clinical practice. Indeed, MRI allows the simultaneous assessment of synovial membranes, articular fluid, cartilage, bones, ligaments, tendons, and tendon sheaths [6]. MRI can detect erosive RA changes such as bone and cartilage damage with greater sensitivity than conventional radiography [7]. An interesting and clear overview of the clinical problem is given in [5] along with indications for the quantification of RA features from MRI.

For the diagnosis and monitoring of disease progression, a longitudinal study is necessary where the physician conducts several observations of the same subject over a period of time. A major objective is then the definition of a simple but effective procedure for volume registration while preserving the original intensity values. To this end, a simplified method is here proposed whose main objective is the achievement of a limited effect of interpolation-related blurring. The registration result allows the detection of changes in the morphological bone shape, the localization of changes with respect to the bone volume and surface, and the direct grey-level comparison after registration.

Repeatability of the wrist tomographic examination is affected by difficulties in perfectly reproducing the location and orientation of the anatomical district. Even if translation or angular rotation differences are expected to be relatively small, changes more strongly affect the coronal plane than the sagittal and axial planes. With the aim of a direct shape comparison, a 3D registration is then required to align two volumes acquired at different times so that they overlap and a quantitative analysis of changes becomes possible.

In this context, this work proposes a segmentation-based registration method allowing the overlapping of two MRI tomographic spatial sequences so that the shape of a bone at two different instant times can be compared. Such a unimodal 3D-3D multitemporal registration makes it possible to monitor disease progression in the rheumatic wrist. A rigid registration was chosen to avoid artificial bone deformations and to allow finding eventual differences in the bone shape due to erosion, disease regression, or other eventual pathological signs.

The automatic segmentation of the bone of interest from both acquired volumes is the starting point for a good estimate of the location and orientation parameters, which is here performed through the computation of the inertia matrix. The subsequent steps of registration are performed on the basis of principal inertial axes of each volume.

In contrast to more classical approaches, the proposed method suggests splitting the 3D rotation into successive basic rotations about one axis at a time (also termed elemental rotations). In such a way, only a single plane is affected at one time (i.e., coronal, sagittal, and axial plane), starting from the most important one, as indicated by the major principal inertia axis. In addition, singularity and ambiguity drawbacks affecting other classical methods (e.g., the Euler angles method) are addressed.

Depending on the carpal bones, stopping at one basic rotation is sometimes better than performing two rotations or the full 3D transform. In fact, the errors introduced by the interpolation process can be bound to one plane at a time and can be avoided in those planes where the error is larger than the expected registration benefit.

Such a result has been achieved by applying and quantitatively evaluating the proposed method to a set of real tomographic sequences from healthy subjects and by choosing as a cost function a direct volumetric comparison of a wrist bone acquired at two different times. Starting from such empirical results, the suggestion for the best rotation steps for each bone is derived and can be applied also in the pathological situation. The method can also be applied to pathological bones because the segmentation process is based on the tomographic acquisition grey-level values and does not make any use of atlas or anatomical models. Furthermore, bone morphological changes do not affect the registration procedure. The direct application of the method to such pathological cases has been compared with the classical Euler method, showing that singularity and ambiguity drawbacks do not affect the proposed solution.

Due to the good registration result, 3D visualization of the registered bone surfaces could be displayed to assess eventual morphological changes related to pathology evolution, or multiplanar reformatted slice cuts of the two registered volumes could be displayed, along with the original grey levels, to visually compare corresponding slices of a given bone.

Some aspects led to better results for a few bones when applying only one elemental rotation instead of all three (or even the classical rotation matrix):(i)Because the first rotation refers to the eigenvector associated with the largest eigenvalue, its effect is the most relevant, it is associated with the most informative content, and it aligns the two bones along their first principal axis.(ii)The second and the third rotations refer to the second and the third principal axes, which are more critical and might be affected by a larger estimation error.(iii)Even though repeatability cannot be assured in two different tomographic acquisitions, the flat shape of the hand and wrist assures that differences in the position and orientation likely affect the coronal sections more than the axial or sagittal ones.(iv)The rotation angle might be so small that the errors introduced by the transformation and interpolation steps are larger than the benefit of the rotation itself.The paper is structured in the following way. In Section 2, the work is placed in the context of the state-of-the-art methods. In Section 3, the various registration steps and the metrics used for a quantitative evaluation are presented, pointing out that the method is automatic, unsupervised, and not based on models or a priori knowledge, yet it adapts to the image content. In Section 4, the results obtained from the application of the method to MRI volumes of carpal bones are shown, and these findings are discussed in Section 5.

#### 2. Previous Works: Image and Volume Registration

The two major registration strategies can be classified into intensity-based and feature-based methods. An intensity-based approach compares pixel/voxel intensity by means of correlation metrics. Feature-based methods find correspondence between significant features such as points, lines, and surfaces, and a geometrical transformation is then determined to map the two volumes [8].

Longitudinal volumetric image analysis is an important topic, especially investigated for brain structure studies. The review paper by Mills and Tamnes [9] gives a comprehensive overview on longitudinal structural imaging for brain development in children and adolescents.

Indeed, from the imaging point of view, brain morphology is very complex and registration turns to be a very crucial task where intensity-based procedures are mainly required. As an example, the research carried out in [10] addresses the analysis of brain atrophy or, more in general, changes in brain size and shape. It is based on an automated linear registration tool, FLIRT (FMRIB’s Linear Image Registration Tool) [11], which makes use of the correlation ratio cost function and of a multiscale optimization strategy.

In [12], the quantification of pointwise changes in surface morphology of the bones of the human wrist is addressed based on CT imaging. The proposed method, referred to as Registration-based Bone Morphometry (RBM), consists of two steps: an atlas selection step and an atlas warping step. Statistical group analysis is reported to demonstrate the application of RBM for the comparison of representative carpal bones based on sex and for the comparison of carpal bones of the left and right wrists for individuals. The utility of RBM is also demonstrated in the context of tracking bone erosion status in RA.

The notable brain size and the structured grey-level content assure a reliable statistical information basis for such intensity-based approaches. In contrast, given their unique shape, hand bones are better suited for applications of feature- and geometrical-based methods. In addition, their small size, the MRI low field, and the low spatial resolution make it harder to use intensity-based algorithms which cannot be based on enough statistically significant data.

As previously reported [13], 2D and 3D feature-based image registration approaches can be grouped according to the four basic steps of the registration procedure:* feature detection, feature matching, mapping function design,* and* image transformation and resampling, *which will be recalled in the next paragraphs.

##### 2.1. Feature-Based Registration Methods

According to previous reports [13], the use of feature-based methods is recommended if the images contain enough distinctive and easily detectable objects. In medical applications, such a problem is usually faced by interactive selections made by an expert or by introducing extrinsic features, rigidly positioned with respect to the patient (skin markers, screw markers, dental adapters, etc.). Recently, the availability of reliable segmentation results has enabled improvement of such a process: the obtained closed-boundary regions can be used to derive surface or volumetric features, such as in [14], where the cervical spine is segmented and then a 3D articulated registration is applied. In [15], a registration of functional data was performed after the segmentation of the left ventricle.

After the identification of the significant features, their correspondence can be estimated using their descriptors, preferably invariant to the expected image deformation. A group of methods devoted to register satellite images uses moment-based invariants [16] for a description of the segmented regions features. Flusser and Suk used the affine transform invariants [17, 18]; Holm [19] extracted closed-boundary regions and proposed to represent them by perimeter, area, compactness, moments, and moment invariants; and Brivio et al. [20] modeled shadows in mountain images by means of their inertia ellipses. For the stereo-correspondence problem, Bhattacharya and Sinha [21] suggested the application of complex moments.

The mapping function design step, subsequent to the feature correspondence step, is devoted to designing a transform algorithm so that the corresponding features of the two volumes would be as close as possible. At this level, the type of mapping function is chosen, and its parameters are estimated.

Volume registration might be a global or a local process. In the present application, we were interested in coregistration of single bones to monitor the eventual disease evolution. To this end, even though the spatial configuration of the various bones strongly depends on the position of the anatomical part, this problem is solved by using a local registration, where the attention is focused on single bones (volumes of interest) separately. Because each bone is a rigid body, affine transformation is sufficient to achieve a good result and allows for a final morphological comparison.

Finally, the mapping functions constructed at the previous step are applied during the transform step so that the volumes are actually registered.

In the next paragraph, a short overview of the most widely used methods is reported, referring to mapping function design and parametric transformation models addressing the rigid registration. They focus on extraction and correspondence analysis of location and orientation descriptors of the objects of interest and propose a few analytical formulations for the registration solution.

##### 2.2. Mapping Functions and Parametric Transformation Models

In space science, the “attitude” of a flying object refers to its orientation with respect to a fixed inertial frame. In robotics, the estimation of orientation and location parameters is a key aspect of rigid-body kinematics. The problem addressed in the present paper has some similarities with these concepts, as the same bone is acquired in a different position and orientation by means of two (nonregistered) tomographic acquisitions performed at some time distance, so that the imaged bone volumes have a different orientation in 3D space.

A very detailed and mathematically sound review of related concepts is given in [22], which was used as a basis for the following considerations.

Let** T** denote an inertial fixed frame whose coordinates are (), and let denote the body-attached frame that is free to rotate. Frame is determined by its orthogonal tern .

As defined in [22], the orientation of a rigid body is defined as the orientation of frame with respect to frame . The attitude of a rigid body can be expressed by a variety of mathematical parameterizations, which can be either constrained with redundant elements or unconstrained with minimal elements. The rotation matrix and the unit quaternion are examples of constrained parameterizations. Euler angles, the Rodrigues parameters, and the modified Rodrigues parameters (MRPs) are examples of minimal parameterizations [22].

In practice, the most common orientation representations are the rotation matrix, Euler angles, and unit quaternion. These constitute the basis of most attitude and orientation estimation techniques, and their properties and group structures are of great importance [22].

The Euler angle parameterization is not global, and angles are not well-defined for *θ* = ±*π*/2. This singularity problem is not unique for the Euler angle parameterization but affects all other 3D parameterizations. Unit quaternion is a 4D parameterization, which allows such singularities to be avoided [22].

Both Euler angles and unit quaternions are proven to have a corresponding rotation matrix. Using Euler’s theorem, which implies that anybody orientation can be specified in terms of a rotation by some angle about an arbitrary axis, unit quaternions can assume a formulation that gives the corresponding rotation matrix through the Rodrigues formula.

The rotation matrix is a global and unique representation of orientation. By using this matrix, corresponding voxels in coordinates of inertial () and body () frames can be mapped to each other as the result of rotations applied to the 3D Euclidean space.

Similar to Euler angles, the rotation matrix can be obtained as a product of three different rotation matrices, each corresponding to an elemental rotation about one of the axes of the fixed coordinate system or about the rotating coordinate system. However, as it was experienced in the present work, for the given application, it is sometimes better to avoid one or even two rotations.

##### 2.3. Medical Applications of Moment-Based Registration Methods

In medical applications, registration methods based on principal axes have been applied in different cases of study such as in [23, 24] for MR and positron emission tomography (PET) brain images and in [25] for computed tomography (CT) and PET images. In [26], this method is used for the registration of CT and transmission single-photon emission computed tomography (SPECT) thorax images and for myocardial SPECT stress and rest scans.

Referring to the present application, several papers exist in the literature about the registration of the whole wrist. Most of these studies correlate the registration with the wrist kinematic [27–29], and they consider both the possible rotations of the whole wrist and the movement limits between the singular bones. When focalizing on the wrist bones, such as in [30], a matching between the object of study and a model previously created under different viewing angles is proposed.

In [31], the movements of the wrist are studied, through the registration of the scaphoid bone. In particular, 10 wrist scans are performed from one subject in 10 different positions, and the orientation of the scaphoid, through the determination of the principal axes, is evaluated. The comparison between the different scaphoid bone positions is carried out by registering the volumes with Euler’s angles method.

Unlike in [31], our method is not based on the Euler angles and a single rotation matrix, but it performs successive basic rotations to achieve a better registration, as will be demonstrated by our experiments. The relationship between the rotation matrix and the three elemental (i.e., basic) rotation matrices is given in Subsection in Supplementary Material available online at https://doi.org/10.1155/2017/7232751.

#### 3. The Proposed Approach

The proposed method is not intended to register the whole wrist where kinematic studies are required, and neither does it address the matching with some idealized model or atlas. On the contrary, the purpose of our work is the registration of individual bones, acquired at different times, to detect eventual changes and support the presence of eventual erosion and evolution of the rheumatic disease in general.

The method is supported by a robust, automatic, and reliable segmentation [32, 33] to apply a local registration only to the bone of interest. The volumetric region referring to the bone of interest is extracted independently from the first and second MRI acquisitions.

Based on absolute and central moments of the obtained 3D bone volumes, the proposed method takes the principal axes as feature descriptors and computes them starting from the inertial matrix. The transformation model approach proposes to align both bone volumes to the inertial frame to give them the same location and orientation in 3D space and to allow the subsequent shape comparison. To this purpose, instead of the direct application of one roto-translation matrix, the procedure proposes successive rotations on each 2D plane (along with appropriate translations). The optimization problem is then related to the possible registration parameter configurations, which describe the execution order of the 2D rotation axis (in the maximum of three rotations), the number of such rotations, and the direction of each rotation. Such geometrical considerations, described at Section 3.3, assure that coherent rotations are applied to the two corresponding bone volumes to have them aligned in orientation and in direction as well, without ambiguity.

As a feature-based approach, the estimate of the transformation parameters and the possible configurations are computed in a deterministic way. The cost function adopted is the volumetric comparison between the two aligned bone volumes, as computed by the confusion matrix, which will be described later. Such quantitative metrics have been applied to the available database of real non-pathological cases to drive the selection of the best parameter configuration for each wrist bone type, as explained in the following section.

The fixed inertial frame with digital coordinates () is given. Let us start from two MRI spatial sequences, acquired in the coronal plane from the same individual at two instant times and . From the tomographic sequences, two digital volumes, and , are generated. Typically, the first (initial) acquisition and the second (follow-up) acquisition are separated by at least a few months in time. These two volumes constitute a mapping whose domains are the voxels coordinate () and codomain is a grey-level value ranging from 0 to 255. For the specific tomograph we used (ESAOTE Artoscan 0.2 Tesla), the image size in the coronal section (represented by [] coordinates) is 256 × 256 pixels. In the case of anisotropic resolution, an interpolation between coronal slices is performed. The coordinate refers to the normal to the coronal plane and ranges from 0 to , where is the number of coronal slices in the tomographic volume, after the eventual interpolation step. As an example, a coronal slice of a wrist is shown in Figure 1(a), together with an axial and a sagittal section in Figures 1(b) and 1(c), respectively.