Automated marker-free longitudinal infrared breast image registration by GA-PSO

The automated marker-free longitudinal Infrared (IR) breast image registration overcomes several challenges like no anatomic fiducial markers on the body surface, blurry boundaries, heat pattern variation by environmental and physiological factors, nonrigid deformation, etc, has the ability of quantitative pixel-wise analysis with the heat energy and patterns change in a time course study. To achieve the goal, scale-invariant feature transform, Harris corner, and Hessian matrix were employed to generate the feature points as anatomic fiducial markers, and hybrid genetic algorithm and particle swarm optimization minimizing the matching errors was used to find the appropriate corresponding pairs between the 1st IR image and the nth IR image. Moreover, the mechanism of the IR spectrogram hardware system has a high level of reproducibility. The performance of the proposed longitudinal image registration system was evaluated by the simulated experiments and the clinical trial. In the simulated experiments, the mean difference of our system is 1.64 mm, which increases 57.58% accuracy than manual determination and makes a 17.4% improvement than the previous study. In the clinical trial, 80 patients were captured several times of IR breast images during chemotherapy. Most of them were well aligned in the spatiotemporal domain. In the few cases with evident heat pattern dissipation and spatial deviation, it still provided a reliable comparison of vascular variation. Therefore, the proposed system is accurate and robust, which could be considered as a reliable tool for longitudinal approaches to breast cancer diagnosis.


Introduction
The number of patients with breast cancer has significantly increased since the 1970s due to lifestyle changes (Laurance 2006).In the last decade, the breast has been the leading site of invasive cancer in women.Breast cancer causes the deaths of several hundreds of thousands of women every year.For breast cancer cases detected early, the survival rate is over 80%.Infrared (IR) imaging features the advantages of being noninvasive, risk-free, and low-cost.IR images reveal the heat distribution on the surface of a human body.Since the 1960s, IR photograms have been used for numerous applications in medicine diagnostics, and this technology was approved by the US Food and Drug Administration as an adjunctive diagnostic instrument for breast cancer in 1982.Thus, IR imaging has long been considered to be a potential medical imaging modality for diagnosing breast cancer.It is suitable for the early detection of breast cancer and frequent follow-up, and it does not incur a substantial financial burden (Wishart et al 2010).
During tumor growth, the metabolic activity and vascular circulation of tumor tissues result in areas with higher temperature than those of the surrounding normal tissues (Guo and Wan 2007, Arora et al 2008, Tan et al 2009).Therefore, tumors may be detected based on observations of the symmetry of vessels and the difference between high-and low-temperature areas (Barnes 1963, Keyserlingk 1997, Jones 1998, Boquete et al 2012, Mashekova et al 2022).Most breast cancer detection methods rely on once-off information, such as Sentinel Breast Scan and No Touch Breast Scan (Mambou et al 2018).However, the results of once-off IR imaging may be too arbitrary to enable diagnosis of breast cancer; IR images are easily influenced by the effects of human physiology and the environment on temperature, and the vessels in the left and right breasts are not always symmetrical in individuals.Therefore, longitudinal approaches based on heat pattern variations over several time points could be more reliable for determining whether cancerous tissues are present in the breast.
In 1974, Stark used the thermogram to screen participants and revealed that all participants had unique and stable heat patterns (Stark andWay 1974a, 1974b).In other words, longitudinal approaches can enable observation of variations in the same regions on the breast surface over time, if longitudinal IR images taken at different times are well aligned in the spatiotemporal domain.Therefore, the registration of longitudinal IR images is an essential step in diagnostic longitudinal IR imaging studies.After longitudinal IR image registration, the pixel-wise tissues on the breast surface can be analyzed according to the change of heat areas (Zylberberg et al 1982, Chung et al 2011).However, longitudinal IR image registration has generally been difficult, mainly due to the following three factors: (1) It is impractical and difficult to keep markers attached to the body surface of participants for several weeks or longer.
(2) X-ray mammograms, magnetic resonance imaging (MRI), and ultrasound imaging employ anatomic fiducial markers as feature points to establish relationships for transformation models (Meyer et al 1997, Kostelec et al 1998, Crum et al 2004), but IR images do not reveal anatomic fiducial markers on the body surface.Furthermore, in IR images, the boundaries and sizes of heat patterns may vary because of environmental and physiological factors, although the patterns themselves are generally preserved.
(3) The manual determination of corresponding pairs of fiducial points between two IR images is not only labor-intensive and time-consuming but also likely to lead to interobserver differences.This renders it impractical for clinical applications.
Additionally, breast image registration, similar to the registration of most medical image types, is classified as a nonrigid problem because of the softness and flexibility of the human skin.Nonrigid registration is more difficult than rigid registration because underlying nonrigid transformations are often imperceptible, complex, and difficult to model (Ma et al 2013).Therefore, biological distinctions between the structures of two images cannot be achieved without localized stretching of the images.
In our previous study, the proposed nonrigid registration method for IR breast image was expected to overcome the aforementioned challenges.The shape context algorithm was used for initial feature point matching and the registration results were optimized by a self-organizing map (SOM) (Lee et al 2017).The mean error of registration results was 1.99 mm, and 85% faster than manual registration.However, the variation of breast thermal images in clinic is usually complex, especially for those subjects treated with chemotherapy.The registration results of those cases with complex vascular structure or evident heat pattern dissipation were poor, because the limitation of the shape context algorithm is that a significant quantity of mismatching pairs may occur between two images with lower global and local similarity, and SOM is easily trapped in local maximum, which may lead to choosing the wrong corresponding pairs.Therefore, the new longitudinal IR image registration proposed in this study was aimed at improving the aforementioned problems and extending the practicability in clinic.
To perform the new longitudinal IR breast image algorithm for clinical application, it was expected to have the following abilities: (1) Feature point descriptor: in contrast to anatomic fiducial markers of other medical images, the feature points in IR images were fewer and difficult to exact.Thus, the feature point descriptor was composed of several algorithms to extract more reliable feature points from the breast surface.
(2) Automatic feature point matching: it is impractical to determine the corresponding pairs manually, and thus the progress of matching was designed to be implemented automatically.
(3) Optimization: the matching results may need to be optimized to find better results because IR images are blurred and easily affected by physiological and environmental influences.Moreover, the optimization would be designed to have the opportunity to prevent trapping in a local maximum.
(4) Deformation model: after the relationship between the two IR images was established successfully, the nonlinear model was used to deform the sequential images to be aligned appropriately in the spatiotemporal domain.
With those aforementioned abilities, the proposed algorithm was expected to have more reliable image registration results for longitudinal approaches to breast cancer diagnosis.
The heat pattern of infrared breast images is like the fingerprint for everyone to describe its content.Thus, the feature extraction from the vascular structure is the primary step in establishing the relationship between the two IR images.In contrast to many feature extractions based on curvature, singular values, and wavelet, corner detection is the most intuitive approach for extracting feature points from the heat pattern.The famous corner detections are the Moravec corner detector (Moravec 1980), the Harris corner detector (Harris and Stephens 1988), and the Fast corner detector (Trajković and Hedley 1998).In 1980, Moravec found that the difference between the adjoining pixels in a uniform region of the image is small, and the difference is significantly high in all directions at the corner.In 1988, Harris improved the problem of the Moravec corner detector.It overcame the restriction of calculating only eight discrete 45-degree angles of the strength value.In 1998, the Fast corner detection approach was proposed and had a better performance in computation time, but it is not as rigorous as Harris corner detection, and lots of useful corners are lost.Moreover, other feature points detections such as scale-invariant feature transform (SIFT) (Lowe 1999) and speeded-up robust features (SURF) (Bay et al 2008) with robustness in resisting interference of changes in scale, rotation, illumination, and local affine distortion, would be also considered.They combine a scale-invariant region detector and a descriptor based on the gradient distribution in the detected region to find the extreme points as the feature points.Furthermore, vessel enhancement is an important preprocessing step in accurate vessel-tree reconstruction in the field of medical images which inspired us that the direction of the vessel could be another feature of the heat pattern.The vascular structure can be enhanced and extracted by the Hessian matrix, because it is based on the fact that a vessel has one principal direction which is indicated a numerical relation between the two eigenvalues of the matrix (Truc et al 2009).However, it could be inadequate to use only one feature point descriptor because of the low contrast and blur of infrared images.Therefore, the feature point descriptor was tried to be composed of several aforementioned approaches in this study.
To find the best feature pairs between two IR images, feature point matching and optimization were the most important steps before the establishment of the deformation model, which are often discussed in the optimal field simultaneously.The most conventional optimization algorithms focus on accurate computation but are difficult to operate in real-world settings.By contrast, evolutionary computation provides a more robust and efficient approach for dealing with complex real-world problems (Goldberg 1989, Yao 1999, Fogel 2006).An evolutionary algorithm known as particle swarm optimization (PSO), which serves as a nature-inspired optimization technique, was originally introduced by Eberhart and Kennedy in 1995 (Eberhart andKennedy 1995, Del Valle et al 2008).PSO is used extensively for optimization problems in numerous areas such as power systems (AlRashidi and El-Hawary 2008), industrial electronics (Ruiz-Cruz et al 2012), wireless sensor networks (Kulkarni and Venayagamoorthy 2010), and feature selection (Xue et al 2012), because its structure is simple, and the algorithm is easy to implement.However, the shortcomings of PSO are similar to SOM.PSO often calculates nearly optimal solutions, which leads to premature and rapid convergence (Van den Bergh andEngelbrecht 2002, 2010).
To improve the aforementioned problems, Alcock proposed that an evolutionary process included not only the behavioral variation of a population but also genotypic variations (Alcock 2009), which enhanced the performance of optimization algorithms.Moreover, genetic algorithms (GA) are the most well-known type of evolutionary algorithm (Juang 2004).They are stochastic search procedures that involve three mechanisms of natural selection, crossover, and mutation.GA is relatively suitable for global searches, and some studies have shown that GA can adjust its mutation step size dynamically to better reflect the granularity of a local search area.Many other evolutionary algorithms, such as memetic algorithms (Radcliffe and Surry 1994), multiobjective evolutionary algorithms (Deb et al 2002), genetic programming (Angeline 1994), evolutionary programming (Fogel 1994), and evolution strategies (Rechenberg et al 1994), which are also heuristic and stochastic, are less likely to exhibit problems related to a local minimum.It indicated the potential of GA used for the enhancement of PSO solution quality.Accordingly, some hybrid GA-PSO methods have been proposed that combine the advantages of PSO and GAs.For example, Juang proposed a hybrid GA and PSO method and applied it to recurrent neural network designs (Kennedy and Eberhart 1995).Subsequently, Moradi and Abedini hybridized PSO with GA to achieve optimal location and sizing in distribution systems (Moradi and Abedini 2012).
Valdezet et al integrated GA and PSO by using a fuzzy logic method for decision-making (Valdez et al 2008).Numerous hybrid methods are based on combinations of GA and PSO (Premalatha and Natarajan 2009).The concept of most hybrid algorithms is similar: the population is separated into two parallel groups controlled by GA and PSO individually and then the groups are recombined at set intervals.However, Gong et al proposed that this type of hybridization mechanism does not operate as well as expected, because the effects of the interactions of GA and PSO are difficult to be recognized (Gong et al 2015).Gong et al hybridized GA and PSO in a cascade loop and then used GA to assist PSO in determining which particles should be eliminated or in generating new particles.Therefore, the optimization in this study was inspired by Gong and used to determine the corresponding pairs of fiducial points on the breast surface between two IR images.
The paper is organized as follows.Section 2 describes the architecture of the IR image registration system, including the IR spectrogram hardware system and the algorithms of feature point descriptor, optimization for matching and deformation model.Section 3 summarizes the experimental results, which demonstrated that the proposed IR image registration achieved subpixel accuracy in the clinical trial.Finally, discussions and conclusions are drawn in sections 4 and 5 respectively.

Materials and methods
The longitudinal IR image registration system is comprised of a photograph hardware module, a feature point detection module, an optimization module and an image deformation module.The pipeline, displayed in figure 1, is described as follows: (1) The hardware part: the IR spectrogram hardware was developed for the longitudinal IR imaging system.
The IR breast images of volunteers were captured and transferred to the computer for image processing.
(2) The feature point detection part: the heat pattern and the branch points of the skeletons of the heat pattern on the breast surface as the initial fiducial points were implemented by the Harris corner detector algorithm and Hessian matrix to find corner points and cross-points of the middle lines in the vessel heat pattern, respectively (Chen et al 2014, Lee et al 2018).Moreover, SIFT has better performance than SURF in resisting interference of changes in scale and rotation (Juan and Gwun 2009).Therefore, the feature point descriptor was composed of SIFT, Harris corner and Hessian matrix, which extracted the feature points as the initial particles in PSO.
(3) The optimization part: the optimization named GA-PSO cascaded GA and PSO to find the best corresponding pairs between two IR images.In this method, particles in PSO were not only guided by the global maximum and best post solution but also by the contribution of the GA.
(4) The image deformation part: the outlier pairs were eliminated by random sample consensus (RANSAC).
After the corresponding pairs were found well, the relationship between the two IR images was established by thin plate spline (TPS), which was employed in image transformation to simulate the surface deformation of nonrigid objects such as the human breast (Lee et al 2018).
Until these IR images of the time course photography were registered, the quantification of the pixel-wise heat pattern variations would be able to be analyzed.

IR spectrogram system and experiment setup
The clinical trial was approved by the Institutional Review Boards of the National Taiwan University Hospital (approval number: 201103117RC), 80 patients between the ages of 25-72 (average age of 50) were involved in this study from July 2011 to December 2015, who suffered from breast cancer and were being treated with chemotherapy.In this study, the experiments divided into simulated and clinical parts are described as follows: (1) In the simulated part, the experiments were not only designed to compare the manual selection of corresponding pairs with automatic matching but also verified the validity of our image registration system.Ten normal subjects were involved in the experiments, and the electric blanket was used to warm up their breasts to simulate the metabolic activity and vascular circulation of tumor tissues, which resulted in areas with higher temperature.Moreover, several makers as the standard were kept on the breasts for evaluating the accuracy of manual selection and automatic matching, and they were not used to be the feature points for registration.
(2) To evaluate the performance of the proposed system in the clinical trial, patients were captured several times of IR breast images during chemotherapy.Due to the physiological and environmental influences, the images would reveal the effects of spatial deviation and heat pattern dissipation, and they were classified into several levels in the experiments.
The features of the IR spectrogram hardware system are a pair of mid-wavelength infrared (MIR, 3-5 μm) and long-wavelength infrared (LIR, 8-9.2 μm) cameras with high-temperature sensitivity and spatial resolution, and the mechanism is displayed in figure 2(a).In figure 2(b), the MIR and LIR cameras employed were manufactured by FLIR Systems (http://www.flir.com/home/),and both had 320 × 256 pixel detectors.Their temperature sensitivity is above 0.02 °C, and the spatial resolution at 1 m is approximately 0.6 mm.The MIR and LIR cameras are positioned on top of a platform, and the system is equipped with seven degrees of freedom.Each camera can move freely along the x-axis and can rotate on the x-y plane, which accounts for four degrees of freedom.The platform supporting the two cameras can move along the x-and z-directions, providing the fifth and sixth degrees of freedom.Additionally, an image-based calibration system is positioned on the top of the chair to ensure that the central optical axes of both cameras are on the same plane.The location of the chair can be adjusted along the rail on the base, which enables the seventh degree of freedom.The IR spectrogram hardware system was designed to provide an imaging environment with favorable reproducibility and flexibility to assist in the development of an IR spectrogram imaging system for practical use.
To achieve a high level of reproducibility, a volunteer was required to sit in the same location to reduce spatial errors (relative to the cameras) when IR images were taken at different times.Moreover, the threedimensional coordinates of the chair and cameras as well as the azimuth angles of both cameras were recorded for volunteers receiving multiple checkups.Before the photograph, each volunteer would sit and wait about 10-20 min for temperature stabilization.After that, they should take off their clothes to expose their breasts, be seated on the chair, and place their arms and hands on the armrests, which is represented in figure 2(c).The process of IR photography takes about 5 min to take several images from different angles, and these breast images are captured by the MIR and LIR cameras simultaneously and harvested for registration.In this study, the spatial relationship between the two cameras and volunteers had been established, and thus MIR and LIR images were deformed by the same nonlinear model, which was calculated by the proposed longitudinal image registration algorithm.

Scale invariant feature transform (SIFT)
SIFT is used for detecting and describing local features in an image because of its robustness in resisting interferences from changes in scale, illumination, and local affine distortion.Therefore, SIFT is suitable for identifying the initial corresponding points between the first and nth IR images.Feature points were detected from two steps: extracting and filtering.In the first step, the IR image was convolved using Gaussian filters at different scales, and the differences of these consecutive Gaussian-blurred images were then determined.The Difference of Gaussian (DoG) of the IR image between different scales, s k i and s k , j is as follows: )is the Gaussian function at scale s k , and I x y , ( ) is the original IR image.Once the DoG images have been obtained, feature points are identified as the local minimum or maximum of the DoG images within different scales.The feature points are generated by comparing each pixel in the DoG images to its eight neighbors at the same scale and nine corresponding neighboring pixels in each of the neighboring scales.However, in this method, too many feature points were produced, and some were unreliable.Therefore, the second step involves filtering out points that exhibited a low contrast or were poorly localized along an edge.
Finally, from the selected feature points, corresponding points were identified through the point description process.A feature point can be represented by the magnitude m x y , ( ) and direction q x y , ( )terms at its scale, which are as follows: The parameters of magnitude and direction are considered for each pixel close to the feature point.An orientation histogram with 36 bins is formed, with each bin covering 10 degrees.Each sample in the neighboring window added to a histogram bin is weighted by its gradient magnitude and by a Gaussian-weighted circular window with a sigma that is 1.5 times that of the scale of the keypoint.The peaks in this histogram correspond to dominant orientations.Once the histogram is filled, the orientations corresponding to the highest peak and local peaks that are within 80% of the highest peaks are assigned to the keypoint.If multiple orientations are assigned, an additional keypoint is created with the same location and scale as the original keypoint for each additional orientation.In figure 3(a), the matching results of corresponding points between the first IR image and the nth IR image were extracted by SIFT.The corresponding pairs with red color were connected by the blue line.

Corner points and branch points of the skeletons of the heat pattern 2.3.1. Harris corner detector
This is used to determine feature points on the corners of a heat pattern.The local autocorrelation function of a signal that measures the local changes of the signal with patches is shifted by a small amount in different directions.The weighted sum of squared differences D W , is the autocorrelation function where I x y , ( ) is the original IR image, ) is approximated by Taylor expansion, and Ga x y , ( )is a Gaussian function.Let I x and I y be the partial derivatives of I with respective to x and y, respectively.Then This leads to the approximation which can be written in matrix form as . 9 x y x x y x y y According to the Harris corner detector, the structure tensor J should have two 'large' eigenvalues for a corner.
In figures 3(b) and (c), the locations of corner points were shown as the small square frames on the first and nth IR images.

Hessian matrix
The skeleton of a heat pattern was generated by evaluating the eigenvalues of the Hessian matrix for each pixel.This can provide the morphological information of heat pattern, and the cross-points of the skeleton were identified as other feature points.The vessel heat pattern was identified using the Hessian matrix where ) is the original IR image, and a VR index is used to quantify the response to the characteristic values of the Hessian matrix in the vessel as follows: where l 1 and l 2 represent the smallest characteristic absolute value and the largest characteristic absolute value, respectively.In this process, b and c are 0.5 and 450, respectively.Subsequently, the threshold of the VR index must be determined to obtain the image represented by l 1 and l 2 in the vessel region of the IR image, which is 0.0005.Finally, MATLAB was used to identify the middle lines of the vessel by applying the command 'bwmorph,' and the cross-points of the skeleton can be easily detected by considering pixels with more than three neighbors.
The corner points and the branch points of the skeletons of the heat pattern were derived, as shown in figures 3(b) and (c), for the first IR image and the deformed nth IR image, respectively.These fiducial points and keypoints of SIFT were defined as the initial corresponding points.

Hybrid genetic algorithm and particle swarm optimization (GA-PSO)
In our proposed methodology, the main approach involved the incorporation of genetic evolution and particle searching into the optimization process of identifying the corresponding pairs between two images.Thus, PSO was conducted before GA to search for the optimal solution, and the GA assisted in preventing premature convergence.Initially, one of the corresponding points on the first IR image was selected and the closest feature point was identified (if several feature points were found within the search area, they competed with each other simultaneously) on the nth IR image as the initial corresponding point, which was marked as the red dot on the 'Initialization' block of figure 4. The GA-PSO mechanism was employed to determine the optimal solution using the following steps, and the pseudo-code was described in algorithm 1.

PSO step
In PSO, searches of the solution space are based on the movement of each particle, which is determined by an assigned vector.Each particle updates its velocity in three parts, namely current velocity, the best previous position of the particle, and the global best position of the population.This involves randomly generating a set of candidate points around the initial corresponding point.At each time step t, a fitness function f x t i ( ( ))is calculated to evaluate the position of particles.Moreover, the velocity and position of particles are consistently updated by using the best previous position of the particle, x t , i *( ) and the global best position of the population, x t , gbest * ( ) until the terminal condition is satisfied.
For all particles, ( ) [ ( ) ( ) ( )] represent velocity and position, respectively, where n is the population size.The new velocity for particle ith is updated by where w denotes the inertia weight, c 1 and c 2 are acceleration constants, and r 1 and r 2 are uniformly distribution numbers in [0, 1].For each particle, TPS is performed together with these corresponding pairs (keypoints pairs and determined pairs) on the nth IR image, which results in a candidate deformed nth IR image, and the mutual information is calculated as the fitness function f x t i ( ( ))between the first IR image and the candidate deformed nth IR image.The best previous position of the particle and the global best position of the population, which are evaluated by the fitness function, are recorded to update the velocity and position.Moreover, one of the terminal conditions is based on the fitness function, and the other one is the number of iterations of particle update.Once one of them is achieved, the corresponding pair is determined by the best particle on the nth IR image and the corresponding point on the first IR image, which is shown in the 'Correspond pair establishment' block of figure 4.

GA step 2.6.1. Selection
To increase the efficiency of PSO, some particles with the least favorable performance are eliminated.After D time steps, the ranking of the particles according to the fitness function reveals which particles are closest to the best solution.In this study, the top 60% of f x t i ( ( ))are saved as the selected competitors, and the rest are removed.The selected competitors compete with the new particles generated from crossover and mutation.

Crossover
Pairs are randomly selected from saved particles, and their offspring are reproduced by the best search experiment and global best position.M offspring o t c ( ) are generated as where G and T respectively denote the genetic weight and the talent weight with Gaussian distribution in [0, 1], and }is the index of offspring.Generally, the offspring exhibit certain characteristics such as the genetic weight of their parents, moreover, they may have individual talents that contribute to the global best position, increasing the efficiency of the search process.

Mutation
To prevent the solution from converging prematurely at the local maximum, the mutation increases the diversity of the searching space and makes the exploration more credible.Whether the mutation has occurred is determined based on the probability of mutation P .mutation In this study, we posited that if the talent weights of offspring were extremely high, their performance was markedly different from that of their parents.Therefore, an offspring with a mutation is reproduced as where b min and b max respectively represent the lower and upper bounds in the search space.Based on the GA- PSO algorithm, the population gradually moves to the best position.

Random sample consensus (RANSAC)
After the corresponding pairs between two images were determined, any incorrect matches between points may substantially influence the registration result.To remove these mismatched points, the RANSAC algorithm was used to estimate a transformation matrix (homography matrix) that suppressed the transformation error of all matched points to a minimum.First, four pairs of matched particles are randomly selected to solve the transformation matrix as follows:

¢ ¢ = x y h h h h h h h h h
x y where x y , ( ) and ¢ ¢ x y , ( )are the coordinates of the pairs.The best transformation matrix is established when the cost function is minimal

x h y h h x h y h y h x h y h h x h y
The transformation model is established based on the estimated matrix to evaluate whether the matching error is below the threshold.Therefore, all of the matched points comply with the model as inliers, and the other points are discarded when they fail to comply with the model as outliers.In figure 5(a), these corresponding pairs between the first IR image and the nth IR image were matched by GA-PSO.The blue lines were identified as inlier pairs by RANSAC.In contrast, the green lines are identified as outliers, which were eliminated before deformation.

Thin plate spline (TPS)
The nonlinear model based on TPS was employed to deform the nth IR image to increase its similarity to the first IR image.For given corresponding pairs of  for = ¼ i n 1, , , the TPS model uses interpolation between these corresponding pairs based on the cubic spline function.Thus, the deformation mapping f p ( ) is defined as where a a a , , x y 1 are the parameters of linear transformation, p is the coordinates (x, y) of the nth IR image, n is the number of corresponding points, i is the ith corresponding point; w i is the nonaffine parameter; and -G p p i ( )   is the shape function for representing the intensity of force between the corresponding points, its function can be described as where r is the distance between two points in the nth IR image.To calculate these parameters of deformation mapping f p , ( ) W used to be represented as them can be defined as where the O is 3 × 3 zero matrix.The U and Q are further defined as followings: Therefore, the unknown parameters a 1 ,a x ,a y , and w i can be calculated by the corresponding pairs between the first and nth IR image.After those parameters are determined, other points in the nth IR image can be deformed by f p .
( ) The pixel-wise heat pattern variations quantification After these IR images were registered, the heat pattern variations can be observed as longitudinal approaches in the same regions on the breast surface over time.As represented in figure 6, the patient with breast cancer underwent thermography during chemotherapy.The region marked by the red circle in the nth IR image denoted the heat pattern of a tumor, the position of which was identified based on the clinical reports of MRI results.In this region, the heat energy decreased and the heat pattern evidently dissipated after treatment.This trend was the same as those reported previously.In other words, the tumor size could be determined based on changes in heat energy and heat pattern.Moreover, the physiological and environmental influences were further minimized through the quantification method (Chen et al 2012), which was one of our previous studies for breast cancer detection.However, the experiment results in this study would still focus on the evaluation of longitudinal image registration.

Error evaluation for the simulated experiments
Ten normal subjects were involved into the simulated experiments and captured 7-16 times IR photograms on the same day.Several makers as the standard (red points) were kept on their breasts only for evaluation of the performance, as shown in figure 5(b).During the IR photography, the heat energy decreased and the heat pattern slightly dissipated on the breast surface, which indicated that most details of the vascular structure were preserved.
The performance comparison is comprised of manual determination, the previous method (Shape Context + SOM) and GA-PSO.The corresponding pairs in manual determination rely on the manual selection of the feature points generated from the Harris corner and Hessian matrix.The comparison results are presented in table 1, the mean distance between the true marker positions and the marker positions on the deformed images was calculated.In the column of manual determination, the mean difference of these cases is 3.88 mm, which is reasonable for image registration.In the GA-PSO column, the mean difference of our proposed system is 1.64 mm, which increases 57.58% accuracy than manual determination.Even if it compares with the previous study, it still makes 17.4% improvement in accuracy.Moreover, the distribution of the distance between the true marker positions and the marker positions on the deformed images within these ten cases is shown in the boxplot.The interquartile range and variation of manual determination are much higher than our methods, and the GA-PSO has the best performance in these trials.The interquartile range of the previous method is close to GA-PSO, but its maximum is obviously bigger than GA-PSO, which indicates that the effectiveness of the proposed algorithm has the ability to search better corresponding points.

The performance of longitudinal ir image registration in the clinical trial
In the clinical trial, 80 patients were captured several times of IR breast images during chemotherapy, and these IR images not only recorded the variation of vascular structure but also may have heat pattern dissipation.According to the classification of the IR images database, there are three levels as follows: (1) slight heat pattern dissipation and spatial deviation, (2) evident heat pattern dissipation and slight spatial deviation, (3) evident heat pattern dissipation and spatial deviation.The parts of image registration results are also represented as follows:

Slight heat pattern dissipation and spatial deviation
As displayed in the first row of figure 7, Case 36 (C36) with the malignant tumor (size 5.5 cm −> 4.4 cm) under the central and upper hemisphere of the right breast was captured 11 times IR photograms before the operation, and the heat pattern of the vascular structure was marked by the red circle, which was located on the same Figure 9.The registration results of a breast cancer subject (C26) were acquired at six different times using IR photograms.In the first row, the heat pattern of the vascular structure marked by the red circle was evidently dissipated during chemotherapy.In the other rows, the Canny edge of the nth IR image was superimposed onto the 1st IR image, both the heat pattern of vascular structure and the contour of breasts still fit the 1st IR image well, even if the IR images had the evident heat pattern dissipation.Moreover, slight spatial deviation marked by the yellow arrows was found in the left breast.In Figure 8(b), the green line was the Canny edge of the target image and the red line was the Canny edge of the deformed Source Image 1, and the overlaid portions were represented as the white lines.Eventually, the quantitative analysis of registration results was presented in table 2. To calculate the target registration error (TRE), about 20 landmark pairs between the Target image and the Source image were manually selected and recognized by the doctor.On average, TRE values were 10.41 ± 3.364 mm before registration and 1.97 ± 0.64 mm after registration.

Evident heat pattern dissipation and slight spatial deviation
Case 26 (C26) with the malignant tumor (size 5.7 cm −>3.6 cm) under the upper hemisphere of the left breast was captured 8 times IR photograms before the operation, and the heat pattern of vascular marked by the red circle was dissipated evidently because of chemotherapy, which was shown as the first row of figure 9.Moreover, the longitudinal image registration results are presented in the fourth row of figure 9. Through the superimposition by Canny edge, both the heat pattern of vascular structure and the contour of breasts still fit the 1st IR image well, even if the IR images had the evident heat pattern dissipation.Comparing with the second row of Figure 9, the spatial deviation of the nth IR image from the 1st IR image decreased after registration.The quantitative analysis is presented in table 3. The mean of TRE values was 8.12 ± 3.88 mm before registration and 2.12 ± 0.65 mm after registration, respectively.The reproducibility of the IR spectrogram hardware system was also demonstrated in this experiment.The position of patients in the nth IR image was similar to the 1st IR image because of the records of threedimensional coordinates of the chair and cameras as well as the azimuth angles of both cameras, as shown in the third row of figure 7 and the second row of figure 9.In the aforementioned cases, primary spatial deviation marked by the yellow arrows was found in the left breast, and it slightly affected the registration results by minimizing the spatial errors.Therefore, the IR spectrogram hardware system provides an imaging environment with high reproducibility and flexibility to assist in the development of the longitudinal IR imaging system for practical use.

Evident heat pattern dissipation and spatial deviation
In a few cases, the spatial deviation of the nth IR image from the 1st IR image was huge, which would lead to the accuracy decrease of the registration results.Moreover, if their heat pattern also evidently dissipated after treatment, the corresponding pairs would be determined difficultly.For example, Case 51 (C51) with the malignant tumor (size 3 cm −> 1.4 cm) under the upper inner quadrant of the left breast was captured 14 times IR photograms before operation, who was a severe case with serious heat pattern change and spatial deviation.As shown in the first row of figure 10, the heat pattern of the vascular structure marked by the red circle was dissipated evidently, especially after the treatment of Chemotherapy 2.Moreover, the spatial deviation of the nth IR image from the 1st IR image is shown in the third row of figure 10, which was found in the whole chest, may be caused by the patient's breath.To increase the accuracy of registration results in the severe case, coefficients of GA-PSO were adjusted to find more appropriate corresponding pairs between the 1st IR image and the nth IR image, which included expanded search area, smaller acceleration constants c 1 and c , 2 higher probability of mutation P mutation and more iterations.The registration results are presented in the second and the fourth row of     Especially the heat pattern of vascular structure, the Canny edge of the nth IR image was superimposed onto the 1st IR image.However, the edge of the images does not fit well because of the insufficiency of precise feature points on the edge, which is marked by the green arrows.It also reveals that the nonlinear model based on TPS still has room for improvement of the deformation of the body.Eventually, the quantitative analysis is presented in table 5, which indicates that the mean of TRE values is 28.36 ± 4.77 mm before registration and 2.99 ± 1.7 mm after registration, respectively.In contrast to the frontal view cases, the lateral view case has the worse TRE results, which may be caused by parallax.

The feature point descriptor for the cold tissue
The abdomen in the infrared images may provide very important information about breast cancer, but it is often ignored because of its' less vascular structure.In our IR photography system, the information of the abdomen was preserved for subsequent studies.However, a large proportion of the abdomen is cold tissue, which results in areas with lower temperature and lower contrast in the IR images, as well as fewer feature points.Thus, the accuracy of the abdomen in our image registration results is lower than the chest.If a feature point descriptor can be employed in the low contrast area, the accuracy may be further increased.Lots of prior studies have been proposed to recognize and track objects with low contrast.For example, Chan et al proposed an active contourbased object detection algorithm that does not require edge information (Chan and Vese 2001).Tsai et al proposed an independent component analysis-based filter for detecting low-contrast defects on liquid crystal display glasses in uniform surface images (Tsai et al 2006).The Kalman filter has been shown to be capable of tracking small low-contrast objects (Davies et al 1998).With the aforementioned methods, the feature point descriptor in this system is going to be a part of our study in the future.

Conclusion
The longitudinal IR image registration is an essential step toward the quantitative pixel-wise analysis of the heat energy and patterns change in a time course study.To achieve the goal of practical use, the proposed longitudinal IR image registration system is comprised of a photograph hardware module, a feature point detection module, an optimization module and an image deformation module.The hardware module, it was developed for the longitudinal IR imaging system with favorable reproducibility and flexibility.In the feature point detection module, SIFT, Harris corner and Hessian matrix were employed to generate the keypoints, the corners of the heat pattern and the cross-points of the skeleton of the heat pattern, which were considered to be anatomic fiducial markers on the breast surface.In the optimization module, the GA-PSO algorithm was used to find the appropriate corresponding pairs from the feature points between the 1st IR image and the nth IR image, and it minimized the registration errors from the heat pattern variation because of the environmental and physiological factors.In the image deformation module, the nth IR image was deformed to align with the 1st IR image by TPS.
The performance of the proposed longitudinal image registration system was evaluated by the simulated experiments and the clinical trial.In the simulated experiments, the ten normal subjects were involved in warming up their breasts with the electric blanket to simulate the metabolic activity and vascular circulation of tumor tissues.The mean difference of our system is 1.64 mm, which increases 57.58% accuracy than manual determination and makes 17.4% improvement than the previous study.In the clinical trial, 80 patients were captured several times of IR breast images during chemotherapy.Most of them were well aligned in the spatiotemporal domain.Even if the few cases with evident heat pattern dissipation and spatial deviation resulted in the accuracy decrease, the registration results still provide a reliable comparison of vascular variation.Therefore, the proposed longitudinal image registration system is accurate and robust, which could be considered as a reliable tool for longitudinal approaches of breast cancer diagnosis.

Figure 1 .
Figure 1.Flowchart of the longitudinal IR image registration algorithm.

Figure 2 .
Figure 2. (a) Mechanism of the IR spectrogram hardware system.(b) The MIR and LIR cameras were manufactured by FLIR Systems and have several degrees of freedom.An image-based calibration system was constructed on the top of the chair, which can be adjusted along the rail on the base.(c) Each volunteer was seated on the chair, took off their clothes and placed their arms on the armrests.

Figure 3 .
Figure 3. (a) Matching results of corresponding points between the first IR image and the nth IR image using SIFT.(b) Feature points on the first IR image.(c) Feature points on the nth IR image.

Figure 5 .
Figure 5. (a) The corresponding pairs between the first IR image and the nth IR image were matched by GA-PSO.The blue lines were identified as inlier pairs by RANSAC, and the green lines were as outliers.(b) The subject with several makers (red points) on her breast on the same day.

Figure 6 .
Figure 6.Regions of the IR images marked by red circles indicate the heat pattern variation related to a tumor during chemotherapy.

Figure 7 .
Figure7.The registration results of a breast cancer subject (C36) were acquired at six different times using IR photograms.In the first row, the heat pattern of the vascular structure marked by the red circle was gradually dissipated during chemotherapy.In the third, fourth and fifth rows, edges of the C36 deformed images (red curves) derived by a Canny edge detector overlayed on the target image, both the heat pattern of vascular structure and the contour of breasts fit the 1st IR image well.Moreover, slight spatial deviation marked by the yellow arrows was found in the left breast.

Figure 8 .
Figure 8.(a) The green line and the red line were the Canny edge of Source Image 1 before registration and after registration, respectively.(b) The green line was the Canny edge of the target image, the red line was the Canny edge of the deformed Source Image 1, and the overlaid portions were represented as white lines.
position of the MRI clinical reports.In contrast to the target image, the heat pattern in the red circle gradually dissipated since the source image 2, which may be affected by chemotherapy and physiological and environmental factors.After image registration by the proposed system, the target image (as the 1st IR image) and those source images (as the nth IR image) were well aligned in the spatiotemporal domain, which is shown in the second row of figure 7.Moreover, the Canny edge of the nth IR image was superimposed onto the 1st IR image in the third, fourth and fifth row of Figure7, which was used to observe spatial deviation of the nth IR image from the 1st IR image.The spatial deviation between the 1st IR image and the nth IR image decreased obviously after registration.Through the superimposition by the Canny edge, both the heat pattern of vascular structure and the contour of breasts fit the 1st IR image well.Moreover, the effectiveness of the proposed algorithm is further shown in figure 8.In figure 8(a), the 'Source Image 1' with a more completed heat pattern was selected as the instance to display the comparison of Canny edge in the same breast image before registration and after registration with different colors.The green Canny edge represented the edge of Source Image 1 before registration, and the red Canny edge represented the edge after registration.It showed the line's movement in the same figure.

Figure 10 .
Figure10.The registration results of a breast cancer subject (C51) were acquired at six different times using IR photograms.In the first row, the heat pattern of the vascular structure marked by the red circle was evidently dissipated after chemotherapy 2. In the third row, obvious spatial deviation of the nth IR image from the 1st IR image was found in the whole chest.The registration results are presented in the second and the fourth-row fourth rows.Even though the accuracy of registration results gradually decreases (marked by the green arrows) after Chemotherapy 2, it still provides the reliable comparison of vascular variation for longitudinal approaches of breast cancer diagnosis.

Figure 11 .
Figure11.The lateral view registration results of a breast cancer subject (C14) were acquired at six different times using IR photograms.In the first row, the heat pattern of the vascular structure is marked by the red circle.In the second row, obvious spatial deviation caused by the patient's breath was found in the whole chest (marked by the yellow arrows).The registration results are presented in the first and the third rows.The heat pattern of vascular structure is aligned in the spatiotemporal domain well.However, the edge of the images does not fit well because of the insufficiency of precise feature points on the edge (marked by the green arrows).
The performance of longitudinal IR image registration in the cases with outer quadrant tumors In the previous studies, the upper outer quadrant of the breast was the most frequent site for incidence of breast cancer, and thus those patients were usually captured in the lateral view of IR images to observe the whole vascular structure in this study.In contrast to the frontal view, the spatial deviation of the lateral view of IR images is usually more severe.The longitudinal image registration results of the lateral view are shown in figure11.Case 14 (C14) with the malignant tumor (size 4.5 cm −> 3.2 cm) under the outer half of the right breast was captured 17 times IR photograms before operation, who was the lateral view case with slight heat pattern change and serious spatial deviation.As shown in the second row of figure 11, the patient's pose in the 1st IR image is very different from others, which is mainly caused by the patient's breath and marked by the yellow arrows.After the image registration, the results are presented in the third row of figure 11.Most of the sequential images are aligned in the spatiotemporal domain well.

Table 1 .
The performance comparison.

Table 2 .
Tre values of the registration results in case c36.

Table 3 .
Tre values of the registration results in case C26.
figure 10, the sequential images are aligned in the spatiotemporal domain.Furthermore, the quantitative analysis was presented in table 4, which indicated that the mean of TRE values was 21.33 ± 3.88 mm before registration and 2.39 ± 1.04 mm after registration, respectively.Even though the accuracy of registration results gradually decreases (marked by the green arrows) after Chemotherapy 2, it still provides a reliable comparison of vascular variation for longitudinal approaches of breast cancer diagnosis.Therefore, the proposed longitudinal image registration system is accurate and robust.

Table 4 .
Tre values of the registration results in case C51.

Table 5 .
The values of the registration results in case C14.