Performance Study of Nonrigid Registration Algorithm for
Investigating Lung Disease on Clusters∗
Fumihiko Ino1
1

Yuya Tanaka1

Hiroko Kitaoka2

Kenichi Hagihara1

Graduate School of Information Science and Technology, Osaka University
1-3 Machikaneyama, Toyonaka, Osaka 560-8531, Japan
2
Graduate School of Medicine, Osaka University
2-2 Yamadaoka, Suita, Osaka 565-0871, Japan
E-mail: ino@ist.osaka-u.ac.jp
Abstract

X-ray computed tomography (CT) scans [8]. Therefore,
some acceleration techniques are needed to utilize registration during surgery, where registration time is strictly restricted within ten minutes to minimize the patient’s strain.
There are many papers [2, 3, 8, 10–12] reporting experiences in using parallel machines to make registration timeshort enough to perform intraoperative registration. Theseprior projects successfully accelerate the registration procedure by reducing computational complexity and by parallelizing performance bottlenecks in the procedure. Although they present that combining some acceleration techniques realizes intraoperative registration, it is still not clear
how efﬁciently each technique reduces registration time.
In this paper, we study the performance of a nonrigid
registration algorithm, aiming at making clear what kindsof acceleration techniques and computing environments are
effective to reduce registration time. To do this, we implement a registration algorithm for investigating lung disease on clusters. Our algorithm is based on an optimization
approach [9], which resolves the registration problem into
an optimization problem. We combine this approach withtwo conventional acceleration techniques: a data-parallel
processing technique for accelerating the registration procedure; and a precomputation technique for reducing thecomputational complexity.
The remainder of the paper is as follows. Section 2
presents a brief overview of the registration algorithm withthe two acceleration techniques mentioned above. Section
3 shows the performance study on three clusters of PCs. Finally, Section 4 summarizes the paper.

This paper presents a performance study of a nonrigid
registration algorithm for investigating lung disease onclusters. Our algorithm combines two conventional acceleration techniques in order to achieve fast registration: adata-parallel processing technique for accelerating the registration procedure; and a precomputation technique for reducing the computational complexity. We perform some experiments on three clusters with different CPU and network
performance in order to make clear what kinds of acceleration techniques and computing environments provide higher
performance. The results show that a cluster with Gigabit
Ethernet (GbE) network is the most cost effective solutionthat reduces registration time from ten hours to ten minuteswith a linear speedup.
Keywords: image registration, medical image processing, cluster computing, MPI, performance evaluation.

1. Introduction
Image registration [1] is a technique for establishingpoint-to-point correspondences between different images
taken at different time, from different viewing points, in different modalities. This technique assists medical doctors invarious diagnoses. For example, it helps them in detecting
cancers by monitoring changes in size, shape, or image intensity over time intervals.
One problem with this technique is that it requires timeconsuming computations in order to realize accurate, robust, and completely automated alignments. For example,
it takes several hours to align three-dimensional (3-D) clinical images obtained by modern imaging systems such as

2. Nonrigid Registration
Let F and R be the ﬂoating and reference images, respectively. The nonrigid registration problem is deﬁned ascomputing a nonrigid transformation T that aligns F to R.

∗ This work was partly supported by JSPS Grant-in-Aid for ScientiﬁcResearch on Priority Areas (17032007).

1

Proceedings of the Sixth International Conference on Parallel and Distributed Computing, Applications and Technologies (PDCAT’05)
0-7695-2405-2/05 $20.00 © 2005 IEEE

(a)

(b)

Figure 2. Hierarchical B-spline free-form deformations (FFDs) [9]. (a) Deformations of theﬂoating image are performed by manipulating an overlaying mesh of control points in acoarse-to-ﬁne manner. δ represents the spatial resolution of control points and φi,j represents a control point. (b) The deformation ofpoint (x, y) is determined by its surrounding
4 × 4 neighborhood of control points.

Figure 1. Optimization-based registration algorithm [5]. It minimizes the cost function C
by the steepest descent optimization [7] in acoarse-to-ﬁne manner. The ﬁnal mesh Φ determines the ﬁnal transformation T.

Figure 2 shows how this interpolation is applied to image
F . The nonrigid transformation T is given by manipulatinga mesh of control points overlayed in the image domain Ω ={(x, y, z) | 0 ≤ x < X, 0 ≤ y < Y, 0 ≤ z < Z}. Givena mesh Φ of control points φi,j,k , the transformation T ofpoint (x, y, z) in image F is deﬁned by

To compute this, our algorithm uses an optimization approach [5], which aligns images through the optimizationof a cost function C associated with a similarity measure
between the images (see Figure 1). This optimization is performed by means of the steepest descent optimization [7] ina coarse-to-ﬁne manner. This hierarchical alignment strategy contributes to reduce the computational complexity.
Our algorithm has four advantages as follows:

T(x, y, z) =
3
3 
3 


Bl (u)Bm (v)Bn (w)φi+l,j+m,k+n ,where u = x/δ − x/δ, v = y/δ − y/δ, w = z/δ −z/δ, i = x/δ−1, j = y/δ−1, k = z/δ−1, and Bl
represents the l-th basis function of cubic B-splines. Thisequation means that the deformation of any point (x, y, z)
is given by interpolating the deformations of its surrounding
4 × 4 × 4 neighborhood of control points (see Figure 2(b)).
Recall here that our algorithm aligns images in a coarseto-ﬁne manner. To do this, the algorithm organizes meshΦ, images F and R in a hierarchy. Then, it decreases thespatial resolutions of images and control points, γ and δ,
respectively, at each level of hierarchy.

• Hierarchical, locally controlled deformations [9] by Bspline functions [4];• Robust similarity measure by information theory [5];• Fast B-spline interpolation by precomputation [8];• Fast optimization by parallelization [11].

2.1. Hierarchical Deformation Model

2.2. Cost Function

Our algorithm employs B-spline functions in order torepresent the nonrigid transformation T at lower complexity. B-spline functions are computationally efﬁcient because they use interpolation to realize hierarchical, locallycontrolled free-form deformations (FFDs).

The cost function for optimization is deﬁned by Rueckert
et al. [9] as follows:
C(Φ) = −Csimilarity (R, T(F )) + αCsmooth (T),2

Proceedings of the Sixth International Conference on Parallel and Distributed Computing, Applications and Technologies (PDCAT’05)
0-7695-2405-2/05 $20.00 © 2005 IEEE(1)

l=0 m=0 n=0(2)

where Csimilarity (R, T(F )) is the similarity measure between the reference image R and the transformed ﬂoating
image T(F ), Csmooth is the penalty term that constrains
the spline-based transformation to be smooth, and α is theweighting parameter which is determined experimentally.
For the similarity measure Csimilarity , our algorithm employs a general measure based on information theory. This
information-based strategy does not depend on speciﬁc features of images, for example, contours and surfaces, sothat it realizes robust registration of multimodality images[5, 8, 9]. Although there are many measures such as mutual information [5] and its normalized version [9], we haveexperimentally determined to use correlation coefﬁcient forlung registration. In summary, the similarity measure between images A and B is given by
Csimilarity (A, B) =
(A(x, y, z) − Ā)(B(x, y, z) − B̄),
(A(x, y, z) − Ā)2 (B(x, y, z) − B̄)2

Figure 3. Precomputation [8] for B-spline interpolation. Gray colored voxels have thesame value φ̂i+l , and thereby this term should
be computed only once.

nique [8] to the algorithm. The key idea of this technique isprecomputation that intends to eliminate redundant computations which are repeatedly for some different points.
Suppose that we have a zyx-loop to compute T(x, y, z)
for all points (x, y, z) in the image domain Ω. Eq. (1)3
can be rewritten as T(x, y, z) =l=0 Bl (u)φ̂i+l , where3
3
φ̂i+l =
m=0
n=0 Bm (v)Bn (w)φi+l,j+m,k+n . Then,term φ̂i+l keeps the same value for all points (x, y, z) in
one row located within the same cell of the mesh (see Figure3). Therefore, placing this computation outside the x-loop
reduces the total computational cost for the zyx-loop.
This modiﬁcation allows us to compute the basic Bspline function only 4n3 +42 n2 +43 n times while the unoptimized implementation computes it 192n3 times, where nrepresents the number of points in the image. Therefore, theexecution time for object deformations will be signiﬁcantly
reduced by 98%.

(3)

where A(x, y, z) is the image intensity at point (x, y, z) inimage A, and Ā is the mean intensity in image A.
The penalty term Csmooth intends to regularize the transformation in order to avoid unsmooth deformations. For
lung registration, we use the following constraint model:
Csmooth =    2 2  2 2  2 2
∂ T
∂ T
∂ T
1 X Y Z
+
+
2
2
V 0
∂x
∂y
∂z 2
0
0
 2 2 2 2 2 2
∂ T
∂ T
∂ T
+2
+2
+2
∂xy
∂xz
∂yz


 2
2
λ ∂ T ∂2T ∂2T+++
dxdydz,
(4)
2 ∂x2
∂y 2
∂z 2

2.4. Parallelization

where V denotes the volume of the image domain Ω, andλ represents the elastic coefﬁcient for lung. From the viewpoint of physics, this term approximates the energy of anisotropic elastic object (namely lung tissue) which is subjected to deformations with volume change.
As we mentioned before, our algorithm employs thesteepest descent optimization [7] to ﬁnd the optimal transformation parameter Φ that minimizes Eq. (2). The algorithm stops this optimization if a local optimum has beenfound. Here, it assumes a local optimum if ||∇C|| ≤ ,
where ∇C = ∂C/∂Φ, and ||∇C|| and  represent the gradient norm of the cost function C and a threshold for minimization, respectively. The gradient ∇C is estimated byusing the ﬁnite-difference approximation [7].

Although the precomputation technique contributes toreduce registration time, parallelization is still required tomake it be compatible with surgical usage.
Because nonrigid registration is not a data-intensive application, we assume that all computing nodes have the entire images in their local memory. Actually, a 512 × 512 ×
512 voxel image consisting of 2-bytes requires 256MB ofmemory, so that this assumption is reasonable for modern
PCs equipped with more than 512MB of main memory.
We ﬁrst analyzed our sequential implementation to locate performance bottlenecks:• Gradient computation. This computation phase is required for the steepest descent optimization of the costfunction. It takes approximately 57% of total registration time.

2.3. Precomputation for B-spline Interpolation
In order to reduce the computational complexity of Bspline FFDs, we applied Rohlﬁng’s precomputation tech-• Similarity computation. Computing the cost function
3

Proceedings of the Sixth International Conference on Parallel and Distributed Computing, Applications and Technologies (PDCAT’05)
0-7695-2405-2/05 $20.00 © 2005 IEEE

(a)

(b)

(c)

Figure 4. Checkerboard visualization of reference and ﬂoating images (a) before and (b) after registration. The registered image in (b) has the same contours as in (c) the reference image.

C is also a performance bottleneck. It takes 42% of thetotal time.

Level
1
2
3
4

• Image sampling. This phase occurs when the algorithm increases the hierarchy of alignment. Although
this phase is a small bottleneck on single CPU machines, it should be parallelized on parallel machines,
because I/O operations usually become a larger bottleneck as the number of CPUs increases.

Table 1. Dataset speciﬁcation.
Image size
Number of control points
33 × 33 × 36
252
65 × 65 × 70
1,848
129 × 129 × 138
13,041257 × 257 × 275
97,785

We performed nonrigid registration using a pair of lung
CT images (see Figure 4). Table 1 shows its speciﬁcation
with the hierarchical organization. The ﬁle size of an image
is 157MB.
Table 2 summarizes our clusters. Each cluster has twointerconnection networks: one is specialized for data communication between computing nodes; and the other, theEthernet-based network, is basically for I/O from/to the centralized storage server. Note here that cluster #1 is a uniprocessor (UP) system while the remaining clusters are 2way symmetric multiprocessor (SMP) systems.

Each of the above computation phases has parallelism in
the image space. Therefore, we can parallelize each phaseby decomposing the image domain into small portions andassigning them to CPUs. After processing these portions
in a data-parallel manner, each CPU has a local computation result, so that a collective communication is needed
to reduce local results into a global result. Thus, the entire workload can be decomposed into many independentpieces.
For this workload distribution, our parallel algorithm
employs a cyclic distribution, because we assume no data
decomposition. This cyclic strategy statically balances
workloads among CPUs. In other words, static load balancing is sufﬁcient for our algorithm whose workload is generated uniformly in the image domain.

3.1. Acceleration Results by Precomputation
Table 3 shows the timing results measured on a single
node in the clusters. The precomputation technique reduces
registration time by 21–39%. These results are reliable, because the execution time for object deformations, which isincluded in registration time, is reduced by 83–86% while
our theoretical analysis predicts a reduction of at most 98%.
Thus, measured reduction is close to predicted reduction.
These timing results also indicate that sufﬁcient tunings
are essential to parallelize sequential implementations in aproper manner. For example, if the precomputation technique was not applied to the sequential implementation, we
might have parallelized the code for B-spline interpolation.

3. Performance Study
To make clear the timing beneﬁts of the acceleration
techniques mentioned above, we now study the performance of the registration algorithm using three Linux clusters with different performance characteristics in terms of
CPU and network. Our registration algorithm is implemented using the C++ language and the Message Passing
Interface (MPI) standard [6].
4Proceedings of the Sixth International Conference on Parallel and Distributed Computing, Applications and Technologies (PDCAT’05)
0-7695-2405-2/05 $20.00 © 2005 IEEE

Table 2. Cluster speciﬁcation. Network bandwidth and latency are measured using a performancemeasurement tool (http://www.scl.ameslab.gov/netpipe/).
Component
CPU
Cache size
Number of nodes
Number of CPUs per node
Network #1
Bandwidth
Latency
Network #2
Bandwidth
Latency

Cluster #1
Xeon 3.2GHz
512KB
32
1
InﬁniBand (InB)
3999Mb/s
10µs
GbE
852Mb/s
22µs

Cluster #2
Xeon 2.8GHz
512KB
16
2
Myrinet-PCIX (MYX)
1848Mb/s
6µs
GbE
897Mb/s
17µs

Cluster #3
Pentium III 1GHz
256KB
64
2
Myrinet-2k (MY)1433Mb/s
9µs
Fast Ethernet (FE)89Mb/s
48µs

smaller part of the image domain, as compared to thesequential implementation. Therefore, as the numberof CPUs increases, we have higher cache utilization.
Note here that, at the same time, CPUs generally spend
longer time for communication. Thus, the superlinear speedup appears only when the program achieves
higher cache utilization with less communication.

Table 3. Precomputation results on 1 CPU.
Cluster Original Optimized Reduction
time (s)time (s)rate (%)#1
18,352
11,20939#221,796
13,53638#359,424
46,99721• When using more than 32 nodes with Fast Ethernet(FE) network (cluster #3), we obtain lower speedups.
This lower acceleration is due to its lower network
bandwidth, because it causes network contention. Thiscan be explained by Table 5, which presents the ratioof communication time. Almost half of the entire timeis spent for communication, mainly due to image loading from the storage server.

However, this parallelization does not make sense if thecode becomes a trivial bottleneck after applying precomputation, as we presented in Table 3.

3.2. Acceleration Results by Parallelization
Table 4 shows the parallelization results on three clusters. Note here that we basically measure the performanceon a UP conﬁguration. SMP conﬁgurations are employed
only if all CPUs are required to measure the performance.
Our implementation successfully reduces registration
time as the number of CPUs increases, so that registration time is reduced from hours to minutes. As a result, it
achieves linear speedups in most cases. Here, the speedup
Sp on p CPUs is given by Sp = T1 /Tp , where Tp represents registration time on p CPUs. In this table, there are
four important points to be mentioned.

• When using more than 32 nodes with Myrinet (MY)
network (cluster #3), we also obtain lower speedups.
In contrast to FE network, the communication ratio on
MY network is kept a lower value even when using
128 nodes (Table 5). Therefore, this lower acceleration is due to computation rather than communication. The reason for this is that SMP systems, whichshare the memory bus with CPUs in a node, have lower
memory bandwidth, as compared to UP systems. Table 6 presents how this lower bandwidth drops the registration performance. In this table, we can see thatthe parallel efﬁciency for a SMP conﬁguration is 10%
lower than that for UP conﬁgurations, though all thethree conﬁgurations use the same number of CPUs.
Thus, when using 128 CPUs on cluster #3, it takes
longer time to fetch the data from the main memory,
as compared with UP conﬁgurations. In summary, the
SMP architecture decreases the parallel efﬁciency formemory-intensive registration applications.

• GbE is most cost effective network in terms of performance. Although some special networks such as
InﬁniBand (InB) and Myrinet-PCIX (MYX) provide
higher bandwidth with shorter latency (see Table 2),
there is no signiﬁcant difference with respect to theperformance.
• When using less than 8 nodes with MYX (cluster #2),
the speedup exceeds the number of CPUs. This superlinear speedup is due to the memory hierarchy, because our parallel scheme allows CPUs to access a
5

Proceedings of the Sixth International Conference on Parallel and Distributed Computing, Applications and Technologies (PDCAT’05)
0-7695-2405-2/05 $20.00 © 2005 IEEE

Table 4. Parallelization results on p CPUs. Times are presented in seconds.
Registration time (speedup)
p=1
p=2
p=4
p=8
p = 16
p = 32
p = 64
p = 128#1 w/ InB
11,209 (1.0) 5,880 (1.9) 2,811 (4.0) 1,491 (7.5)
767 (15)413 (27)——#1 w/ GbE 11,209 (1.0) 5,682 (1.9) 2,843 (3.9) 1,513 (7.4)
811 (14)457 (25)——#2 w/ MYX 13,536 (1.0) 6,617 (2.0) 3,423 (4.0) 1,681 (8.1)
872 (16)507 (27)——#2 w/ GbE 13,536 (1.0) 7,107 (1.9) 3,685 (3.7) 1,692 (8.0)
889 (15)539 (25)——#3 w/ MY46,997 (1.0) 24,896 (1.9) 12,556 (3.7) 6,369 (7.4) 3,250 (14) 1,756 (27)965 (49)760 (62)#3 w/ FE
46,997 (1.0) 24,915 (1.9) 12,793 (3.7) 6,726 (7.0) 3,705 (13) 2,332 (20) 1,661 (28) 1,547 (30)
ClusterTable 5. Communication-to-computation ratio
presented in percentages.
Cluster#1 w/ InB
#1 w/ GbE
#2 w/ MYX
#2 w/ GbE
#3 w/ MY
#3 w/ FE

2
0.02
0.1
0.04
0.07
0.02
0.4

4
0.1
0.5
0.1
0.2
0.1
1.6

p: number of CPUs
8
16
32
64
0.2
0.5
0.8 —
2.2
5.6 10.8 —
0.3
0.6
1.6 —
0.6
1.2
3.7 —
0.1
0.3
0.9 2.4
5.0 12.7 25.4 42.8

Table 6. Performance comparison between
UP and SMP systems running on a 2-CPU
conﬁguration (cluster #3). Parallel efﬁciency
is given by the speedup divided by 2 CPUs.

128
—
—
—
—
2.6
52.0

Breakdown
SMP UP w/ MY
Image loading (s)3433Gradient Computation (s)16,072
14,052
Similarity Computation (s) 11,598
10,503
Others (s)314308
Total (s)28,01824,896
Parallel efﬁciency (%)
84
94

UP w/ FE113
13,926
10,554322
24,91594

4. Conclusions[4] S. Lee, G. Wolberg, and S. Y. Shin. Scattered data interpolation with multilevel B-splines. IEEE Trans. Visualization
and Computer Graphics, 3(3):228–244, July 1997.
[5] F. Maes, A. Collignon, D. Vandermeulen, G. Marchal, and
P. Suetens. Multimodality image registration by maximization of mutual information. IEEE Trans. Medical Imaging,16(2):187–198, Apr. 1997.
[6] Message Passing Interface Forum. MPI: A message-passing
interface standard. Int’l J. Supercomputer Applications andHigh Performance Computing, 8(3/4):159–416, 1994.
[7] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P.
Flannery. NUMERICAL RECIPES in C: The Art of Scientiﬁc
Computing. Cambridge University Press, Cambridge, UK,
1988.
[8] T. Rohlﬁng and C. R. Maurer. Nonrigid image registrationin shared-memory multiprocessor environments with application to brains, breasts, and bees. IEEE Trans. Information
Technology in Biomedicine, 7(1):16–25, Mar. 2003.
[9] D. Rueckert, L. I. Sonoda, C. Hayes, D. L. G. Hill, M. O.
Leach, and D. J. Hawkes. Nonrigid registration using freeform deformations: Application to breast MR images. IEEE
Trans. Medical Imaging, 18(8):712–721, Aug. 1999.
[10] R. Stefanescu, X. Pennec, and N. Ayache. Grid-enabled nonrigid registration of medical images. Parallel Processing Letters, 14(2):197–216, June 2004.
[11] S. Warﬁeld, F. Jolesz, and R. Kikinis. A high performancecomputing approach to the registration of medical imaging
data. Parallel Computing, 24(9/10):1345–1368, Sept. 1998.
[12] S. K. Warﬁeld, M. Ferrant, X. Gallez, A. Nabavi, F. A.
Jolesz, and R. Kikinis. Real-time biomechanical simulation of volumetric brain deformation for image guided neurosurgery. In Proc. SC’00, pages 1–16, Nov. 2000.

We have presented a performance study of a nonrigid
registration algorithm using three clusters. Our algorithm isbased on an optimization approach with two conventional
acceleration techniques: a data-parallel processing technique; and a precomputation technique.
The experimental results show that (1) the precomputation technique realizes efﬁcient B-spline deformations andreduces registration time by 39%; (2) the data-parallel processing technique is necessary to reduce registration timefrom hours to minutes; (3) it also achieves a linear speedup
on clusters with GbE network; and (4) faster I/O and network infrastructures are essential to achieve higher parallel efﬁciency for large-scale clusters with more than 32
nodes. We also have indicated that sequential implementations should be optimized well enough in advance of parallelization in order to avoid unnecessary parallelization.

References[1] J. V. Hajnal, D. L. Hill, and D. J. Hawkes, editors. Medical
Image Registration. CRC Press, Boca Raton, FL, 2001.
[2] S. Hastings, T. Kurc, S. Langella, U. Catalyurek, T. Pan, and
J. Saltz. Image processing for the grid: A toolkit for building
grid-enabled image processing applications. In Proc. CCGrid’03, pages 36–43, May 2003.
[3] F. Ino, K. Ooyama, and K. Hagihara. A data distributed
parallel algorithm for nonrigid image registration. Parallel
Computing, 31(1):19–43, Jan. 2005.

6

Proceedings of the Sixth International Conference on Parallel and Distributed Computing, Applications and Technologies (PDCAT’05)
0-7695-2405-2/05 $20.00 © 2005 IEEE