Producing Themes from Image Data: Registration and Warping

Simple Image Registration

Scope
Registration of image data against geo-positioned vector shapefiles in a GIS environment

Software
ESRI ArcView GIS v3.2, Win32
TIFF 6.0 Image Support

Data Inputs
Region of Interest Carp Village, Ontario, Canada
Projection UTM, zone 18, NAD 27, units meters
Image Data carp1_obm.tif
Vector Coverages villaoi Area of Interest
  villbldg Buildings
  villroads Roads
  villstream Streams
  villwater Water

Method

Output Registered image against vector data. (See Figure 1).

Figure 1: Simple Image Registration
[Simple Image Registration]

Analysis / Summary of Results

The world file must be present parallel to the image file that it corresponds to, with the identical base filename. When ArcView loads an image, it checks the directory of the image for a corresponding world file (.tfw, .jpgw, etc.). If the software does not find a world file, ArcView places the file in the upper left corner, starting form 0,0 in image co-ordinates, with no geographic relation to other themes in the view window. If a world file is found, the software registers the image into the view as per the world file contents.

The contents of the world file are the driver for correct geo-location. Some testing was done with inputting rogue values for pixel size. When setting the x pixel size to an arbitrary value such as 1440, the image is registered but is shown stretched horizontally, a result of the large x pixel spacing. Similar vertical errors occurred for the y pixel spacing values.

The input parameters are taken and an affine transformation is systematically performed for every row and column in the image.

One could prototype Avenue code for this functionality: 

' assuming other variables declared (A,B,C,D,E,F)
my_img = ISrc.Make(SrcName.Make("carp1_obm.tif"))

numrows  = my_img.GetNumRows
numcols  = my_img.GetNumColumns

for each x in 1..numrows
    for each y in 1..numcols
      x1 = (A * x) + (B * y) + C
      y1 = (D * x) + (E * y) + F
    ' other code goes here…
  end
end
This transformation is valuable and useful for web mapping navigation. One could write code to perform affine transformation calculations on an entire image with the six parameters and image rows and columns as arguments, outputting a text .map file to geographically reference/hyperlink an image on a web page/site. One could also write web applications to navigate through map images with calculations performed based on a user click of the map (through CGI programming).

Another formula, which could be implemented, given no rotation: 

' assuming other variables declared (A,B,C,D,E,F)
for each x in 1..numrows
  for each y in 1..numcols
    x1 =  C + (x * A)
    y1 =  F -  (y * A)
   'other code goes here
  end
end
As a result, ArcView georeferences the image based on the top left co-ordinates in a carpet-like fashion, pinning each pixel down to the 2-d surface view. One must be consider the accuracy of the raw image, taking into account distortions, areas of high relief, etc.

From professional experience, this method is best for geometrically corrected data/imagery. A project this author worked on delivered over 40000 aerial photo products of Canada, geometrically corrected, in which a module was written to create world files for each product. This method (in theory) requires no correction on the image data from the end-user's perspective, serving as the ground control reference data.

Warping and Image Registration

Scope Image warping / resampling and referencing using ground control points instead of world file references

Software
ESRI ArcView GIS v3.2, Win32
ImageWarp 2.0 extension
TIFF 6.0 Image Support
JPEG (JFIF) Image Support

Data Inputs
Region of Interest Kanata, Ontario, Canada
Projection UTM, zone 18, MTM NAD 83, units meters
Image Data 3050518warped.tif
Vector Coverages kanaoi Area of Interest
  kanroads Roads
  kanwater Water

Method

Selecting GCPs
GCPs were selected on the image based on locations that can be easily found and precisely located. These may include road and stream intersections, building corners, bridges, tips of lakes, and other features that are less likely to change position over time. Each point chosen corresponds to the selected point on the coverage that is already in a geographic co-ordinate system.

Figure 2: ImageWarp GUI Environment
[ImageWarp GUI Registration]

RMS Errors
After selecting GCPs, the RMS report function was run, which transforms all cells in the image according to a polynomial transformation, trying to achieve the best overall fit. The suggested rate of error is less than one map unit of the output image.

The results (RMS report) helped determine which points selected caused the most error, or skew among the distribution of GCP errors. These points were deleted, moved, reselected and the transformation was run again until the RMS error rate was acceptable.

The more control points used, assuming they are of equal quality, the better the polynomial will fit the points. A higher number may provide a more accurate general formula since the influence of any single control point is spread out over multiple points.

The RMS error was not as low as desired, and could always be lessened with reselecting GCPs, but the concept was understood, given the scope of the task.

Considering: n = (p + 1) * (p + 2) / 2

A higher polynomial application should be used only if the actual distortion being modelled requires it. The first order polynomial is usually adequate for scale change, rotation.

Figure 3 illustrates the relationship between the level of polynomial transformation and number of GCPs.

Figure 3: GCPs Required proportional to the level of polynomial fitting (n = (p + 1) * (p + 2) / 2)
[GCPs Required proportional to the level of polynomial fitting]
Table 1 shows the RMS reports based on the first order of polynomial fitting.
Table 2 shows the RMS reports based on the second order of polynomial fitting.
Table 3 shows the RMS reports based on the third order of polynomial fitting.

Image Warping The image was then resampled with Nearest Neighbour, Bi-Linear, and Cubic Convolution parameters for comparison.

The image was outputted to the same cell size as the input parameters. During some testing, it was found that the images were blockier when resampling to nearest neighbour at coarser output resolutions. Since bi-linear and cubic functions take into account the weighted average of 4 and 16 input cell centres and their values, the output image appears smoother. However, since the output resolution was identical to the input resolution, these differences were minimal. The differences between the three approaches were minimal during this procedure.

Figure 4 illustrates a flowchart of events involved in this procedure. Figure 5 illustrates the warped output image (nearest neighbour) as a basemap to the Kanata shapefiles.

Figure 4: Image Warping Procedure Flowchart
[Image Warping Procedure Flowchart]

Output Figure 5 illustrates the resampled image atop the Kanata area coverages.

Figure 5: Resampled Imagery over shapefiles
[Resampled Imagery over shapefiles]

A Note about pixel is point vs. pixel is area
It should be noted that the upper-left co-ordinates should refer to the centre of the pixel to minimize spatial error to ½ pixel at all times.

I've written up some sample ANSI C code to handle polynomial level vs. GCPs needed, and affine transformations here.

September 2000
Geospatial Analysis Home