Multi-source image registration

GTC-Geocoding with DEM in map coordinates

Geocoding of SAR images in ground-range / azimuth coordinates.

The GAMMA SAR Geocoding and Image Registration Software (GEO) is a collection of programs designed to support the SAR geocoding between range-Doppler coordinates and map projections and the multi-source image registration. Geocoding of satellite SAR images from all the current available platforms (ERS, JERS, RADARSAT, ENVISAT) is supported by the software. Multi-source image registration is used to improve the geocoding accuracy and facilitate the combined interpretation of different image sources (e.g. geocoded SAR images at different frequencies, incidence angles and orbits, and geocoded satellite SAR and optical images).

Geocoding is the coordinate transformation between the coordinates of an imaging system, in this case range-Doppler coordinates of the SAR, and orthonormal map coordinates. Geocoding is necessary to combine information retrieved by the imaging system (e.g. the SAR image and products derived from it) with information in map coordinates (e.g. a digital elevation model, a land-use inventory, geocoded information from optical remote sensing, etc.). The position of a SAR image point (pixel) is on a specific circle around the velocity vector of the sensor (intersection of Doppler cone and sphere for constant range). The intersection of this circle with the ground surface together with the knowledge that the SAR is either right- or left-looking determines the target location and the look vector.

We distinguish between ellipsoid corrected and terrain corrected geocoding. The corresponding results transformed from SAR to map coordinates are Geocoded Ellipsoid Corrected (GEC) and Geocoded Terrain Corrected (GTC) products. While the generation of GEC products does not require a digital elevation model (DEM), this is necessary for the generation of GTC products. The DEM can either be provided in map projection or in SAR coordinates (i.e. as derived from interferometry). Surface topography has indeed an effect on the look vector.

The three main steps required in geocoding are (i) the initial determination of the geometric transformation, (ii) the refinement of the geometric transformation, and (iii) the resampling of image data sets from one coordinate system to the other.

In order to describe the geometric transformation a lookup table is used. For each image point defined in one coordinate system the lookup table contains the corresponding coordinates in the other coordinate system. Real valued numbers (as compared to integer numbers) are used to optimize the geocoding accuracy. For each pixel the lookup table contains a complex number (fcomplex format) with the real part corresponding to the x (resp. column / range / easting / longitude) coordinate and the imaginary part corresponding to the y (resp. row / azimuth / northing / latitude) coordinate.

Two types of lookup tables are used. The first type of lookup table has the dimension of an image in SAR coordinates. For each pixel (of the SAR image) the lookup table contains a pair of real valued numbers which are the corresponding map coordinates (i.e. the real valued row and column numbers of the corresponding target/location in the map projection). The second type of lookup table has the dimension of an image in map coordinates. For each pixel (of the map) the lookup table contains a pair of real valued numbers which are the corresponding SAR range-Doppler coordinates (i.e. the real valued SAR image row and column numbers of the corresponding target/location).

We call the geometric transformation from the coordinates in which the lookup table is given to the coordinates of the lookup table values forward transformation. In the forward transformation the input data is of the same dimension (number of columns and rows) as the lookup table. The interpolation necessary to determine the transformed values at the output pixel locations is based on irregularly spaced data. The program geocode supports this forward transformation.

We call the geometric transformation from the coordinates of the lookup table values to the coordinates in which the lookup table is given backward transformation. In this case the output data is of the same dimension (number of columns and rows) as the lookup table. The interpolation necessary to determine the transformed values at the output pixel locations is based on regularly spaced data; as a consequence backward geocoding is computationally more efficient than the forward transformation. The backward transformation is generally preferred for its simpler, and therefore better controlled and more efficient, resampling. The program geocode_back supports this backward transformation.

For the resampling step it is recommended that the spatial resolution of the image remains about the same. A practical reason for this is ensure that the sampling is adequate for the interpolator. But there are also other reasons. In resampling SAR images from a higher resolution to a lower resolution the signal noise level is kept unnecessarily high. It is better to first apply multi-looking (or filtering) to reduce the spatial resolution to a similar value as required for the resampled output image. To select an output spatial resolution much higher than the input spatial resolution strongly reduces the geocoding efficiency without much information gain. In such a case the geocoding can be conducted to a lower resolution followed by an oversampling to the required higher resolution.

Due to the special SAR imaging geometry effects such as layover and shadowing occur. In the case of layover the ground surface interacts more than once with the circle of possible target locations, i.e. several locations contribute to the signal. There are various possibilities to deal with this problem. With the GAMMA software the users has the possibility to select one of the following options: (a) set areas of layover and shadow as missing values; (b) keep the actual value, i.e. assuming the same SAR values to different positions of the map; (c) interpolate across regions of layover and shadow.

The typical geocoding processing sequence is summarized in the flow chart given below.

A wide variety of projections is used by the mapping agencies. In order to be compatible with many of the existing map projections different projections are supported, as listed in the table below.

EQA |
Equiangular coordinates (latitude/longitude) |

UTM |
Universal Transverse Mercator Projection |

TM |
Transverse Mercator Projection |

OMCH |
Oblique Mercator Projection (Switzerland) |

LCC |
Lambert Conformal Conic (France) |

PS |
Polar Stereographic Projection |

SCH |
SAR Coordinates used at JPL |

PC |
Polyconic Projection |

AEAC |
Albers Equal Areas Conic |

The specific coordinate system used for an application is specified by what we call the DEM/MAP parameter file. Different projections require different parameters. The exact definitions and formats of the individual parameters are specified in dem_par.

Supported formats for Digital Elevation Models are (4-byte) floating point and (2-byte/short) integer.

Programs to interpolate a DEM to a user defined spatial resolution (dem_trans) and to display it as a shaded relief (disshd, disdem_par, rasshd) are included. To display a DEM and assess its quality a shaded relief can be generated. The example shows the UTM DEM (shaded relief) for Death Valley.

No digital elevation model is required for ellipsoid corrected geocoding.

After the definition of the map geometry through the creation of a corresponding DEM/MAP parameter file (using the program create_dem_par) the geocoding lookup table (*.map_to_rdc) is determined using the program gec_map with a constant terrain height. The resulting lookup table is used for the (backward) gecoding from SAR to map coordinates (geocode_back) and (forward) geocoding from map to SAR coordinates (geocode).

It is possible to refine the geocoding lookup table based on a list of manually selected control points. Typically control points are selected in the SAR image and on a map with with the program gcp_ras. Such a list can be used to derive an offset registration polynomial which is then used to update (refine) the geocoding lookup table (create_diff_par followed by offset_list_fitm, and gc_map_fine).

Ellipsoid corrected geocoding of SAR images in ground-range / azimuth coordinates is supported by the program gec_map_grd.

For terrain corrected geocoding a digital elevation model has to be available. If a DEM in map coordinates is available geocoding requires the following steps:

DEM/MAP parameter file creation using the program create_dem_par. Derivation of initial geocoding lookup table (map geometry to range-Doppler geometry, *.map_to_rdc) and simulation of a SAR intensity image in map geometry based on the DEM and the SAR imaging geometry using the program gc_map.

Resampling of the simulated SAR intensity image from map to range-Doppler SAR coordinates using the program geocode to conduct a forward transforamtion based on the initial geocoding lookup table.

Fine registration of the simulated intensity image with the real SAR image intensity (create_diff_par followed by init_offsetm, offset_pwrm, and offset_fitm. The offset registration polynomials in the DIFF parameter file are then used to refine the geocoding lookup table (gc_map_fine).

As an alternative approach, a list of control points may be manually selected using gcp_ras from the SAR image and a map and used to determine the offset registration polynomials in the DIFF parameter file (create_diff_par followed by offset_list_fitm).

The resulting refined lookup table (*.MAP_to_RDC) is used for the (backward) geocoding from SAR to map coordinates (geocode_back) and (forward) geocoding from map to SAR coordinates (geocode).

Terrain corrected geocoding of SAR images in ground-range / azimuth coordinates is supported by the program gc_map_grd.

The initial geocoding lookup table is calculated based on the heights in MAP coordinates and parameters describing the SAR imaging geometry. Data in non-deskewed or deskewed range Doppler geometry may be geocoded. The lookup table is complex valued and has the dimension of data files in MAP coordinates with each value pair corresponding to the corresponding SAR coordinates (real part = column, imaginary part = row).

The initial lookup table is not completely accurate because of limitations in the knowledge of the SAR imaging geometry (the orbit positions are only known with a certain accuracy). As a result geocoding with the initial lookup table leads to geocoding errors of the order of a few pixels, depending mainly on the accuracy of the orbit data. In a refinement step registration offsets are determined between corresponding locations in an original data set and a data set brought into the same coordinate system as the original one. One way to determine such registration offsets (i.e. geocoding errors of the initial lookup table) is visual interpretation of the geocoded SAR image and a topographic map.

In order to avoid the often tedious and labor intensive extraction of ground control points by an operator the following more efficient (and probably more accurate) method was included. A SAR intensity image is simulated based on the DEM and the SAR imaging geometry. The registration offsets between the simulated SAR image and the true SAR image are then determined automatically using cross-correlation analysis.

Either approach (manual control point selection or automated cross-correlation analysis) results in registration offsets which are then modeled as a bi-linear function of range and azimuth which is then applied to the initial geocoding lookup table to result in the refined geocoding lookup table (representing the complete geometric transformation, i.e. initial geocoding and fine registration). The refined geocoding lookup table is then used to geocode data from one coordinate system into the other coordinate system in a one step, with the advantage of just a single resampling and interpolation. The initial lookup table is no longer used.

If no DEM is available, a constant height value may be used instead of a DEM. In other words the generation of Geocoded Ellipsoid Corrected (GEC) products is included.

Forward transformation, i.e. the geometric transformation from the coordinates in which the lookup table is given to the coordinates of the lookup table values is included. In this case the input data is of the same dimension (number of columns and rows) as the lookup table. The interpolation necessary to determine the transformed values at the output pixel locations is based on irregularly spaced data.

Nearest neighbor or various other interpolation algorithms can be selected for data resampling.

The forward transformation of real valued data (4 byte float format, examples are intensity, coherence map, height, etc.), complex valued data (pairs of 4 byte float format, examples are complex interferograms and complex differential interferograms), data of unsigned char format (examples are the phase unwrapping flag file and the layover and shadow map), and image data in SUN rasterfile format (examples are visualization products, classification results).

Forward transformation example: DEM resampled from UTM coordinates to SAR geometry for Death Valley. One color cycle corresponds to 512m of altitude.

Backward transformation, i.e. the geometric transformation from the coordinates of the lookup table values to the coordinates in which the lookup table is given, is included. In this case the output data is of the same dimension (number of columns and rows) as the lookup table. The interpolation necessary to determine the transformed values at the output pixel locations is based on regularly spaced data. Therefore, backward transformation is much faster than forward transformation. Nearest neighbor, spline, spline-log, bilinear, and bilinear-log interpolation algorithms can be selected for data resampling.

Backward transformation is supporte for real valued data (4 byte float format, examples are intensity, coherence map, height, etc.), complex valued data (pairs of 4 byte float format, examples are complex interferograms and complex differential interferograms), data of unsigned char format (examples are the phase unwrapping flag file and the layover and shadow map), and image data in SUN rasterfile format (examples are visualization products, classification results).

Backward transformation example: geocoded ERS SAR intensity image for Death Valley (in UTM coordinates)

In order to provide maximum flexibility a program to invert a geocoding lookup table is included. The "inverted" lookup table has the dimension of the target coordinates of the input lookup table. If a DEM in MAP coordinates is used to generate a lookup table *.MAP_to_SAR then the inverted lookup table SAR_to_UTM has the dimension of the SAR data files with each complex value containing the corresponding map coordinates.

This program allows to have full flexibility in geocoding. Backward transformation is much faster than forward transformation; therefore inversion of the lookup table may help to save processing time. Nearest neighbor resampling (mandatory for SUN rasterfile and unsigned char formats) may perform better using the inverted lookup table because the resampling and interpolation necessary in the lookup table inversion may have removed noise.

The DEM products listed in the table below are calculated from the DEM in map coordinates and parameters describing the SAR imaging geometry.

sim_sar |
simulated SAR intensity image in map geometry |

u,v |
spherical angles of local surface normal |

&theta |
local incidence angle |

&psi |
projection angle |

pix |
pixel size normalization factor |

ls_map |
layover & shadow map |

The local incidence angle, &theta, is defined as the angle between the local surface normal, n, and the negative of the radar look vector &rho. The projection angle, &psi, is defined as the angle between the local surface normal, n, and the image plane normal, m, the later defined by the (&rho x v) with v the sensor velocity vector. As an effect of surface topography the true pixel size, that is the area contributing to the signal represented in one pixel, is not constant. The pixel size normalization factor is the factor used to normalize the backscattering for the terrain induced variation of the pixel size. The spherical angles u and v of the local surface normal describe the orientation of the local normal vector. We distinguish between true layover, defined as areas with &psi > 90.0, and layover consisting of true layover plus the other areas with slant ranges within the slant range interval of the true layover (at the same azimuth).

The imaging geometry with the different angles and vectors used is shown below.

An example for the pixel size normalization factor (in SAR coordinates) is the shown in the image below for the Death Valley test site. Death Valley: pixel size normalization factor in SAR coordinates

A simulated SAR intensity image can be generated based on a DEM in map coordinates and used to automate fine registration as described above.

Death Valley: SAR image intensity simulated from DEM

Co-registration of two data sets facilitates the combined interpretation of the data. As a part of the software, tools to co-register data sets of similar but not equal geometry to identical geometry are included. This includes SAR data of different frequencies, incidence angles and orbits or geocoded satellite SAR and optical data of similar resolution (e.g. LANDSAT TM, ASTER). Automated fine registration based on real valued or complex valued data and manual fine registration based on operator selected tie points are supported.

In the geocoding sequence, this registration step is required for the fine registration of an image which was geocoded based on the initial geocoding lookup table but does contain slight offsets due to imperfect orbit data.

The co-registration of real or complex valued data to identical geometry consists of two main steps. In the first step the registration function is determined. In the second step the registration function is applied in the resampling of the data to the common geometry.

In the automated fine registration approach the registration function is determined in an initial offset estimation followed by a cross correlation analysis for a large number of image segments, and the determination of bilinear registration offset polynomials.

The manual fine registration approach is more flexible in the way that an "analog" reference data set such as a topographic map may be used. In a first step an operator selects tie points and determines registration offsets for the tie points. These values are then used to calculate the bilinear registration offset polynomials.

Once the registration function is known both real and complex valued data sets can be resampled to the reference geometry. For the interpolation of real valued data the user may select between a sinc interpolation algorithm and a nearest neighbor algorithm. The nearest neighbor algorithm is preferred for the resampling of real valued data sets which include gaps, such as the unwrapped phase image which usually includes areas which are set to the NULL value because no unwrapped phase was obtained in the phase unwrapping. For most other cases the sinc interpolator resampling algorithm is preferred. For the interpolation of complex valued data nearest-neighbor, sinc, and biquadratic spline interpolation algorithms are included.

© Copyrights for Documentation, Users Guide and Reference Manual by
Gamma Remote Sensing, 2002.

UW,
CW,
TS,
last change 5-Nov-2002.