3 D OBJECT COORDINATES EXTRACTION BY RADARGRAMMETRY AND MULTI STEP IMAGE MATCHING

Nowadays by high resolution SAR imaging systems as Radarsat-2, TerraSAR-X and COSMO-skyMed, three-dimensional terrain data extraction using SAR images is growing. InSAR and Radargrammetry are two most common approaches for removing 3D object coordinate from SAR images. Research has shown that extraction of terrain elevation data using satellite repeat pass interferometry SAR technique due to atmospheric factors and the lack of coherence between the images in areas with dense vegetation cover is a problematic. So the use of Radargrammetry technique can be effective. Generally height derived method by Radargrammetry consists of two stages: Images matching and space intersection. In this paper we propose a multi-stage algorithm founded on the combination of feature based and area based image matching. Then the RPCs that calculate for each images use for extracting 3D coordinate in matched points. At the end, the coordinates calculating that compare with coordinates extracted from 1 meters DEM. The results show root mean square errors for 360 points are 3.09 meters. We use a pair of spotlight TerraSAR-X images from JAM (IRAN) in this article. * E-Mail:akrameftekhary@gmail.com. 1. INTRODUCTIO Generally, there are two main techniques to extract 3D object coordinate from SAR images which in both of them two or more images are taken from different points and different look angels in the space. At the first method that is interferometry SAR, different phase between two images use far extracting height information. Unlike Radargrammetric technique analyzes the signal amplitude and exploits the stereoscopy as optical photogrammetric methods. Radargrammetric theory was introduced in 1950. It was the first method used to derive DSMs from airborne radar data in 1986(Leberl, Domik et al. 1986). Because of low resolution in amplitude of radar image, it was less used than InSAR technique. Toutin research on the use of Radargrammetry techniques by RADARSAT-1 and ERS1 / 2 images and model was offered as a three-dimensional Toutin Model(Toutin and Gray 2000). Since 2007 with the advancement of technology and the ability to obtain SAR images with a high spatial resolution of about one meter by TerraSAR-X, Radarsat-2, the development took place and where not possible use optical images or interferometry techniques, Radargrammetry methods were used and good results were obtained (Raggam, Gutjahr et al. 2010) and (Toutin and Chénier 2009). Generally height derived method by Radargrammetry involves of two stages: Images matching and space intersection. By considering not very good radiometric quality for amplitude SAR images and the angles between these two images, matching that used in optical images does not provide good results. So we propose a new method for images matching in this paper. The multi-stage algorithm originated on the combination of feature and area based images matching. Range and Doppler (R&D) are physical equations that commonly used for the SAR images space intersection. Recent researches shown the generic Rational Function Models (RFM) are appropriate replace for R&D equations and the RFM are fitted with accuracy better than 10 pixels on the physical model(Zhang, He et al. 2011). Because of simplicity and high speed processing of RFMs, we use of RPCs instead of R&D for space intersection.


INTRODUCTIO
Generally, there are two main techniques to extract 3D object coordinate from SAR images which in both of them two or more images are taken from different points and different look angels in the space.At the first method that is interferometry SAR, different phase between two images use far extracting height information.Unlike Radargrammetric technique analyzes the signal amplitude and exploits the stereoscopy as optical photogrammetric methods.Radargrammetric theory was introduced in 1950.It was the first method used to derive DSMs from airborne radar data in 1986 (Leberl, Domik et al. 1986).Because of low resolution in amplitude of radar image, it was less used than InSAR technique.Toutin research on the use of Radargrammetry techniques by RADARSAT-1 and ERS1 / 2 images and model was offered as a three-dimensional Toutin Model (Toutin and Gray 2000).Since 2007 with the advancement of technology and the ability to obtain SAR images with a high spatial resolution of about one meter by TerraSAR-X, Radarsat-2, the development took place and where not possible use optical images or interferometry techniques, Radargrammetry methods were used and good results were obtained (Raggam, Gutjahr et al. 2010) and (Toutin and Chénier 2009).Generally height derived method by Radargrammetry involves of two stages: Images matching and space intersection.By considering not very good radiometric quality for amplitude SAR images and the angles between these two images, matching that used in optical images does not provide good results.So we propose a new method for images matching in this paper.The multi-stage algorithm originated on the combination of feature and area based images matching.Range and Doppler (R&D) are physical equations that commonly used for the SAR images space intersection.Recent researches shown the generic Rational Function Models (RFM) are appropriate replace for R&D equations and the RFM are fitted with accuracy better than 10 -3 pixels on the physical model (Zhang, He et al. 2011).Because of simplicity and high speed processing of RFMs, we use of RPCs instead of R&D for space intersection.

IMAGE MATCHING
Due to the presence of intrinsic speckle noise in SAR images, first use a Wiener filter to reduce noise and then the matching is done.For matching process, at first the SIFT feature extraction method is used Then the OGM method that is used to assign attributes to the features is perform.Finally, the local matching on the least squares matching technique to achieve sub-pixel accuracy is done.Details of this process are described at the follow.

Feature extraction
Generally feature extraction methods used in optical images can take in SAR images.Because of speckle noise, some of these methods are not suitable for SAR images.Therefore, if the speckle noise is reduced, some of feature extraction methods can be used.Among these algorithms Harris and Forstner are performed to point features and gradient based algorithm are used on the SAR images edge detection (Pan, Chen et al. 2010).
Scale Invariant Feature Transform (SIFT) is one of the most popular algorithm that is applied to SAR images for feature extraction.A strategy based on SIFT matching algorithm that initially feature based process for optical pattern recognition in images is presented.SIFT algorithm is independent from rotation, illumination changes, 3D viewpoint changes and affine distortion.Also it is stable against noise and illumination.Therefore it applied to feature extraction of SAR images.This algorithm consists of four major stages: Scale-space extreme detection, key point localization, Orientation assignment and key point descriptor (Lowe 2004).

Feature allocation
In the feature based matching, after the feature extraction, image feature should be allocated.These features can be a window of a neighboring pixel gray level, scale and direction of each feature and the bit set relationship between neighboring and central pixels.In general, each descriptor can be scalar value, vector, matrix and etc. that will be assigned to a specific feature.These descriptors are used in adopting features in the matching.Speckle noise and different look angle between two SAR images are main problems that cause SAR images matching difficultly.Also layover, foreshortening and shadow disturb the correlation process.Therefore the OGM is used as descriptor to allocate features for image.In the one-dimensional case, if F (x) is a filter: Then Gradient Amplitude of I(x) function is obtained from convolution I(x) and F(x).The edges occur at the maximum of the convolution: These relations can be extended to two-dimensional case.X(i,j) and Y(i,j) are masks that convolve with input image and I x (i,j) and I y (i,j) get.
If A(i,j) shows Gradient Amplitude Image and D(i,j) shows Gradient: This process reduces speckle noise and improves image structure.In this stage, although the edges are extracted, the image is smoothed and the noise is reduced (Paillou and Gelautz 1999).Therefore, if Linear Optimal Gradient is performed on the window's pixels around features, it can use as good descriptor for features.

Local Matching
Two strategies is possible for matched candidate points, they are local and global matching.In the local matching, a candidate point feature on master image matches with one of the candidate point feature extracted in slave image.This is done to other points of the two images so that each candidate point for the master image with the slave image points of similarity are considered separately.At this stage, based on the correlation of similarity between the features of locally assessed.According to equation (5) around the image window candidate features, the similarity criterion is calculated.The highest value obtained for the correlation coefficient determines the maximum similarity.It is worth noting that the image windows in A and B are the same as the original image windows which Linear Optimal Gradient operator is applied on them.Thus, the values obtained for the coefficient of linear correlation between the various permutations of the feature windows image belong to the master image; a number of points on the master are corresponding by a number of points on the slave image.

Least Square Matching
Least Square Matching similar to any other Area Based Matching (ABM) followed by an area of template in the searched image.Cost function using gray level values of the pixels between the two areas, such as the same size of the image is defined.The sum of square differences in the same area will be cost function.The cost function must be minimized in an iterative process (Förstner 1982).Because of changing in imaging condition between two images, geometric and radiometric distortion is occurred in the corresponding area in slave image (Gruen 1985).Equation ( 6) expresses this relationship.
( , ) ( , ) 1( ( , )) 1( 2( , ), 3( , )) f x y e x y T g u v T T x y T x y Where, f(x,y) is the image window and g(x,y) is the searched window.Also T1 is a linear radiometric transformation, T2 and T3 are geometric transformation.These transformations are shown in the equations 7 to 9.
In these equations, r 1 and r 2 are coefficients of the linear radiometric transformation, a i and b i are the affine coefficients of geometric transformation.By performing matching processing that is described, matching accuracy improved and it reached to sub-pixel.

SPACE INTERSECTION
In this step Rational Coefficient Polynomials are used to extract height information from SAR images pair.According to research performed by Zhang et al, Rational Function Models are a good alternative for Range-Doppler equation (Zhang, He et al. 2011).Here's a short explanation of how to calculate RPCs and their improvement will be express.The Rational Function Model (RFM) is defined as the ratio of two bi-cubic polynomials involving 39 parameters in any direction in a 3D form: (10) In this equation, a ijk represents the polynomial coefficients called rational polynomial coefficients (RPC).The variables r and c are the normalised image coordinates.We use third-order polynomials with an unequal denominator, as they are capable of modelling most of the distortion in the image space.
Consequently, 39 terms are used, 20 terms for the numerator and 19 for the denominator.Thus, there are 78 unknown RPCs to be solved in the general RF model (Zhang, Li et al. 2011).

RPC Generation
The general workflow of developing the RF model is shown in Figure 1.The procedure consists of three main processing steps.First, a grid of Control Points (CP) is established.These points can be either in real or virtual form, depending on the availability of SAR orbital information.In the second step, the RF model is solved by transforming it into a linear parametric model and determining the unknown RPCs.In this stage, the challenge is to select the appropriate regularisation method (Zhang, He et al. 2011).In the third processing step, RPCs are refined, the accuracy of the RPCs is checked and suitable models are fit to the data to improve precision.
Figure 1.General workflow of developing an RF model

RPC Model Refinement with an Affine Model
Because the RF model is calculated from the physical imaging model without the aid of ground control points, errors in the direct measurement of sensor orientation can cause biases in the RPC mapping.The biases can be taken into account by using a bias-corrected RPC model, expressed as follows: (11) In these equations, and represent the differences between the measured and the nominal line, respectively, and the pixel coordinates.These deviations can generally be described as polynomials of the image line and sample coordinates.( 12) Here, A i and B i (i = 1, 2, 3 …) are the correction parameters in the bias-correction model.In this study, four sets of correction parameters are tested, and the results are compared: (1) Case one.A 0 and B 0 , which model only the shift bias.
(2) Case two.A 0 , A 1 , B 0 and B 1 , which model the shift and time-dependent, drift bias.
(3) Case three.A 0 , A 1 , A 2 , B 0 , B 1 and B 2 , which model the bias using an affine transformation model.For more information about the physical meanings of the parameters, please see the relevant reference (Tong, Liu et al. 2010).

Datasets
In this study, we use a pair of TerraSAR-X Single-Look Slant-Range Complex (SSC) images with long baseline.The images were acquired over the city of Jam, southern Iran, in spotlight mode and descending orbit.We use a part of these images by 1000×1000 pixels.These images are shown in Figure2.The ground elevation in the study area is between 620 and 680 m.In the experiment, a total of 4 distinct ground control points in the four corners of images were measured by GPS.The accuracy of these points is less than 10 cm.These points are used to refine RPCs.The image coordinates of these 4 points were carefully measured up to a nominal accuracy of 1 pixel.Also DEM by 1 meters resolution and 0.5 meter accuracy are used to evaluate the results in this research.

Results
By applying matching method that described in the second part, 360 matched points are extracted from images.Matching accuracy is about 0.2 pixels.
In the space intersection step, firstly, the independent RPCs are calculated as method that described in Section 3.1 for both Master and Slave images.The results show that using the RF model with the independent method provides a good fit to the rigorous Range-Doppler model, with an accuracy of greater than 10 −3 pixels in image space.Then, we use an affine model for RPC refinement.In this case, the coefficients of the affine model (A 0 -A 2 and B 0 -B 2 in Equation ( 12)) are calculated by using four GCPs in the corners of images, Δl and Δs are applied to all images points and the RPCs are regenerated.At the end, refined RPCs are used to calculate 3D coordinates of matched points by Equation 11.The height coordinates of these points are shown in Figure 3.

CONCLUTION
We did Radargrammetry process to extract elevation information from a pair of long baseline TerraSAR-X images.Image matching and space intersection are two important stage of Radargrammetry.In this paper we propose a multi-step algorithm founded on the combination of feature and area based image matching.Also we used RPCs calculated for two images.Since RPCs are the best instead of R&D equations according to the previous investigation, using RPCs in terms of computational time is affordable.It is considerable that the RPCs must be refined then affine equation was used for modelling errors that there were in the RPCs.In a part with 1000 × 1000 pixels of these images 360 points are extracted with the proposed matching algorithm and 3D object coordinates are calculated for all matched points.At the end, the calculating 3D coordinates by offering method compare with true coordinates extracted DEM by 1 meters resolution.The results showed mean error and standard deviation is 0.59 and 3.05 meters for 360 points.
, ρ is the correlation coefficient between A and B windows in the master and slave images, A  and B  are the mean of gray values in the A and B windows.
Master (A) and Slave (B) images by size 1000×1000 pixels that used in this study

Figure 3 .
Figure 3. Matched Points Height extracted by proposed method in this paperAlso, we use of DEM by 1 meters resolution to evaluate calculation height in this study.By comparing computational height by proposed method and extracted height from DEM, the errors were found.Table1show statistical parameters of errors such as Mean, standard deviation (SD) and root mean square error (RMSE).
Table 1 show statistical parameters of errors such as Mean, standard deviation (SD) and root mean square error (RMSE).