An Overview on the Latest Nature-Inspired and Metaheuristics-Based Image Registration Algorithms

: The development of automated image registration (IR) methods is a well-known issue within the computer vision (CV) ﬁeld and it has been largely addressed from multiple viewpoints. IR has been applied to a high number of real-world scenarios ranging from remote sensing to medical imaging, artiﬁcial vision, and computer-aided design. In the last two decades, there has been an outstanding interest in the application of new optimization approaches for dealing with the main drawbacks present in the early IR methods, e.g., the Iterative Closest Point (ICP) algorithm. In particular, nature-inspired computation, e.g., evolutionary computation (EC), provides computational models that have their origin in evolution theories of nature. Moreover, other general purpose algorithms known as metaheuristics are also considered in this category of methods. Both nature-inspired and metaheuristic algorithms have been extensively adopted for tackling the IR problem, thus becoming a reliable alternative for optimization purposes. In this contribution, we aim to perform a comprehensive overview of the last decade (2009–2019) regarding the successful usage of this family of optimization approaches when facing the IR problem. Speciﬁcally, twenty-four methods (around 16 percent) of more than one hundred and ﬁfty different contributions in the state-of-the-art have been selected. Several enhancements have been accordingly provided based on the promising outcomes shown by speciﬁc algorithmic designs. Finally, our research has shown that the ﬁeld of nature-inspired and metaheuristic algorithms has increased its interest in the last decade to address the IR problem, and it has been highlighted that there is still room for improvement.


Introduction
It has been proven that Image registration (IR) [1] plays a vital role within the computer vision (CV) field. Specifically, IR is one of the most important CV tasks within image processing. It is concerned with the process of accurately aligning two or more images (taken at different times, from different sensors or viewpoints) into a common coordinate system. Then, IR is aimed at finding either a geometric transformation or a correspondence (matching of feature image points) which provides the best possible overlapping of the considered images. It has been applied in a huge amount of real-world applications, medical imaging and remote sensing bwing two of those most addressed in the State-of-The-Art (SoTA).
Over the years, different algorithms facing the IR problem have been proposed resulting in a large body of research [2][3][4]. To sum up, the IR problem can be formulated from an optimization view-point, in which the IR procedure involves finding the optimal transformation achieving the best fit between the considered images. Such estimation/optimization of the IR transformation is typically addressed by means of an iterative process properly searching in the space of promising IR solutions.
In the early contributions to IR, the optimization process carried out by these methods was negatively influenced by the image noise, the image discretization, and orders of magnitude in the scale of the parameters of the IR transformation, among other factors. In particular, it was the case of the canonical Iterative Closest Point (ICP) algorithm [5], in which the process was highly prone to get trapped in local optima IR solutions [6]. Moreover, the latter pitfall is due to that kind of IR method assuming that the alignment procedure starts from accurate pre-alignments of images usually provided by experts.
Soft computing (SC) [7] is a term applied to a field within computer science which is characterized by the use of inexact solutions to computationally hard tasks such as the NP-complete problems, in which a global optima solution cannot be derived in (desirable) polynomial time. In particular, SC differs from conventional (hard) computing in that it is tolerant regarding imprecision, uncertainty, partial truth, and approximation. Then, the guiding principle of SC is exploiting the tolerance for such an imprecision, uncertainty, partial truth, and approximation in order to achieve tractability, robustness and low solution cost. The main SC paradigms include fuzzy systems, nature-inspired computation, and artificial neural computing, among others.
In the last two decades, nature-inspired computation, e.g., evolutionary computation (EC) [8,9], has demonstrated its ability to deal with complex real-world problems within the field of CV. Additionally, metaheuristics are other general purpose optimization algorithms successfully applied in this field. As an example, several special issues and books on the topic have been published in international forums in the last few years [10][11][12].
Specifically, both the nature-inspired and the metaheuristic (NI&M) approaches have been extensively applied to tackle IR problems without requiring a good initial estimation of the alignment of the images. That advantage is mainly motivated by the global optimization nature of these optimization strategies, which allows them to perform a robust search in complex and ill-defined search spaces.
Since the first attempts facing the IR problem using NI&M, this topic has become a very active niche of research and a significant number of proposals have been contributed. Nevertheless, an updated revision of the SoTA is needed due to new and novel enhancements being proposed in the last decade. In this work, we focused on introducing an extensive survey of those most relevant NI&M-based IR methods in the SoTA. Subsequently, it allowed us to provide a more comprehensive analysis of the inherent nature of these new IR methods, in which several quality parameters have been considered for their proper assesment. Moreover, some suggestions of improvement have been accordingly introduced as concluding remarks.
The structure of this contribution is as follows. Section 2 describes the IR problem and introduces the most important strengths and pitfalls of the early IR methods. Next, Section 3 describes the key concepts of NI&M and presents a quantitative analysis of the increasing interest of the community in the application of NI&M to IR. Section 4 is devoted to carrying out an extensive overview of the current SoTA dealing with these methods. Finally, Section 5 provides a comprehensive analysis of the overviewed methods.

Image Registration
It has been proven that there is not a universal design for a hypothetical IR method that could be applicable to every real-world application [2]. However, IR methods consist of the following four components: • Two input images named scene I s = { p 1 , p 2 , . . . , p n } and model I m = { p 1 , p 2 , . . . , p m }, with p i and p j being image points.

•
A registration transformation, named f . It corresponds to a geometric/parametric function that relates the coordinate systems of both images.
• A similarity metric function, named F. It is aimed at assessing the level of resemblance (i.e., the degree of overlapping) between the transformed scene image (i.e., f (I s )) and the model one.

•
An optimization procedure. The optimizer regards a method seeking the optimal registration transformation ( f ). Likewise, it is necessary to provide a search space for the suitable representation of the IR solutions.
As stated, due to the non-linear nature of the inherent optimization procedure of the IR problem, an iterative procedure is often followed. As depicted in Figure 1, such a procedure usually stops once convergence is achieved, specifically, after the F value (i.e., the similarity metric) is lower than a given tolerance threshold. Then, we focused our attention on the optimization task, because of its vital importance in the success of the whole performance of the IR method. Moreover, two different search space schemes were adopted in the SoTA: (i) The matching-based approach, in which the optimization process is aimed at looking for an optimal set of pairs of similar image points; and (ii) the parameter-based one, which directly explores the search space of the values of IR solutions, i.e., f .
The early proposals tackling the IR problem were based on the first approach (i.e., the matching-based). This is the case of the well-known ICP method [5]. Its main pitfall is its sensitiveness to the initial IR transformation ( f ) given by the user [13]. Then, ICP is prone to get trapped in local optima solutions. A detailed description of the IR framework is out of the scope of this contribution. We refer the interested reader to [2]. Likewise, the formulation of the IR problem is dependent on the particular environment in which it is involved (remote sensing, medical imaging, CAD, etc.).
In our case, the 3D model reconstruction application can be considered in order to formulate the IR problem. Range scanners are 3D image acquisition devices able to capture images of the surface of the sensed object (see Figure 2), named range images. Each range image is acquired by considering a specific view-point different of the previously sensed ones. The full set of range images models the whole geometry of the scanned object. Nevertheless, each of them is in a different coordinate system. Thus, the reconstruction procedure carries out the accurate integration of all the 3D images in order to achieve a complete and reliable model of the physical object. This framework is commonly named 3D modeling (see Figure 2) and it concerns the application of, first, successive pair-wise IR steps of adjacent range images. The final stage considers all the range images at once in order to achieve the optimal integration of the whole set ( [14]). Notice how the IR problem is present in two stages of the reconstruction procedure, i.e., in the pair-wise alignment and in the final integration of the whole set of range images. Thus, the IR problem can be formulated as follows. Every pair-wise IR step is aimed at finding the optimal Euclidean transformation ( f ) that brings the scene view (I s ) into the best possible overlapping onto the model view (I m ). Usually, it is considered a 3D rigid transformation composed by seven real-coded parameters: A rotation R = (θ, Axis x , Axis y , Axis z ) and a translation t = (t x , t y , t z ), with θ and Axis being the angle and axis of rotation, respectively. Then, the transformed image points of the scene view are noted by Then, the main goal of the Optimizer is to search for the optimal IR transformation f * which provides the best alignment according to the Similarity metric F(): Specifically, the median square error (MedSE) is the Similarity metric considered in 3D modeling [4]: where MedSE() corresponds to the median d i value. We define d i = f ( p i ) − q cl 2 as the squared Euclidean distance between the transformed scene point, f ( p i ), and its corresponding closest point, q cl , in the model view I m .
Moreover, it has been used ( [3,15]) indexing structures such as the kd-tree [16] and the grid closest point (GCP [17]) in order to speed-up the computation of the closest point q cl of I m .

Nature-Inspired and Metaheuritcs-Based Image Registration
As stated, NI&M is a specific computational approach within SC which has been extensively applied in addressing many real-world optimization problems as those within CV. The reader can find an extensive overview on the topic in references [11,12,18].
Briefly, it is necessary to describe some key concepts of NI&M to achieve a suitable understanding of the inherent foundations of such an optimization approach. Specifically, EC proposes models of computation inspired in biological processes to evolve populations of solutions, thus providing key elements in the development of computer-based problem solving systems. On the other hand, Metaheuristics [19,20] is a more general family of enhanced search optimization algorithms and, as EC, they constitute a very compelling alternative due to their capability to obtain good quality solutions when global optima solutions cannot be found within a reasonable time of computation.
In particular, it has been proposed a variety of EC models referred as evolutionary algorithms (EAs) [21]. EAs is a family of Nature-inspired algorithms widely used for solving complex optimization problems. Among them, GAs are probably the most used EAs in the literature to face real-world optimization problems. Some other EAs have been proposed in the last few years improving the SoTA of this field by adopting more suitable optimization strategies: Evolution strategies (ES) [22], scatter search (SS) [23], differential evolution (DE) [24,25], memetic algorithms (MAs) [26], particle swarm optimization (PSO) [27,28], estimation distribution algorithms (EDAs) [29], membrane computing (MC) [30], and cellular automata (CA) [31], among others. Moreover, other EC models inpired by similar principles have also been recently proposed, e.g., bacterial foraging optimization algorithm (BFOA) [32] and artificial bee colony (ABC) [33].
Specifically, the application of NI&M to the IR optimization process has caused an outstanding interest in the last decade. Unlike traditional ICP-based IR methods, those supported by NI&M need neither rough nor near-optimal pre-alignment of the images to proceed. The first attempt to tackle IR using NI&M can be found in the eighties, in which a binary-coded canonical GA was proposed for addressing medical scenarios of 2D angiography images [34]. Figure 3 depicts the evolution of the interest of the scientific community in the application of NI&M to IR in the last decade (The bar-graph in Figure 3 was directly obtained from Thomson Reuter's Web of Science using the query TI = (image AND (registration OR alignment OR matching) AND (evolution* OR swarm OR optimization)).), in which more than one hundred and fifty papers appeared.

Revision of the State-of-The-Art
As shown in Figure 3, an important number of NI&M-based IR methods have been proposed in the last decade. Nevertheless, an updated revision of the SoTA is needed due to new and novel enhancements having been contributed, and none of them having been reviewed in previous surveys on the topic [3,35].
In this contribution, our revision of the SoTA regarding the application of this kind of optimization strategy is based on the selection of twenty-four papers coming from more than one hundred and fifty papers plotted in Figure 3. The considered selection criteria was three-fold and inspired by the following: • Novelty: Does the contribution make use of a new NI&M-based optimization algorithm which has not been tested for the IR problem yet? • Technology: Is parallel computing taken into account to improve optimization in any way? • Continuity: Does the contribution provide a new advancement based on any previous research?
Next, we analyze in deep detail and chronologically the resultant twenty-four selected contributions.

Santamaría et al.'s Memetic-Based Proposal
This contribution [36] deals with real-world forensic tasks by using 3D models of forensic objects, e.g., skulls, bones, corpses, etc. The 3D data was acquired by means of 3D range scanners (see Figure 2) and the optimization challenge addressed by the authors was to provide a near optimal solution to the pair-wise IR problem. Regarding the nature-inspired approach adopted in this contribution, the authors provide a new way of facing the problem, usually centred in a sequential hybridization of optimization techniques. In their work, they proposed an opposite scheme based on the memetic (MA) framework [26], and aimed at analysing the optimization capabilities of MAs for tackling really complex and challenging 3D reconstructions of real-world forensic objects. The contribution provides an extensive setup of different MAs adapted to face this specific IR scenario. The conducted experiments demonstrated the better performance of these new methods, i.e., MAs, when compared with others in the SoTA.

Queirolo et al.'s SA-Based Proposal
This work [37] extends their previous research on using nature-inspired algorithms for addressing the IR problem [38]. In this case, authors focused on dealing with complex 3D face recognition instances using the FRGC v2 database simulating both verification and identification systems. They used a simulated annealing (SA) algorithm for tackling optimization tasks present in face recognition. Moreover, as in [38] they newly made use of the Surface Interpenetration Measure (SIM) for the corresponding assessment of the quality of the matching of pairs of range images. From the conducted experiments, their method achieved outstanding results, at a False Acceptance Rate (FAR) of 0.1% and, regarding identification scenarios, they obtained a rank-one accuracy of 99.5%.

Rusu and Birmanns' GA+TS-Based Proposal
In this case, authors provided a more elaborated design based on the GAs and the Tabu Search (TS) algorithms [39]. Then, their novel IR method deals with the simultaneous registration of multiple atomic structures into cryo-EM envelopes, named MOSAEC (Multi-Object Simultaneous Alignment by Evolutionary Computing). In MOSAEC, authors adapted the classic scheme of GAs to enhance the exploration of the search space. In particular, new genetic operators were introduced to preserve the genetic diversity of the population and it was used in combination with parallel evolution of subpopulations. Moreover, the exploration of the complex search space was improved by including tabu regions, i.e., areas of the search space which are marked as local optima and thereby should not be further sampled. In the conducted experiments, the performance of MOSAEC was assessed by carrying out simultaneous IRs of the bacterial ribosome and of the chaperonin GroEL. The reported results indicated that their proposal was able to simultaneously register multiple component structures, and identifying their placement and orientation with accuracies within one order of magnitude of the nominal resolution of the cryo-EM maps. Finally, the optimization procedure was enhanced with parallel computing strategies accompanied by horizontal gene transfer.

Maia et al.'s EvSOM-Based Proposal
In this work, a variant of the Evolutionary Self-Organizing Map (EvSOM) algorithm was proposed in order to fully exploit the neighbourhood preservation property of the topographic maps induced by it when tackling the IR problem [40], in particular, using affine transformations. EvSO is a well-known robust approach for topographic map formation that is based on EAs [41], which has proved to be strongly resilient to outliers, less sensitive to control parameter selection, and less susceptible to the effects of multimodality and local optima. Preliminary experiments with black & white retinal blood vessel images from the DRIVE repository were discussed and their IR proposal was compared against two well-known IR methods: ICP and template matching. Despite the reduced experimentation carried out, EvSOM showed promising results when dealing with near-complex 2D IR scenarios.

Das and Bhattacharya's Proposal
This work [42] presents a non-linear 2D/2D affine IR technique for MR and CT modality images of sections of human brains, in which a Correlation-based similarity metric is considered for the assessment of the results. The authors focused their study on the comparison of the performance of two nature-inspired optimization techniques, GAs and PSO, for the maximization of the considered similarity. Moreover, the proposed study considered a multiresolution strategy for tackling the IR problem in order to decrease the sensitivity of the IR optimization algorithm to local optima, in the sense that the images were downscaled to a number of resolution levels. The conducted experiments showed that the IR methods based on PSO achieved better performance than those using GAs, thus becoming the more suitable optimization approach facing this specific IR scenario.

Santamaría et al.'s GRASP&PR-Based Proposal
In this paper [43], the authors focused on the point matching-based approach of the IR problem using 3D medical images. They contributed with improved hybridizations of two well-known metaheuristics, named greedy randomized adaptive search procedure (GRASP) and PathRelinking (PR) [44]. For instance, for the static GRASP and PR (noted by Stc-G&PR), the hybridization variant is an static one in which, first, GRASP is applied to build an elite set (ES) of solutions (see steps 1 to 14 in Figure 4 [43]) and next, as a second step, PR is performed in order to generate solutions among all the pairs of solutions in ES. Thus, the optimization procedure aimed at finding the best matching of pairs of feature 3D points and the corresponding registration transformation is numerically derived by means of least squares methods. The best the matching, the more accurate the numerical estimation. Their hybrid proposals were based on specific designs focused on achieving a good trade-off between the intensification and the diversification components of the search space in order to obtain high quality IR solutions. Finally, the authors tested their new designs addressing real-world problems of medical imaging. In particular, they considered six medical images from two different image datasets for testing purposes. The first dataset consists of MR images obtained from the well-known database at McGill University, named BrainWeb [45], and the second dataset provides CT images from the Rhode Island Hospital [46]. Several IR methods of the SoTA were considered for computational comparison against their proposals. The reported results demonstrated the better trade-off between intensification and diversification offered by these new optimization designs when facing the point-matching-based IR problem.

Yang et al.'s SDE-Based Proposal
This contribution [47] proposed a self-adaptive differential evolution (SDE) algorithm for 2D multi-modal IR using normalized mutual information (NMI) as the similarity metric, cubic B-splines for image voxel interpolation and image down-sampling for image multiresolution. Usually, the control parameters F and CR are hand-matched to the problem, but for the SDE algorithm [48] they adapt during successive iterations. For instance, for the SDE/rand/1/bin algorithm, F j and CR j are vectors additionally coded in the solution together with the transformation parameters. Then, the searching space is augmented with both two control parameters. From the reported experimental results, authors concluded that the IR transformation could be obtained significantly faster than if no down-sampled images are used for seeding their SDE method.

Santamaría et al.'s MA+AIS-Based Proposal
This work [49] presented a new self-adaptive evolution model to deal with IR problems facing pair-wise range IR problem instances. The specific design of their self-adaptive evolutionary optimization method, named SaEvO, takes advantage of the synergy between two different EAs: A memetic algorithm (MA) [26] based on differential evolution (DE) [25] and the variable neighbourhood search (VNS) [50], and combined with an artificial immune system (AIS) algorithm [51]. In similar manner to DE, VNS is easy to implement and it requires few control parameters. Moreover, VNS extends the capabilities of DE by means of carrying out a local search step. While the optimization stage of MA (first stage of SaEvO) is aimed at searching for IR solutions, the latter stage (based on AIS) is focused on tuning the control parameters, i.e., for both the DE and the VNS algorithms. The conducted experiments considered several range image datasets obtained from the well-known public repositories, e.g., the Signal Analysis and Machine Perception Laboratory (SAMPL), and others from a Konica Minolta c Laser Range Scanner. Their conducted experiments reported outstanding results when comparing SaEvO against several IR algorithms in the SoTA.

Alderliesten et al.'s EDA-Based Proposal
In [52], the authors developed, by considering a previously introduced multi-objective optimization framework, a dual-dynamic transformation model to be able to tackle prone-supine deformable IR scenarios using 2D axial slices of breast MRI scans. According to optimization purposes, they used a specific variant of multi-objective EA based on EDAs, which aims to exploit features of the structure of the problem by means of probabilistic modeling [29]. The method is known as iMAMaLGaM-X+ (i.e., incremental Multi-objective Adapted Maximum-Likelihood Gaussian Model miXture) and it makes use of a population of solutions. In particular, iMAMaLGaM-X selects 35% of the best solutions according to a domination-rank ordering, next it estimates an l-dimensional normal mixture distribution from these selected solutions (where l is the number of real-valued variables to be optimized), and generates new solutions by sampling the estimated distribution. Their method reduces the risk of premature convergence by means of using adaptive techniques for adjusting the parameters of covariance matrices of the normal distributions. Their proposed variant is capable of obtaining an ever better spread of solutions by maintaining m additional components in the mixture distribution, one for each objective. The selection of these components is done independently on the basis of each respective individual objective, thereby specifically targeting convergence at the extreme regions of the Pareto. As the authors suggested, their preliminary results indicated that the outlined methodology is also highly likely capable of paving the road to an elegant solution to some of the hardest deformable IR problems.

Bermejo et al.'s BFOA-Based proposal
In this contribution [53], the authors faced the IR problem for new 3D image acquisition devices based on the time-of-flight technology, which obtains range images from real-time 3D video sequences. They contributed with a broad study of the performance of variants of IR based on the bacterial foraging optimization algorithm (BFOA) [32]. Briefly, BFOA is based on allowing cells to stochastically and collectively swarm toward optima. This is achieved through a series of three processes on a population of simulated bacteria: (i) Chemotaxis, where the cost of bacteria is reduced by the proximity to other bacteria, and they move along the manipulated cost surface one at a time; (ii) reproduction, in which only those bacteria that performed well over their lifetime may contribute to the next generation; and (iii) elimination-dispersal, where each individual/bacteria is discarded and new random samples are inserted with a low probability. In particular, this contribution added new components to the previous recipe. For instance, they introduced a local search (LS) strategy as a hybridization with the Dasgupta et al.'s BFOA variant. Specifically, they considered the crossover-based LS (XLS) technique which obtained promising results in their previous works. XLS is applied after the reproduction step.
In the experimental setting, they used range images obtained from video sequences acquired from the PMD[Vision] c Camcube 2.0 camera. The reported results revealed the suitability of BFOA to address this kind of IR problem.

Ma et al.'s OLDE-Based Proposal
Ma et al. [54] proposed an intensity-based IR method for remote sensing. This IR variant made use of the hybridization of both the orthogonal learning (OL) and the DE algorithms (OLDE). It has been demonstrated that the OL strategy provides an efficient alternative when searching for in complex problem spaces. In OLDE, the OL strategy guides the DE algorithm to select promising search directions towards the global optimum. In particular, they provide a crossover component for DE based on OL, named orthogonal crossover. They tested the performance of OLDE against the one from other EA-based IR methods. The experimental section considered 2D images of RADARSAT SAR sensors. The reported results showed that OLDE was the best IR algorithm which achieved robust and efficient outcomes with respect to the initial values of the IR transformation.

De Falco et al.'s AIM-dDE-Based Proposal
De Falco et al. [55] proposed an Adaptive Invasion-based Model (AIM) to tackle range IR problems. Their AIM-dDE method made use of a distributed DE (dDE) algorithm which considers a migration model. The latter is inspired in the natural phenomenon known as biological invasion. Briefly, at each generation t, the subpopulation P p (t) of each node p performs a sequential DE until t equals t max generations, and at every T generations, neighbouring subpopulations exchange individuals. Then, at a given invasion time, t = kT, the pool of individuals in each subpopulation P p (t) sends to its propagule (i.e., neighbours) according to M P p (t) and it is carried out by gathering the individuals of P p (t) which are better than its average fitness. Next, a founding subpopulation formed by both native and exogenous individuals is built in each node by means of adding this huge number of invading individuals to the local subpopulation. In this manner, the authors provided a novel mechanism based on the founding subpopulation as a source of diversity exploitable by their algorithm to properly enhance its capability to search for new niches of solutions. Moreover, their AIM-dDE made use of an adaptive strategy for searching for the optimum control parameters. Specifically, they proposed two updating schemes: RandAvg and ChAvg. In the experimental section, they collected several range image datasets from a well-known source in the field, i.e., the SAMPL at the Ohio State University. They provided parallel versions of the proposed algorithms by making use of the Message Passing Interface (MPI) paradigm. Their proposed range IR method was compared against several algorithms in the SoTA. The reported statistical results suggested that the AIM-dDE-ChAvg variant turn out to be the best one in all the three considered tests.

Pirpinia et al.'s HybGA-Based Proposal
In Pirpinia et al. [56], the authors faced a hard optimization problem based on a multi-objective approach for tackling deformable image registration (DIR) problems. In this multi-objective DIR algorithm, which is a version of gradient descent method, the set of all pareto-optimal improving directions for a vector function f is initially calculated. The directions that are of interest, i.e., that lead to maximum improvement of all objectives, correspond to the negative non-dominated directional derivatives. One of these promising directions is randomly chosen and a multi-objective line search is performed along that direction. This scheme is named Combined Objectives Repeated Line-search (CORL). The authors provided several hybridizations based on the CORL and the Genetic Local Search (GLS) algorithms. Specifically, the Adaptive GLS variant integrates the CORL method in an adaptive hybridization scheme taken from the literature, which exploits gradient information in three different ways. For achieving maximum efficiency, AGLS decides during the run which exploitation method to use at the end of each evolutionary cycle based on the number of improved solutions found. Moreover, instead of applying conjugate gradients to a randomly chosen objective for a solution, AGLS applies the procedure to only members of the population following again a rank-based scheme, i.e., RL-CORL. Therefore, conjugate gradient is accordingly applied to the solution that ranks best in one objective to further improve it. The proposed multi-objective DIR approach was tested in medical IR scenarios considering 2D MRIs. The reported results showed the improvement of their proposal over the existing adaptive hybridization schemes. Actually, the unbiased improvement of all objectives in a multi-objective gradient technique as part of a hybridization scheme was indeed successful at addressing the tested hard DIR cases, thus providing the best balance between exploitation and exploration capabilities.

Bermejo et al.'s BFOA-Based Proposal
In this contribution [57], the authors carried out an extended research of their previous work [53] using BFOA and mainly aimed at taking a step ahead in the development of the discipline by introducing several enhanced variants of this nature-inspired algorithm. In this case, they carried out the combination of some advanced BFOA designs proposed in [58]. Among this variants, they proposed an alternative solution for enhancing the social behaviour of the bacteria simulated in the original swarming stage. Moreover, their new designs are aimed to ease the negative effects of the tumble step by means of the promotion of a bacterium following the progress of another successful bacterium in the population, thus mimicking its direction. Besides, it was reported that the usage of the best bacterium found so far to guide the dispersal step of BFOA allowed to minimize the negative impact of killing the best bacterium proposed in the canonical approach. Another major refinement of their method concerns the inclusion of an elite set of the M = 5 best solutions, and after each chemotactic step the recombination is applied by using the BLX − α crossover. The proposed variants were mainly focused on addressing a lack of exploitative capability present in their previous canonical BFOA. In the experimental setup, the authors tackled a new IR scenario from medical imaging by means of considering 3D images of feature points extracted from MRI datasets of human brains. Specifically, fifty-two different BFOA variants were tested and accordingly compared against several of the best performing EA-based IR methods in the SoTA. Many of the proposed variants obtained better results than those methods in the field, which proved that their new designs achieved improved and more robust IR outcomes.

Yang et al.'s HLCSO-Based Proposal
In [59], the authors tackled non-rigid multi-modal IR scenarios using the NMI metric and proposed a hybrid method that combines both the L-BFGS-B and the cat swarm optimization (CSO) algorithms, named HLCSO. On the one hand, CSO makes use of cooperative coevolving and block grouping as the grouping strategy to capture the interdependency among variables. On the other hand, L-BFGS-B provides a way of achieving faster convergence and higher accuracy of the final solution. Moreover, the roulette wheel method is introduced into both the seeking and the tracing modes of HLCSO in order to improve the performance of the algorithm. Extensive experiments on 3D CT, PET, T1, T2 and PD weighted MR images demonstrated the superiority of the proposed hybrid approach, HLCSO, against separately using the L-BFGS-B and the CSO methods.

Li et al.'s HFFO-Based Proposal
Li et al. [60] contributed to the field with an IR algorithm using 3D range images and the fruit fly optimization (FFO) algorithm. Moreover, they proposed the computation of the normal angle histogram (NAH), thus the maximum NAH can be used as a criterion to direct the step size of the FFO-based optimization algorithm. In FFO, when searching for the food, fruit flies smell all kinds of scents in the air based on their osphresis organs. Then, the route is obtained with the location of food. When getting close to the food, the visions are used to accurately discriminate the food. FFO has a simple structure and few control parameters to tune. Specifically, HFFO self-adaptively adjust the step size of fruit fly based on NAH in order to reduce the search space. In the experimental setup, the authors considered several 3D modeling scenarios for testing their proposal. The HFFO-based IR algorithm was compared against other methods of the field, in particular, ICP and a variant of FFO, named improved FFO (IFFO), in which a new parameter was introduced in FFO for dynamically changing the search radius along each iteration number. The reported results showed that HFFO played the best performance in all considered scenarios. Then, the authors demonstrated the suitability of considering the NAH information for better guiding the exploration capabilities of FFO.

Qin et al.'s ABChDE-Based Proposal
In [61], the authors contributed with an hybridization of an improved version of the ABC and the DE algorithms (ABChDE) applied to remote sensing of 2D images. In the canonical version of the ABC algorithm, onlooker bees use the roulette-wheel selection scheme to select a nectar source, and they explore around a vicinity. A probabilistic scheme was considered by ABChDE for the onlooker step, in which either the canonical version or the one based on DE can be considered. At the initial stage of ABChDE, onlooker bees have a large probability to choose the updating equation of DE, which makes them search in a wide rough space, thus promoting scape from local optima. Several 2D image datasets were considered in the experimental section and the canonical versions of ABC, DE, and PSO algorithms were accordingly tested against their proposal. ABChDE showed a stronger global search ability than the traditional ABC algorithm and it became the best method compared with the other algorithms.

De Falco et al.'s AsAMP-dDE-Based Proposal
In this work [15], the authors provided an extended version of their previous research [55]. Specifically, their improved IR algorithm is characterized by following an asynchronous migration mechanism and a multi-population recombination information exchange (AsAMP). Moreover, it was also supplied with adaptive updating schemes for automatically setting the control parameters of DE (dDE). The experimental setup considered several 3D range image datasets, and AsAMP-dDE was tested against some of the best EA-based IR methods in the SoTA. While the statistical results do not revealed a clear evidence of the pre-eminence of AsAMP-dDE over the other methods, their proposal exhibited a remarkable improvement in terms of both accuracy and CPU-time due to its distributed approach.

Costin et al.'s BFOA-Based Proposal
Costin et al. [62] contributed with the application of the canonical BFOA to PET-CT multi-modal and MR image rigid IR scenarios. The novelty of this work is the development of a parallel version of the BFOA algorithm. Several medical IR scenarios considering 2D rigid transformations were addressed. Both the Normalized Cross Correlation (NCC) and the NMI metrics were also used. The parallel version of their BFOA-based IR method was tested on a large set of DICOM medical images from an available database. Regarding the efficiency of the proposal, the parallel version of BFOA showed a greater speed-up of CPU-time computation.

Bouter et al.'s GOMEA-Based Proposal
In [63], the authors addressed the DIR problem by means of the proposal of a real-coded version of the recently contributed Gene-pool Optimal Mixing Evolutionary Algorithm (GOMEA). Their GOMEA variant tackled the IR problem from a multi-objective view-point. As in their previous contributions ( [52,56]), each IR solution in the Pareto set corresponds to a different trade-off between the objectives of interest. Specifically, the solution that is estimated the most suitable can manually be selected a-posteriori. The main difference between their new GOMEA-based IR method and the ones in the SoTA [56] was that the former has a much higher selection pressure as a result of increasingly enhancing sections of existing Parent solutions. In this new version, the authors separately proposed the estimation of such distribution for each subset of variables in the dependency model. Thus, the solutions are then independently improved by means of sampling new values for these subsets of variables in the linkage model. Then, only those sampled values that lead to an improved solution are accepted, otherwise, the main solution returns to its previous state. Because of the solutions are partially disturbed, it is not necessary to completely evaluate them to assess the contribution of a partial alteration. This strategy allowed their proposal to enhance its efficiency. The authors performed several experiments on three different IR scenarios: (i) An artificial problem with disappearing structures; (ii) a pair of pre-and post-operative breast CT scans; and (iii) a pair of breast MRI scans that were acquired in prone and supine positions. The reported results showed a speed-up of up to a factor of ∼1600 compared to the previously used algorithms.

Panda et al.'s ERBD-Based Proposal
In this contribution [64], the authors introduced a novel evolutionary rigid body docking (ERBD) algorithm for addressing medical IR problems. They considered different image data formats like MR, CTs and PET images of the brain of the same patient, all of them from the well-known Retrospective Image Registration Evaluation (RIRE) project. Regarding their proposal, i.e., ERBD, the docking process is successfully used in the design of drugs and it has never been used in IR. In particular, the docking task predicts the optimal configuration and energy between the protein and ligand. It changes the orientation of molecules and maximizes their contact. When the interaction energy between the protein and the ligand is minimized, they achieve best binding. Minimization of the energy also preserves the discontinuities in the low texture region. In their contribution, the authors used a canonical GA within ERBD to minimize the total energy. From the reported results, it is showed that the proposed method leads to higher IR accuracy since the interaction of both the energy and the MI were considered as the similarity metrics.

Li et al.'s PDE-Based Proposal
Li et al. [65] proposed a design model-based inspection method to address IR problems of 3D range images. Their method combined the DE algorithm and an EFT-based 3D point descriptor to improve the IR results. This point property-based DE algorithm was named PDE, in which the control parameters F i and CR i were set adaptively as proposed in [66]. Initially, the contributed algorithm randomly selects three point pairs n (n < n p ) times to form the candidate solution pool (CSP). The solutions in CSP are approximate and provide pre-registration for the two tested image models. In practice, with the increasing complexity of the geometry of the image the rate of false solutions becomes higher. Thus, PDE generates the population from two different sets. One set, X 1 , is randomly generated, and the other one, X 2 , is generated based on the CSP. Additionally, their method made use of the k best individuals of X 2 to replace the k worst individuals of X 1 as the current population. This strategy was aimed at speeding-up the convergence of the algorithm. Moreover, the mutation strategy proposed in [66] was also included in PDE. In the experimental section, the authors conducted several IR tests and it was considered both synthetic and range image datasets. PDE was compared against three deterministic (non-coarse) algorithms and other five heuristic algorithms. They validated the robustness of the EFT-based point descriptor, and the resultant outcomes demonstrated it suitability to achieve reliable IR results. Moreover, according to the analysis of the statistical results, in most of the cases, PDE showed significant better results compared to the other five heuristic-based IR algorithms. Finally, the authors concluded that PDE was insensitive (robust) to outliers, and PDE can obtain the global optima solutions successfully for all tested situations. However, the performance of PDE is dependent on the image content because the lower the similarity between the images, the more incorrect the correspondence of the image points.

Bermejo et al.'s CRO-SL-Based Proposal
In [67], the authors addressed medical IR scenarios by means of the Coral Reef Optimization (CRO) [68] algorithm with Substrate Layers (CRO-SL). The conducted experimental studies focused on both mono-modal and inter-modal medical IRs in which the scenarios belong to different kinds of 3D affine transformations. Besides, the authors adapted CRO-SL to tackle both the intensity and the feature-based IR approaches. The (CRO) algorithm is a recent contribution to the field of EAs which provides a global-search strategy based on the biology of the formation of coral reefs. Specifically, the algorithm is based on a reef of solutions to a given optimization problem (i.e., corals), that reproduce in the same manner as corals do in nature, by means of using the broadcast spawning, brooding and budding operators. In its canonical version, CRO lacks of exploitation capabilities. Then, the authors adapted CRO in order to address the IR problem by cleverly providing an enhanced version, named CRO-SL [68]. The new design of CRO-SL refers to a well-balanced combination of two crossover operators (i.e., BLX-α and SBX), two evolutionary-based operators (i.e., HS and DE), and one mutation operator (i.e., GM). The proposed IR algorithm was tested against two gradient-based IR algorithms and other tree heuristic-based methods in the SoTA. From the experimental results, the authors concluded that CRO-SL provided a suitable trade-off between exploration and exploitation capabilities when dealing with medical IR scenarios.

Cocianu and Stan's ES-APSO-Based Proposal
Last, in [69] the authors proposed a new IR method based on the hybridization of an self-adaptive ES algorithm and an accelerated PSO algorithm (ES-APSO). This recent contribution has been developed for signature recognition purposes. Specifically, their PSO strategy is based on the Firefly Optimization (FiFO) [70] algorithm. FiFO is based on the flashing behaviour of fireflies. In essence, it uses the following three idealized rules. First, fireflies are unisex, so that one firefly will be attracted to other fireflies regardless of their sex. The second rule states that the attractiveness is proportional to the brightness, and both fireflies decrease as their distance increases. The third rule concerns the brightness of a firefly, which is determined by the landscape of the objective function (i.e., F). ES-APSO was tested along a broad range of experiments on several 2D images representing signatures, and it was accordingly compared against several canonical versions of FiFO. The experimental results revealed the suitable application of ES-APSO for tackling signature recognition instances.

Analysis and Discussion
We can see how the number of publications facing the IR problem by means of NI&M keeps raising in the last decade with more than one hundred and fifty in total (see Figure 3). Nevertheless, what about the quality of such amount of contributions?
The overview carried out in this contribution can be considered a comprehensive one, in a way that several aspects covering broadly the quality of the proposal can be followed in order to measure this parameter. We schematized the (around 16 percent) twenty-four contributions revised in Table 1 according to relevant design parameters which can be considered for quality measurement, such as coding scheme of solutions (i.e., real, integer, or binary), search space (i.e., parameters vs. matching), and tuning of control parameters (i.e., hand-made vs. self-adaptive), among others.
Specifically, the first column in Table 1 identifies the publication year of the algorithm. Next, the second column refers to the coding scheme of solutions, i.e., R (real coding), B (binary coding), and I (integer coding). The third column concerns the IR approach followed (see Section 2): P (parameter-based) and M (matching-based). The NI&M-based optimization technique considered is shown in the fourth column. Moreover, the next two columns are devoted to revealing both those methods recently addressing the optimization problem from a multi-objective view-point and also considering self-tuning of the control parameters. According to the latter, S stands for those methods using a self-tuning approach, and M the canonical hand-made tuning one. Finally, the main application, the image modality, and the computer architecture (S and P stands for serial and parallel, respectively) are comprised in columns seven, eight, and nine, respectively. Figure 4 depicts the statistics regarding the level of impact caused by a subset of these quality parameters extracted from Table 1, in which QP#? Refers to the level of usage of the real-coding scheme (QP#1), the parameters-based search space (QP#2), the multi-objective approach (QP#3), the self-tuning strategy (QP#4), and the parallel computing paradigm (QP#5). The considered subset (QP#1 to QP#5) allows us to provide a more comprehensive analysis of the maturity of the field, as well as how healthy it has become as we argue next. It has been proven by the early NI&M-based IR methods that the usage of real-coding schemes provided a better performance [3], and this trend can also be shown in the last decade due to more than ninety percent (see QP#1 in Figure 4) of the overviewed contributions that make use of this kind of coding scheme (Other contributions also echoing this fact [71][72][73].). Only two of them did not adopt such design. In similar manner, it is remarkable the preference (87.5%) for the usage of the search space based on parameters (QP#2). The latter fact confirms the previous results shown in [35], in which those methods using such approach outperformed the results of those based on the point-matching scheme. Moreover, the last decade is also distinguished by the increased diversity of NI&M-based algorithms proposed for addressing the IR problem. This is a feature which exhibits the good health of the field. Moreover, several of the selected contributions provided extensive experiments comparing several of the best NI&M-based IR methods in the SoTA, which brings the opportunity of deepen on the enhancement of both the understanding of the behavior of the new NI&M-based optimization techniques addressing the IR problem and the performance of the best methods in the field. Additionally, it is clear that the number of contributions using novel designs based on self-tuning of the control parameters (QP#4) represent the twenty-nine percent of the methods, which depicts the rise in the interest in adding this feature. Specifically, two of the best methods dealing with the 3D reconstruction of range images made use of this particular design [15,49]. However, more contributions dealing with this feature are needed in the near future in order to achieve a better understanding of the mechanisms and inherent benefits behind it. Another aspect of quality is the consideration of a multi-objective approach (QP#3). Only three contributions, i.e., 12.5% (one of them being an enhancement of previous research), have been proposed and none of them compared their performance with other NI&M-based IR methods in the SoTA. In similar manner, it is necessary to increase the research in this aspect of design in order to reveal the degree of improvement provided by this new approach tackling the IR problem. Finally, only sixteen percent of the revised papers contributed with a parallel implementation (QP#5) of the proposal. Actually, in [37] the authors avoided using their previous IR algorithm based on GAs [16] and they adopted a different one such as the simulated annealing (SA) algorithm, due to it being considered less CPU-time demanding when facing recognition tasks. Then, it would be desirable to increase the ratio of contributions dealing with efficient IR methods by considering the current availability of High Performance Computing (HPC) resources. In fact, this trend would notably assist in the rapid expansion of new and promising approaches based on the application of deep learning techniques to medical IR [74,75].
Finally, it is necessary to introduce a brief discussion on the level of advancement of the field regarding how the new contributions built their research from previous results published in the SoTA. As stated in [4], "There are no benchmarks in this area to confirm the efficiency and universality of each method". The latter claim is clearly remarkable considering the high amount of contributions in the field in the last decade. It is necessary to share efforts in order to establish a standard benchmarking for testing IR proposals; then it will easily and feasibly allow the community to test the performance of the new proposals against those old ones in the SoTA. More importantly, an inherent principle of publication is that other authors should be able to replicate and build their research upon previously published claims. Thus, additional efforts need to be made to achieve the goal of making available the code of the proposed methods, which undoubtedly will reinforce the previous request regarding the achievement of the suitable benchmarking of future NI&M-based IR contributions.