1. Introduction

The orthorectification process assumes that each image pixel has an accurate ground coordinates represented in the digital surface model (DSM). Given the ground X, Y, and Z coordinates of the object point and the sensor parameters collinearity condition equations can be applied to determine the image coordinates of the object point. After reaching the image position the intensity value of this image point is assigned to its corresponding location in the orthoimage.

Most of orthoimage generation problems originates from the non-corresponding representation between image and DSM. The first problem originates from the non-complete digital surface modeling. Traditional orthoimage production uses digital terrain models (DTM) representing the topography of terrain without representing the detailed man made features on the ground. These man-made features cannot be rectified correctly and the top and bottom of these objects are projected as two points still displaced from each other as shown in figure 1 on the left hand where a DTM is considered. Starting from point A on the ground and applying collinearity condition we will reach to the image point of the building wall and then assigning its intensity value g(a) to the position of point A on orthoimage plane. All intensity values that would be assigned to points from A to B represent the wall of the building. Applying back projection of terrain points from B to C lead to image points of the roof of the building. Thus, the wall appears in the resulting orthoimage and the building roof does not appear in its right location. To overcome this problem a complete representation of the surface is required, in which all surface features such as trees and man-made features are represented. Referring to the right hand side of the figure 1 the terrain and building are represented in the DSM. Points before E are terrain points that when projected to the image plane will correctly lead to image points related to ground. Also, points from E to F are roof points that when projected will lead to image points related to the roof. As noticed there are no walls represented in the orthoimage generated with DSM, which means that all surface features are rectified correctly.

Figure 1. The conventional DTM on left hand side and the right handed is DSM with building represented

Ghost effect in the form of double mapped areas at areas with sudden relief is the second problem. This serious problem is the major obstacle to the interpretability of the orthoimage generated. Orthorectification process performs back projection of DSM through collinearity condition to determine the image position for the object ground point. Due to perspective projection, at areas of sudden relief there would be a hidden area by that relief. This leads to incorrectly doubling of intensity values in the vacancy hidden by this feature appearing in the image as shown in figure 2. Starting from the object point A on the roof of the building and applying collinearity condition leads to image point a. The intensity value g(a) will be assigned to its corresponding location on the orthoimage plane. Also, object point B on the roof of the building will be projected to image point b. The intensity value g(b) will be assigned to the corresponding location in the orthoimage. Object point C is a terrain point on the same ray of object point B and hidden by point B as it is nearer to the perspective center. Projecting point C to image plane leads to image point b. Thus, the intensity value g(b) of point b will incorrectly assigned to the object coordinate of point C. This leads to a duplication of the assigned value g(b) that is referred to as double mapping effect or ghosting effect. The ghost of the building appears in the orthoimage within the area from B to C, the area obscured by the building itself, in the form of a double mapped roof. All true-orthoimage strategies focus mainly on identification of these hidden areas where the double mapping effect occurs, marking them as blank and finally filling them from neighboring images.

Figure 2. Double Mapping (Ghost Effect)

To conclude, all the problems for true-orthoimage generation originates from the non-correspondence between the information sampled in the image with these represented in the DSM. If there are features represented in the image but do not exist in the DSM, they would not be rectified correctly. On the other hand, if there is information about the surface recorded in the DSM and hidden in the image by an object’s relief the double mapping effect would appear in the areas obscured by the object.

2. Literature Review

As discussed above the main artifact in orthoimage generation in large scale urban areas where high buildings obscure the ground or hide each other is the double mapping or ghost effect. True orthoimage generation methods focus on how to find the areas hidden on the image by other higher objects and mark them as blank and try to fill them from other views if possible. Amhar et al. [10] proposed a methodology for making true orthoimages using DTM and digital building models (DBM) in which the building heights are only represented. In this methodology, two orthoimages, one corresponds to the terrain while the other corresponds to the buildings, are independently generated. The final true orthoimage is created by combining the terrain and building orthoimages.

The Z-buffer method [11-14] resolves the ambiguity of which object point should be assigned to the image location by considering the distances between the perspective center and the object points in question. Among the competing object points, the closest point to the perspective center is considered visible while the other points are judged to be invisible in that image. This method was sensitive to the sampling interval of the digital elevation model (DSM) relative to the ground sampling distance (GSD) of imaging sensor.

A polygon-based method was proposed by Kuzmin et al. [15] for the detection of obscured areas for true orthoimage generation. In this method, conventional differential rectification is first applied. Afterwards, hidden areas are detected by projecting polygonal surfaces, which are generated from a DBM, onto the image plane.

Habib et al. [17] proposed two new methodologies for occlusion detection in which the presence of occlusions can be discerned by sequentially checking the off-nadir angles to the lines of sight connecting the perspective center to the DSM points along a radial direction starting from the object space nadir point. Another method for occlusion detection is the Multi-Visibility Analysis of the object [18]. This method performs an analysis of the radial angular visibility from the terrain point, in the nadiral direction of each image in which this point appears.

Projective methodology for occlusion detection was discussed in Volotão [20] and Wang and Xie [21]. This method adopts the way of projecting roof and wall polygons to the ground horizontal planes using direct collinearity equations, and obtains the occlusions by iteration strategy. Oliveira and Galo [22] presented a method for occlusion detection based on height gradient computations applied to a DSM of the region. These height gradients computed in radial directions are important for the identification of the beginning of the occlusions in these directions. The final limits of the occlusions are obtained from the projection of these initial points in the DSM.

Yong Hu et al. [27] proposed an elevation buffer (e-buffer) technique for occlusion detection employing the plain elevations instead of the distances from perspective center by z-buffer, This method does not depend on the perspective center so it is applicable for all types of imaging sensor models. The drawback of this method is the sensitivity to the direction of the building roofs with respect to the imaging line of sight.

After occlusion detection and ghost effect elimination these methods combine multiple images captured from different points of view so that occluded areas in one image can be filled by images from other views [1-5]. In general, the quality of any orthoimage depends on image resolution, accuracy of the estimated camera calibration and orientation parameters, and density of the applied DSM [6].

It is worth mentioning that elimination of the ghost effect in these methods is done by finding the parts of orthoimage that is double mapped. This is performed by comparing the DSM candidate cells that have the same intensity values by measuring distances or angles from perspective center or by measuring gradients within the DSM itself. In this research, the main idea of the proposed method for true orthoimage generation is to use a pixelwise DSM generated from images itself filtered from occlusions and hidden areas, i.e., has empty cells in these regions and then projecting it to image plane. The empty or no data DSM cells will be skipped in orthorectification process and there is no intensity values would be assigned to these cells, leading to no double mapping. To generate such DSM a very accurate and dense matching algorithm should be applied. In computer vision community there is a suitable matching algorithm referred to as Semi-Global Matching (SGM) developed by Hirschmüller [9]. This matching algorithm tries to match all pixels and imposes restrict conditions that filters out occlusions and outliers. Hirschmüller [16] used the pixel-wise SGM stereo technique for occlusion detection, subpixel refinement and multi-baseline matching. The mutual information-based matching cost has been applied to compensate for radiometric differences of images. The results showed that SGM is suited for reconstruction from huge aerial and pushbroom images, in case of subpixel accuracy.

Gehrke et al. [19] compared the SGM method against LiDAR in terms of DSM/DEM resolution and accuracy. The results showed that SGM can be used as an alternative to LiDAR for high resolution DSM generation and orthoimage generation. Haala [8] used the SGM stereo method for a pixel-wise matching of highly overlapping aerial imagery to generate dense 3D point clouds. The results confirmed that SGM has efficiently increased the accuracy and reliability of the obtained results.

Michael et al. [23] applied a stereo vision system based on non-parametric rank transform and SGM. The system covers the entire stereo vision process and includes noise reduction, rectification, disparity estimation and visualization. The results confirmed that the system is applicable in terms of frame rate and high resolution. Dall'Asta and Roncella [24] compared the SGM and the Local Dense Matching (LSM) algorithms for DSM generation based on real and simulated stereo. The comparison showed that SGM provided much better results than LSM which produced noisy data especially in low-contrast and blurred areas. Gong and Fritsch [25] presented a comprehensive study about the DSM generation from WorldView-2 stereo satellite imagery. The Bundle Block Adjustment based on bias-compensated Rational Polynomial Coefficients (RPCs) has been applied for image orientation and rectification. In terms of robustness and precision, the results have verified that very dense DSMs can be generated from high resolution satellite imagery.

d’Angelo [26] applied the more global matching (MGM) technique to improve the regularization step of the SGM based on a Pleiades Triplet and a Cartosat-1 Stereo pair. The results showed that the MGM cost aggregation slightly increased the accuracy of the obtained results.

3. Study Area and Data Used

The area to be processed has a size of 1.2 x 1.3km in the city of München. The imagery captured by the DMC II 230 camera. The aerial images were collected at height above datum of 2225m and a ground sampling distance (GSD) of 10 cm. The block consists of three strips with 5 images each as shown in figure 3. The available overlap of 80% in flight and 80% cross flight, results in a considerable redundancy of up to fifteen-folded object points. Table 1 presents the exterior orientation parameters of the block images. The rotation angles between image and object space is defined as follows: Omega: X as Primary; Phi: Y as Secondary; and Kappa: Z as Tertiary. Table 2 represents the interior orientation parameters for the camera used as indicated in the camera calibration file.

Figure 3. A block containing three strips of imagery

Table 1. The exterior orientation parameters of the block images

Table 2. The interior orientation parameters of the camera

The image data set is a collection of fifteen 16-bit RGB images with a pixel size of 10 cm. Figure 4 is a typical view of these image data with building walls and the roofs are shifted from the bases due to the reliefs of the buildings.

Figure 4. Typical view of the raw image data

4. Methodology

Orientation of the Images within the block is a fundamental prerequisite for any metric reconstruction from these images. Thus, the operational procedure for true-orthoimage generation involves three phases; aerial triangulation starting from the known exterior orientation and camera calibration parameters, automated dense image matching for DSM generation, and true-orthoimage generation from the generated DSM.

4.1. Aerial Triangulation

2D affine transformation has been applied to transform the pixel coordinates (the origin is at the upper left corner of the image) to film coordinates (the origin is the principal point) using the following equations:

(1)

(2)

Six transformation parameters (a₀, a₁, a₂, b₀, b₁ and b₂) define the scale, shift and rotation between the two coordinate systems. Once the transformation parameters are calculated, the pixel coordinates can be transformed to film coordinates.

Aerial triangulation is then performed using a bundle block adjustment based on collinearity equations (3), (4).

(3)

(4)

Where f is the focal length, x_o and y_o are the principal point coordinates The elements m₁₁ through m₃₃ define a 3x3 rotation matrix derived by applying sequential rotations of: omega about the x-axis; phi about the y-axis; and kappa about the z-axis. X₀, Y₀, Z₀ are the perspective center coordinates. And, X_p , Y_P, Z_P are the object point coordinates.

The input exterior orientation parameters (X₀, Y₀, Z₀, omega, phi and kappa) are applied as fixed values for each image in the block.

Once the observation equations are formulated, they can be solved utilizing an iterative least squares solution for determination of ground point coordinates.

4.2. DSM Generation by High Density Image Matching

The orthorectification process assumes that each image pixel has 3D ground coordinates. All problems related to large scale orthoimages in urban areas originate from the non-correspondence between the information recorded on the image with these on the DSM. If there are more details recorded on the image than these represented in the DSM, these features will not rectified correctly. Also, if there is more information about the surface represented in the DSM and hidden in image by reliefs of tall features, the double mapping effect will appear in the resulted image. Thus, to overcome all the above problem we need a detailed digital surface model with the same cell size as the ground sampling distance (GSD) of the image pixel and produced from these images themselves (i.e. pixel wise DSM).

With the advances in digital airborne camera technology automatic image based 3D surface construction is highly motivated. Thus, high density image matching algorithms have been introduced. Dense image matching algorithms obtain the corresponding points for almost all pixels in the image rather than defining interest points or features and trying to match them. They compare the overlapping image row by row which assumes that epipolar geometry is determined first. Then, for each pixel they search in the corresponding row for the pixel with best match. The disparity (or parallax) is stored for each pixel in a disparity map. Finally, the 3D position is computed by the traditional photogrammetric techniques.

Semi-global matching is anew matching algorithm developed by Hirschmüller [28, 29]. The main idea of the algorithm is a computation of per-pixel matching costs utilizing the radiometrically robust mutual information and disparity map. Mutual information is the measure of similarities between two random variables based on their probability distribution. It is a measure of how much information is communicated in one random variable with another. If one random variable indicates something about the other then these two variables share mutual information. If the mutual information between the two variables is zero then these variables are independent. The advantage of using mutual information over intensity based matching algorithms in image matching is that they are not sensitive to illumination and recording changes between the two images. Also, the use of a correlation window in intensity based matching algorithms causes blurred effects on object boundaries.

Mutual information is defined from the entropy H of two images as well as their jointed entropy

(5)

The entropy of a random variable is a measure of the "unpredictability" of random variable. It is defined by a probability distribution P of the intensities of the two images.

(6)

(7)

For each image, entropy is computed from the probability distribution of its intensities, i.e. from a histogram. Similarly, the joint entropy is computed from the 2D histogram plot that indicates how the DNs of image I₁ relate to DNs of image I₂ (An ideal match of identical images is indicated by a diagonal straight line of all histogram bins in the 2D plot). According to equation (1) this leads to a pixel based MI (mi) and the matching cost CMI defined by

(8)

Semi-global matching minimizes the global energy in horizontal, vertical, and diagonal directions by aggregation of costs along 16 directions as indicated in figure 5.

Figure 5. Paths from all directions r

The costs-path Lr(p,d) of the pixel p at disparity d in direction r is computed recursively as follows:

(9)

Where P₁ is a penalty, which is added if the disparities differ by one, a larger constant penalty P₂ is added to all pixels in the neighbourhood of p, if the disparities differ by more than one. D is a set of all possible disparities. Then the costs S are summed up over paths in all 16 directions r.

(10)

A disparity map (parallax map) D_b for the base image (I_b) is computed for each pixel p which corresponds to the minimum matching cost, i.e. min S(p,d). A quadratic curve is fitted through the neighbouring costs for sub pixel determination. Let it be the pixel q in the match image (I_m). Another disparity D_m is computed by the same costs for the pixel q of the match image. The determination of occlusions and outliers is determined by comparing each D_b with its corresponding disparity of Dm. the disparity is set to be invalid (D_inv) if both differ by more than 1 pixel.

(11)

This consistency check enforces the uniqueness constraint, by permitting one to one mapping only.

4.3. True-Ortho Image Mapping

The orthorectification process takes the raw aerial imagery and applies the created DSM refined from occlusions and outliers as indicated in equation (11) and the obtained aerial triangulation results to create a true-orthoimage. Relief displacement has been corrected by determining the equivalent position in the image for each pixel in the DSM. By resampling of the surrounding pixels, intensity value can be determined. Figure 6 illustrates how to find the intensity values of a true-orthoimage. In the figure, P is the ground point; P₁ is the image point; O is the perspective center (origin); X and Z are the ground coordinates in the DSM; and f is the focal length.

Figure 6. Finding the intensity values of a true-orthoimage

5. Results and Discussions

At the beginning, the 2D affine transformation has been applied to transform the pixel coordinates to film coordinates. The obtained six transformation parameters (a₀, a₁, a₂, b₀, b₁ and b₂) defining the scale, shift and rotations between the two coordinate systems are given in table 3. Once the transformation parameters are determined, the pixel coordinates can be transformed to film coordinates. After that, aerial triangulation has been performed using a bundle block adjustment based on collinearity equations; yielding a fraction of a pixel as a measure of global quality.

Table 3. Affine coefficients from file (pixels) to film (millimeters)

The SGM penalties, P1 and P2, have been determined based on the statistics of the input data. Images were matched in both directions, from left to right and from right to left, and matching results have been compared. The disparity maps have been reprojected into a DSM in UTM Zone 32 North with a grid spacing of 10cm. Disparity segments smaller than 1 pixel have been removed. Figure 7 is a typical example showing part of the generated DSM using SGM algorithm where the black regions with no data value are the pixels that voided in the reverse matching. The dense built up area result in occlusions especially for surface parts close to the façades. Thus, visibility was limited and indicates mismatches, and respective results are voided according to SGM processes during DSM generation. The black regions indicate null values, the mismatches near building walls where the results are voided.

Figure 7. The generated DSM with the 10cm cell size

At first, occluded areas have been detected and then the true-orthoimages have been produced. The pixel size of the obtained orthoimage is selected to be similar to that of the original image. A smaller pixel size can oversample the original image. On the other hand, any further enlargement from the original image will not improve the image details. In this regard, the final obtained GSD is 10cm which is sufficient for mapping with scale 1:500. Figure 8 indicates the classical orthoimage generated by traditional methods using DEM in which we can see leans and unrectified objects and double mapped areas (ghost effect). Figure 9 shows the true orthoimage with vertical up walls and without lean and double mapped areas. We can see straight sharp edges of building roofs. The roads and sidewalks are completely visible. The areas with null value in the DSM are filled with blank image content.

Figure 8. Orthoimage generated by traditional methods with leans, walls, unrectified objects and serious ghost effects

Figure 9. True orthoimage generated using DSM generated with semiglobal matching algorithm

In order to get a robust decision about the performance of the proposed approach, experiments were carried out using the entire dataset. Figure 10 represents samples of the obtained true-orthoimage, evenly distributed through the block. The samples have been selected in complex built-up areas. No double mapped areas or walls can be observed. On the other hand, edges are very sharp due to the complete correspondence between the generated DSM and the images. These results confirm that the relief displacements have been largely removed, and the occluded areas as well as the ghost effect have been avoided in the obtained true-orthoimage.

Figure 10. Different views from the generated true-orthoimage. Left column: the position of the image in the block marked in red colour; right column: the obtained true-orthoimage

6. Conclusions

In order to generate a true orthoimage for a large scale urban imagery we presumes that each pixel in the imagery has a correct 3D ground coordinates. Thus, a very precise and high detailed DSM is required. Implementing high density image matching algorithms such as semi-global matching for DSM generation produces high accurate digital surface models with cell size in the same size of digital image ground sampling distance. Due to restricted matching criterion the matching results of occluded areas are voided. Therefore, almost all occluded areas and shadows disappear from the orthoimage generated using these digital surface models. This saves the hard and time consuming tasks of occlusion detection and elimination procedures. The drawbacks of this method are the processing time and capabilities for generating these high accurate digital surface models especially in dense built up urban areas. As a future work, a global optimization algorithm has to be tested in order to improve the regularization performance of the SGM algorithm.

ACKNOWLEDGEMENTS

The authors would like to thank the data providers for their generous support.