| GISdevelopment.net ---> Technology ---> Image Processing |
|
Matching of IKONOS Stereo and Multitemporal GEO Images for DSM Generation
Li Zhang, Maria Pateraki, Emmanuel Baltsavias Institute of Geodesy and Photogrammetry, ETH-Hoenggerberg, CH-8093, Zurich, Switzerland Email : zhangl@geod.baug.ethz.ch, maria@geod.baug.ethz.ch, manos@geod.baug.ethz.ch
1. Introduction
Although IKONOS imagery has been commercially available since early 2000, the use of this imagery for DSM generation has been restricted due to various reasons. Firstly, up until beginning of 2002 stereo imagery was available only for governments and national mapping organisations. Secondly, Space Imaging follows a closed policy and does not release the sensor model but provides rational polynomial coefficients instead. These are provided for stereo images and special versions of Geo imagery (so called OrthoKit) but at significantly higher prices. To these reasons, one should add the general problems of IKONOS imagery like cost, availability, lack of user control on imaging parameters and conditions etc. Since summer 2001, some commercial systems (Erdas Imagine, LHS Socet Set, Z/I Imaging Image Station, PCI Geomatics OrthoEngine) support to one or the other extent IKONOS imagery for import, stereo viewing and processing, orthoimage and DSM generation by using the rational polynomial coefficients (RPCs) provided by SI, estimating RPCs from ground control or using alternative (bundle adjustment, DLT) and partly proprietary mathematical models. While estimating RPCs from GCPs is both expensive (high number of GCPs required) and sub-optimal, use of bundle adjustment due to the undisclosed IKONOS sensor model is impossible. DLT, as shown in Fraser et al. (2002) could be used with GCPs, while proprietary sensor models (PCI) are of unknown quality, while raw IKONOS data which would be most appropriate for use with a strict sensor model are not available. As shown in Baltsavias et al. (2001) and Fraser et al. (2002), a very simple alternative geometric model, needing at least 3 GCPs, can achieve better positioning accuracy than the RPCs and at sub-meter level. Investigations on DSM generation from high-resolution imagery include the following. Ridley et al. (1997) evaluated the potential of generating a national mapping database of maximum building heights for buildings at least 5 x 10 m in planimetry by using DSMs extracted by matching of 1m aerial imagery. They report that matching has a potential to provide the requested information if the DSM had a spacing of 1-3 m but with lower accuracy (1.5 – 3 m RMS) and completeness compared to manual measurements. Muller et al. (2001) use simulated 1m-resolution and IKONOS data for DSM generation and landuse determination to estimate effective aerodynamic roughness for air pollution modelling and determine position of trees close to buildings that may cause soil subsidence for insurance risk assessment. The most extensive results on DSM generation from IKONOS up to now have been reported by Toutin et al. (2001). Dial (2000) also reports the expected mapping accuracy of IKONOS. There are no published investigations on the other two operational high-resolution spaceborne systems (EROS-A1 and Quickbird-2). As far as the authors know, this is the first investigation deriving DSMs from multitemporal IKONOS images. This data are much more available than stereo data, they do not include epipolar resampling which possibly degrades geometric accuracy and are cheaper than stereo IKONOS. Naturally, multitemporal differences make matching more difficult, but as it is shown in Baltsavias and Stallmann (1992) it is possible with intelligent techniques to achieve good results even in case of large image differences. 2. Test Data 2.1 Melbourne Dataset The first dataset comprised a stereopair of epipolar-resampled Geo PAN images of Melbourne. As indicated in Table 1, the sensor and sun elevation angles for the stereopair (imaged in winter) were less than optimal. The azimuths of sensor and sun of the left image differed considerably, leading to strong shadows in non-occluded areas. The base/height ratio of the stereopair was 1.2. GCPS measured with an accuracy of 1-2 dm in both image and object space existed and could be used for our IKONOS sensor model and aerial image orientation. More details on the dataset and the GCPs can be found in Fraser et al. (2002). An existing stereopair of 1:15,000 scale wide-angle colour aerial photography was used to measure reference data. 2.2 Nisyros Dataset This dataset consisted of two Pan-Sharpened Geo images of a Greek island with little man-made objects and sufficient vegetation. It is an non-active volcano with elevation range 0-700 m. The parameters of the images are listed in Table 2. As it can be seen, they were similar except of the sensor azimuth. The different viewing angle let to illumination differences between the two images. The fact that shadows were mostly in occluded areas was an advantage. The major multitemporal differences were due to different clouds and their shadows. 38 GCPs measured with GPS were available with an object accuracy of less than 1 dm. About half of them were not well defined in the image. The reference DTM had a grid spacing of 2 m and was interpolated from 1:5,000 map contours and additionally measured break lines. Its accuracy, estimated using the GCPs, was ca. 3.3 m RMS. More details on the data can be found in Vassilopoulou et al. (2002).
3. Own Matching Method and Test Results with Nisyros Dataset 3.1 Quasi-Epipolar Image Generation Unlike frame-based imagery, where all pixels in the image are exposed simultaneously, each scan line of the IKONOS image is collected in a pushbroom fashion at different instant of time. Thus, epipolar lines with linear CCDs become curves. Using Orun and Natarajan' sensor model (Orun and Natarajan, 1994), which models the position of the sensor and its yaw angle variation as second-order polynomials of scan line (or time), and assuming the pitch and roll angles to be constant, the derived epipolar curve for a certain point shows a hyperbola-like shape. It can be approximated by a straight line for a small length but not for the entire image. An epipolar curve for the entire image can be approximated only by piecewise linear segments. In our approach, the epipolar curve for point (xl, yl) in the left IKONOS image is approximated by a quadratic polynomial and has the following form: Where, (xr, yr) are the pixel coordinates in the right IKONOS image and A1-9 are unknown parameters. Using Eq. (1), the quasi-epipolar image pair can be generated by re-arranging of the original IKONOS image pair. However, the computation of the unknown parameters of the epipolar geometry needs a certain number of well distributed conjugate points. These conjugate points are extracted by the following feature-point based matching:
3.2 Derivation of Approximations by Using Grid Point Matching based on Relaxation The matching procedure can be treated as a labeling problem and solved by the relaxation technique. That is, we regard the template image (in our case, the left quasi-epipolar image) as the model and another image as the scene, and then use the feature points in the model as a set of labels to label the feature points extracted from the scene. Relaxation is one of the efficient methods to solve the labeling problem. The work of Hancock and Kittler (1990) that uses Bayesian probability theory has provided a theoretical framework and rigorous basis for the relaxation method. The important aspect of the relaxation matching algorithm that distinguishes it from the single point matching is its compatible coefficient function and its smoothness constraint satisfaction scheme (Baltsavias, 1991; Zhang et al., 1992). This is most important for areas with homogeneous or only little texture. In such areas, the single point matching is unable to match images because of the lack of information. With the smoothness constraint, such areas can be bridged over, assuming that the terrain surface varies smoothly over the area. Firstly, the match points are selected and distributed in the form of a regular grid in the template image. Given a grid point in the template image, a search window in the search image (the right quasi-epipolar image, in our case) can be determined. The correct match of this point should lie in this search window. However, due to repetitive texture or poor texture information, there could be several candidate matches appearing in the search window. These candidate matches are located along the epipolar curve. They can be derived by traditional cross-correlation technique, and the candidate matches are selected if their correlation coefficient lies above a certain user-defined threshold (we choose this threshold value as 0.7). The approximate centers of the search windows can be derived from the matching results of the previous pyramid level. Let Ii be one of the grid points on the template image and Ij (j=1,…,m) its candidate matches in the search image. P(i,j) is the probability of match Ik¬Ii. Moreover, let Ik be one of the points located in the neighborhood of point Ii and Il (l=1,…,m) be its corresponding candidate matches. In order to link the matching results of the neighboring grid points to each other, we define the following compatible coefficient function C(i,j;k,l) which quantifies the compatibility 
between the match Ii?Ij and a neighboring match Ik¬Ii: In equation (2), Dp expresses the difference of the x-parallaxes in point Ii and its neighboring point Ik. The bigger the Dp, the smaller the compatibility. This corresponds to a smoothness constraint on the image matching results, and it provides an ability to bridge over the poor texture areas. T is a value quantified by the texture information and it is defined as inversely proportional to the minimum of four gray value variances (horizontal, vertical, and along the two main diagonals) at the image window around the point Ii. Normally, if one point is located in rich texture areas or at linear features, this value is small. This value can be treated as the weight of the smoothness constraint and it provides the possibility to control the continuity of the terrain surface. ß is a constant and its value was set to 400 empirically. In the relaxation scheme, the so-called global consistency of matching can be achieved by an iterative scheme where the probabilities P(i,j) in iteration n+1are updated by the following rule: 
C(i,j;k,l) is the compatible coefficient function defined as above, W(Ii) is the neighbourhood of point Ii (can be its 8 or 24 neighboring points), and n is the iteration number. The quantity Q(n)(i,j) expresses the support the match Ik¬Ii receives at the nth iteration step from the matches Ik?Il in its neighbourhood W(Ii). The iteration scheme can be initialized by assigning the normalized correlation coefficient to P(0)(i,j) and, ideally the process will terminate when an unambiguous match result is reached, that is when each point Ii is matched with one candidate with probability 1, the probabilities for all other candidate matches for this point being zero. In practice, we terminate the process if any one of the following two conditions holds:
3. 3 Modified MPGC MPGC is described in detail in Baltsavias (1991). It combines least squares matching and geometric constraints formulated either in image or object space. The constraints lead to a 1-D search space along an epipolar line, thus to an increase of success rate, accuracy and reliability, and permit a simultaneous determination of pixel and object coordinates. Any number of images (more than two) can be used simultaneously. The achieved accuracy is in the sub-pixel range. The algorithm also provides criteria for the detection of observation errors and blunders, and adaptation of the matching parameters to the image and scene content. Usually, the sensor model used for IKONOS is the rational function model, other models like polynomial mapping functions (Baltsavias and Stallmann, 1992) or DLT can be treated as a special case of the rational polynomial functions. For description of the RPCs see Grodecki (2000). Their form is given by: 
Where, (xn, Yn) and (Xn, Yn, Zn) are normalised image and object coordinates respectively. To compute the normalised coordinates, the following equations are used: 
Where, (x, y) and (X, Y, Z) are image and object coordinates; x0, y0 and xs, ys are offset and scale values for the image coordinates respectively, and similarly X0,Y0, Z0 and Xs, Ys, Zs are offset and scale values for the object coordinates. In the RPCs, the maximum power of each object coordinate and the total power of all object coordinates are limited to 3. In such a case and following the SI definition of coefficient sequence, each polynomial has the following form (for convenience, the subscripts are omitted), which leads to the 80 RPCs per IKONOS image: 
The geometric constraints were derived from equation (4) and have the following form: 
Equations (5) can be treated as weighted observation equations in MPGC, where (?x, ?y) are the unknown shifts in pixels, which are the common unknowns, appearing also in the affine parameters of the MPGC. Weighted geometric constraint force the matching to search for a conjugate point only along the epipolar curve. If the initial match of the point in the search image does not lie on this epipolar curve, at the first iteration of MPGC matching, it jumps onto this curve. The grid point matching results based on relaxation provide quite good approximations for the MPGC procedure and increase the convergence rate. In our implementation of the modified MPGC, the adjustment starts only with the two shift parameters and after the first convergence all affine transformation parameters are used. Also, at the first two iterations, the weight value for the geometric constraints takes a large value in order to speed up the convergence and then it decreases to consider remaining errors in the orientations/epipolar geometry and allow matching slightly off the epipolar curve. The modified MPGC matching procedure has the potential of obtaining sub-pixel matching results and exploiting more than 2 IKONOS images simultaneously. 3.4 Test Results For the Nisyros dataset, the DSM was generated for the entire island, including a part on the sea. The approach starts with a feature-point based matching to obtain conjugate points for epipolar geometry recovery and approximation values for relaxation matching. This stage provides 114 conjugate points that are refined by least squares matching without geometric constraints. The refined matching results are used for generation of quasi-epipolar images by using bilinear interpolation of the grey values. Then, the grid relaxation matching is performed on the quasi-epipolar image pair. The grid interval was 7´7 pixels and the patch size also 7´7 pixels. This step results in 852,615 grid points (493,916 left after cutting out the results of cloud covered areas). These points were then transformed to the original IKONOS images and 3D coordinates needed by MPGC were computed by spatial resection. The modified MPGC procedure is performed by taking the relaxation matching results as approximations. For testing the potential accuracy and performance of our matching approach the following 5 different versions were executed: Version 1: relaxation matching only, by cutting out the cloud covered areas manually. Version 2: modified MPGC 2D matching, patch size is 21´21 pixels and no geometric constraints, by cutting out the cloud covered areas manually. Version 3: modified MPGC matching, patch size is 21´21 pixels and with geometric constraints, by cutting out the cloud covered areas manually. Version 4: modified MPGC matching, patch size is 11´11 pixels and with geometric constraints, by cutting out the cloud covered areas manually. Version 5: modified MPGC matching, patch size is 21´21 pixels and with geometric constraints, without cutting out the cloud covered areas manually. Versions 2-5 were applied to the raw images, but after noise reduction, contrast enhancement and radiometric equalisation. In all versions, the same match area and manually selected cloud areas were used. 3D coordinates were computed using the terrain-corrected affine transformation (Baltsavias et al., 2001) as sensor model. The affine parameters for each image were derived from just 4 GCPs. Table 3 shows the accuracy test results of different versions (matching – reference DTM). The results are according to the expectations with version 3 delivering the highest accuracy of ca. 3 m RMS, a remarkable result if one considers the multitemporal nature of the images. The positive bias is probably because matching measures on top of buildings, trees etc. In Fig. 1, an image area with contours derived from the matched results of Version 3 is shown. 
 Figure 1 An image region overlaid with contours derived from Version 3 matching results and 10 m interval 4. DSM Generation and Test Results using Commercial Digital Photogrammetric Systems A DSM was generated with matching using the LHS DPW770 and VirtuoZo digital photogrammetric systems. The IKONOS images used were pre-processed with the third approach mentioned in Baltsavias et al. (2001), but without Wallis filter and for VirtuoZo with reduction to 8-bit (VirtuoZo can not process 11-bit images for DSM generation). Wallis filter removes the low frequencies, and in the image pyramid generation process, which is used by both systems, this leads to homogeneous images with little texture. The proper procedure is to generate the pyramid levels first and then apply Wallis. This could be done with LHS but not with VirtuoZo, which does not provide access to the pyramid levels, but to allow a more objective comparison of the two systems it was not used with LHS. LHS can read the RPCs of IKONOS images, and we assume that after matching these RPCs are used to compute a spatial resection, although no information is provided on how this is done. With VirtuoZo, the pixel coordinates were used to compute object coordinates using the RPCs and spatial resection. For both systems, the geographic coordinates (RPCs use geographic coordinates) were transformed to UTM, and then the RPC bias, computed from 28 GCPs, was subtracted and a regular 2m grid was interpolated. With LHS a so-called Adaptive Automatic Terrain Extraction (AATE) method was used, associated with a modified steep_1 matching strategy, which is suitable for rough terrain. Nonadaptive strategies perform an epipolar resampling of the images first. Although the IKONOS images were already epipolarly resampled, the system could not recognise this, and thus was trying to transform the images, but unsuccessfully, and thus matching with nonadaptive strategies was failing (it is unknown whether AATE performed on-the-fly epipolar resampling or not, and if yes how). With AATE, 6 pyramid levels in 8 computational passes and a patch size of 9x9 pixels at the lowest level were used to generate a regular grid with 2m spacing. LHS uses a simple 2x2 averaging to generate the pyramid levels. This fast, but qualitatively poor, method leads to a smoothing of details in higher pyramid levels. For the Melbourne city scenes, this led to images that had very repetitive signal due to the repetitive scene content with buildings and city blocks. Another filter was used for pyramid generation (7x7 approximation of the ideal sinc filter, see Baltsavias, 1991) but this, although finally used in the matching, led only to minor improvement. With VirtuoZo, the number of pyramid levels is fixed (4, including the original image, whereby pyramid levels are generated with 3x3 averaging and a reduction factor of 3). A regular image grid with 7 pixels spacing was matched, using a patch size of 11x11 pixels. The matching results were evaluated visually by superimposing contours on the stereopair (that was extremely slow and cumbersome with LHS) and DSM-shading (see Figs. 2 and 3), but no manual editing was performed. Visually, the results of VirtuoZo looked initially better (not necessarily when comparing Fig. 2 and 3, because of the different visualisation and the lower match point density of the VirtuoZo DSM), and the contours were following the objects (e.g. buildings) in most cases nicely. It should be noted that the Melbourne scene was a very difficult case for matching. Although the terrain is relatively flat, high buildings created parallaxes of over 200 pixels. Their occlusions and shadows created additional problems. The scene objects were highly discontinuous, leading to smoothed surface results, when using area-based matching as is the case with both systems. The repetitive image content in higher pyramid levels could easily lead to multiple solutions and wrong matching. Although the images had practically no y-parallax, with LHS a 1D matching in x-direction could not be enforced, while with VirtuoZo this was possible, but only because we had access to a separate, commercially unavailable, matching module.  Figure 2 Shaded 2m IKONOS-DSM generated by the LHS system with AATE. The oval University of Melbourne campus is at the top left part.  Figure 3 Shaded 2m IKONOS-DSM (a part only) generated by VirtuoZo. The oval University of Melbourne campus is at the top left. For a quantitative analysis, we used the aerial imagery, which was covering only a part of the IKONOS stereopair. A DSM was generated using matching with VirtuoZo. The colour images were scanned at a Vexcel UltraScan 5000 with a pixel size of 25 microns (ca. 0.40 m pixel footprint) and were first transformed to B/W using a weighted combination of the RGB channels, whereby the weights were a linear function of the noise of each channel, determined by the method described in Baltsavias et al. (2001). Then, they were pre-processed using the third approach in Baltsavias et al. (2001), including Wallis filtering (the pyramid levels had sufficient signal in this case). The images were oriented using 8 GCPs from the Melbourne test-field in UTM, although their distribution was not optimal, and the RMS of the orientation residuals in X, Y and Z were 0.31, 0.70 and 0.35 m respectively. GCPs were measured in one image manually and transferred to the second one semi-automatically using matching. A regular image grid with 7 pixels spacing was matched, using a patch size of 11x11 pixels, and then a 2m regular grid was interpolated. DSM matching with the aerial images had the same problems as IKONOS, with the exception of reduced repetitive patterns, due to the higher resolution of the aerial imagery. Thus, the use of this DSM as a reference dataset would not be optimal. Another problem in comparing the DSMs from aerial and IKONOS imagery were the multitemporal differences. The aerial imagery was acquired on 11/1/1993. Thus, there were many differences in buildings, roads and vegetation, while the occlusions and shadows were also different. This problem could be avoided by manually measuring in the aerial images only points that were visible in both image datasets and had remained stable (72 points). In addition, points in open areas matched with VirtuoZo were also used leading to a total number of 8627 reference points. With LHS, a Figure Of Merit (FOM) is provided to indicate suspicious points. 22% of the matched points had poor FOM (7% had high elevation slope, 7% low (less than 0.45) correlation coefficient, ca. 3% each had large signal power difference between the two images and excessive correlation shift, and ca. 2% were elevation spikes). A dataset without poor-FOM points was also used (marked “cleaned” in Table 4). As with previous DSM tests with LHS, the results without poor-FOM points are better but large blunders still remain in the presumably good points. A quantitative analysis of the DSMs gave the results listed in Table 4, whereby the reference data were interpolated in the raw matching data of both systems. 
The results of Table 4 are to a certain extent deteriorated by some remaining errors in the reference DSM or differences between IKONOS and aerial images (e.g. due to some isolated trees existing only in IKONOS imagery). Due to the latter, the bias, as also most of the blunders, is positive since the IKONOS DSMs was higher. The statistics for VirtuoZo are worse, contrary to the visual checks, but this may be due to aerial DSM errors and the partly different DSM coverage for the two systems. 5. Conclusions This investigation has shown that IKONOS Geo stereo imagery has high geometric integrity and the potential to yield accurate DSM results, even from the cheapest Geo version, without RPCs and even with multi-temporal images. Investigations on other linear-CCD space-borne sensors have shown that under ideal circumstances raw matched DSM data can be ca. 1 pixel accurate or even less for good B/H ratio. In these tests, the conditions were far from ideal and thus a lower accuracy was achieved, which is however still remarkable. An additional problem for Nisyros was the insufficient quality of GCPs and reference DTM that could not allow achievement of accuracies at 1-m level. GCPs and reference data in Melbourne would allow that, and the fact that in these results the standard deviation was close to 1-m may show the full potential of DSM generation from IKONOS. Acknowledgements The provision of imagery and ground control information by the Department of Geomatics, University of Melbourne and the co-operation with C. Fraser, H. Hanley and T. Yamakawa are gratefully acknowledged. The IKONOS image of Nisyros was made available by Prof. E. Lagios, University of Athens, Greece within the EU project Geowarn (www.geowarn.org). References
|