Most of the image processing techniques have been first proposed and developed on small size images and progressively applied to larger and larger data sets resulting from new sensors and application requirements. In geosciences, digital cameras and remote sensing images can be used to monitor glaciers and to measure their surface velocity by different techniques. However, the image size and the number of acquisitions to be processed to analyze time series become a critical issue to derive displacement fields by the conventional correlation technique. In this paper, a mathematical optimization of the classical normalized cross-correlation and its implementation are described to overcome the computation time and window size limitations. The proposed implementation is performed with a specific memory management to avoid most of the temporary result re-computations. The performances of the software resulting from this optimization are assessed by computing the correlation between optical images of a serac fall, and between Synthetic Aperture Radar (SAR) images of Alpine glaciers. The optical images are acquired by a digital camera installed near the Argentière glacier (Chamonix, France) and the SAR images are acquired by the high resolution TerraSAR-X satellite over the Mont-Blanc area. The results illustrate the potential of this implementation to derive dense displacement fields with a computational time compatible with the camera images acquired every 2 h and with the size of the TerraSAR-X scenes covering 30 × 50 km2.
In the last decades, the warmer climate, together with less precipitation in the glacial accumulation areas, has resulted in a spectacular retreat of most of the monitored temperate glaciers . If confirmed in the coming years, this evolution will have important consequences in terms of water resources, economical development and risk management in the surrounding areas. To monitor glacier displacements and surface evolutions, two main complementary sources of information are available:
• in-situ data collected for instance using accumulation/ablation stakes, Global Positioning System (GPS) stations, or digital cameras installed near the glaciers to acquire regular images of specific areas such as serac falls, unstable moraines ...
• remote sensing data acquired by air-borne or space-borne sensors such as multispectral optical images or Synthetic Aperture Radar (SAR) images.
Optical data sets are often used to observe changes and allow the computation of high resolution (HR) information such as the surface elevation or glacier displacement fields during the summer [2-4], but they cannot be regularly acquired along the year and efficiently used because of clouds or snow cover uniformity. Space-borne SAR data, especially the recently lunched HR satellites such as TerraSAR-X, COSMO-SkyMed or Radarsat-2, are a new source of information which may allow global evolution monitoring and provide regular measurements thanks to the all-weather capabilities of SAR imagery. They are used to derive surface changes and velocity fields , or to detect and track rocks and crevasses .
With the increase of the sensor spatial resolution, the data transmission and storage possibilities, the use of image time series for Earth observation is facing computational challenges which can be separated into two groups: the need to develop new signal/image processing methods to extract information from huge amount of data, but also the need to improve existing robust techniques applied at the early processing stages to be able to apply them in a reasonable computation time on very large images and on large number of images to explore temporal evolution. Image co-registration is one of the first tasks to be performed to handle time series of images acquired by a sensor in similar conditions. When motion-free areas and moving features can be distinguished, this co-registration stage also provides displacement information which is useful to derive surface displacement fields. This task is often performed by the well-known correlation technique, which can be applied in different ways.
Several tools have been developed to solve the classical correlation problem. For optical imagery, a software like Co-registration of Optically Sensed Images and Correlation (COSI-Corr) [7,8] is widely used in the geoscience community. Due to its integration to ENVI, COSI-Corr is easy to use and offers classical and Fast Fourier Transform (FFT) techniques to compute correlation. However, its use for large images is limited by computation time. For SAR imagery, the well-known software called Repeat Orbit Interferometry Package (ROI-PAC)  is dedicated to SAR interferometry, but it also includes tools to solve the amplitude correlation problem. A two steps strategy has been adopted: a first global co-registration of the two images on a sparse grid, followed by the refined computation of the correlation on a regular grid. A disadvantage of ROI-PAC is that the computation time can dramatically increase with the image size and the number of correlation points in the image.
There are many different techniques developed for image co-registration [10,11]. Those based on sub-image correlation operate either in the temporal domain (the spatial domain for the 2 dimensional (2D) signal images) by directly computing the values of the cross-correlation function and searching for its peak, or in the spectral domain after the computation of the discrete Fourier transform of the two sub-images. The methods developed in the spectral domain are meant to speedup the computation using the FFT algorithm proposed with optimized implementation in signal/image processing libraries . They derive the sub-image shift either from the phase of the cross-spectrum , or by computing its inverse Fourier transform and identifying the correlation peak in the spatial domain . A basic computation of the cross-correlation in the spatial domain requires a number of operations proportional to N2, whereas with an implementation in the spectral domain, it is proportional to N log N. A speedup of the process is expected when the window size increases, with the constraint of being a power of 2 in both directions to benefit from the FFT optimizations.
Compared with the conventional implementation of the correlation in the spatial domain, the benefit of the spectral approach depends on the window sizes. An efficient implementation in the spatial domain also presents some advantages. It is more flexible since there is no constraint on window sizes, which allows the limitation of the local stationarity hypothesis to be taken into account. It has also the advantage of being more generic since it allows the choice of different similarity criteria according to the statistics of the images. Several alternatives to the conventional "cross-correlation function" have been proposed for image co-registration , especially in the case of SAR images which are affected by the speckle effect for distributed targets. The properties of the "true correlation function" in the Fourier domain cannot be transposed for more complex criteria derived for instance from a maximum likelihood approach [16,17].
In this paper, an implementation strategy of the correlation function in the spatial domain is proposed. The objective is to preserve the flexibility and the genericness of the spatial domain approach, and to benefit from the computation efficiency of parallel or distributed processing architectures which become more and more common on conventional computers. The originality of this approach is to be able to efficiently compute the disparity measure at the initial resolution and to derive a dense displacement field. To our knowledge, it is difficult to find efficient tools for such fast computation over large remote sensing images, whereas they are essential to manage the new data sets from HR sensors, time series and large scenes. The potential and the performances of this approach are illustrated on two kinds of data: remote sensing data with repeated pass acquisitions of HR TerraSAR-X images over fast moving glaciers in the Alps, and proximal sensing image time series from a digital camera installed in front of a serac fall of the Argentière glacier in the Mont-Blanc area.
This paper is organized as follows: Section 2 details the Normalized Cross-Correlation (NCC) algorithm, its optimization and its implementation, so as to obtain an efficient correlation software. In the next sections, Sections 3 and 4, the correlation software is applied to a realistic problem. Section 3 is dedicated to the computation of the displacement of serac falls in front of the Argentière glacier. The results show a set of serac displacements and highlight the impacts of the optimized software. Section 4 illustrates the computation of glacier flow by correlation of SAR images. This section confirms the results obtained with optical images and shows the impact of the master window size on the computation time. Finally, Section 5 concludes this paper and projects future work.
2 Implementation techniques for fast correlation
2.1 Similarity function
The objective of the correlation consists in finding the best match between sub-images, a slave image I' compared with a master image I. To simplify the algorithms in this paper, the sizes of images I and I' are the same and given by the number of rows Ir and the number of columns Ic. Figure 1 illustrates the algorithm and the chosen notations.
Figure 1. Schematic illustration of the correlation algorithm with used notations.
The master window M and the search window S are, respectively, defined by their numbers of rows Mr and Sr and their numbers of columns Mc and Sc. To simplify the notations and to make the presentation easier, Mr, Mc, Sr and Sc are odd. In this manner, the correlation objective is finding for each point (k, l) of the master image, the best position of the window Mk,l centered on (k, l) in Sk,l, according to a similarity function D(p, q), where p and q are the displacements of Mk,l in Sk,l. The search window definition implies that Sr ≥ Mr and Sc ≥ Mc. The best position is defined by the maximum of the similarity function for a given couple (Mk,l, Sk,l):
As Mk,l and Sk,l always depend on the position (k, l), they will be denoted by M and S, respectively, in the rest of this paper. The similarity function D(p, q) is not fixed and depends on the user needs. In this paper, the classical NCC defined by Equations 2, 3 and 4 is used.
The correlation result is the computation of for all (k, l) such that and . Thus, for each point, the result is defined by and . The values of are, respectively, the displacement on the lines and the displacement on the columns of the point (k, l), and is the cross-correlation level for these displacements, which varies between 0 and 1.
2.2 Optimized algorithm
To optimize the algorithm and to reduce the computation time, the correlation algorithm must be rewritten to highlight the computation dependencies. The first objective is to avoid re-computing an already computed value. The second one is to introduce a flow computation technique to reduce the number of operations of the algorithm. These two techniques are the well-known summed-area tables algorithms . They have been more recently used in the method proposed by Viola and Jones  for object fast detection. According to these points, the correlation equations given in Section 2.1 can be rewritten as follows:
For a given master point (k, l):
where the computation of U, V and W can depend on their previous computation. For the first master point (k0, l0) given by , no optimization can be used, thus U, V and W are computed according to Equations 2, 3 and 4:
For the points (k, l) such that k ≠ k0 or l ≠ l0, the values of U , V and W can be expressed depending on the previous point (k - 1, l) or (k, l - 1). If the point (k, l) is not on the first column--(l ≠ l0)--U can be computed like in Equation 5,
Let us note that the indices of the first column j0 and the last column jn of the current master window are, respectively, given by and . The values of Vk,l and Wk,l(p, q) can be computed in the same way. If l ≠ l0:
Respectively, for Wk,l(p, q):
If k ≠ k0 the same optimizations--Equations 5-11--can be performed using the line dependencies.
Let us note that this optimization strongly reduces the number of operations compared with a naive implementation. As the number of operations is one of the most critical criteria for the efficiency, the correlation algorithm must be implemented according to this optimization.
For the implementation, one of the main problems is the memory to be used. The input and output images can be too big to be stored in the memory, and hard drive access can be very time consuming. Moreover, the optimizations presented in Section 2.2 need memory to store the precomputed values. Thus, an important point is to manage the required memory according to the available memory to execute the correlation algorithm as fast as possible.
The common point about the optimization given by Equations 7-11 is the use of rolling vectors or matrices. A rolling vector is a vector of N + k elements where N is the common size of the vector and k is the number of slide steps. At each slide step, the head of the vector is increased by 1 and only the last element is re-computed (see Figure 2). To create a rolling matrix, a data vector of the size of the matrix plus the number of slide step is allocated. In the following example, a 3 × 3 matrix is allocated and three slide steps are planed (see Figure 3(a)). Next, each start point of lines is correctly placed to have the following matrix (see Figure 3(b)). To slide the matrix on the right, each start point of lines is incremented by 1 (see Figure 3(c)). After that, the last column can be recomputed (see Figure 3(d)). With a rolling vector or matrix, only the last column is recomputed, which is equivalent to a condition of the re-computation of Equations 7-11. Thus, these equations are rolling matrices.
Figure 2. Initial step and first slide of a rolling vector.
Figure 3. Vector model and first slide of a rolling matrix, the data illustrate the line and column numbers.
The optimizations presented in Section 2.2 can be applied using line dependencies or column dependencies. Both are necessary. In our case, a point that is not on the first column is computed depending on the point on the previous column. A point that is on the first column, except on the first line, is computed depending on the previous line. In this way, the memory corresponding to the pre-computation of two points must be allocated, one for the next point on the same line and one to start the next line.
The required memory to compute the correlation can be greater than the available memory. That is why the implementation of the algorithm must manage the computation lines block. This kind of implementation has two advantages. First, it allows the distribution of the algorithm. If N Central Processing Units (CPU) are available, the images can be split in N blocks and each CPU computes the correlation on its block. Second, if on a machine there is not enough memory to compute the correlation, the implementation computes on a first block that can be stored in the memory, saves the results and then computes the next block, and so on.
This approach can be realized due to the fact that the needed memory for each part of the algorithm can be predicted according to the previous optimizations.
This optimized implementation is available in the Extraction and Fusion of Information for ground Displacement measurements with Radar Imagery (EFIDIR) Tools under GNU General Public License (GPL). These tools can be downloaded from the EFIDIR web site (see Acknowledgments).
3 Experiments and results on digital video camera images
In this section, the performances of the implementation proposed in Section 2 are assessed and illustrated on the processing of optical images from a digital camera installed for glacier monitoring. In the literature, two types of cameras have been used to measure glacier flow: the analog and the digital cameras. Initially, the traditional analog technology has been used in [20-23]. At the beginning of the twenty-first century, digital photography development has made the glacier flow monitoring with HR digital cameras possible. Up to now, only few experiments have been reported with HR digital cameras, as for example in Greenland polar glaciers [24,25]. To our knowledge, no experiment on an Alpine temperate glacier has been performed.
3.1 Digital camera data set
Since 2007, in the framework of the ANR (French National Research Agency) Hydro-Sensor-FLOWS project, HR-automated digital cameras have been developed and installed around the Mont-Blanc massif (see Table 1). In this paper, one of the Argentière cameras is used: the camera installed at 2,300 m Above Sea Level (ASL) which is focused on "Lognan serac falls" (see Figure 4).
Table 1. Cameras installed around the Mont-Blanc massif
Figure 4. Automated digital camera installed near the Argentière glacier in the summer of 2008.
The HR-automated digital cameras installed around the Mont-Blanc massif are based on customized Leica DLux 3 and DLux 4 units. They have been heavily modified to allow a custom low-power microcontroller-based board to control any basic function, including switching on and off the camera, focusing and triggering the shutter. When the user-defined alarm condition is met, the camera triggering sequence is started and a pre-defined amount of time is provided for the camera to focus and grab the picture before power is switched off to save battery life. All functions provided by the camera manufacturer for operator handling are simulated through analog switches. A custom software allows the user to define on the field the wake up hour, time interval between images and number of images taken every day. The default configuration is to wake up at 8 a.m. local time and grab six images every day, with 2 h intervals between images.
The system grabs 16:9 HR images of 10 Mega pixels (4, 224 × 2, 376 pixels) with the same field of view over time. The angle of view of the camera is calibrated to 65°, with a width of 4,224 pixels, the angle of view of a single pixel is 0.015° (aperture angle) .
All the images are stored as HR JPEG images: this format was selected as a compromise between storage efficiency (since the cameras are running autonomously for up to 6 months without supervision) and data quality. However, the JPEG format is not compatible with the mono-band fast correlation approach presented in this paper. Moreover, the weather conditions are often extreme above 2, 000 m ASL in mountain areas such as the Alps. Wind and strong temperature variations might move the camera, as observed previously on a similar setup . In such a case, a translation, up to 4 pixels in both directions, can be observed between two images.
According to these two previous points, the digital images acquired over "Lognan serac falls" are processed in three steps:
1. The initial RGB JPEG images Ijpeg are converted in grayscale images Igray to obtain mono-band images. This conversion is processed according to the following formula:
The resulting image is typically called luminance in the digital image processing domain .
2. An initial co-registration between the images is made on the motion-free part of the images. In practice, the motion-free parts, i.e. mountains on the background, are used to perform it. This initial image co-registration on motion-free areas is realized by a translation without applying sub-pixel offsets.
3. The proposed fast correlation technique is applied on the image pair with 31 × 31 pixels master window (i.e. Mr × Mc) and 51 × 51 pixels slave window (i.e. Sr × Sc), corresponding to a maximum offset of 10 pixels in each direction. On motion-free areas, the sub-pixel offsets provide an accurate estimation of the remaining offset due to the camera instability. On the glacier, the measured offset is the sum of the displacement offset and the geometrical offset which has not been compensated for at step 2.
The correlation results obtained are illustrated by the magnitude and the orientation of the pixel offset vector in Figure 5. The values close to zero (in black) on magnitude correspond to the motion-free parts around the Argentière glacier which are well co-registered by the initial translation (step 2). The areas in purple correspond to the motion-free parts where there are remaining offset variations due to the camera rotation . The glacier displacement appears with stronger magnitude in blue, green and yellow colors. The heterogeneity is due to either the glacier flow physics or the scene configuration: the nearest parts of the glacier appear to be flowing faster than the farthest parts. The displacement map of Figure 5(b) highlights the differences between the ice blocks in the foreground (mostly in yellow), in the middle distance (mostly in green) and in the background (mostly in blue). One can notice a large ice block in green color on the right part of the blue background, corresponding to a larger displacement: this ice block is about to fall. There are also a few parts where the magnitude and the orientation maps look like noise and inconsistence. These parts correspond to ice falls which happened between the two image acquisitions.
Figure 5. Displacement field over the Argentière glacier derived by NCC between digital camera image pair (2008-10-09/2008-10-10).
3.3 Computation speedup
To highlight the effect of the optimization, the correlation is executed with and without optimization, using 1-8 CPU. The objectives of this execution are to illustrate the speedup given by the optimization and the number of used CPU. This experiment, and the following one, are computed on an octo-core Intel(R) Core(TM) i7 3 GHz with 24 Go memory. In the experiment, this machine is considered as eight independent CPU with 3 Go of memory per CPU. Figure 6 shows these results: the computation time without optimization (T) and with optimization (Topt) depending on the optimization and the number of used CPU. A first observation of Figure 6 shows that the benefit of the optimization is very important. According to the delay between two image captures (2 h in our case), it is interesting to observe that the correlation can be computed between two captures only if the optimized method is used. Moreover, the impact of the number of used CPU gives an almost linear speedup. When the number of used CPU doubles, the computation time is divided by two. Thus, the combination of optimization and distribution reduces the computation time, in our context from more than 36 h to 10 min, or even less time if more than 8 CPU are used.
Figure 6. Computation time depending on the optimization and CPU number (displayed with y-time-log scale).
Figure 7 illustrates the gain obtained by the optimization. The first curve named "absolute gain" shows the difference of computation time without and with optimization (T - Topt), for each number of used CPU. The second curve named "relative gain" shows the ratio between the previous curve and the computation time without optimization .
Figure 7. Gain given by the optimization depending on the CPU number.
From Figure 7, the relative gain can be considered constant for our experiment and it is very significant: more than 96%. Since the computation without optimization can be very long--more than 1 day--the absolute gain can change the work habits. The prospects with many days of computation are not the same as with a few hours. The computation time and the absolute benefit decrease when the number of used CPU increases, but even with 8 CPU, several hours are saved thanks to the optimization.
This first experiment highlights the benefit of the optimization and the distribution of the correlation algorithm for optical images. It is important to note that this benefit enables to decrease the interval between two image captures. Consequently, a real time glacier flow monitoring becomes feasible. With the appropriate computation system, an acceleration of the glacier and an important loss of correlation corresponding to serac falls can be quickly detected.
4 Experiments and results on SAR images
Despite improved acquisition, transmission and processing performances, the proximal sensing by ground-based optical cameras, as illustrated in Section 3, is limited to specific parts of a few glaciers. In this section, the proposed fast correlation technique is applied to remote sensing data which can cover large areas: space-borne images allow the whole glacier surface, and even all the glaciers of a mountain area, to be observed simultaneously. The feasibility in a reasonable computation time and the interest of the dense correlation measurements of this fast correlation technique are illustrated on HR SAR images which can be regularly acquired by repeated satellite passes.
4.1 TerraSAR-X data set
In the framework of the TerraSAR-X science project MTH0232 , 35 stripmap TerraSAR-X images have been acquired on the Mont-Blanc test site (see Table 2). There are three time series in descending configuration (orbit 25) and 1 in ascending configuration (orbit 154). Most images have been acquired in a single polarization mode (HH), except a winter series in the dual polarization mode (HH/HV) for the analysis of the snow backscattering properties. Ascending and descending measurements provide four different projections of the surface 3D displacement field. The combination of these projections allows retrieving the 3 components (East, North, Up) of the surface 3D displacement field .
Table 2. Temporal series of TerraSAR-X images acquired on the Cha-monix Mont-Blanc test site.
In this paper, the cross-correlation results are illustrated on the single polarization (HH) descending images which are acquired with an incidence angle of 37° and a spacing of 1.36 m per pixel in range and 2.04 m per pixel in azimuth direction. The range and azimuth image axes correspond to the radar line of sight (LOS) and the sensor displacement direction, respectively. The stripmap mode has been chosen because it supplies a large scene coverage (about 30 × 50 km2) and HR at the same time. Figure 8 shows a whole strip-map image which covers almost the whole Mont-Blanc area, i.e. French, Italian and Swiss parts.
Figure 8. TerraSAR-X amplitude strip-map image (15, 790 × 24, 183 pixels), Mont-Blanc area, 2008-09-29 (averaged by 10 × 10 for display). The blue rectangle corresponds to the sub-image chosen to illustrate glacier displacements in Figure 10.
In the mountainous areas where most of the Alpine glaciers are located, the "range sampling" of SAR images introduces strong geometrical distortions. To avoid geocoding artifacts, the SAR images of the Mont-Blanc test site have been ordered in their initial geometry. The offsets measured in range direction between two images are sensitive to the position along the swath (near rangea/far rangeb), to the topography, as well as to the surface displacement occurred between the two acquisition dates. The offsets measured in azimuth direction mainly depend on the surface displacement (a linear correction is sufficient to remove along-track registration variations over long scenes). The range variations due to the topography depend on the perpendicular baseline between the two orbits as in a stereo configuration. These variations can be predicted using a Digital Elevation Model (DEM) of the area and the orbital data (antenna state vectors) which are provided together with the images.
In the studied area, the altitude varies between 1,000 m ASL (in the Chamonix valley) up to 4,800 m ASL (on the Mont-Blanc). For the image pair (2008-09-29/2008-10-10) whose perpendicular baseline is around 138 m, the range registration offsets due to this baseline vary between 28.9 and 82.4 pixels in near and far range, respectively. The glaciers of this test site might move up to 1.5 m per day in the fastest areas, according to in situ measurements. The glacier displacements vary between 0 and 16 m in 11 days, hence 0-8 pixels with the resolution of the TerraSAR-X images used in this paper.
According to this a priori displacement information, the TerraSAR-X data acquired over the Mont-Blanc test site are processed in three steps:
1. An initial co-registration by a simple translation (without resampling) is applied by matching an area of the image located at an intermediate elevation of about 2,000 m ASL.
2. The proposed fast correlation technique is applied to the whole image with 61 × 61 pixels master window (i.e. Mr × Mc) and 77 × 77 pixels slave window (i.e. Sr × Sc), corresponding to an offset of ±16 m in each direction. On motion-free areas, the sub-pixel offsets provide an accurate estimation of the remaining offset due to the SAR geometry. On the moving glaciers, the measured offset is the sum of the displacement offset and the geometrical offset which has not been compensated for at step 1.
3. Depending on the variations of the geometrical offset along the glaciers, a post-processing step can be necessary to deduce the offsets only due to the glacier movement. The remaining geometrical offset can be subtracted using either the predictions from the DEM and the orbits, or the results of the sub-pixel correlation around the glaciers.
The correlation results obtained on the whole TerraSAR-X image presented in Figure 8 are illustrated by the magnitude of the offset vector in Figure 9. The values close to zero (in black) correspond to the motion-free areas which are well co-registered by the initial translation (step 1). The remaining offset variation due to the SAR geometry can be observed in the dark and light blue areas. The shapes of moving glaciers (Mer de Glace, Argentière, Les Bossons ...) appear with stronger magnitude in green. Some of the stronger magnitudes are due to misregistration when the correlation technique fails, in areas with strong discorrelation between the two images because of the surface changes.
Figure 9. Magnitude of the offset vector derived by NCC between two TerraSAR-X strip-map images (2008-09-29 and 2008-10-10) over the Mont-Blanc area (averaged by 10 × 10 for display).
The results obtained on moving glaciers are illustrated with the Taconnaz, the Bossons and the Bionnassay glaciers in Figure 10. The displacement magnitude and orientation show that the motion is not uniform: the velocity is higher in the center of the glacier, and two acceleration areas appear on the Bossons glacier. These results are consistent with the glacier behaviors, but there is no ground truth available since crevasses and seracs make in situ measurements too dangerous. On the higher part of these glaciers, the magnitude and orientation are very noisy: the correlation technique does not provide reliable results. A larger window size could improve the results on poorly correlated areas, but the window size cannot be very large since the displacement field is not homogeneous over the glaciers and the border discontinuity should be preserved. A flexible choice of window size is useful to find a good trade-off between reducing the "false alarms" (wrong match) and preserving the spatial resolution (displacement field heterogeneity).
Figure 10. Displacement field over the Bossons, Taconnaz and Bionnas-say glaciers derived by NCC between TerraSAR-X image pair (2008-09-29/2008-10-10).
4.3 Computation speedup
This second experiment on SAR images is realized in the same context as the experiment presented in Section 3: the same algorithms and the same CPU configuration are used. Figure 11 shows the speedup given by the optimization and the number of used CPU. These results confirm those obtained with optical images. The only difference is that the computation time is longer with SAR images. This is mainly due to the master and search window sizes that are larger for SAR images. For the whole SAR image of Figure 8, the computation of Figure 9 takes 15 h with optimization and 8 CPU against 18 days without optimization (30 times longer). As in Section 3, the gain given by the optimization is very important. Figure 12 illustrates these results.
Figure 11. Computation time depending on the optimization and CPU number (displayed with y-time-log scale).
Figure 12. Gain given by the optimization depending on the CPU number.
The relative gain is close to that obtained with the optical images: more than 96%. As the computation time without optimization is very long--many days--in the case of SAR images, the benefit can be expressed in computation days. Thus, the impacts of the optimization and distribution for SAR images are more important than for the smaller images of the digital camera.
To extend these results, another experiment is realized. The objective is to highlight the impact of the master window size on the computation time and the optimization. For this experiment the master window size is increased from 41 to 81 pixels with a step of 4 pixels. The search window size is 16 pixels larger than the master window size. The computation time depending on the optimization and the master window size is given by Figure 13.
Figure 13. Computation time depending on the optimization and the master window size.
On one hand, this figure shows that the computation time dramatically increases when the master window size increases. In our case, when the size doubles, the computation time is multiplied by more than 3. Despite this impact, this time stays reasonable with the optimized implementation. Thus, the "best" master window size can be searched by experiments. On the other hand, the impact of the master window size is quantified in an absolute and a relative point of view as shown in Figure 14.
Figure 14. Gain given by the optimization depending on the master window size.
Let us note that the absolute gain increases with the master window size and several hours are saved. It is also important to note that the relative gain increases with the size of the master window. In other words, the larger the master window size, the more efficient the optimization.
5 Conclusions and future work
This paper details an optimized implementation of the NCC algorithm. The objective is to reduce the computation time of the correlation technique to handle large data set for Earth change monitoring. The saved time induced by the optimization has multiple impacts. The computation on each point of the image can be achieved in a reasonable time: 0.02 min/mega pixel instead of 0.4 min/mega pixel with a conventional approach. High resolution remote sensing images covering large scenes can be processed in few hours. This fast correlation technique is very useful to extend experimental researches. For example, it allows researchers to experiment different processing parameters and to analyze large data sets.
Two experiments illustrate the benefits of the proposed approach. The evolution of serac falls is studied with optical images and the whole glacier surface evolution can be observed with SAR images. On the Mont-Blanc area, the correlation reveals particular areas like glaciers, lakes or other changing features that can be studied. These experimental results highlight the potential of proximally and remotely sensed images to monitor the glacier flow and to contribute to risk assessment: the Taconnaz glacier is for instance an important source of risk for the access road to the Mont-Blanc tunnel.
Future work includes a comparison between this optimization and different implementations of the FFT approach to illustrate the advantages and limitations of those techniques. Regarding the optical images, a stereo camera will be installed near Argentière glacier to measure simultaneously the topography and the displacement of the serac fall. Regarding the SAR images, as the NCC is only one of the available similarity functions, the study and the optimization of new criteria, different from the NCC, will also be investigated.
The authors declare that they have no competing interests.
aThe portion of a radar image which is the closest to the satellite flight path.
bThe portion of a radar image which is the farthest from the satellite flight path.
The authors wish to thank the French Research Agency (ANR) for supporting this work through the Hydro-Sensor-FLOWS project and the EFIDIR project (ANR-2007-MCDC0-04, http://www.efidir.fr). They also wish to acknowledge the German Aerospace Agency (DLR) for the Terra-SAR-X images (project MTH0232) and Électricité Emosson SA for their support.
C Vincent, A Soruco, D Six, E Le Meur, Glacier thickening and decay analysis from 50 years of glaciological observations performed on glacier d'Argentière, Mont Blanc area, France. Ann Glaciol 50, 73–79 (2009). Publisher Full Text
E Berthier, H Vadon, D Baratoux, Y Arnaud, C Vincent, KL Feigl, F Rémy, B Legrésy, Mountain glaciers surface motion derived from satellite optical imagery. Remote Sens Environ 95(1), 14–28 (2005). Publisher Full Text
D Scherler, S Leprince, MR Strecker, Glacier-surface velocities in alpine terrain from optical satellite imagery-accuracy improvement and quality assessment. Remote Sens Environ 112(10), 3806–3819 (2008). Publisher Full Text
UC Herzfeld, GKC Clarke, H Mayer, R Greve, Derivation of deformation characteristics in fast-moving glaciers. Comput Geosci 30(3), 291–302 (2004). Publisher Full Text
E Trouvé, G Vasile, M Gay, L Bombrun, P Grussenmeyer, T Landes, JM Nicolas, P Bolon, I Petillot, A Julea, L Valet, J Chanussot, M Koehl, Combining airborne photographs and spaceborne SAR data to monitor temperate glaciers. Potentials and limits. IEEE Trans Geosci Remote Sens 45(4), 905–923 (2007)
R Fallourd, O Harant, E Trouvé, J-M Nicolas, F Tupin, M Gay, G Vasile, L Bombrun, A Walpersdorf, J Serafini, N Cotte, F Vernier, L Moreau, Ph Bolonm, Monitoring temperate glacier by multi-temporal TerraSAR-X images and continuous GPS measurements. IEEE J Sel Top Appl Earth Observ Remote Sens (2011) (to appear)
S Leprince, S Barbot, F Ayoub, JP Avouac, Automatic and precise ortho-rectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements. IEEE Trans Geosci Remote Sens 45(6), 1529–1558 (2007)
S Leprince, F Ayoub, Y Klinger, JP Avouac, in Co-registration of optically sensed images and correlation (COSI-Corr): an operational methodology for ground deformation measurements. IEEE International Geoscience and Remote Sensing Symposium (IGARSS 2007) (Barcelona, Spain, 2007), pp. 1943–1946
B Zitová, J Flusser, Image registration methods: a survey. Image Vis Comput 21(11), 977–1000 (2003). Publisher Full Text
J Gao, MB Lythe, The maximum cross-correlation approach to detecting translational motions from sequential remote-sensing images. Comput Geosci 22(5), 525–534 (1996). Publisher Full Text
H Stone, M Orchard, C Ee-Chien, S Martucci, Fourier-based algorithm for subpixel registration of images. IEEE Trans Geosci Remote Sens 39(10), 2235–2243 (2001). Publisher Full Text
O Harant, L Bombrun, G Vasile, M Gay, L Ferro-Famil, R Fallourd, E Trouvé, J-M Nicolas, F Tupin, in Fisher pdf for maximum likelihood texture tracking with high resolution polsar data. EUSAR 2010 Proceedings, Aachen, Germany, 418–421 (2010)
FC Crow, in Summed-area tables for texture mapping. SIGGRAPH '84: Proceedings of the 11th annual conference on Computer graphics and interactive techniques, New York, NY, USA (ACM, USA, 1984), pp. 207–212 PubMed Abstract | Publisher Full Text
P Viola, M Jones, Robust real-time object detection.. Int J Comput Vis (2001) [http://www.cs.cmu.edu/~efros/courses/AP06/Papers/viola-IJCV-01.pdf] webcite
AN Evans, Glacier surface motion computation from digital image sequences. IEEE Trans Geosci Remote Sens 38(2), 1064–1071 (2000). Publisher Full Text
H-G Maas, E Schwalbe, R Dietrich, M Bässler, H Ewert, in Determination of spatio-temporal velocity fields on glaciers in West-Greenland by terrestrial image sequence analysis. IAPRS, Beijing, China, XXXVII, Part B8, 1419–1424 (2008)
J-M Friedt, C Ferrandez, G Martin, L Moreau, M Griselin, E Bernard, D Laffly, C Marlin, in Automated high resolution image acquisition in polar regions. European Geosciences Union, Vienna, Austria (2008)
R Fallourd, F Vernier, J-M Friedt, G Martin, E Trouvé, L Moreau, J-M Nicolas, in Monitoring temperate glacier with high resolution automated digital cameras--application to the Argentière glacier. PCV 2010, ISPRS Commission III Symposium (Paris, France, 2010)
R Fallourd, F Vernier, Y Yan, E Trouvé, Ph Bolon, J-M Nicolas, F Tupin, O Harant, M Gay, G Vasile, L Moreau, A Walpersdorf, N Cotte, J-L Mugnier, in Alpine glacier 3D displacement derived from ascending and descending TerraSAR-X images on Mont-Blanc test site. EUSAR 2010 Proceedings, Aachen, Germany, 556–559 (2010)