FFTMPOLY

Description

FFTMPOLY extracts ground control points (GCPs) using image to polygon matching. Note that this function is designed to work with high-resolution data; using low-resolution data such as Landsat may not produce the desired results.

Parameters

Name	Type	Caption	Length	Value range
FILE *	String	Input file name	1 - 192
DBIC	Integer	Input raster channels or layers	0 -
DBINVC	Integer	Inverted input raster channels or layers	0 -
DBIW	Integer	Input image window	0 -
FILV *	String	Input vector file	1 - 192
DBVS *	Integer	Input vector segment or layer	1 - 1
DBGC	Integer	Output GCP segment or layer	0 -
SARIMAGE *	String	SAR image	2 - 3	YES, NO Default: NO
LASC	Integer	Last segment created	0 - 1
REPORT	String	Report mode	0 - 192	Quick links

* Required parameter

Parameter descriptions

FILE

Specifies the name of the file containing the image. The image in FILE must be in a format supported by GDB and contain enough georeferencing information that it can be transformed to the projection of the polygon data specified by FILV. The file must be writable, and must support the creation and update of GCP segments. To allow sufficient room for matching windows, the image must have at least 512 pixels and 512 lines.

DBIC

Specifies the list of image channels used in the matching process. If DBIC is defaulted, at least one DBINVC channel must be specified. The polygons should appear darker than their surroundings in the DBIC channels.

The same channel cannot appear in the DBIC and DBINV lists.

DBINVC

Specifies the list of inverted image channels used in the matching process. Each channel is inverted relative to its maximum. If DBINVC is defaulted, at least one DBIC channel must be specified. The polygons should appear brighter than their surroundings in the DBINVC channels.

The same channel cannot appear in the DBIC and DBINV lists.

DBIW

Specifies the raster window (Xoffset, Yoffset, Xsize, Ysize) in the image within which the GCPs are collected. If DBIW is not specified, the entire layer is used by default. Xoffset, Yoffset define the upper-left starting pixel coordinates of the window. Xsize is the number of pixels that define the window width. Ysize is the number of lines that define the window height.

FILV

Specifies the name of the file containing the polygon data. The polygon data must be in a format supported by GDB and contain enough georeferencing information that it can be transformed to the projection of the image in FILE.

DBVS

Specifies the vector segment or layer containing the polygon data used for matching.

DBGC

Specifies the existing GCP segment from the image file (FILE) that receives the newly collected GCPs. If DBGC is specified but does not exist, is not a valid GCP segment, or cannot be updated, FFTMPOLY throws an exception. If no valid GCPs are collected, the DBGC segment is not modified.

If DBGC is not specified, a new GCP segment is created and its number is returned in the output parameter LASC. The new segment name is 'GCR GCPs', and its description is 'GCPs from FFTMPOLY'. If no valid GCPs are collected, the new segment is not created.

GCP identifiers are composed of hash values of the raw image and reference vector files, respectively, and the consecutive GCP numbers. Each hash value is derived from the corresponding file base name and its creation date. Hash values are not guaranteed to be unique, but are likely to be distinct in most practical situations. The hash values for the raw and reference files are printed in the report.

FFTMPOLY produces GCP IDs in the following format:

  XXXXXXXX_YYYY_N

where:

XXXXXXXX = the hash code of the input image
YYYY = the hash code of the reference image
N = the number of the GCP, 1 to any number

For example:

1BE489CD_F23A_001
1BE489CD_F23A_102

The collected GCPs are always written to the standard report file or device.

SARIMAGE

Specifies whether the input channel or layer contains SAR image data. The two supported values are YES and NO. If YES is specified, the internal processing parameters are tuned for the radar image characteristics. If NO is specified, the parameters are tuned for the optical image characteristics. The default value is NO.

LASC

The existing value of LASC is always overwritten. If a new GCP segment is created, the segment number for the new segment is returned in the output parameter LASC. If the segment cannot be created, FFTMPOLY throws an exception.

LASC is set to DBGC if valid GCPs are collected when FFTMPOLY is run, and successfully stored in an existing segment.

LASC is set to 0 if no valid GCPs are collected when FFTMPOLY is run, regardless of the setting of the DBGC parameter.

REPORT

Specifies where to direct the generated report.

Available options are:

TERM: generates a report on the terminal (default)
DISK: generates a report on file "IMPRPT.LST"
OFF: turns off report generation
<filename>: appends a report to the specified file

Details

FFTMPOLY extracts ground control points (GCPs) by matching an optical or SAR image with rasterized polygons.

FFTMPOLY automatically collects GCPs by matching patches in the image with polygons, usually lakes. The matching is performed in the spatial frequency domain using the Fast Fourier Transform (FFT) to transform and match image patches and rasterized polygons. The matching algorithm is based on the Kuglin and Hines (1975) paper in the "References" section. The number of extracted GCPs depends on the number of polygons in the area used for GCP extraction. The GCPs can then be used in automated image orthorectification workflows.

The input image can be in the raw (image acquisition) geometry, but it must contain nominal georeferencing or a sensor-specific math model. All formats supported by GDB can be used in FFTMPOLY. The file must be writable, and its format must support the creation or update of GCP segments. Note also that FFTMPOLY is designed to work with high-resolution images; low-resolution data such as Landsat may not generate desired output.

The DBIC parameter (InputRaster port) specifies the list of image channels used as is in the matching process. The DBINVC parameter (InvertedInputRaster port) specifies the list of inverted channels used for matching. At least one channel or layer must be specified between the DBIC parameter (InputRaster port) and DBINVC parameter (InvertedInputRaster port). The polygons should appear darker than their surroundings in the DBIC channels (InputRaster layers), and brighter in the DBINVC channels (InvertedInputRaster layers).

If only one DBIC channel (InputRaster layer) is specified, its pixels are used directly for matching. If only one DBINVC channel (InvertedInputRaster layer) is specified, its pixels are inverted with respect to the maximum value in the patch before matching.

If two or more channels are combined, the pixels in channels from DBIC (InputRaster) and the inverted pixels from channels in DBINVC (InvertedInputRaster) are averaged before matching.

When extracting GCPs from lakes in multi-channel (polarimetric) SAR images, it is recommended that you derive the total intensity at each pixel, and specify it as a single channel or layer. When extracting GCPs from lakes in multi-channel optical images, and the infrared channel is available, specify it as a single DBIC channel (InputRaster layer). Vegetation is bright in IR channels and provides good contrast for the lakes. Otherwise specify in the DBIC parameter (InputRaster port) the combination of channels or layers in which the lakes appear the darkest. The DBINVC parameter (InvertedInputRaster port) should be used for lakes only in special circumstances, for example for scenes with frozen lakes.

If DBIW is specified, GCPs are only extracted in the window area. The specified window should not be too small since the algorithm needs to shift between the predicted and actual position of the polygons in the image.

The FILV parameter specifies the source of the reference polygons used in the matching process. The polygons in FILV and the image in FILE must be reprojectable to the other's projection. The GCPs are extracted using FILV's coordinate system. All polygons in the vector segment specified by DBVS are examined and used if they intersect the bounding box of the image or of the window specified by DBIW and their size is suitable.

The extracted GCPs are stored in the GCP segment specified by the DBGC parameter. If the DBGC segment cannot be updated, FFTMPOLY throws an exception. If DBGC is defaulted, a new GCP segment is created, and its number returned in the output parameter LASC. If the new segment cannot be created, FFTMPOLY throws an exception. The details of the extracted GCPs are written to the destination specified by the REPORT parameter.

The SARIMAGE (SAR Image) parameter specifies whether the input image is SAR or optical. This information is required to set some of the internal processing parameters differently for the two image types.

All of the polygon features in the image must have the same contrast with their surroundings. If some features appear bright and some appear dark when compared against their surroundings, split the reference data into two segments: one containing polygons that represent the bright features in the image, and the other the dark. Process each segment in a separate run of FFTMPOLY and combine the resulting GCPs before further processing.

The LASC parameter is an output parameter, and its existing value is always overwritten. If valid GCPs are extracted in the run, LASC provides the number of the GCP segment containing them.

If no valid GCPs are collected in the run, a message is written to the REPORT destination, and the return parameter LASC is set to 0 (defaulted). If DBGC is specified, the corresponding GCP segment is not modified. If DBGC is defaulted, a new GCP segment is not created.

The REPORT (Report) parameter specifies the destination of the FFTMPOLY report. The report provides some details of the matching process and the list of all extracted GCPs.

The entry for each GCP has the following format:

GCP 1     2957.79    9661.26  53.99900767 -122.98764174   -34.21    -2.74  50.70

where:

GCP: keyword and space; indicates that the line contains a GCP
1: assessed GCP status: 1 = good, 0 = failed
2957.79 9661.26: the actual pixel and line position of the GCP in the image
53.99900767 -122.98764174: the ground Y, X position of the GCP in the image (in this example, geographic coordinates).
-34.21 -2.74: pixel and line shifts between the nominal and actual GCP positions
50.70: GCP matching weight (the higher the better)

FFTMPOLY assigns GCP elevations extracted from the global DEM. Please see point 2 in the Algorithm section for further details.

The matching algorithm is robust and relatively insensitive to spurious or incorrect data such as inaccurate or seasonally varying lake shorelines. It properly handles polygons with islands provided they have sufficient contrast. The match may fail if the polygons are not uniform or their contrast differs from the setting in the DBIC parameter (InputRaster port) and the DBINVC parameter (InvertedInputRaster port). For polygons representing lakes, this may be due to sun glint in optical images, for example, or waves or aquatic vegetation in SAR images.

The number of extracted and accepted GCPs depends on the sizes and distribution of the polygons, and on the contrast of the image areas corresponding to individual polygons. For well-matching images, typically fewer than 10 percent of incorrect matches are classified as valid, and this is acceptable for the intended use of FFTMPOLY.

Example

You have already imported a SAR image to a PCIDSK format and have lake polygons in another PCIDSK file. You want to collect GCPs for automatic image orthorectification. The GCPs for the full image area are extracted using FFTMPOLY.

EASI>FILE="Line_1003_SIR_C.pix"
EASI>DBIC=1
EASI>DBINVC=
EASI>DBIW=
EASI>FILV="lakes_088w59n_069w48n.pix"
EASI>DBVS=2
EASI>DBGC=
EASI>SARIMAGE="yes"
EASI>REPORT="disk"
EASI>r FFTMPOLY

In a typical scenario, GCPELEV is then used to add GCP elevations from an existing DEM, followed by GCPREFN to refine the GCPs and reject any surviving outliers. If GCPs from other sources are available, they can be combined with the extracted GCPs prior to running GCPREFN.

Algorithm

FFTMPOLY extracts GCPs by matching the image with rasterized polygon patches. Matching is performed in the Fast Fourier Transform (FFT) spatial frequency domain. The GCP extraction follows these main steps:

Parameter extraction and process initialization:

In this step, the run parameters are extracted and verified, and the FFT software is set up. It is critical that it is possible to convert coordinates between the "reference" polygon data in FILV (InputVector) and the "object" image data in FILE (InputRaster and InvertedInputRaster). The conversion does not need to be very accurate, as long as the positions agree to within about 200 image lines and pixels.
GCP candidate extraction:

The locations and sizes of the polygons determine their use as potential GCPs (GCP candidates). For each polygon, its bounding box in the image coordinates is determined. If the center of the bounding box lies outside the window selected for processing, the polygon is rejected. For the remaining polygons, their sizes are examined. A polygon is rejected if the area of its bounding box is less than 1,600 image pixels or if the line or pixel dimension of the bounding box exceeds 1,920. All remaining polygons are accepted as GCP candidates. FFTMPOLY throws an exception if there are no GCP candidates in the specified window.

Ground coordinates of each GCP candidate are set at the center of its bounding box. The elevation of each GCP candidate is retrieved automatically at this location from the PCI-supplied global DEM. If a more detailed DEM is available it can be used, by setting its full path in an environment variable PCI_GLOBAL_DEM just before running FFTMPOLY, and unsetting it after the module completion.
Patch contrast threshold derivation:

Not all image patches corresponding to GCP candidate polygons may have sufficient contrast for the matching to be successful, for example, lakes surrounded by dark forests in optical images. Therefore, a threshold for the acceptable patch contrast is derived from an approximate histogram of the image.
GCP patch extraction:

The resolution adjustment factor is derived for each GCP candidate polygon to fit it into a matching patch. The requirement is that both bounding box dimensions must not exceed 256 lines or pixels of the adjusted-resolution patch. Then, a co-located adjusted-resolution image patch (1024 by 1024) is extracted, and its contrast verified against the threshold derived in the previous step. If the contrast of the image patch is greater than the threshold, the polygon is smoothly rasterized into a 256 by 256 patch. During the rasterization process the GCP candidate elevation is used to reproject each GCP polygon vertex to image coordinates. The image patch is derived by averaging all channels or layers in the DBIC parameter (InputRaster port) with inverted channels or layers specified in the DBINVC parameter (InvertedInputRaster port). The pixels are inverted by subtracting them from the maximum value of all pixels of the corresponding channel or layer in the patch. The edges of the polygons are smoothed out to avoid step-like discontinuities.
Matching:

For each GCP candidate, the image and rasterized polygon patches are matched iteratively. A 256 by 256 image subpatch is extracted and centered on the predicted GCP position. The spatial displacement (X, Y shift) between the two patches is derived by image to image matching based on phase correlation in the 2D Fourier transform domain. The shifts are used to adjust the GCP position before the next iteration. The iterations continue until the shifts between the two iterations are sufficiently small or the iteration limit is reached. The magnitude of the last iteration phase peak provides the GCP weight.
Matching evaluation:

The matches are accepted or rejected based on the GCP weight. The acceptable weight is reduced for SAR images. Every match is assessed on its own without reference to the other matches.

Acknowledgements

FFTMPOLY is a PCI Geomatics implementation of the patchMatcher module from the GEOCOR package developed by Dr. Jack Gibson from the Canada Centre for Remote Sensing (CCRS). Produced under license from Her Majesty the Queen in Right of Canada, represented by the Minister of Natural Resources.

The implementation was carried out under the Canadian Space Agency EOADP contract 28/7003796.

The Fast Fourier Transforms required by the algorithm are computed by the code adapted from www.kurims.kyoto-u.ac.jp/~ooura/fft.html (Copyright Takuya OOURA, 1996-2001).

References

Kuglin, C.D. and D.C. Hines (1975). "The phase correlation image alignment method". Proceedings of the IEEE 1975 International Conference on Cybernetics and Society, San Francisco, California, pp. 163-165.

Guo, X. and Y. He (2008). "Mismatch of band sequences between an image and header file: a potential error in SPOT L1A products". Canadian Journal of Remote Sensing, Volume 34, Number 1, February 2008, pp. 1-4.

Environments	PYTHON :: EASI :: MODELER
Quick links	Description :: Parameters :: Parameter descriptions :: Details :: Example :: Algorithm :: Acknowledgements :: References :: Related