# Appendix Scatter Search for the Point Matching Problem in 3D Image Registration

Appendix "Scatter Search for the Point Matching Problem in 3D Image Registration"

Oscar Cordón

Departamento de Ciencias de la Computación e Inteligencia Artificial, Universidad de Granada, Daniel Saucedo Aranda s/n, 18071 Granada, Spain, ocordon@decsai.ugr.es European Centre for Soft Computing. Edificio Científico-Tecnológico, C/ Gonzalo Gutiérrez Quirós, s/n, 33600 Mieres (Asturias), Spain, oscar.cordon@softcomputing.es

Sergio Damas

Departamento de Lenguajes y Sistemas Informáticos, Universidad de Granada, Daniel Saucedo Aranda s/n, 18071 Granada, Spain, sdamas@ugr.es European Centre for Soft Computing. Edificio Científico-Tecnológico, C/ Gonzalo Gutiérrez Quirós, s/n, 33600 Mieres (Asturias), Spain, sergio.damas@softcomputing.es

Jose Santamaría

Departamento de Lenguajes y Sistemas Informáticos, Universidad de Cádiz, Chile 1, 11002 Cádiz, Spain, jose.santamarialopez@uca.es

Rafael Martí

Departamento de Estadística e I.O., Universidad de Valencia, Dr. Moliner 50, 46100 Burjassot, Valencia, Spain, rafael.marti@uv.es

1. Scatter-Search Methodology

The fact that the mechanisms within SS are not restricted to a single uniform design allows exploration of strategic possibilities that may prove effective in a particular implementation. These observations and principles lead to the following “five-method template” for implementing SS: 1. A diversification-generation method to generate a collection of diverse trial solutions, using an arbitrary trial solution (or seed solution) as an input. 2. An improvement method to transform a trial solution into one or more enhanced trial solutions. Neither the input nor the output solutions are required to be feasible, though the output solutions will more usually be expected to be so. If no improvement of the input trial solution results, the “enhanced” solution is considered to be the same as the input solution. 3. A reference-set update method to build and maintain a reference set consisting of the b “best” solutions found (where the value of b is typically small, e.g. no more than

1

20), organized to provide efficient accessing by other parts of the method. Solutions gain membership to the reference set according to their quality or their diversity. 4. A subset-generation method to operate on the reference set, to produce several subsets of its solutions as a basis for creating combined solutions. 5. A solution-combination method to transform a given subset of solutions produced by the Subset Generation Method into one or more combined solution vectors. Figure 1 shows the interaction among these five methods and highlights the central role of the reference set. This basic design starts with the creation of an initial set of solutions P, and then extracts the reference set (RefSet) of solutions from it. The darker circles represent improved solutions resulting from the application of the improvement method.

P Diversification Generation Method Repeat until |P| = PSize

Improvement Method

Improvement Method

Reference Set Update Method

Solution Combination Method

Subset Generation Method RefSet Stop if no more new solutions

Figure 1: Schematic Representation of a Basic SS Design The diversification-generation method is used to build a large set P of diverse solutions. The size of P (PSize) is typically at least 10 times the size of RefSet. The initial reference set is built according to the reference-set-update method, which can take the b better solutions (with regard to their quality in the problem solving) from P to

2

compose the RefSet. However, diversity can be considered instead of, or in addition to, quality for the updating. For example, the reference-set update method could consist of selecting b distinct and maximally diverse solutions from P. Regardless of the rules used to select the reference solutions, the solutions in RefSet are ordered according to quality, where the best solution is the first one in the list. The search is then initiated by applying the subset-generation method that, in its simplest form, involves generating all pairs of reference solutions. The pairs of solutions in RefSet are selected one at a time and the solution-combination method is applied to generate one or more trial solutions. These trial solutions are subjected to the improvement method. The reference-set-update method is applied once again to build the new RefSet with the best solutions, according to the objective-function value, from the current RefSet and the set of trial solutions. The basic procedure terminates after all the subsets generated are subjected to the combination method and none of the improved trial solutions are admitted to RefSet under the rules of the reference-set-update method. The reference set, RefSet, is a collection of both high quality solutions and diverse solutions that are used to generate new solutions by applying the combination method. We can use a simple mechanism to construct an initial reference set and then update it during the search. The size of the reference set is denoted by |RefSet| = b1 + b2 = b. The construction of the initial reference set starts with the selection of the best b1 solutions from P. These solutions are added to RefSet and deleted from P. For each solution in PRefSet, the minimum of the distances to the solutions in RefSet is computed. The solution with the maximum of these minimum distances is then selected. This solution is added to the RefSet and deleted from P, and the minimum distances are updated. The process is repeated b2 times, where b2 = b – b1. The resulting reference set has b1 highquality solutions and b2 diverse solutions. Of the five methods in SS methodology, only four are strictly required. The improvement method is usually needed if high-quality outcomes are desired, but an SS procedure can be implemented without it. On the other hand, hybrid SS designs could incorporate a short-term tabu search or another complex metaheuristic such as the improvement method (usually demanding more running time).

3

2. 3D Images Considered in the Computational Study

a) The MRI dataset of human brains Our results for this dataset correspond to a number of registration problems with four different simulated real-world MRIs. These images have been obtained from the BrainWeb database at McGill University (Collins et al. 1998, Kwan et al. 1999). The purpose of this simulator is to provide researchers with ground truth data for imageanalysis techniques and algorithms. BrainWeb has been widely used by the IR research community (Rogelj and Kovacic 2002, Held et al. 2004, Wachowiak et al. 2004). Moreover, we have added different levels of noise to three of the four images used. The reason is to model noisy conditions related to the images acquired by some devices. Likewise, we cannot avoid one of the most important goals of IR: supporting critical decisions concerning the evolution of a patient’s lesion. To do so, two of our images will include a multiple sclerosis lesion. The influence of these two factors (noise intensity and the presence or absence of a lesion) will allow us to design a set of experiments with different complexity levels. Note that these decisions on the introduction of noise and the consideration of a lesion are related to the capabilities provided by the BrainWeb database where only some choices can be performed, i.e. we are constrained to adding only Gaussian noise to our testing images, although other noise models would be more suitable when dealing with MRIs. Nevertheless, as stated by Yetik and Nehorai (2006): “It is not possible to model the noise as an additive Gaussian model when different techniques are used to acquire images. In certain cases, e.g., in different modalities of medical imaging, the images are related by a nonlinear and nonrandom function. Therefore, the results we obtain here are valid for the same modality where additive Gaussian noise is a reasonable assumption.” Images from the same modality are also used in our contribution and the Gaussian noise model is a reasonable assumption as well. A pre-processing step has been carried out on all these 3D images to obtain problem-dependent information to guide the IR process, as well as to reduce the huge amount of data stored in the initial instances of the images (see Section 3.1 in the paper). Therefore, from every original image, we extract the isosurface and select crestlines points with relevant curvature information. The first image (“I1”, see Figure 2) corresponds to an MRI of a healthy person obtained with an ideal scanner, i.e. no lesion is present and it is a noise free scenario.

4

583 points have been selected after performing both the isosurface extraction to identify the brain and the study of the crest line points to choose those features with relevant curvature information

a)

b)

c)

Figure 2: Image I1. a) Original MRI with Three Views (Transverse, Sagital, and Coronal). Different Organs (Skull, Brain, Eyes, etc.) can be Clearly Identified. b) Isosurface Corresponding to the Brain from the MRI. c) Crest Line Points with Relevant Curvature Information. Image “I2” (Figure 3) corresponds to a low level of noise scenario (1% of Gaussian noise) of a healthy person. After the isosurface extraction, 393 crest line points have been chosen.

a)

b)

c)

Figure 3: Image I2. a) MRI with a 1% of Gaussian Noise. b) Isosurface Corresponding to the Brain from the MRI. c) Crest Line Points with Relevant Curvature Information. The third image (“I3”, see Figure 4) includes a multiple sclerosis lesion and the same level of noise than I2. As in the previous images, we extract 348 points with relevant curvature information.

5

a)

b)

c)

Figure 4: Image I3. a) MRI with a 1% of Gaussian Noise. The Sclerosis Multiple Lesion is Located using a Circle. b) Isosurface Corresponding to the Brain from the MRI. c) Crest Line Points with Relevant Curvature Information. Image “I4” (Figure 5) corresponds to a multiple sclerosis patient but the MRI has been acquired using a poor device (5% of Gaussian noise is introduced). Blurring can be easily observed in the extracted isosurface. Finally, 284 points with relevant curvature have been identified.

a

b

c

Figure 5: Image I4. a) Original MRI with a 5% of Gaussian Noise. Multiple Sclerosis Lesion is Located using a Circle. b) Isosurface Corresponding to the Brain from the MRI. Blurring is Easily Identified. c) Crest Line Points with Relevant Curvature Information. b) The CT dataset of human wrists The Bioengineering Laboratory of the Brown Medical School/Rhode Island Hospital has provided us with an in-vivo dataset. Specifically, it is composed of two real-world CT scans (3D images) of human wrists acquired from the same patient at different points in time, thus becoming an intra-object registration problem.

6

We want to underline the complexity of the new IR problems to be tackled. The difficulty is due to different reasons: the images are not realistic but real-world images (that do not come from a simulator but from a real CT scanner), the bones composing the wrist are prone to performing a local movement that does not fit into the typical profile of a similarity transformation (unlike the brain, they are not protected by the skull and the registration might not be perfectly achieved by a similarity transformation), and the degree of occlusion between the images is higher than in the previously described dataset. With the same goal as in the MRI dataset, a preprocessing step is performed for the two CT images of this dataset to extract the isosurface and select their crest-lines points. Figures 6 and 7 show the pipelining of this image-preprocessing step. After the isosurface extraction to identify the bones of the wrist and the study of the crest line points to choose those features with relevant curvature information, 575 and 412 points have been selected, respectively.

a

b

c

Figure 6: Image I5. a) Original CT Image. b) Isosurface Corresponding to the First CT Scan of the Patient Wrist. c) 575 Crest Line Points Extracted.

a

B

c

Figure 7: Image I6. a) Original CT Image. b) Isosurface Corresponding to the Second CT Scan of the Patient Wrist. c) 412 Crest Line Points Extracted.

7

3. Parameter Settings in the SS algorithm

In this section we present our preliminary study on the more suitable parameter values for the different IR algorithms. Both the preliminary test and the subsequent experimentation have been made on a platform with an Intel Pentium IV 2.6 MHz processor and all the algorithms have been implemented in C++ and compiled with GNU/g++. The parameters considered for each method follow: ? The I-ICP algorithm is run for a maximum of 40 iterations with the parameter value used in Liu (2004). The initial conditions (rotation, translation and scaling) are estimated as in Horn (1987). ? The ICP+SA algorithm is run for one complete iteration with a maximum of 40 iterations for the wrapped I-ICP algorithm. The annealing process has 20 iterations and 50 trial movements around each annealing iteration, with an initial temperature value estimated with T0 = [?/-ln(φ)]C(S0), where C(S0) is the cost of the given solution generated by the previous run of I-ICP, and both the ? and the φ factors take value 0.3. ? In the dynamic genetic algorithm (Dyn-GA), the size of the initial population has been established at 100 individuals and the remainder of the specific parameters have been kept with their original values (Chow et al. 2004). ? Both SS proposals (the previous one in Cordón et al. 2004, SS, and the current one, SS*) deal with an initial set P comprising 80 diverse solutions, and a RefSet with b1=7 quality solutions and b2=3 diverse solutions. The local search algorithm is run a maximum of 80 iterations at each execution, updating the δ matching vector every k=10 local search iterations (see Section 3.3 in the paper). In this way, δ will be updated seven times during each local search run. ? For both the Dyn-GA and the two SS algorithms, we have established a maximum CPU time of 20 seconds at each run. Furthermore, for each one of the latter algorithms as well as for the ICP+SA one, we have performed a total of 15 runs (with different random seeds) for each of the twenty problem instances (sixteen related to the BrainWeb MRI dataset and four corresponding to the CT wrist one) in order to avoid the usual random bias of probabilistic algorithms.

8

References

Chow, C.K., H.T. Tsui, T. Lee. 2004. Surface registration using a dynamic genetic algorithm. Pattern Recognition 37 105–117. Collins, D.L., A.P. Zijdenbos, V. Kollkian, J.G. Sled, N.J. Kabani, C.J. Holmes, A.C. Evans. 1998. Design and construction of a realistic digital brain phantom. IEEE Transactions on Medical Imaging 17 463–468. Cordón, O., S. Damas, J. Santamaría. 2004. A scatter search algorithm for the 3D image registration problem. X. Yao, E.K. Burke, J.A. Lozano, J. Smith, J.J. Merelo Guervós, J.A. Bullinaria, J.E. Rowe, P. Ti?o, A. Kabán, H.P. Schwefel, eds. Parallel Problem Solving from Nature - PPSN VIII, 8th International Conference, Birmingham, UK, September 18-22. Lecture Notes in Computer Science 3242, Springer, 471–480. Held, M., W. Weiser, F. Wilhelmst?tter. 2004. Fully automatic elastic registration of MR images with statistical feature extraction. Journal of WSCG 12 153–160. Horn, B.K.P. 1987. Closed-form solution of absolute orientation using unit quaternions. Journal of the Optical Society of America A 4 629–642. Kwan, R.K.S., A.C. Evans, G.B., Pike. 1999. MRI simulation-based evaluation of image-processing and classification methods. IEEE Transactions on Medical Imaging 18 1085–1097. Liu, Y. 2004. Improving ICP with easy implementation for free-form surface matching. Pattern Recognition 37 211–226. Rogelj, P., S. Kova?i?. 2002. Validation of a Non-Rigid Registration Algorithm for Multi-Modal Data. M. Sonka, J.M. Fitzpatrick, eds. SPIE in Medical Imaging 4684 299–307. Wachowiak, M.P., R. Smolíkova, Y. Zheng, J.M. Zurada, A.S. Elmaghraby. 2004. An Approach to Multimodal Biomedical Image Registration Utilizing Particle Swarm Optimization, IEEE Transactions on Evolutionary Computation 8:3 289–301. Yetik, I.S., A. Nehorai. 2006. Performance Bounds on Image Registration. IEEE Transactions on Signal Processing 54 1737–1749.

9