ILLUMINATION INVARIANT CHANGE DETECTION ( IICD ) : FROM EARTH TO MARS

Multi-temporal Earth Observation and Mars orbital imagery data with frequent repeat coverage provide great capability for planetary surface change detection. When comparing two images taken at different times of day or in different seasons for change detection, the variation of topographic shades and shadows caused by the change of sunlight angle can be so significant that it overwhelms the real object and environmental changes, making automatic detection unreliable. An effective change detection algorithm therefore has to be robust to the illumination variation. This paper presents our research on developing and testing an Illumination Invariant Change Detection (IICD) method based on the robustness of phase correlation (PC) to the variation of solar illumination for image matching. The IICD is based on two key functions: i) initial change detection based on a saliency map derived from pixel-wise dense PC matching and ii) change quantization which combines change type identification, motion estimation and precise appearance change identification. Experiment using multi-temporal Landsat 7 ETM+ satellite images, Rapid eye satellite images and Mars HiRiSE images demonstrate that our frequency based image matching method can reach sub-pixel accuracy and thus the proposed IICD method can effectively detect and precisely segment large scale change such as landslide as well as small object change such as Mars rover, under daily and seasonal sunlight changes.


INTRODUCTION
The detection of dynamic changes on the planet's surface is not only important for the monitoring of agriculture, deforestation and urbanization on Earth but also essential for the understanding of geology structure, climate and environment on Mars.Some of these changes are big and significant, such as geo-hazards (Cheng, Wei, and Chang 2004;Lu et al. 2011), glacier motion (Berthier et al. 2005;Kaufmann and Ladstä dter 2003) and sand dune migration, and some are small, for example a newly built house or a recently landed Mars rover.The growing archive of repeat coverage and the increasing spatial resolution of Earth Observation and Mars orbital images, enable the detection and quantitative measurement of these changes.
When comparing two images taken at different times of day or in different seasons for change detection, the variation of topographic shades and shadows caused by the change of sunlight angle can be so significant that it overwhelms the real object and environmental changes.This effect caused by the local illumination variation is more problematic on Mars than that on Earth as the shadows and shades appearing on Mars orbital images are darker than those on Earth Observation images because of the thinner density of atmosphere.As this illumination variation not only alters image contrast and brightness but also imagery spatial patterns/textures, making image co-registration inaccurate and automatic detection unreliable.

RELATED WORK
The state-of-the-art change detection algorithms remove the illumination effects by radiometric correction (Schroeder et al. 2006), which can be divided into three main categories: normalization (Lillestrand 1972;Dai and Khorram 1998), filtering-based correction (Toth, Aach, and Metzler 2000) and shadow removal (Duguay and Ledrew 1992).The normalization and filtering-based radiometric correction methods, though robust to global illumination variation caused by solar strength variation, are not robust to local illumination variation caused by the sunlight angle variation.The shadow removal method may result in over compensation in areas of deep shadows of steep slopes (Tan et al. 2013).
While the global illumination variation issue has been well studied and effectively dealt with, to the best of our knowledge, there are no change detection methods truly robust to both global and local illumination variation caused by the sunlight angle variation.Thus, a Phase Correlation based IICD (Illumination Invariant Change Detection) is proposed to robustly detect and precisely segment the changes on planet's surface regardless of illumination variation and imaging geometry.

Illumination-insensitive property of Phase Correlation
Phase Correlation (Kuglin 1975) is an image matching algorithm based on Fourier shift property, which states that a translation shift between two similar images generates a linear phase difference in the Fourier frequency domain.For two images, (, ) and (, ), with a horizontal and vertical translation (, ) , the phase correlation is defined as the cross power spectrum of the Fourier Transform (FT) of the two images, which is a complex conjugation of the two frequency spectrums.
In Equation (1), (, ) is a 2π wrapped complex matrix and the density and orientation of phase correlation fringes are √ 2 +  2 and / .As long as the two images are largely overlapped (usually more than1/4 image size), the shifts (, ) can be resolved at integer level via IFT (Inverse Fourier Transform) to convert (, ) to a Dirac delta function .
Dirac delta function  is a unit impulse function, which the coordinate of the impulse peak (, ) is the image shift.The translation (, ) can also be solved directly in frequency domain by unwrapping and fitting the fringe patterns in the cross power spectrum (, ) (Foroosh, Zerubia, and Berthod 2002;Hoge 2003;Liu and Yan 2008;Morgan, Liu, and Yan 2010).The peak value of Dirac delta function , ranging from 0 to 1, indicates the quality of PC matching.If the two images are exactly the same, the peak value of  equals 1.
Equation ( 1) is established based on the assumption that the two images for matching are radiometrically similar however, the two images are not very similar in radiometry if they are taken under different lighting conditions.According to our research (Wan, Liu, and Yan 2014), the PC matrix (, ) of two images with a translational shift and under different illumination conditions can be decomposed into two independent phase correlation matrices: the illumination impact matrix  1 (, ) corresponding to texture difference caused by illumination variation and the translation matrix  2 (, ) =  −(+) resulted from the 2D shifting.
The illumination-invariant matching can be achieved by suppressing or eliminating the impact of  1 (, ) from (, ).
Our study (Wan, Liu et al. 2015) reveals how  1 is determined by sun illumination conditions.An auxiliary 3D space named as SAI (Slope-Aspect-Intensity) was proposed to present the intensity of hill shading image and terrain slope/aspect angles in a given illumination conditions.Based on the SAI space, the relationship between PC cross power spectra and local illumination conditions in terms of azimuth angle and zenith angle is investigated via mathematical derivation and experiments.The azimuth angle variation divides the PC fringes into positive correlation and negative correlation corresponding to the azimuth angle difference, because partial image texture is reversed by azimuth angle variation.The change of zenith angle alters image brightness and contrast but does not cause reverse of image texture, which is similar to adding random noise to the image phase correlation and the larger difference the zenith angles are, the more noisy the  1 matrix will be.The detailed derivation of impact of azimuth and zenith change on phase correlation matrix can be found in (Wan, Liu et al. 2015).Recent study demonstrates that the zenith angle will also modulate the strength of negative correlation; the negative fringe pattern becomes clearer with the increase of zenith angles.
Consequently, the inversed fringes caused by azimuth angle difference and increased noise caused by zenith angle variation appear in the PC matrix but the fringe density √ 2 +  2 and fringe orientation a/b remain unchanged and therefore PC algorithm is not sensitive to illumination variation.

Change type identification based on a Dirac delta function of Phase Correlation
As proved in §3.1, PC is robust to local illumination change, so the illumination change will not alter the distinctiveness of Dirac delta function peak which indicate image shift.Thus the alteration of Dirac delta function performance is only related to change information.In this section, we will further investigate the ability of PC to detect the motion and appearance change under illumination change conditions is investigated.For change detection, all the pixels in images can be generalised into the three scenarios, no change, appearance change and motion.The corresponding performances of Dirac delta function , according to Equation (2), are presented.An approach to differentiate the change information into motion and appearance change from the PC change detection result is thus introduced.

Case 1 No change
As PC is robust to illumination change proved by §3.1, when there is no change in two images taken under different illumination conditions, the calculated Dirac delta function  has a clear the peak value at (0,0) position reflecting 0 image shift in x and y direction.

Case 2 Large appearance change
When there a large appearance change (> ½ size of the image matching window) in two images,  function becomes very noisy indicating that the two regions in the matching window have a low correspondence rate.In this case, the  peak value is too low to be distinctive and consequently, the image shift calculated from  is a random number rather than 0.

Case 3 Small appearance change
When the appearance change is smaller than ½ size of the matching window, the remaining part in the matching windows can still be correlated by Phase Correlation.In this case, the peak value of  drops considerably owing to partial de-correlation, caused by small appearance change, but the peak is distinctive enough at the (0,0) location reflecting 0 image shift in x and y direction.

Case 4 Motion
When there is object motion within the matching window of PC, the related  will have a clear peak because the appearance of the moving object remains the same and the location of the maximum peak value indicates the object shift in x and y direction.In summary,  from PC matching will have a clear peak in (0,0) position when there is no change; the peak value of  decreases for appearance change depending on the areal size of the change; and the object motion shifts the location of  peak away from (0,0) position depending on the motion value.The combination of  peak values and peak locations can therefore detect the change and classify different types of change as characterised in Figure 1.A pixel will be labelled as 'no change' for a very small disparity value with a high  peak; as 'motion' for a large disparity value with a high  peak; as 'appearance variation' for a large disparity value with a very low  peak.
The 6 intensity vector maps and 24 orientation vector maps are combined together to generate a single normalised saliency map.
The GBVS method (Harel, Koch, and Perona 2006) uses a bottom-up approach based on graph cut, which includes two steps: i) activation map generation based on finding unusual values in feature vector maps and ii) generation of saliency map by normalisation of activation map.The final saliency map ranges from [0,1].Large saliency values indicate high probability of changes.Thus, an initial change detection map, a binary mask, is produced by extracting the large values in the saliency map using a global threshold   .The suggested threshold value for   is 0.5, which has been proved effective through all the tested image pairs.This threshold can be used a filter to remove the unwanted changes, if a large   is given, for example 0.8, only the significant changes will be detected, and some small changes may be neglected.where  ℎ =  3 , and  3 is mean value of  peak probability distribution map (), this area will be labelled as 'motion'.

Change
The reason for motion identification in the first step of the approach, is that object motion leads to false appearance change around the moved object and to the so called 'hallo effect' (Wang et al. 2014).It is assumed that every candidate region of change only contains one change type, and thus if the candidate change region is assigned to 'motion', it cannot proceed to the next stage and has no chance of being labelled as 'appearance change'.
For the candidate areas which do not belong to 'object motion', the following judgment is carried out to check whether they belong to 'appearance change' or not.
If the value   in the candidate area meet the following requirement,   <   , where   =  3 −  3 2 , this candidate area will be labelled as 'Appearance change'.For each blob within the neighbourhood region defined by the bounding box, three binary masks,  1 ,  2 and  3 are firstly generated based on   ,   and P. The criterion is similar to the generation of  0 , but Bayes decision is made within the local probability distribution rather than the whole image.
Then, we added two more binary masks,   To balance the robustness and the precision of the IICD algorithm, a binary mask integration function  is defined to stack the five binary mask together with various weighing factor   , The determination of   is carried out as follows: for each binary mask   () , the percentage of change area   within the bounding box () is calculated.
where  and ℎ are width and height of the bounding box ().
If the five  have similar values, this means that the final result has a high confidence, since five binary mask   () have similar change detection results.If one of the binary mask   () has a completely different value to other four masks, the relating weighing factor   will be assigned to a very small value.If we assume that the probability density function of all  obey a normal distribution, the binary mask   () which are anomalous to others will have a significant lower probability.Thus, the weighing factor   is determined as follows, If the related weighing factor  3 and  4 are large, this means that the change detection results determined by  4 and  5 agree with the initial change detection results.As algebraic manipulation based masks,  4 and  5 , are better at localising the changing boundaries of object whereas less robust to texture change, a bonus factor  = 0.3 is added to  3 and  4 , and the area is labelled as 'object change'.This means that the refined change mask  1 in this area is largely depended on the local algebraic manipulation based segmentation.Otherwise, the candidate area will be marked as 'texture change', and refined change mask  1 will be largely determined by the PC based dense matching.
The refined change mask  1 is generated by the weighted binary map stacking according to Equation ( 12) and then transform in to a binary map by the morphological based segmentation algorithm as demonstrated in Figure 5.

Snowdonia model under Daily Illumination Variation
In this section, one datasets of Snowdonia terrain shading images under daily illumination conditions were acquired.The daily illumination conditions for the simple hill and Snowdonia dataset was shown in Figure 6 in terms of azimuth and zenith angles.It can be concluded that the major change of light source within a day is azimuth angle variation, and the azimuth angle difference between morning and afternoon can be as large as 115°.Thus, both of the dataset contain reversed image patterns caused by the severe azimuth angel variation.The largest zenith angle variation is 24°.where  is the number of true positives: the ground truth image says it is a 'change' and our algorithm assigned it to a 'change';  is the number of false positive: the ground truth image says it is a 'no change' and the algorithm assigned it to 'change';  is false negative: the ground truth image says it is a 'no change' and the algorithm assigned it to 'no change'.
A good change detection algorithm needs to have a high precision as well as recall, and thus a  1 score is calculated using the combination of  and The higher value  1 score the IICD algorithm has, the better change detection result the algorithm performs.
For the Snowdonia model dataset, the ground truth image is generated by the subtraction of images with and without changes under the same illumination conditions, and then manually rectified the shape of each change area.The change detection accuracies of the Snowdonia model is presented in Figure 7   In this experiment, the vegetation seasonal colour change will not be detected by the IICD software.This is because the vegetation seasonal change alters the global grey value change of the vegetation covered areas, and thus it can be regard as global illumination change within certain areas.As the proposed method is robust to both global and local illumination change, the seasonal vegetation colour change will not be a problem for the identification of real surface change.The comparison results using state-of-the-art methods and our method are demonstrated in Figure 10.
As shown in Figure 10, filtering based CD methods perform better than differencing and normalisation based CD method in change detection, which means that the radiometric correction using homographic filter can partially remove the seasonal illumination effect, however, the average  1 of 0.1-0.4 is still unsatisfactory for effective change detection.In contrast, the proposed IICD method is able to achieve  1 of 0.87 stating that the subtle change in mountainous area can be detected and extracted while suppressing the effect of long dark shadow effect and snow coverage variation.
Figure 9 Change detection and results using the proposed IICD method.
Figure 10 Comparison change detection accuracies using stateof-the-art CD and our IICD methods for change detection

Mars rover detection
In this experiment, three HiRISE (High Resolution Imaging Science Experiment) images captured the Opportunity rover in different positions were used.The Opportunity rover is 2.3m wide and 1.6m long and in the 25cm resolution HiRISE images, the size of the rover is approximately 6×9 pixels.
According to the observation statics provided by NASA/JPL/University of Arizona, the illumination conditions vary considerably in the three HiRISE images, with the largest azimuth angle difference of 56.1° between ESP_012820_1780 and ESP_016644_1780, and the largest zenith angle difference of 8° between ESP_011765_1780 and ESP_016644_1780.
The proposed IICD algorithm was applied then to the three image pairs.For illustration, the rover detection using the image pair of  A comprehensive package of experiments using images of terrain models acquired under real illumination conditions, as well as satellite images of Mars, has been conducted to rigorously assess the capability of the algorithm.The experimental results demonstrate that the proposed PC based change detection approach is able to robustly detect real changes in object space regardless of illumination variation.Moreover, the proposed approach can also differentiate appearance change from object motion, and can correctly calculate the motion with sub-pixel accuracy.The proposed IICD algorithm the detection of subtle changes in short period and quantitative measurement from multi-temporal airborne and satellite imagery data, which is essential for fully exploit the big data information capacity.Future work will be carried out for automatically change object recognition by machine learning.For the study area, the database of multi-temporal EO images will firstly be trained through a classifier using machine learning algorithms.Based on the trained classifier, the change information detected from our PC based IICD algorithm will be automatically recognised and assigned to a more specific object type.

Figure 1
Figure 1 Basic criterion for PC based change detection Figure 2. The system is composed of two main steps: initial IICD and change quantization.The two multi-temporal optical images, T1 and T2, are firstly co-registered pixel by pixel by PC based dense matching ( §4.1.1),which results in three outputs: disparity maps in x and y direction,   and   , and  peak value map P.Meanwhile, a subtraction map S and a ratio map R are generated by the co-registered T2 images and T1 image.Based on the values of P,   and   , an initial change detection map  0 is generated, in which no change areas are regions which have small disparity values and large  peak values, and change areas are the rest regions.

Figure 2
Figure 2 Flow chart of and development path of the proposed IICD.

Figure 3
Figure 3 Pixel-wise PC based dense matching Meanwhile the corresponding Dirac delta peak value (, ) is also calculated through the PC based dense matching.Finally, disparity maps   and   are generated to rectify image T1 to the reference image T2 pixel-to-pixel.As   and   are calculated in sub-pixel accuracy, the pixel-wise image geometric distortion has been rectified in a new co-registered T2 image.A subtraction image S and a ratio map R is generated by the initial image T1 and the co-registered T2 image.
type identification: From the previous step, an initial change detection mask  0 has been generated with several candidate areas.The next step is to classify the candidate change areas into two different change types: motion and appearance change, according to the criterion proposed in §3.2.For each candidate change region, a high-low thresholding scheme is proposed to differentiate 'motion' and 'appearance change'.According to the analysis in §3.2, 'motion' produces high  peaks while 'appearance change' produces low  peaks.Two thresholds,  ℎ as high threshold and   as low threshold are introduced for the change type identification as demonstrated in Figure4.

Figure 4
Figure 4 Flowchart of change type determinationMotion identification is carried out as follows: in the current candidate area k, the mean value of  peak   is calculated.If the value   is larger than the high threshold  ℎ ,   >  ℎ , where  ℎ =  3 , and  3 is mean value of  peak probability distribution map (), this area will be labelled as 'motion'.
appearance change identification: For the appearance change areas, it is essential to extract the precise boundaries of the changing targets.As PC is a window based correlation operation, it introduces blurring effects(Aach, Kaup, and Mester 1993) to the change area boundaries.Based on the three features,   ,   and P, from PC based dense matching, the change areas can be roughly detected but the boundaries of change areas cannot be precisely determined.Thus, we added two more features based on subtraction map and ratio map for precise change area boundary estimation.As the change areas can be coarsely determined, the refinement of initial change mask  0 is based on local regions rather than the whole image.The binary mask  0 is processed by a blob analysis, which the adjacent '1' region in  0 are clustered into one blob and the statistics (size and shape) of every blob in the change map  0 are calculated.For each blob, a bounding rectangle box is calculated and the precise boundary is determined within the bounding box.Take the blurring effect caused by window based correlation into consideration; the bounding box expanded  pixels (  = 1 2 window size) to the up, down, right and left directions to the minimum bounding box subjects to the outlines of the blob.
4 and  5 , based on local ratio and subtraction map for precise change area boundary estimation.The change mask  4 is generated based on the extraction of large values in local ration map ().Then, the anomalous large values in local ratio map () which subject to object change or texture change are extracted according to the following criterion to generate change mask  4 and variance values of ().The binary mask  5 is generated using a morphological based segmentation algorithm to provide precise shape and contour information of the change areas.The local subtraction image (), as shown in Figure 5(b), is firstly transferred into a binary gradient mask, shown in Figure 5(c), by Canny edge detection algorithm (Canny 1986).The small gaps within the binary gradient mask is filled by dilation (Van Den Boomgaard and Van Balen 1992) and connectivity check (Soille 2013), presented in Figure 5(d) and (e).The change areas are segmented (Figure 5(f)) using a morphological erosion (Van Den Boomgaard and Van Balen 1992), and the largest area is detected to produce the final binary mask  5 , as shown in Figure 5(g).The boundary tracing result based on  5 is shown in Figure 5(h).

Figure 5
Figure 5 The process of generation binary mask  5 is based on a local subtraction image Although the local ratio and subtraction are algebraic manipulation which better preserve the boundary information, there are some limitations.Firstly, they are not truly local illumination invariant.If edges of long casting shadows remains in a bounding box (), algebraic manipulation based methods are prune to have wrong estimations.Moreover, algebraic manipulation is better at boundary identification for object change rather than texture change; the latter type of change does not have a clear boundary to be extracted from the background.

Figure 6
Figure 6 Daily sun illumination variation in azimuth and zenith angles of the two dataset For the Snowdonia dataset, test objects were placed on the model to simulate change on terrain surface.There are eight areas with appearance difference between the two images of three types: linear object, 3D object and camouflage.The change detection accuracy is assessed by both precision and recall  =   +   =   + and the change detection results of one image pair taken at 11:45 and 16:00 are demonstrated in Figure 8(b).

Figure 7 Figure 8
Figure 7 Change detection accuracies by the proposed IICD algorithm using the Snowdonia dataset ESP_012820_1780 and ESP_016644_1780 is presented here in detail.The two images were taken under quite different illumination conditions and have large SNR difference.The final detection results locating the Opportunity Rover in different positions in the three image pairs are shown in Figure 11.The final rover positions estimated by this approach are marked with red circles.In the three cases, although under different illumination conditions, irregular geometric distortions and different SNR levels, the positions of the Opportunity Rover have been correctly detected.The experimental results demonstrate the effectiveness and robustness of the proposed subtle object change detection approach under substantial illumination variation and geometric distortion.The comparison results using state-of-the-art change detection methods and our method for Mars rover localization are demonstrated in Figure12.As can be seen in Figure12that only the proposed IICD method is able to achieve high values both in precision and recall.

Figure 11
Figure 11 Rover Detection results using three stereo pairs

:
In §4.1.1,the motion shift of every pixel has been determined by PC based pixel-wise dense matching, and thus if the candidate area is assigned to 'motion' tag, the related object motion can be solved at sub-pixel accuracy.As indicated in §3.1, the value of  peak presents the matching quality of PC, thus the motion value calculated with large  peak value usually is more accurate in sub-pixel precision.Thus the estimated motion in x and y directions,   and   , are calculated by weighted average of disparity values in disparity maps   and   , and the weights are determined by the  peak value.
where (, ) is the value of the pixel (, ) in  peak value map , and   (, ) and   (, ) are calculated pixel motion values of the pixel (, ) in disparity maps   and   .