Registration of Volumetric Prostate Scans using Curvature Flow

Registration of Volumetric Prostate Scans using Curvature Flow


Radiological imaging of the prostate is becoming more popular among researchers and clinicians in searching for diseases, primarily cancer. Scans might be acquired with different equipment or at different times for prognosis monitoring, with patient movement between scans, resulting in multiple datasets that need to be registered. For these cases, we introduce a method for volumetric registration using curvature flow. Multiple prostate datasets are mapped to canonical solid spheres, which are in turn aligned and registered through the use of identified landmarks on or within the gland. Theoretical proof and experimental results show that our method produces homeomorphisms with feature constraints. We provide thorough validation of our method by registering prostate scans of the same patient in different orientations, from different days and using different modes of MRI. Our method also provides the foundation for a general group-wise registration using a standard reference, defined on the complex plane, for any input. In the present context, this can be used for registering as many scans as needed for a single patient or different patients on the basis of age, weight or even malignant and non-malignant attributes to study the differences in general population. Though we present this technique with a specific application to the prostate, it is generally applicable for volumetric registration problems.

Shape registration, geometry-based techniques, medical visualization, mathematical foundations for visualization

000 \vgtccategoryResearch \vgtcinsertpkg \authorfooter Saad Nadeem, Rui Shi, Joseph Marino, Xianfeng Gu, and Arie Kaufman are with Stony Brook University, Stony Brook, NY 11794-4400.
E-mail: {sanadeem, ruishi, jmarino, gu, ari} Wei Zeng is with Florida International University.
Email: \CCScatlist \CCScatComputer GraphicsI.3.3Data RegistrationGeometry-based Techniques, Medical Visualization, Mathematical Foundations for Visualization \teaser Registration of and axial MR scans for the same patient on the same day. MR scans are taken from (a) and (b) axial view and segmented (c & d). These segments are reconstructed into a surface (e & f) and mapped onto respective spheres and tetrahedralized (g & h). These spheres are then registered using surface feature points onto another sphere (i). On this registered sphere, we can mark regions (j) and since our mapping is diffeomorphic, these regions are mapped back to the original tetrahedralized prostate shapes (k & l). \nocopyrightspace
\firstsection Introduction

Cancer of the prostate is the second leading cause of cancer-related mortality among males in the United States, and is the most commonly diagnosed cancer in general [18]. Although it is such a common cancer, diagnosis methods remain primitive and inexact. Detection relies primarily on the use of a simple blood test to check the level of prostate specific antigen (PSA) and on the digital rectal examination (DRE). If an elevated PSA level is found, or if a physical abnormality is felt by the physician during a DRE, then biopsies will be performed. Though guided by transrectal ultrasound (TRUS), these biopsies are inexact, and large numbers (143 in one case) are often necessary to try and retrieve a sample from a cancerous area [15]. More recently, it has been noted that Magnetic Resonance Imaging (MRI) can be used for the detection of prostate cancer [10].

Multiple MR images obtained with different settings are necessary for the detection of prostate cancer. Most commonly used is a combination of and image sequences. images are generally of higher quality than images, and prostate cancer will appear as areas of reduced intensity. The acquisition of these MR image sequences is often done with varying orientations and resolutions per sequence. In cases where the image sequences are acquired during a single session without the patient moving, the resulting volumes will be naturally registered in world space. However, there are many scenarios, where registration methods are necessary for multiple sequences. One such reason would be the situation where a patient shifts during a session, such that the patient’s body position is not constant, and thus the volumes will not orient themselves properly with respect to each other in 3D space. Registration is also required when the sequences are acquired at different times. This could be a common occurrence to monitor the progression of the tumor(s) during watchful waiting or brachytherapy. Thirdly, it is possible that one might want to correlate the MR images to histology information in order to confirm results during development and testing of a system.

In this paper, we proposed a novel volumetric registration framework, which can reduce the work of radiologists, and help to register the datasets, to monitor the progression of abnormalities, and to facilitate the diagnoses of prostate diseases. The pipeline consists of reconstruction, mapping, and registration. We first reconstruct the volume of the prostate from or MR images. Then we utilize our volumetric mapping algorithm to map the prostate to a solid ball. Finally we register the solid balls obtained from the mapping, resulting in a registration of the original prostates.

Our volumetric registration framework is based on a volumetric parameterization algorithm using discrete volumetric Ricci flow, by which a ball-shape volume is mapped to a canonical domain, i.e., a unit ball. This parameterization technique is a homeomorphism, which means one-to-one and continuous. This property is important to medical imaging and visualization because we do not want to lose the fidelity of the data, or miss any part which could possibly be an abnormality, or a sign of the diseases. This mapping also has a low stretching energy in a physical sense, so that the shape distortion is limited for small regions after the mapping. After mapping two prostates to the unit ball, we use anatomic feature points to further align and register them. We have thoroughly tested our registration method on data from ACRIN 6659 study, Prostate MR Image Database [1] and Stony Brook UHMC Database for scans acquired in different orientations, different modes of MRI and different days. The results show that we can easily handle all these cases and give good registration.

The main contributions of this work are as follows:

  1. We propose a novel volumetric registration framework

  2. We show registration of prostate MR scans in different orientations, different modes and different days (for prognosis monitoring)

The rest of the paper is organized as follows. In Section 2 we review the previous work on volume registration. Then we introduce the theoretical foundation of our algorithm in Section 3, followed by the details of our framework in Section 4. Evaluation methods and experimental results are shown in Section 5 and 6, respectively. Finally, we conclude our discussion and sketch the future research plan in Section 7.

1 Related Works

There have been some related works on medical volumetric registration, but most of them are based, which cannot ensure the accuracy, or even the one-to-one property. In particular, registration of prostate volumes from MRI has often relied on the use of a similarity matrix and the application of a deformation field to the original MR image slices [4]. More recently, it has been suggested to perform the nonrigid registration by extracting tetrahedral grids from the volume datasets and deforming these grids by using the finite element method [6]. We also use the tetrahedral grid structure, but perform the deformation through the use of the volumetric curvature flow algorithm.

Ricci flow [9] is a powerful tool for geometric analysis. It has been applied to prove the Poincaré conjecture successfully. Surface Ricci flow presents an efficient approach to computing the uniformization theorem [8]. Ricci flow refers to the process of deforming Riemannian metric proportional to the curvature, such that the curvature evolves according to a heat diffusion process, eventually the curvature becomes constant everywhere [5]. Discrete surface Ricci flow [11, 22] generalizes the curvature flow method from smooth surface to discrete triangular meshes. The key insight to discrete Ricci flow is based on the following observation: conformal mappings transform infinitesimal circle fields to infinitesimal circle fields. Furthermore, Ricci flow can be used to construct the unique conformal Riemannian metric with prescribed Gaussian curvatures. The resulting conformal surface mappings are free of angle distortion.

Volumetric mapping is a generalization of the surface mapping problem. In theory, three dimensional manifolds (volumetric) also admit uniformization metrics, which induce constant sectional curvature. The only existing approach for generalizing of surface Ricci curvature flow to volume was presented by Yin et al. [23]. It generalizes the discrete curvature flow from surfaces to hyperbolic 3-manifolds with complete geodesic boundaries. The metric deforms according to the curvature, until the curvature is constant everywhere. However, it can only support the hyperbolic 3-manifolds and cannot handle the cases of topological ball. In this paper, we present the discrete volumetric Ricci flow to compute the volumetric mappings for topological balls, such as the volumetric prostates. To the best of our knowledge, the proposed method is the first work to generalize the discrete surface curvature flow to 3-manifolds with solid ball topology.

The uniformization metrics on 3-manifolds have great potential for volumetric parameterization, volumetric shape matching and registration, and volumetric shape analysis. In this work, we propose to use this discrete volumetric curvature flow for registration of topological ball volumes with large deformations. Compared to the previous methods, our proposed method can generate diffeomorphisms (one-to-one and onto mappings). It is global and steady, and can handle large deformations.

There are some related works on volumetric mapping. The volumetric harmonic map depends on volumetric Laplacian; Zhou et al. [24] have applied volumetric graph Laplacian to large mesh deformation. Harmonicity in volumes can be similarly defined via the vanishing Laplacian, which governs the smoothness of the mapping function. Wang et al. [21] have studied the formula of harmonic energy defined on tetrahedral meshes and computed the discrete volumetric harmonic maps by a variational procedure. Volumetric parameterization using fundamental solution method [14] was applied to volumetric deformation and morphing. Other than that, harmonic volumetric parameterization for cylinder volumes is applied for constructing tri-variate spline fitting  [17]. All the above approaches rely on volumetric harmonic maps. Unfortunately, these volumetric harmonic maps cannot guarantee bijective mappings even though the target domain is convex. Besides the volumetric harmonic map, there is another approach of mean value coordinates for closed triangular meshes [12, 7]. Mean value coordinates are a powerful and flexible tool to define a map between two volumes. However, there is no guarantee that the computed map is a diffeomorphism.

Our approach differs intrinsically from these existing approaches in two ways. First, we solve the curvature flow for arbitrary topological balls (prostate volumes are of this case). As a result, the curvature flow induced map is guaranteed to be a diffeomorphism. Second, we use the theory of Ricci flow rather than the conventional volumetric harmonic map [21].

2 Theory

This section briefly introduces the theoretic foundation of our discrete curvature flow method.

2.1 Surface Ricci Flow

Let be a smooth surface with a Riemannian metric . One can always find isothermal coordinates for locally, which satisfies


where is a conformal factor. The Gaussian curvature of the surface is given by

where is the Laplace-Beltrami operator induced by . Although the Gaussian curvature is intrinsic to the Riemannian metric, the total Gaussian curvature is a topological invariant: the total Gaussian curvature of a closed metric surface is

where is the Euler number of the surface. Suppose and are two Riemannian metrics on the smooth surface . If there is a differential function , such that

then the two metrics are conformal equivalent. Let the Gaussian curvatures of and be and respectively. Then they satisfy the following Yamabe equation


Suppose the metric in local coordinate. Hamilton introduced the Ricci flow as


During the flow, the Gaussian curvature will evolve according to a heat diffusion process.

Theorem 1 (Hamilton and Chow)

[5] Suppose is a closed surface with a Riemannian metric. If the total area is preserved, the surface Ricci flow will converge to a Riemannian metric of constant Gaussian curvature.

2.2 General Ricci Flow

This section uses the Einstein summation convention. Given coordinate functions , the tangent vector can be described by its components in the basis . The Riemannian metric tensor is given by

Let vector fields and . The Riemannian metric defines an inner product

The so-called Christoffel symbols are given by


where is the inverse matrix of .

The covariant derivative of a basis vector along a basis vector is a vector, and is given by

We get

The curvature tensor of the manifold is given by

Curvature tensor measures non-commutativity of the covariant derivative. The sectional curvature is given by


where is the wedge product. Intuitively, vector and span a tangent plane in the tangent space at a point . All the geodesics emanating from and tangent to the plane form a surface . The Gaussian curvature of at is .

Assume is an orthonormal basis of the tangent space at , the scalar curvature is the trace of the curvature tensor,

Ricci curvature is a linear operator on tangent space at a point,

For general Riemannian manifold, Ricci flow is given by


The evolution behavior of Ricci flow on 3-manifolds is much more complicated, because singularities will emerge in the flow. We refer readers to [5] for detailed explanation.

2.3 Discrete Surface Ricci Flow

Discrete surface Ricci flow generalizes the curvature flow method from smooth surface to discrete triangular meshes. The key insight to discrete Ricci flow is based on the following observation: conformal mappings transform infinitesimal circle fields to infinitesimal circle fields. Discrete Ricci flow replaces infinitesimal circles by circles with finite radii, and modifies the circle radii to deform the discrete metric, to achieve the desired curvature.

In engineering fields, surfaces are approximated by their triangulations, the triangle meshes. A triangle mesh is a dimensional simplical complex embedded in .

The discrete Riemannian metric is defined as the set of edge lengths,

as long as, for each face , the edge lengths satisfy the triangle inequality: . The edge lengths determine the corner angles by cosine law. The discrete Gaussian curvature at a vertex can be computed as the angle deficit,


where represents the corner angle attached to vertex in the face , and represents the boundary of the mesh.

The Gauss-Bonnet theorem holds on meshes as follows.


where M is a mesh.

In the smooth case, Conformal deformation of a Riemannian metric is defined as

In the discrete case, there are many ways to define conformal metric deformation. Generally, we associate each vertex with a circle centered at with radius .

On an edge , two circles are separated, the edge length is given by,


where is the so-called inversive distance. During the conformal deformation, the radii of circles can be modified, but the inversive distances are preserved. There exists a unique circle, the so called radial circle, that is orthogonal to three vertex circles. The radial circle center is denoted as (see Figure 1(a)).

(a) (b)
Figure 1: Configuration of discrete Ricci flow (a) Surface case, (b) Volume case. (a) Inversive circle packing. (b) Discrete edge Ricci Curvature.
(a) (b)
Figure 2: Geometric interpretation of discrete surface Ricci energy. The inversive distance circle packing on each face (a) truncates a hyperbolic tetrahedron (b). The Ricci energy is the total hyperbolic volumes of the truncated tetrahedra.

Let be the logarithm of , the discrete Ricci flow is defined as follows:


where is the user defined target curvature and is the curvature induced by the current metric. The discrete Ricci flow has exactly the same form as the smooth Ricci flow, which conformally deforms the discrete metric according to the Gaussian curvature.

The discrete Ricci flow can be formulated in the variational setting, namely, it is a negative gradient flow of a special energy form, the so-called entropy energy. The energy is given by


where is an arbitrary initial metric. Figure 2 shows the geometric interpretation of the discrete Ricci energy. Each triangle with a circle pattern corresponds to a truncated hyperbolic tetrahedron. As shown in frame (a), the triangle is laid on the xy-plane, a prism is built based on the triangle. Through each circle, one hemisphere is constructed. The intersection between the inside of the prism and the outside of all hemispheres is a truncated tetrahedron. The upper half space is treated as three dimensional hyperbolic space with the Riemannian metric

The truncated tetrahedron is treated as a hyperbolic tetrahedron, with edge lengths . The dihedral angles on edges are (see Figure 2(b)). Denote the hyperbolic volume of the tetrahedron as , then according to Schlafli formula:


Then the Ricci energy is the Legendre dual to the volume,


therefore . The Ricci energy for the whole mesh is the summation of those of all faces,


Computing the desired metric with user-defined curvature is equivalent to minimizing the discrete entropy energy. Because the discrete entropy energy is locally convex, it has global minima.

The Hessian matrices for discrete entropy are positive definite. The energy can be optimized using Newton’s method. The Hessian matrix can be computed using the following formula. Given a triangle , let be the distance from the radial circle center to the edge (see Figure 1(a)), then


We define the edge weight for edge , which is adjacent to and as


The Hessian matrix is given by the discrete Laplace form


Algorithmic details for discrete surface Ricci flow can be found in [11] and [22].

2.4 Discrete Volumetric Ricci Flow

Volumetric data are represented as tetrahedral meshes. Suppose is a tetrahedron, the dihedral angle on edge is denoted as (see Figure 1(b)). If the edge is an interior edge, (i.e. is not on the boundary surface), the discrete Ricci curvature is defined as


If is on the boundary surface, its curvature is given by


The discrete curvature flow for volume is defined by


where is the edge length of .

Similar to the discrete surface Ricci flow, volumetric curvature flow is also variational, which is the gradient flow of the following energy:


where is the volume of the tetrahedron. Therefore, the total energy for the whole tetrahedron mesh is the summation of those of each tetrahedron,


In the procedure of the curvature flow, the image of each tetrahedra is non-degenerated. Therefore, the mapping is a homeomorphism.

Unlike surface case, the volumetric Ricci energy is not convex. The const curvature metric corresponds to a critical point of the total energy, which is a saddle point. Therefore, in the computational process, the choice of initial Riemannian metric is important. In practice, we use homotopy method explained as follows.

In our current project, we deform the prostate volume to the unit solid ball . The initial metric on is denoted as . The boundary surface of is denoted as . First we use discrete surface Ricci flow method to map to the unit sphere . This gives us the final metric on the boundary surface , denoted as . Then we define a one parameter family of Riemannian metrics on as


For each , we fix the metric on the boundary, and then use volumetric curvature flow to compute a metric, such that all interior edge curvatures are zeros. Then we use as the initial metric, to compute . This homotopy method greatly improves the stability and robustness of volumetric curvature flow.

3 Algorithm

Our volume registration pipeline consists of three steps: reconstruction, parameterization and registration. First the tetrahedral mesh is reconstructed after we segment the prostate from the scans. Then each scan is parameterized to a canonical unit ball respectively. After that, the resulting balls are registered using marked feature points, which could give us a registration of the prostates. Here we explain in details the algorithms in each step of our pipeline.

3.1 Volume Reconstruction

The segmentation of the prostate from the surrounding anatomy is the first step. Though there has been research into the automatic segmentation of the prostate from axial MR images [25], it is not general to all possible types of images. Since our focus is on the registration of the prostate, segmentation is outside the scope of this paper, and therefore, we apply manual segmentation to the datasets. After we extract the prostate from the raw image slices, we apply the Marching Cubes algorithm [16] to build the corresponding boundary surface of the prostate, where tricubic interpolation [13] is used for resampling in case of possible different resolutions in the three dimensions.

The triangle mesh which we get from the previous step acts as the boundary constraint of the volume tetrahedralization. Then we use Tetgen [20] to generate a tetrahedral mesh for a given surface mesh to meet the boundary constraints. To improve the meshing quality, we employ the variational tetrahedral meshing technique [3] which can significantly reduce the slivers and produce well-shaped tetrahedral meshes, as shown in Figure 3(a).

3.2 Volumetric Parameterization

Here we assume the volume is represented as a tetrahedral mesh as follows. is the vertex set; is the edge set; is the face set; and is the tetrahedron set. Suppose is a tetrahedron with vertices , and is the edge connecting vertices and .

Surface Parameterization

We first use our surface curvature flow method (Algorithm 3.2) to calculate a parameterization from the prostate boundary surface to the unit sphere.


[!] Surface parameterization of a sphere. \SetKwInOutRequireRequire \RequireThe prostate volume

  1. Extract the boundary surface of the prostate volume .

  2. Remove one triangle from , isometrically map the triangle onto the plane, denote the image triangle on the plane as

  3. Using discrete surface Ricci flow introduced in Section 2.3 to compute a conformal map from by setting the target curvature

  4. Fill in the removed face . Scale the image of .

  5. Use stereo-graphic projection to map to the unit sphere. The north pole is inside the image of .

Metric Homotopy Method for Volumetric Curvature Flow

After we get the boundary surface mapping, we can compute the discrete volume Ricci flow using the surface mapping as boundary constraint. In order to improve the stability and the robustness of the volumetric curvature flow, we use metric homotopy method to map the the prostate volume to the solid unit ball, the mapping is denoted as .

(a) (b)
Figure 3: Examples of the tetrahedral mesh structure for (a) a prostate dataset and (b) the resulting ball from the volumetric curvature flow. The mesh is cut through the center to show the interior structure.

The boundary mesh has the initial induced Euclidean metric and the final metric on the unit sphere. The discrete metrics are represented as edge lengths. For each edge , we use to denote its initial length (metric), and its final length on the sphere. Then we define the metric at time , as

At each time , we fix the boundary edge lengths and use Ricci flow to adjust the interior edge lengths, such that all interior edge curvatures become zeros. For each interior edge , we deform its length as

After we compute the flat metric for the whole volume at time , we step further to compute that at time , using the resulting metric at as the initial input, and boundary mesh metric as the boundary condition. During the evolvement, if we encounter a degenerated tetrahedron, we perform local re-triangulation. See Algorithm 3 for details.


[!] \SetKwInOutRequireRequire \SetKwInOutEnsureEnsure Volume parameterization of a ball. \RequireThe prostate volume , and boundary mapping

  1. Set , the boundary metric to be the initial Euclidean

  2. For all interior edges , compute the edge curvature.

  3. Evolve the interior edge length .

  4. If the tetrahedron is close to be degenerated, locally remesh
    the neighborhood by swapping the shared face between two adjacent tetrahedra.

  5. Repeat step 2,3, until all the interior edge curvature are zeros.

  6. Update the boundary metric. For each edge on the boundary, .

  7. Update time .

  8. Repeat step 2 through 5, until .

Figure 3(b) shows the structure of the resulting tetrahedral mesh.

3.3 Volumetric Registration

After mapping different prostates to the unit solid ball, we register them using anatomical feature points.

Prostate Feature Detection

The prostate, a gland like a walnut in size and shape, does not contain a complicated geometric structure. The prostate gland, which surrounds the urethra, is located in front of the rectum, and just below the bladder.

For volumetric registration, we need to match at least three identical anatomical features within the MRI images of different directions to obtain the accurate and reliable registration result. A pair of glands called the seminal vesicles are tucked between the rectum and the bladder, and attached to the prostate as shown in Figure 4(a). The urethra goes through prostate and joins with two seminal vesicles at the ejaculatory ducts. Therefore, some distinctive anatomical structures, such as the prostatic capsule and seminal vesicle contours, ejaculatory ducts, urethra, and dilated glands as represented in Figure 4(c), can be applied for the registration between different scan directions of one dataset or between MR slices and histology maps.

(a) (b)
Figure 4: Anatomical prostate points. (a) The anatomical position of the prostate [2]. (b) Multi-view of a prostate MR image along three viewing directions. (c) The prostate structures [2]. Each feature point with the pre-defined index number is highlighted by the red circle. S_No’s are the surface feature points while {I_No} is the interior feature point set.

MRI provides images with excellent anatomical detail and soft-tissue contrast, and MRI sequences are displayed in a serial order. We analyze -weighted datasets along the axial, sagittal and coronal view, as shown in Figure 4(b). On each MRI prostate view direction, the exact outline of prostate boundary is traced and all corresponding feature points are manually marked with predefined index numbers as shown in Figure 5. For registration, we use three exterior feature points and a set of interior ones based on the structure information of urethra and seminal vesicles: two surface feature points are exactly the entrance and exit points of the urethra because the urethra goes through the entire prostate, while the third surface point is marked at the intersection between each seminal vesicle and prostate with respect to the fact that two seminal vesicles attach to the surface of prostate and merge with urethra at the ejaculatory ducts. A set of interior feature points are marked along the urethra, beginning with the entrance surface feature point and ending with the exit one.

(a) (b)
Figure 5: Manually marked features on the MRI slices using (a) axial, (b) coronal and (c) sagittal view of the same prostate dataset. Red contour highlights the boundary of prostate and green points show the corresponding features. The urethra is marked using yellow contour.

Registration Framework

Let and denote the volumes of the two prostates. The computational algorithm for registration is as shown in the following diagram.

We first compute two volumetric maps , and . Then, we compute an automorphism of , , which aligns all the feature points. Then the final registration is given by

The key component is to construct . Suppose and are the corresponding feature points on and . Then should align them .

Figure 6: Registration of and coronal view MR scans for the same patient at (a) day 102 and (b) day 1801. MR scans are taken in coronal modes and segmented. These segments are then reconstructed into a surface and mapped onto respective spheres and tetrahedralized. These spheres are then registered using surface feature points on another sphere (c). On this registered sphere, we mark regions (d) and since our mapping is diffeomorphic, these regions are mapped back to the original tetrahedralized prostate shapes (e & f).

First we choose three feature points on the boundaries (see e.g. Figure 6(a)). Assume and are the three points on the two volumes. We use stereo-graphic projection to map the boundary of the ball onto the complex plane. Then, we construct a Möbius transformation to align these three feature points. By abusing the symbols, we still use for their images of the stereo-graphic projection on the complex plane. Define

then maps to . Similarly, we define which maps to . Then aligns there feature points on the complex plane, align the feature points on the sphere.

In order to align interior feature points, we add position and target curvature constraints to volumetric curvature flow. Then in the parameter ball, the interior feature points will be placed at exactly the target position. To align the feature points, we first map model A to a solid sphere with no constrains, and get the result coordinates of the features; then map model B to the solid sphere with target feature position the same as these feature coordinates in the solid sphere. By doing this, we can perfectly align these feature points, which leads to a good alignment of the models. Furthermore, our method does not create flipped tetrahedrons as in the case of harmonic mapping (as explained in the following section), so the registration result is good, both locally and globally.

4 Evaluation

(a) (b)
Figure 7: The comparison of (a) our ricci flow method with (b) harmonic mapping. The edges are flipped at the boundary in harmonic mapping, making in unsuitable for our application.

We test the algorithm using and -weighted MR postate images. In it is not possible to identify interior feature points, so wherever we consider for registration, we make use of only the surface feature points. In , we can identify interior feature points so we try to identify as many of these as possible on the corresponding prostates during the registration. We also observe that the urethra is easily identified in the coronal view (in high quality MR Images) so to take advantage of the feature points along urethra, whenever possible, we use coronal views. We also compare our method against harmonic maps.

Harmonic map method has been broadly applied for surface and volume registration in medical imaging. One of the major merits for surface harmonic map is that it produces diffeomorphism, as stated in Rado’s theorem [19]: for a harmonic map from a surface with the disk topology to a convex planar domain, if the restriction of the map on the boundary is a homeomorphism, then the interior harmonic map is a diffeomorphism. Unfortunately, this theorem doesn’t hold for volumetric case. As shown in Fig. 7, the volumetric harmonic map introduces flipped tetrehedra (the orientation of the image tetrhedra is reversed) near the boundary area. In contrast, our current method guarantees the mapping to be one-to-one and flipping-free.

Figure 8: Registration of -weighted MR scans of the same patient in different orientations. (a) Axial view, (b) Coronal view and (c) Sagittal view. These are registered and tetrahedralized (d). A region is marked on the registered sphere (e) and the corresponding regions on the registered prostate models are identified (f,g & h).

5 Results

The data we have used for this work is from ACRIN 6659 study, the Prostate MR Image Database [1] and Stony Brook UHMC Database. Our work focuses on the registration, so for this project we have used manual segmentation (done by an expert/radiologist) to obtain the prostate gland. We will make the data from Stony Brook UHMC Database publicly available as well.

The Prostate MR Image Database contains MR prostate scans for 231 patients with multiple days and in different modes. The Prostate MR Image Database also contains expert segmentations for most of the datasets, but there were instances where we were not able to find the segmentation corresponding to a specific dataset. For these instances, we employed the help of a radiologist to do the segmentation. The Stony Brook UHMC Database contains MR scans for more than 200 patients. In this case, there were no prior segmentations available so a radiologist segmented out the datasets tested.

MRI scans Vertices Tetrahedra Edges Faces Running time
a 96536 564478 678062 1146005 170 secs
b 98898 580013 695976 1177092 191 secs
c 115946 689982 822853 1396890 210 secs
Table 1: Running time of the volumetric curvature flow algorithm for the three prostate MR scans, corresponding to the models reconstructed in Figure 8(a), (b) & (c).

This work was performed on an Intel Xeon E2540 2.5GHz machine with 16GB of RAM. The running time to map the mesh from its original shape to the sphere is around 3 minutes for the 3 MRI scans for Figure 8 (see Table 1).

Feature# Reg Error Reg Error
w/ seg distortion
1 0.0035 0.0119
2 0.0111 0.0148
3 0.0096 0.0108
4 0.0089 0.0091
5 0.0053 0.0078
6 0.0032 0.0051
7 0.0196 0.0201
8 0.0187 0.0221
9 0.0117 0.0175
10 0.0112 0.0129
11 0.0148 0.0174
12 0.0106 0.0181
13 0.0179 0.0209
Table 2: Registration error computed with respect to the Euclidean distances between the registered feature points on different prostates (Figure 9), then divided by the diagonal length of the model’s bounding box (in this case, the length is 100). The Feature# are the labels shown in Figure 9(a). The second column shows the error distances computed between feature points in Figure 9(a) and (b), whereas the third column shows the distances between Figure 9(a) and (c).

Different Days. We have tested our method with different modes of MRI, i.e. and -weighted scans. For prognosis monitoring of a tumor, we register coronal views from multiple days (as shown in Figure 6) since the urethra is easily visible in this view, as discussed in Section 4 above. Then we can identify the reduced intensity regions, if found, and mark the corresponding region on the registered sphere. Since our mapping is diffeomorphic, we can get corresponding regions on the registered prostate models, as shown in Figure 6. We see that the suspected region is larger at Day 1801 than at Day 102. This allows us to measure the shape of the tumor, if marked accurately. Furthermore, we can register as many scans as possible, preferably -weighted in coronal view, to get a continuous progression of the abnormality. In case of Figure 6, we have registered these prostate models on the basis of four feature points, 3 surface feature points and one interior feature point (marked at the intersection of three lines of urethra). We then mark out the other feature points along urethra to study the error.

Figure 9: (a) The feature points marked out for the urethra on the prostate after registration in Figure 6(a). (b) The corresponding urethra feature points on the prostate (b) in Figure 6(b) after registration. (c) The corresponding urethra feature points after a segment was removed from the original MR scans.

For the alignment of these feature points along urethra, we actually create a one-to-one map between prostate models A (Figure 6(a)) and B (Figure 6(a)) by mapping these two models to the same parameter sphere with the original 4 features aligned. Hence, we map the urethra feature points in model B to model A by first mapping them to the parameter sphere and then to model A. As shown in Figure 9, the original urethra feature points in model A are well aligned with that from model B, which gives a strong evidence of our good registration quality. We further provide the evidence of the quality of our registration using the euclidean distances computed between these points as shown in Table 2. Since we are using manual segmentations, we also provide the error if the segment boundary is distorted and a complete segment slice is removed. The results for this combined distortion and removal of a segment slice are shown as the column of Table 2 as well. As seen, the Euclidean distances does not change drastically giving the strong foundations for our registration.

Different Modes of MRI. Moreover, we also register different modes of MR scans for a patient on the same day. The result of this registration is shown in Figure Registration of Volumetric Prostate Scans using Curvature Flow. We register and scans in axial view on the basis of 3 surface feature points, since interior feature points can not be identified in . We observe that if the interior feature points are marked on the registered sphere from scans, we can get a good estimate for the location of these interior feature points on the registered prostate model, reconstructed from scans (given the good quality of our registration regardless of the minor distortions in the initial segmentation). To prove this is the case, we registered 3 -weighted MR scans in axial view from the same day for a patient on the basis of just the surface feature points and then computed the euclidean distances (as done above) to study the error if a feature point is marked at the intersection of the three lines of urethra on one of the prostate models and the rest are registered with reference to this. The error computed did not exceed in this case (again the diagonal length of the bounding box was 100 as in Figure 9), showing strong evidence that given good segmentation, the feature points from one prostate model can be accurately marked out on the other registered models.

Different Orientations. We also register MR scans acquired in different orientations, as shown in Figure 8. The axial, coronal, and sagittal images are all segmented, reconstructed, and parameterized. The resulting models are then registered using the feature points stated in Section 3.3. Table 3 shows the error registration for three extra geometric identical points used to evaluate the quality of registration. We use similar method as discussed above for error evaluation by keeping one model as a reference and mapping the identical points from the other to the reference and compute the euclidean distances.

From the results, we observe that the error for feature points in sagittal and coronal views is less than the axial view and the other views. This is due to the error of segmentation and feature points location. In the coronal or sagittal direction, all the feature points are within one or two slices, and we are able to precisely distinguish their locations. However, in the axial direction, features are scattered in multiple slices. Since the interslice resolution is much lower than the intraslice one, the euclidean distance error for the extra geometrical feature points is larger in the axial and sagittal/coronal combination than the coronal and sagittal combination. However, still, our results show good quality of our registration even in this scenario.

Feature# Reg Error Reg Error Reg Error
(a) & (b) (a) & (c) (b) & (c)
1 0.0025 0.0029 0.0014
2 0.0021 0.0024 0.0018
3 0.0022 0.0027 0.0019
Table 3: Registration error computed with respect to the Euclidean distances between the registered feature points on different prostates (Figure 8), then divided by the diagonal length of the model’s bounding box (in this case, the length is 100 again). The feature# are the labels of the identical geometric feature points to evaluate the method. The second column shows the error distances computed between feature points in Figure 8(a) Axial and (b) Coronal, the third column shows the distances between Figure 8(a) Axial and (c) Coronal and the fourth column shows the error between feature points in Figure 8(b) Coronal and (c) Sagittal

6 Conclusion

In this work, we performed a registration of volumetric meshes of the prostate gland using volumetric curvature flow. We have tested our method by registering the MR scans of the same patient in different orientations, from different days and with different modes of MRI. Our method gives good results on data from the ACRIN 6659 study, Prostate MR Image Database and Stony Brook UHMC Database. Moreover, analysis of the resulting registration shows accurate demarcation of the same region can be achieved across different registered prostates (from different scans), even with distortions in the initial segmentation.

In the future, we plan to integrate our registration method with automatic/semi-automatic segmentation to create a complete framework for the prostate and eventually, extend this framework for other volumetric objects, such as the brain. We also plan to validate our method across different modalities, such as CT, MRI and ultrasound.

The data we have use for this work is from ACRIN 6659 study, the Prostate MR Image Database [1] and Stony Brook UHMC Database. This work has been supported by NSF grants IIS0916235, CCF-0702699 and CNS0959979 and NIH grant R01EB7530.


  1. Prostate MR Image Database:
  3. P. Alliez, D. Cohen-Steiner, M. Yvinec, and M. Desbrun. Variational tetrahedral meshing. ACM Trans. Graph., 24(3):617–625, 2005.
  4. A. Bharatha, M. Hirose, N. Hata, S. Warfield, M. Ferrant, K. Zou, E. Suarez-Santana, J. Ruiz-Alzola, A. D’Amico, R. Cormack, R. Kikinis, F. Jolesz, and C. Tempany. Evaluation of three-dimensional finite element-based deformable registration of pre- and intra-operative prostate imaging. Medical Physics, 28:2551–2560, Dec. 2001.
  5. B. Chow, P. Lu, and L. Ni. Hamilton’s Ricci Flow. American Mathematical Society, 2006.
  6. A. D. d’Aische, M. De Craene, S. Haker, N. I. Weisenfeld, C. Tempany, B. Macq, and S. K. Warfield. Improved non-rigid registration of prostate MRI. Proc. of MICCAI, pages 845–852, 2004.
  7. M. S. Floater, G. Kós, and M. Reimers. Mean value coordinates in 3d. Comput. Aided Geom. Des., 22(7):623–631, 2005.
  8. O. Forster. Lectures on Riemann Surfaces, volume 81 of Graduate texts in mathematics. Springer, 1991.
  9. R. S. Hamilton. The ricci flow on surfaces. Mathematics and general relativity, 71:237–262, 1988.
  10. M. R. Hersh, E. L. Knapp, and J. Choi. Newer imaging modalities to assess tumor in the prostate. Cancer Control, 11(6):353–357, Nov. 2004.
  11. M. Jin, J. Kim, F. Luo, and X. Gu. Discrete surface Ricci flow. IEEE TVCG, 14(5):1030–1043, 2008.
  12. T. Ju, S. Schaefer, and J. Warren. Mean value coordinates for closed triangular meshes. In SIGGRAPH ’05: ACM SIGGRAPH 2005 Papers, pages 561–566, New York, NY, USA, 2005. ACM.
  13. F. Lekien and J. Marsden. Tricubic interpolation in three dimensions. Journal of Numerical Methods and Engineering, 63:455–471, 2005.
  14. X. Li, X. Guo, H. Wang, Y. He, X. Gu, and H. Qin. Harmonic volumetric mapping for solid modeling applications. In SPM ’07: Proceedings of the 2007 ACM symposium on Solid and physical modeling, pages 109–120, New York, NY, USA, 2007. ACM.
  15. T. Loch. Urologic imaging for localized prostate cancer in 2007. World Journal of Urology, 25:121–129, 2007.
  16. W. E. Lorensen and H. E. Cline. Marching cubes: a high resolution 3d surface construction algorithm. Computer Graphics, 21(4):163–169, 1987.
  17. T. Martin, E. Cohen, and R. M. Kirby. Volumetric parameterization and trivariate b-spline fitting using harmonic functions. Comput. Aided Geom. Des., 26(6):648–664, 2009.
  18. L. A. G. Ries, D. Melbert, M. Krapcho, D. G. Stinchcomb, N. Howlader, M. J. Horner, A. Mariotto, B. A. Miller, E. J. Feuer, S. F. Altekruse, D. R. Leqis, L. Clegg, M. P. Eisner, M. Reichman, and B. K. Edwards (eds). SEER cancer statistics review, 1975-2005, national cancer institute. Bethesda, MD, based on November 2007 SEER data submission, posted to the SEER web site, 2008.
  19. R. Schoen and S.-T. Yau. Lectures on harmonic maps. International Press Cambridge, 1997.
  20. H. Si. Tetgen: A quality tetrahedral mesh generator and three dimensional delaunay triangulator.
  21. Y. Wang, X. Gu, P. M. Thompson, and S.-T. Yau. 3d harmonic mapping and tetrahedral meshing of brain imaging data. Medical Image Computing and Computer Assisted Intervention, 2, 2004.
  22. Y.-L. Yang, R. Guo, F. Luo, S.-M. Hu, and X. Gu. Generalized discrete ricci flow. Computer Graphics Forum, 28(7):2005–2014.
  23. X. Yin, M. Jin, F. Luo, and X. Gu. Discrete curvature flow for hyperbolic 3-manifolds with complete geodesic boundaries. In International Symposium on Visual Computing (ISVC2008), Las Vegas, Nevada, USA, December 2008.
  24. K. Zhou, J. Huang, J. Snyder, X. Liu, H. Bao, B. Guo, and H.-Y. Shum. Large mesh deformation using the volumetric graph laplacian. ACM Trans. Graph., 24(3):496–503, 2005.
  25. R. Zwiggelaar, Y. Zhu, and S. Williams. Semi-automatic segmentation of the prostate. Proc. of IbPRIA, pages 1108–1116, 2003.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description