Biomedical Image Registration by Means of Bacterial Foraging [601786]

Biomedical Image Registration by Means of Bacterial Foraging
Paradigm
Hariton Costin Silviu Bejinariu Diana Costin
Hariton Costin
1Faculty of Medical Bioengineering Grigore T. Popa University of Medicine and Pharmacy, Ia» si,
Romania
2Institute of Computer Science of Romanian Academy Ia» si Branch, Romania
hncostin@um¯asi.ro
Silviu Bejinariu
2Institute of Computer Science of Romanian Academy Ia» si Branch, Romania
[anonimizat]
Diana Costin*
3Faculty of Medicine, Grigore T. Popa University of Medicine and Pharmacy, Ia» si, Romania
diana.costin@um¯asi.ro
*Corresponding author
Abstract
Image registration (IR) is the process of geometric overlaying or alignment of two or
more 2D/3D images of the same scene (unimodal registration), taken or not at di®erent
time slots, from di®erent angles, and/or by di®erent image acquisition systems (multimodal
registration). Technically, image registration implies a complex optimization of di®erent
parameters, performed at local or/and global level. Local optimization methods often fail
because functions of the involved metrics with respect to transformation parameters are ge-
nerally nonconvex and irregular, and global methods are required, at least at the beginning
of the procedure. This paper presents a new evolutionary and bio-inspired robust approach
for IR, Bacterial Foraging Optimization Algorithm (BFOA), which is adapted for PET-CT
multimodal and magnetic resonance image rigid registration. Results of optimizing the nor-
malized mutual information and normalized cross correlation similarity metrics validated
the e±cacy and precision of the proposed method by using a freely available medical image
database.
Keywords: medical imaging, image registration, soft computing, evolutionary strategies,
bacterial foraging algorithm, global optimization.
1 About Multimodal Image Registration
Image registration (IR) is a fundamental task in computer vision used to ¯nd either a spatial
transformation (e.g., rotation, translation, etc.) or a correspondence (matching of similar image
entities) among two (or more) images taken under di®erent conditions (at di®erent times, using
di®erent sensors, from di®erent viewpoints, or a combination of them), with the aim of overlaying
such images into a common one [1], [2], [3], [4]. Over the years, IR has been applied to a broad
range of situations from remote sensing to medical images or arti¯cial vision and CAD systems,
and di®erent techniques have been independently studied resulting in a large body of research.

INTERNATIONAL JOURNAL OF COMPUTERS COMMUNICATIONS & CONTROL,
ISSN 1841-9836
IR methods can be classi¯ed in two groups according to the nature of images: pixel/voxel –
based IR methods (also called intensity -based), where the whole image is considered for the
registration process; and, on the other side, feature -based methods, which consider prominent
information extracted from the images, being a reduced subset of them. The latter methods
take advantage of the lesser amount of information managed in order to overcome the problems
found in the former when the images present some inconsistences to deal with, for example,
regardless of changes in the geometry of the images, radiometric conditions, and appearance of
noise and occlusion. These features correspond to geometric primitives (points, lines, surfaces,
etc.) which are invariant to the transformation to be considered between the input images.
Moreover, the latter methods perform faster than the former ones due to the reduced amount
of data they take into account, at the expense of achieving coarse results. Likewise, IR is the
process of ¯nding the optimal spatial transformation (e.g., rigid, similarity, a±ne, etc.) achieving
the best overlaying between two (or more) di®erent images named scene and model images
(Figure 1). They both are related with the speci¯c transformation, measured by a similarity
metric function. Such transformation estimation is interpreted into an iterative optimization
procedure in order to properly explore the search space. Two search approaches have been
considered in the IR literature: matching-based , where the optimization problem is intended
to look for a set of correspondences of pairs of those more similar image entities in both the
scene and the model images, from which the registration transformation is derived; and the
transformation parameter-based , where the strategy is to directly explore inside each range of
the transformation parameters. Both strategies can be used with either a voxel-based or a
feature-based approach.
Speci¯c aspects such as the presence of noise, image discretization, di®erent amplitudes
in the scale of the IR transformation parameters, the magnitude of the transformation to be
estimated cause di±culties for traditional local optimizers (gradient- and nongradient-based)
and they become prone to be trapped in local minima. As a consequence, global methods are
preferred, at least at the beginning of the IR process. As for image segmentation procedure,
there is not a universal design for an IR method that could be applicable to all registration
tasks, since various considerations on the particular application must be taken into account.
In recent years a lot of studies and papers were dedicated to medical IR, with more or less
good results [5], [6], [7], [8], [9], [10], [11], [12], [13], [14], [15], [16], [17], [18], [19]. Thus, e.g.
rigid 3D transformations were performed, e.g., by Alpert [20] using the images principal axes
and center of gravity. A±ne registration was obtained by Wahl [21], employing user identi¯ed
anatomical landmarks and external markers, and Maguire et al. [22], who optimized cross-
correlation around such user identi¯ed anatomical landmarks and external markers. In [11]
a robust surface registration using a Gaussian-weighted distance map (GWDM) for PET-CT
brain fusion was proposed. A similarity measure was evaluated repeatedly by weighted cross-
correlation (WCC).
1.1 Transformations
The IR methods can also be classi¯ed according to the registration transformation used to re-
late both the scene and the model images. The ¯rst category of transformation models includes
linear transformations , which preserves the operations of vector addition and scalar multiplica-
tion, being a combination of translation, rotation, global scaling, and shear components. The
most common linear transformations are rigid, similarity, a±ne, projective, and curved. Linear
transformations are global in nature, thus not being able to model local deformations. The
second category of transformation models includes elastic andnonrigid transformations, which
allow local warping of image features, thus providing support for local deformations.

1.2 Similarity Metric
The similarity metric is a function Fthat measures the goodness of a given registration solution,
that is, of a registration transformation f. The ¯nal performance of any IR method strongly de-
pends on its accurate estimation. Each solution is evaluated by Fapplying such transformation
f on one of the two images, usually to the scene image ( f(Is)). Next, the degree of closeness or
¯tting between the transformed scene and the model images, ă( ¢) must be determined
F(Is; Im; f) = ă( f(Is); Im) (1)
The main approaches trying to estimate the function ă( ¢) depend on the dimensionality (2D
or 3D) and the nature of the considered images. There are: (a) voxel-based approach: sum
of squared di®erences, normalized cross-correlation (i.e., correlation coe±cient or phase corre-
lation), and mutual information; (b) feature-based approach: feature values-based metrics (i.e.,
registration based on the curvature) and distance between corresponding geometric primitives.
Unfortunately, the Ffunction is a®ected by both the discretization of images and the presence
of noise, yielding worse estimations and favoring the IR to get trapped in local minima.
1.3 Search Space Strategies
The IR process performs an iterative exploration to obtain that optimal transformation f. So,
the closer fto the unknown global optimum, the better the ¯tting (measured by the similarity
metric F) between scene and model. The optimization process considered to obtain those
solutions can be deterministic or stochastic (either a global or a local one).
Although the ¯nal registration problem solution consists of the right values for the parameters
which determine f, we can distinguish two di®erent strategies to solve the problem, each of them
working in a di®erent solution space: (i) the ¯rst searches in the matching space to obtain a
set of correspondences of pairs of the most similar image entities in both the scene and the
model images, from which the registration transformation is derived; (ii) the second directly
makes a search in the space of the f parameters guided by the F function, called transformation
parameters space. The matching-based search space exploration usually consists of the two
following stages: ¯rst, a set of correspondences with those more similar regions of pixels (voxel-
based) or geometric primitives (feature-based) in both the scene and the model images must
be computed; second, the transformation f is assessed by numerical methods considering the
previous matching.
On the contrary, transformation parameters-based search space involves directly searching
for the solution in the space of parameters of the transformation f. In this respect, each solution
to the IR problem is encoded as a vector composed of the values for the parameters of f, and
the IR method generates possible vectors of parameter values, that is, possible registration
transformations. As a consequence, the search space exploration is guided by the similarity
metric F. In this way, each solution vector is evaluated by the chosen metric, and the IR problem
becomes a parameter optimization procedure of ¯nding the best values of fthat maximize
the similarity metric F. Other classi¯cation divides search strategies in local andglobal ones.
Local optimization techniques frequently fail because functions of these metrics with respect to
transformation parameters are generally nonconvex and irregular and, therefore, global methods
{ such as those based on evolutionary algorithms { are often required.
2 Optimization using BFOA
In recent years, the application of several well-known evolutionary algorithms (EAs) [23] to
the IR optimization process has introduced an outstanding interest in order to solve those

INTERNATIONAL JOURNAL OF COMPUTERS COMMUNICATIONS & CONTROL,
ISSN 1841-9836
problems due to their global optimization techniques nature. The ¯rst attempts to solve IR
using evolutionary computation can be found in the early eighties, when Fitzpatrick et al.
[24] proposed such approach based on a genetic algorithm for the 2D case and applied it to
angiographic images. Since then, several evolutionary approaches have been proposed to solve
the IR problem, mainly in connection with the transformation parameters-based search space,
as shown e.g. in [25], [1], [26], [27], [28], [29], [30]. The main reason of using global optimization
techniques, such as EAs-based algorithms for IR, is that they may not require an optimum
solution to achieve high accuracy of registration.
Introduced by Passino [31], [32], bacterial foraging paradigm is a bio-inspired optimization
method based on the foraging model. This paradigm belongs to the broader class of distributed
nongradient global optimization. A foraging animal takes actions to maximize the energy ob-
tained per unit time spent foraging, E/T, in the face of constraints presented by its own physi-
ology (e.g., sensing and cognitive capabilities) and environment (e.g., density of prey, risks from
predators, physical characteristics of the search area). In other words, these social animals, like
E. coli { a bacterium, try to maximize their long-term average rate of energy intake.
Prominent applications in medical image processing are related to edge detection [33], image
segmentation [34], [35], but also for image registration [36], [37], knowing that these procedures
may be viewed as optimization tasks.
2.1 BFO Algorithm
The bacterial foraging paradigm is suitable as model for optimization algorithms because
animals/bacteria behavior is to search for nutrients and avoid noxious substances to maximize
their energy. As in all evolutionary models, individuals with a good strategy to ¯nd nutrients
are replicated and those having poor foraging strategy are eliminated. In contrast to genetic al-
gorithms and evolutionary strategies, which exploit the competitive characteristics of biological
evolution (e.g., survival of the ¯ttest), bacterial foraging (BF) exploits cooperative and social
aspects of animal colonies (like E. coli bacterium) in their attempts to obtain nutrients that
maximizes energy intake per unit time spent for foraging.
Each member of the bacteria colony is characterized by its position in the n-dimensional
space which is a possible solution of the optimization problem. The solution is computed as
the position in which a bacterium is in the best healthy state or has the lowest cost value.
During foraging, the bacteria colony (swarm) proceeds through four foraging steps: chemotaxis,
swarming, reproduction and elimination-dispersal.
Let's consider a bacteria colony with Sindividuals; P(j; k; l ) =fµt(j; k; l ); i= 1:::Sgis the
position of colony members in the jthchemotactic step, kthis the k-threproduction step and lth
{ the l-thelimination-dispersal step; J(i; j; k; l ) denotes the cost of the ithbacterium in position
µt(j; k; l ).
-Chemotaxis : E. Coli bacteria have two types of movements: tumble and swim. The
chemotactic step is de¯ned as a tumble followed by a tumble or a tumble followed by
a run. In the chemotactic step each bacterium changes its position to: µt(j+ 1; k; l) =
µi(j; k; l ) +C(i)'(i), where C(i) is the size of the chemotactic step and is a unit length
random generated direction [4]. If the cost computed in the new position is lower than in
the previous position, then the swim is continued in the same direction as long as the cost
is reduced but not more than a maximum number of steps.
-Swarming : In case the bacteria have the ability to signal to others the existence of a
favorable or poisonous environment, they will tend to swarm together in the direction of
nutrients. The cell to cell attraction or rejection is modeled by adding to the cost function
J(i; j; k; l ) computed for a speci¯c bacterium, components computed as function of the
status of all the other bacteria in the colony.

-Reproduction : All bacteria reach the reproduction state after a number of chemotactic
steps. The healthy state is computed for all bacteria and it may be expressed as the total
quantity of accumulated nutrients or simply by the value of the cost function in the current
position. The least healthy bacteria die while and to keep constant the size of the colony,
an equal number of healthier bacteria split into two bacteria without mutation.
-Elimination and Dispersal : After a number of reproduction steps, some bacteria are dis-
persed into the environment (moved in a random position) with a speci¯ed probability,
without taking into account their healthy state.
BFOA starts with a colony of Sbacteria placed in randomly generated positions. The evo-
lutionary process goes through Nedelimination/dispersal steps, each of these consists of Nre
reproduction steps. Each reproduction step consists of Ncchemotactic steps. In each chemo-
tactic step a bacterium may do at most Nsswimming steps if the value of the cost function
decreases. The position in which a bacterium has the greatest healthy status is the solution of
the optimization problem. In case of image registration, the size of the search space is equal to
the number of parameters of the geometric transform and as healthy status is used the value of
a measure that evaluates the similarity between the transformed source image and the model
image.
2.2 Parallel Version of BFO Algorithm
A closer look at BFOA reveals that it contains 4 nested loops: elimination/dispersal, reproduc-
tion and chemotaxis for each bacterium in the colony. The body of the inner loop is executed
Ned£Nre£NC£Stimes, which may be a fairly large number. In the examples presented
in the next sections, it is executed 256000 times, but the cost function evaluation is performed
about 600000 times due to the fact that each bacterium may perform more swim steps in a
single chemotactic step.
We propose a parallelization based on the shared memory model that is suitable for multi-core
processor based systems. It must be noticed that in this case the number of available processors
is reduced (2 ;4 or 8). An excessive tasks partitioning obviously leads to poor performances due
to the large number of synchronization operations.
If we consider to not use the attractant/repellant e®ect in the optimization algorithm, then
the calculations performed for each individual bacterium in the inner loop are independent
excepting the test in which the best value of the cost function is checked. So, we can execute in
parallel a chemotactic step for all bacteria, taking care to not simultaneously call the function
to check for the best value of the cost function.
3 Pixel Based Image Registration
The proposed IR procedures use the Normalized Mutual Information and Normalized Cross
Correlation as measures to evaluate the quality of the registration process.
Our study approaches the rigid body image registration , which initially determines global
alignment, followed by local elastic registration. Let Tdenote the spatial transformation that
maps features or coordinates (spatial locations) from one image or coordinate space to another
image or coordinate space. Let pAandpBdenote coordinate points (pixel locations) in images
AandB, respectively. The image registration problem is to determine Tso that the mapping
T:pA!pB,T(pA) =pBresults in the "best" alignment of AandB. For 3-D rigid body
registration, the mapping of coordinates p= [x; y; z ]Tintop0= [x0; y0; z0]Tcan be formulated
as a matrix multiplication in homogeneous coordinates, as shown in equation (2) in an explicit
manner. That is, the goal of the optimization is to determine the parameters tx; ty; tz; ®; ¯; and

INTERNATIONAL JOURNAL OF COMPUTERS COMMUNICATIONS & CONTROL,
ISSN 1841-9836
°in (2). Usually, optimization in image registration means to maximize similarity. Similarity
metric values, as functions of transformation parameters, refer to the objective function , denoted
asf(x). Alternatively, one may formulate the image registration as a minimization problem and,
thus, the goal is to minimize ¡f(x). Although there is yet no proof for its optimality, because
of its robustness (usually it attains its maximum at correct alignment) and good results in
previous works, normalized mutual information was selected as the similarity measure in our
study. Moreover, it is still generally non-smooth and prone to local optima. For this reason,
global optimization approaches are preferred.
"
p0
1#
=T"
p
1#
,
2
6666664×0
y0
z0
13
7777775=2
66664cos¯cos°cos®sin°+ sin ®sin¯cos°sin®sin°¡cos®sin¯cos° t x
¡cos¯sin°cos®cos°¡sin®sin¯sin°sin®cos°¡cos®sin¯sin° t y
sin¯ ¡sin®cos¯ cos®cos¯ t z
0 0 0 13
777752
66664x
y
z
13
77775(2)
3.1 Normalized Mutual Information
Mutual Information ( MI) and Normalized Mutual Information ( NMI ) evaluate the relative
independence of two images and do not depend on the speci¯c dynamic range or intensity
scaling of the images [1], [10]:
MI(A; B) =H(A) +H(B)¡H(A; B) (3)
NMI (A; B) = (H(A) +H(B))=H(A; B); (4)
where H(A); H(B) are the image entropies and H(A; B) is the joint entropy of the two images.
High values of mutual information indicate high dependence between images. Because the goal
of the optimization algorithms is to minimize a cost function, the value of ( ¡1)¤NMI will be
used to evaluate the quality of a certain solution. In the cost function evaluation, the geometric
transform corresponding to the current solution is applied to the source image and then the NMI
value is computed for the model image and the transformed source image. The area based IR
implementations are time consuming because each cost evaluation requires a geometric transform
to be applied and also image and matrix operations to compute NMI .
3.2 Normalized Cross Correlation
Cross correlation is used for estimating the degree to which two series are correlated. One of the
most encountered applications of the normalized cross correlation is to determine the position
of a template sub-image Bin a source image A. The normalized cross correlation (NCC) is
computed by
NCC (i; j) =§³§´A(i+³; j+´)¢B(³;´)p
§³§´B(³;´)2¢p
§³§´A(i+³; j+´)2: (5)
The problem is to determine the position of a given pattern in a two dimensional image f.
Letf(x; y) denote the intensity value of the image fof size Mx£Myat the point ( x; y);
x2 f0; : : : ; M x¡1g; y2 f0; : : : ; M y¡1g. The pattern is represented by a given template tof
sizeNx£Ny. A common way to calculate the position ( upos; vpos) of the pattern in the image
fis to evaluate the normalized cross correlation value °at each point ( u; v) for fand template

t, which has been shifted by usteps in the xdirection and by vsteps in the ydirection. The
following equation gives a basic de¯nition for the normalized cross correlation coe±cient
°=§x;y(f(x; y)¡fu;v)(t(x¡u; y¡v)¡t)q
§x;y(f(x; y)¡fu;v)2§x;y(t(x¡u; y¡v)¡t)2; (6)
where ( fu;v) denotes the mean value of f(x; y) within the area of the template t shifted to ( u; v),
which is computed by:
(fu;v) =1
NxNyu+Nx¡1X
x=uv+Ny¡1X
y=vf(x; y): (7)
With similar notation, tis the mean value of the template t. The denominator in (6) is the
variance of the zero mean image function f(x; y)¡fu; v and the shifted zero mean template
function t(x¡u; y¡v)¡t. Due to this normalization, °(u; v) is independent to changes in
brightness or contrast of the image, which are related to the mean value and the standard
deviation.
4 Experiments and Results
The BFOA based image registration procedure was tested on a large set of DICOM medical
images from a database at the address http://www.osirix-viewer.com/datasets/ [39]. In
Figure 1 below, information about some test images are shown.
Image Description:
File : img 1.tif
Size : 256 £256£8b
Name : PETCETIX
Modality : PET-CT
Description: Whole body FDG PET-CT study
in a patient with abdominal lymphoma :
File: img 2.tif
Size: 256 x 256 x 8b
Name: BRAINIX
Modality: MR
Description: Brain tumor.
File: img 3.tif
Size: 256 x 256 x 8b
Name: WRIX
Modality: MRI
Description: Scaphoid fracture. T1 / STIR fusion.
Figure 1: Test images used in the experiment
In the experiment, 4 images were used: img 1;img2;img2HistEq (that was obtained by
applying the histogram equalization to img 2 having low contrast), and img 3.
The source images were obtained by applying a rotation (angle µ= 10±) against the rotation
center ( cx=¡20 and cy= 20) followed by an isotropic scaling ( scale = 1:2). The transform
matrix is:
T=2
64® ¯ (1¡®)cx¡¯cy
¡¯ ® ¯c x+ (1¡®)cy
0 0 13
75 (8)

INTERNATIONAL JOURNAL OF COMPUTERS COMMUNICATIONS & CONTROL,
ISSN 1841-9836
where ®=scale¢cosµand¯=scale¢sinµ. The actual value of the transform matrix is
T=2
641:1818 0 :2084 ¡0:5322
¡0:2084 1 :1818 ¡7:8029
0 0 13
75 (9)
The inverse transform matrix is T¡1=2
640:8207 ¡1:1477 ¡0:6924
1:1447 0 :8207 6 :4807
0 0 13
75that corresponds to an
a±ne transform with the following parameters: µ0=¡10±,c0
x=¡20,c0
y= 20 and scale0=
0:8333:The search space for the optimization problem is R4.
The BFO parameters values used in the experiment are:
Colony size, S= 400; Number of chemotactic steps, Nc= 20; Maximum number of swim
steps, Ns= 10; Number of reproduction steps, Nre= 16; Number of elimination/dispersal steps,
Ned= 2; Probability of dispersal, Ped= 0:25; Length of the move step= 0:001:
As it can be observed in Figure 1 the test images have the same size but are di®erent in
terms of contrast. It was expected that the image contrast a®ects the quality of the registration
process, and this assumption was found to be true.
In vision, contrast is the di®erence in luminance that makes an object distinguishable. The
test images we used in this paper were not analyzed for their content (i.e., it is a context-free
registration), so the contrast and his evaluation was performed by means of the histogram of
the images.
Root mean square contrast is computed as standard deviation of pixel values. It does not
depend on the spatial frequency or spatial distribution in the image.
Contrast RMS =q
1
MNPN
i=1PM
j=1(I(i; j)¡Iavg)2; (10)
where M; N are the image dimensions, I(i; j) is the image pixel at ( i; j) coordinates and Iavgis
the average of pixel values in image I:
Visibility (Michelson contrast) is represented by formula
Contrast Michelson =Imax¡Imin
Imax+Imin; (11)
where IminandImaxare the lowest and highest pixel values in image I, respectively.
Contrast values computed for the tested model images are described in Table 1.
Table 1: Contrast values computed for the test images
For image 2:tif with initial low contrast a histogram equalization procedure was performed
and image named img 2HistEq was obtained. The source images were obtained by applying
the a±ne transform speci¯ed above followed by salt-and-pepper and Gaussian noise alteration.

The Signal-To-Noise Ratio was computed using the formula below and is shown in Table 2.
SNR = 10 logPN
i=1PM
j=1I(i; j)2
PN
i=1PM
j=1(I(i; j)¡Inoise(i; j))2dB; (12)
Table 2: Signal-to-Noise Ratio (SNR) values for testimages with noise
The registered images for test image img 1 are shown in Figure 2.
a:Model image,
img1b:Source image,
no noise addedc:Source image,
salt and pepper
noise added,
SNR=1.63 dBd:Source image,
Gaussian noise
added
SNR=20.14 dB
e. Registered image f. Registered image g. Registered image
Figure 2: Registered images for test image img 1
In Table 3, the columns NMI/Expected and NCC/Expected contain respectively the va-
lues of Mutual Information and Normalized Cross Correlation computed for the model image
and the image obtained by applying the inverse a±ne transform on the source images. The
NMI/Computed and NCC/Computed contain respectively the Mutual Information and Nor-
malized Cross Correlation computed for the model image and the registered image obtained
by applying the approximated a±ne transform. The ¯rst 3 rows in Table 3 contain the values
obtained using MI as similarity measure while the last 3 rows contain the results obtained by
using NCC as similarity measure in the registration evaluation.

INTERNATIONAL JOURNAL OF COMPUTERS COMMUNICATIONS & CONTROL,
ISSN 1841-9836
Table 3: Image registration results obtained for test image img 1
In the case of the image drawn in Figure 2, it must be noticed that in all cases the value of the
similarity measure computed for the registered image does not exceed the expected value. This
happens because the image histogram is uniform (the gray levels are more uniformly distributed
in the image), even if there does exist a considerable number of black pixels. Accordingly, the
results obtained using the two similarity measures are quite similar.
Similarly, registered images for test image img 2 are depicted in Figure 3.
a:Model image,
img2b:Source image,
no noise addedc:Source image,
salt and pepper
noise added,
SNR=-15.04 dBd:Source image,
Gaussian noise
added
SNR=5.76
e. Registered image f. Registered image g. Registered image
Figure 3: Registered images for test image img 2
In the case of test image img 2, even if the image has a low contrast, the IR result is also good.
The image contrast is Contrast RMS = 0:05 (a very low value) and in all cases the similarity
measure error is less than 1/1000, so the proposed BFOA-based IR method works well in case
of low contrast images. Moreover, for very noisy images (SNR as low as ¡15:04dB) this method
also works very well, as shown in Figures 3 ( f; g):

Table 4: Image registration results obtained for test image img 2
a:Model image,
img3b:Source image,
no noise addedc:Source image,
salt and pepper
noise added,
SNR=-2.28 dBd:Source image,
Gaussian noise
added
SNR=15.80
e. Registered image f. Registered image g. Registered image
Figure 4: Registered images for test image img 3
Table 5: Image registration results obtained for test image img 3
In case of the test image img 3(Table 5), it must be noticed that in all experimental situations the
computed NCC value is very slightly greater than the expected value. This may be explained
by: (1) the low contrast of model and source images, (2) the images have a high portion of
dark background and (3) the black areas inserted (as a necessity, from experimental reasons)

INTERNATIONAL JOURNAL OF COMPUTERS COMMUNICATIONS & CONTROL,
ISSN 1841-9836
in the transformed and registered images after the a±ne transforms application. But also in
these conditions the values of MI seem to be more relevant and give a real evaluation of the
registration process.
a:Model image b:Source image,
no noise addedc:Source image,
salt and pepper
noise added,
SNR=4.12 dBd:Source image,
Gaussian noise
added
SNR=22.49
e. Registered image f. Registered image g. Registered image
Figure 5: Registered images for test image img 2HistEq
The proposed IR procedure was also applied to the low contrast image img 2, processed by
applying the histogram equalization (Figure 5). In this way the contrast value was increased
from 0.05 to 0.35. It must be noticed that in this case the registration results are similar to the
case of the original image, img 2 (Table 6).
Table 6: Image registration results obtained for test image img 2HistEq
4.1 Execution Time and Parallel Approach of the Algorithm
The following tables contain for each test image the execution time in seconds, the number of
cost function evaluations for sequential and parallel executions and also the parallel e±ciency
(Eff). The most common evaluation of parallel algorithms is performed using the parallel
e±ciency Eff =ts
tp£n, where tsis the time used by the sequential version of the algorithm, tp
is the processing time for the parallel version and nis the number of used processors [38].
When comparing these values, the following issues must be considered:

{The parallel implementation was evaluated on Intel Core i5 3.10 GHz processor, with 4
cores. The system has 4 GB RAM and uses Windows 7 (64 bits) as operating system.
The application was compiled as a 32 bits application;
{The execution times have to be considered as approximations, because the application
was executed in Windows OS without obtaining the exclusive access of processor. The
processor allocation during parallel execution was made by the operating system;
{In all executions of BFOA, the random number string was the same (the generator was
initialized using the same value).
In Table 7 and Table 8 below the execution time (in seconds), the number of cost function
evaluations and parallel e±ciency are presented.
Table 7: Parallel e±ciency obtained for images img 1 and img 3
Table 8: Parallel e±ciency obtained for images img2 and img 2HistEq
Also, Figures 6 and 7 show parallel versus sequential execution time for img 2 and img 2HistEq
respectively, for clean and di®erent added noise.
Figure 6: Parallel and sequential execution time in case of image img 2HistEq

INTERNATIONAL JOURNAL OF COMPUTERS COMMUNICATIONS & CONTROL,
ISSN 1841-9836
Figure 7: Parallel and sequential execution time in case of image img 2
As is depicted in Table 7 and Table 8, the parallel e±ciency value is between 0.78 and 0.86.
Higher values of e±ciency are obtained for longer tasks, when normalized cross correlation is
used for image similarity evaluation. It must be noticed that in this case, the parallel version
execution time is similar to the sequential version execution time when NMI is used for similarity
evaluation. Considering that also the registration precision is higher in case of NMI usage, it can
be concluded that NMI is recommended as evaluation measure. The obtained parallel e±ciency
value is high, considering the time required for synchronization tasks and critical sections which
cannot be parallelized. It must be noticed that in case of images img 2 and img 2HistEq, for
the sequential image registration, the cost function is called the same number of times. The
execution time is di®erent because in Windows operating system the processor is also used for
other tasks. In case of parallel image registration, the number of cost function evaluation is
di®erent for the di®erent allocation of the processor cores. Figures 6 and 7 show processing
time for di®erent added noise and the two similarity measures considered, for sequential and
parallel versions of the proposed algorithm, when registering images img 2HistEqandimg 2,
respectively.
5 Conclusions
This study proved the feasibility of rigid and mono-modal image registration by using a new
optimization approach – Bacterial Foraging Optimization Algorithm (BFOA), a bio-inspired
technique that belongs to the large family of evolutionary computing and metaheuristic meth-
ods. As similarity measures between two images performing registration process, we used the
normalized mutual information and normalized cross correlation which had to be optimized by
BFOA. The obtained results are encouraging, i.e. the accuracy of registration process as high
even in the case of noisy images, with very low signal-to-noise ratios. In this way our method
might be considered as a robust registration technique. Yet, the maximum expected value of the
IR evaluation measure was not reached because when processing images from the used database
two technical details have limited the overall registration:
{after the a±ne transform applied to obtain the source image and inverse transform is
applied to obtain the registered image, some black areas are inserted in some images due
to geometrical experimental reasons;
{the pixels values are changed because during the direct and inverse a±ne transform the
pixels values are approximated using interpolation methods. In this experiment the bilinear
interpolation is used.
Concerning the possibility to accelerate the execution time of the registration process, it is
noteworthy that the proposed algorithm is suited to parallelization, as shown above. In this
respect runtimes at least 3 times lower of the parallel version than of sequential approach were
obtained, so the parameter named "speedup" equals 3. When computing parallel e±ciency, it is
worth mentioning that values of 0.78{0.86 were obtained when using the Intel Core i5 processor
(4 cores), that meaning a good e±ciency of the parallel approach.

References
[1]Cordon, O.; Damas, S.; Santamaria, J. (2006); Feature-based Image Registration by Means
of the CHC Evolutionary Algorithm, Image and Vision Computing , 24(5): 525-533.
[2] Gonzales, R.C.; Woods, R.E. (2002); Digital Image Processing (2nd ed.), Prentice Hall,
New Jersey.
[3]Pratt, W.K. (2001); Digital Image Processing , John Wiley & Sons, New York.
[4]Rangayyan, R.M. (2005); Biomedical Image Analysis , CRC Press, Boca Raton, 2005.
[5]Alterovitz, R., et al. (2006); Registration of MR Prostate Images with Biomechanical
Modeling and Nonlinear Parameter Estimation, Medical Physics , 33(2): 446-454.
[6]Cooper, J. (2003); Optical Flow for Validating Medical Image Registration, Proc. of the
9th IASTED Int. Conference on Signal and Image Processing , IASTED/ACTA Press:
502-506.
[7]He, R.; Narayana, P.A. (2002); Global Optimization of Mutual Information: Application
to Three Dimensional Retrospective Registration of Magnetic Resonance Images, Com-
puterized Medical Imaging and Graphics , 26(4): 277-292.
[8]Hill, D.; Studholme, C.; Hawkes, D. (1994); Voxel Similarity Measures for Automated Im-
age Registration, Proceedings of the Third SPIE Conference on Visualization in Biomedical
Computing : 205-216.
[9]Lavavely, W.C.; Scarfone, C. et al. (2004); Phantom Validation of Co-Registration of PET
and CT for Image-Guided Radiotherapy, Medical Physics , 31(5): 1083-92.
[10]Levin, D.N.; Pelizzari, C.A. et al. (1988); Retrospective Geometric Correlation of MR,
CT, and PET Images, Radiology , 169(3): 817-823.
[11]Lee, H.; Hong, H. (2006); Robust Surface Registration for Brain PET-CT Fusion, Medical
Imaging 2006: Visualization, Image-Guided Procedures, and Display , Ed. by Cleary, Kevin
R.; Galloway, Robert L., Jr., Proceedings of the SPIE , Vol. 6141: 684-693.
[12]Maintz, B.A.; van den Elsen, P.A.; Viergever, M.A. (1996); Registration of SPECT and
MR Brain Images Using a Fuzzy Surface, in Loew, M.H. and Hanson, K.M. (eds), Medical
Imaging: Image processing , Bellingham, WA. SPIE, Vol. 2710: 821-829.
[13]Maintz, B.A.; Viergever, M.A. (1998); A Survey of Medical Image Registration, Medical
Image Analysis , Oxford University Press, 2(1): 1-37.
[14]Pluim, J.P.; Maintz, J.B.A.; Viergever, M.A. (2003); Mutual Information-based Regis-
tration of Medical Images: A Survey, IEEE Transactions on Medical Imaging , 22(8):
986-1004.
[15]Pietrzyk, U., Herholz, K. et al., (1996); Clinical Applications of Registration and Fusion
of Multimodality Brain Images from PET, SPECT, CT, and MRI, European Journal of
Radiology , 21: 174-182.
[16]Chen, Qin-Sheng (1993); Image Registration and Its Applications in Medical Imaging , PhD
thesis, Vrije Universiteit Brussel.
[17]Zibaeifard, M.; Rahmati, M. (2001); An Improved Multi-stage Method for Medical Image
Registration Based on Mutual Information, Proceedings of the 8-th International Confer-
ence on Computer Vision : 718-725.
[18]Xuan, J.; Wang,Y. et al.(2006); Nonrigid Medical Image Registration by Finite-element
Deformable Sheet-curve Models, Int. Journal of Biomedical Imaging , 2006, Article ID
73430: 1-9, Hindawi Publishing Corporation.
[19]Zitova, B., Flusser, J. (2003); Image Registration Methods: A Survey, Image and Vision
Computing , 21, Elsevier, 977-1000.
[20]Alpert, N.M.; Bradshaw, J.F.; Kennedy, D.; Correia, J.A. (1990); The Principal Axis
Transformation { A Method for Image Registration, Journal of Nuclear Medicine , 31,
1717-1722.

INTERNATIONAL JOURNAL OF COMPUTERS COMMUNICATIONS & CONTROL,
ISSN 1841-9836
[21]Wahl, R.L.; Quint, L.E.; et al. (1993); "Anametabolic" Tumor Imaging: Fusion of FDG
PET with CT or MRI to Localize Foci of Increased Activity. Journal of Nuclear Medicine ,
34: 1190-1197.
[22]Maguire, G.Q. Noz, M. et al. (1991); Graphics Applied to Medical Image Registration.
IEEE Computer Graphics and Applications , 11(2): 20-28.
[23]Back, T.; Fogel, D.B.; Michalewicz, Z. (Eds.) (1997); Handbook of Evolutionary Compu-
tation , IOP, Bristol, UK; Oxford University Press, Oxford, UK.
[24]Fitzpatrick, J.M.; Grefenstette, J.J.; Van Gucht, D. (1984); Image Registration by Genetic
Search, Proc. of the IEEE Southeast Conference , Louisville, USA, 460-464.
[25]Chow, C. K.; Tsui, H.T.; Lee, T. (2004); Surface Registration Using a Dynamic Genetic
Algorithm, Pattern Recognition , 37(1): 105-117.
[26]Cordon, O.; Damas, S.; Santamaria, J. (2007); A Practical Review on the Applicability
of Di®erent Evolutionary Algorithms to 3D Feature-based Image Registration, in Genetic
and Evolutionary Computation for Image Processing and Analysis , Eds.: S. Cagnoni, E.
Lutton, G. Olague, Hindawi Publishing Corporation, 241-264.
[27]Etienne, E.K.; Nachtegael, M. (2000) (eds.). Fuzzy Techniques in Image Processing ,
Physica-Verlag, N.Y.
[28]Reyes-Sierra, M.; Coello C.A. (2006); Multi-objective Particle Swarm Optimizers: A Sur-
vey of the State-of-the-art, International Journal of Computational Intelligence Research ,
2(3): 287-308.
[29]Rouet, J.-M.; Jacq, J.-J.; Roux, C. (2000); Genetic Algorithms for a Robust 3-D MR-CT
Registration, IEEE Trans. on Information Technology in Biomedicine , 4(2): 126-136.
[30]Wachowiak, M. P.; Smolkov, R. et al. (2004); An Approach to Multimodal Biomedical
Image Registration Utilizing Particle Swarm Optimization, IEEE Transactions on Evolu-
tionary Computation , 8(3): 289-301.
[31]Passino, K.M. (2002); Biomimicry of Bacterial Foraging for Distributed Optimization and
Control, IEEE Control Systems Magazine , 52-67.
[32]Liu, Y.; Passino, K.M. (2002); Biomimicry of Social Foraging Bacteria for Distributed Op-
timization: Models, Principles, and Emergent Behaviors, Journal of Optimization Theory
and Applications , 115(3): 603-628.
[33]Verma, O. P.; Hanmandlu, M.; Kumar, P.; Chhabra, S.; Jindal, A. (2011); A Novel
Bacterial Foraging Technique for Edge Detection, Pattern Recognition Letters , Elsevier,
32: 1187-1196.
[34]N. Sanyal, A. Chatterjee, S. Munshi, "An Adaptive Bacterial Foraging Algorithm for Fuzzy
Entropy Based Image Segmentation", Expert Systems with Applications , 38, Elsevier, 2011,
pp. 15489-15498.
[35]Sathya, P.D.; Kayalvizhi, R. (2011); Modi¯ed Bacterial Foraging Algorithm Based Multi-
level Thresholding for Image Segmentation, Engineering Applications of Arti¯cial Intelli-
gence , Elsevier, 24: 595-615.
[36]Yudong, Z.; Lenan, W. (2008); Multi-resolution Rigid Image Registration Using Bacterial
Multiple Colony Chemotaxis, 5th Int. Conf. on Visual Information Engineering , VIE 2008:
528-532.
[37]Bejinariu, S. I.; Costin, H., Rotaru, F.; Nit »¸ a, C.; Luca, R.; Laz¸ ar, C. (2014); Parallel
Processing and Bioinspired Computing for Biomedical Image Registration, The Computer
Science Journal of Moldova , 22(2): 253-277.
[38]Campbell, C.; Miller, A. (2012); Parallel Programming with Microsoft Visual C++, Mi-
crosoft Corporation.
[39]DICOM Sample Image Sets (http://www.osirix-viewer.com/datasets/), accessed on
1.09.2014