Geometric transformations and registration of images, orthoimage generation and mosaicing

Research Collection

Educational Material

Geometric transformations and registration of images, orthoimage generation and mosaicing

Author(s):

Baltsavias, Emmanuel P.

Publication Date:

2000

Permanent Link:

https://doi.org/10.3929/ethz-a-004659029

Rights / License:

In Copyright - Non-Commercial Use Permitted

Geometric transformations and registration of images, orthoimage generation and mosaicing

E.P. Baltsavias

Institute of Geodesy and Photogrammetry, ETHZ Zurich

Geometric Transformations Consist of 2 steps:

1.Transformation of pixel coordinates 2.Interpolation of grey values

Both steps together often called “resampling”. Grey value interpolation not always necessary, e.g. rotation by 90

.

Applications:

Change of the images (zoom, rotation etc.), correction of geometric errors,

sensor calibration, transformations from one system to another (image to

image/map/vector tranformation), rectification, transformation in epipolar

images (for stereo display, matching), orthoimage generation, matching, ge-

neration of synthetic 3-D scenes (animation, simulation), ...

Examples of geometric operations

• Mirroring, transposing, rotation, shear, zoom-in/zoom-out (with or without grey level interpolation)

• Different polynomial transformations (normallly of 1st - 3rd degree). Number of terms: (n+1)

, with n ...degree of polynomial. Important ones:

- Similarity transformation (Helmert)

X = V

+ M cos ( φ ) x + sin ( φ ) y ; Y = V

- sin ( φ ) x + M cos( φ ) y V

, V

...x-, y-shift ; M ...x-, y-scale ; φ ...rotation angle

- Affine transformation

X = a

+ a

x + a

y ; Y = b

+ b

x + b

y

a

, b

...x-, y-shift; a

, b

...x-, y-scale; a

, b

...x-, y-shear - Bilinear transformation

X = a

+ a

x + a

y + a

xy ; Y = b

+ b

x + b

y + b

xy

Other transformations:

- projective (8 parameters, plane to plane transformation), e.g. for rectificati- on of fassades

X = (a

x + a

y + a

a

x + a

y + 1) Y = (a

x + a

y + a

a

x + a

y + 1)

- perspective (s. collinearity equation), e.g. for orthoimage generation

Representation of transformations with homogeneous coordinates

and X = u / w ; Y = v/w

a

, a

...x-, y-scale ; a

, a

...x-, y-shear ; a

, a

...x-, y-shift a

, a

...x-,y- “principal point” coordinates (normally = 0) a

... “focal length”, i.e. global scale (normally = 1)

Advantage -> combination of transformations by matrix multiplication, e.g.

Transformation T (= T

x T

x T

) instead of each transformation T

separately u

v w

a

a

a

a

a

a

a

a

a

x y z

×

=

H ...Grey value of background (user defined, or fixed in programme)

In geometric transformations with grey level interpolation the indirect me- thod is almost always used

(transformation T) direct method

(inverse transformation T

) indirect method

H

With large data sets with a regular grid structure an update philosophy can be used -> Transformation is much faster. Example:

1 2 Affine Transformation

X

= a

+ a

x

+ a

y

X

= a

+ a

x

+ a

y

; Mit x

= x

, y

= y

+ D

X

= a

+ a

x

+ a

(y

+ D

) = X

+ a

D

= X

+ A

D

with A = constant -> instead of 2 additions/2 multiplications only 1 addition

Grey Level Interpolation

• Polynomial of nth degree -> (n+1)

terms (unknown coefficients) ->

(n+1)

observations (known grey values of neighbours) are required n = 0 -> nearest neighbour interpolation (NINT operation), 1 neighbour n = 1 -> bilinear interpolation, 2 x 2 neighbours

n = 3 -> bicubic interpolation, 4 x 4 neighbours

• Other interpolation methods: Lagrange polynomial

Lagrange: 4 x 4 neighbours, similar quality as bicubic but faster

• In interpolation methods with negative coefficients (e.g. bicubic) the result

can be outside the range [0 , 255] -> Normalisation needed -> additional

• The indirect method is used

direct method

indirect method

Example of computation of the coefficients

• Bilinear interpolation

g = a + b x + c y + d xy with

a = f(0,0) , b = (f1,0) - f(0,0) , c = f(0,1) - f(0,0) , d = f(1,1) + f(0,0) - f(0,1) - f(1,0) g

1,0

0,1 1,1 0,0

local coordinate system of grey values f(i,j) y

x

• Computing time aspects

Bilinear needs more time than nearest neighbour, bicubic significantly more than bilinear

• Reduction of computing time by

- Update philosophy (when interpolated points are a densification of the known points and on a regular grid, e.g. DTM densification)

- Use of precomputed and saved Look-Up Tables (LUTs) for computation of the

coefficients (especially for the bicubic interpolation)

• Quality aspects

Bilinear is clearly better than nearest neighbour (staircase effect at straight ed- ges), bicubic a bit better than bilinear (image is more sharp, but also artifacts ->

grey level oscillations (“ringing”) at high contrast edges). In most cases bilinear interpolation is used.

Nearest neighbour is sometimes used in Remote Sensing -> original grey values

are kept-> important for multispectral classification

Registration of Image to Image/Map/Vectors

• Sometimes also called overlaying or rectification

• Other definitions

Georeferencing: the geometric relation between image and object space (ge- nerally a map coordinate system) is known, but the image has not been trans- formed in the object space system

Geocoding (rectification): the image is transformed in a map coordinate sy- stem (normally with a polynomial transformation)

Orthoimage generation (differential rectification, orthorectification): strict transformation in a map coordinate system with help of a DTM (relief correc- tion)

• Rectification without DTM used often with satellite imagery

• Classes of rectification

- parametric, i.e. flight/orbit, sensor etc. parameters are used

- nonparametric, i.e. normally 2-D polynomial transformations

Nonparametric rectification

• Measurement of control points in image and the reference (image, map, vectors, GPS)

Requirement: well identifiable and measurable in image and reference

• Measurement of x,y coordinates in the reference (easier than below) - manually (cursor, digitiser)

- GPS

- use of a database with known identifiable points

• Identification and measurement of x,y pixel coordinates in image (major problem) - manually

- semiautomatically, e.g. using matching between two images (second one is refe-

Difficulty of identifying and measuring control points in images

- independently whether they come from topographic maps, GPS etc.

- difficulty depends on image resolution. See example below (left: Landsat TM (bands 1,

2 and 3 combined as B/W) 30 m ground resolution, middle: SPOT Pan orthoimage 8.33

m, right: Ikonos GEO 1m, all in the region Luzern, Switzerland). Problems increase as

ground resolution becomes coarser.

Examples of GCP definition (identification quality) in Ikonos GEO images (Nisyros volca-

nic isalnd, Greece). Top:good definition ; bottom: poor. GCPs came from GPS field mea-

surements.

Example of good GCP definition in Ikonos GEO images (Melbourne). The round-abouts

were measured indirectly through measurements along the round-about perimeter. In the

field with GPS, and in the images with manual measurements and ellipse fit with least

squares adjustment, or least squares template matching (see next two pages).

Measurement of GCP (round-about) with best-fit ellipse.

• Estimation of a polynomial transformation from the reference in the image (indirect method), generally with Least Squares Adjustment (LSA)

• Degree of polynomial (generally 1-3) depends on terrain form and terrain area co- vered by image, e.g. smaller degree when terrain flat and area small

• After LSA the residuals are examined and eventually some points remeasured or deleted

• Grey level interpolation

• Important cases:

- Image to map (or orthoimage), image to image (comparison, operations between the images, multitemporal analysis)

- Often used with satellite images

• Factors influencing accuracy: number and quality of control points, degree of transformation, terrain roughness

Image is transformed to the average plane (height) of the control points, i.e. image

regions with heights far away from this plane have large image errors

Orthoimage Generation

• Orthoimage -> orthogonal projection, image with map geometry

• Input data: DTM, digital image, sensor orientation (interior, exterior), transformation from photo- in pixel-coordinate system (only with analog images), optionally other correction values (lens distortion etc.)

• Generation method (coarse description)

- Choice of the orthoimage grid spacing W (usually as W = D/n), D ...DTM grid spacing, n ...integer > 1)

“Rule of thumb”: P/2 < W < 2P, P ...pixel size of original image in object space called also pixel “footprint”

If W >> P, subsample the original image (what other input data must be changed?)

- Choice of the orthoimage region (normally in object space)

digital image

transformation from photo- to pixel- coordinate system (normally affine) projection center

orthoimage regular grid

• With non area sensors (e.g. linear CCDs) the respective geometric sensor mo- dels are used

• With satellite imagery the strict geometric transformation is often approxima- ted by polynomials -> 3-D to 2-D transformation (save computing time!)

• Orthoimage grid spacing is usually smaller than the DTM spacing. Two options exist:

- DTM densified by bilinear interpolation, then all orthoimage grid points trans- formed in pixel coordinate system (strict solution)

- DTM points treated as anchor points and transformed in pixel coordinate sy- stem. The remaining orthoimage grid points get pixel coordinates through bi- linear interpolation using the known pixel coordinates of the DTM points

Anchor point method is faster and results in similar accuracy as the strict solu-

Left: digital input image ; right: orthoimage. Black squares represent

anchor point locations. Note that the regular grid on object space is distor-

ted in image plane. All white squares represent orthoimage grid points

whose pixel coordinates are interpolated using the known pixel coordina-

tes of the anchor points.

• Geometric corrections (usually not made)

- Lens distortion (very small in modern aerial cameras) - Atmospheric refraction (never used with satellite images) - Earth curvature (only when the imaged area is large)

Estimation of the influence of these corrections on the accuracy by using approximative formulas -> decision whether necessary or not

• Other important factors

- Scale of the input imagery (often 2 - 3 smaller than scale of imagery)

- Pixel size used in scanning the imagery (depends on desired data size, re-

solution of the sensor and the film, application). Often in range 20 - 50

• Planimetric accuracy of the orthoimage - under good conditions -> less than 1 pixel

- important influence factors: DTM errors (especially at the image borders), exte- rior orientation accuracy (depends on control point quality and algorithm for its computation, e.g. by resection or bundle block adjustment).

Geometric scanner errors are also possible.

• Non-DTM objects (buildings, bridges, trees etc.) - are radially displaced in orthoimages

- can be corrected a posteriori, or

combination of DTM with 3-D description of the visible surface of these objects and then orthoimage generation

- reduction of radial displacement by using camera with 30 cm lens and/or crea-

ting orthoimage only in the central image part

image

terrain

m

...scale factor of

orthoimage enlargement

Above formulas valid for vertical imagery

∆ρ ∆R c h ---

× ∆Z ρ

--- h

×

= =

∆r ∆ Z r h ---

×

=

• Computing time

- very little in comparison to the time needed for acquisition of data required for orthoimage generation, batch run possible

- Example: 14 min CPU time, image and orthoimage 100 MB, Sun Sparc 20

Planimetric Accuracy Control of Orthoimages

• Use of well defined control points. Possible control sources:

- maps and large scale plans (are they accurate enough?) - GPS (are the point well identifiable?)

- photogrammetric determination using stereo pair

- derivation of control points using ortherectified stereo pairs (treated later) - use of existing higher accuracy orthoimages

• Control should cover whole image format and all existing terrain types

• If orthoimage produced using a DTM (not DSM), control must lie on the ground

• Analysis of the residuals and plot can reveal systematic errors (e.g. in sensor orien-

Mosaicking

• Prerequisite: images must be in the same geometric reference system

• Errors in mosaicking

- geometric -> sometimes covered up by “rubber-sheeting”

- radiometric

- colour balance of the 3 channels in colour orthoimages

• Radiometric correction

- based on equalisation of global statistical functions of the grey values (cu- mulative histogram, mean and standard deviation etc.) in the overlap regi- ons

- Grey level transformation for achieving this equalisation determines grey

level corrections. Corrections in all overlap regions of an image are combi-

ned and applied via LUTs

• Overlap regions used in the radiometric correction can be determined auto- matically (e.g. a rectangle) or manually (e.g. a closed polygon line)

In latter case aim -> exclude regions with differences between the images, e.g. because of clouds

• User determines the priority of the images -> which is up (visible) ? or

gives manually the seam line (polygon line) between the two images with aim to minimise geometric and radiometric differences, e.g. choice of the seam line in regions without texture or along sharp edges

• A linearly weighted averaging of the grey values of the overlapping images

vertically to the seam line reduces radiometric and covers up geometric dif-

ferences

Example of radiometric correction for image A using the cumulative histo- gram (see following figure). Procedure:

• Find the overlapping regions of image A with all other images (in this case region 1 with image B, region 2 with image C, and region 3 with images B and C). For each region, do:

A B

C

1

2 3

• Compute the cumulative histogram of all images in the overlap region. Defi- ne a target cumulative histogram (e.g. the histogram of one image, or the average of the three histograms)

• Compute transformations (corrections) to the histogram of each image, so that the actual histogram is transform to the target histogram. These correc- tions are saved in a LUT with 256 entries and output for 8-bit grey values.

• For each image do:

combine the LUTs to compute the final correction LUT. The combination can be a weighted average using as weights a function of the area of each over- lap region from which the LUTs were computed.

• Read the images from disk in the appropriate sequence, apply the correction

LUT and save them.

Problems in radiometric correction of colour images

• If each R,G,B channel is corrected separately, then when combining the cor- rected channels new unnatural colours appear (colour shift)

• A possible solution is to transform RGB into another colour space, e.g. IHS

(intensity, hue, saturation), then do radiometric correction to the intensity

channel (or if necessary to all), and then transform back to RGB.

Example of mosaicking of false colour IR images (Maggia valley, CH) without

(left) and with (right) radiometric correction

Geometric correction of mosaics

• Errors mainly due to sensor orientation, but locally also due to DTM, scanner

• Generally only aesthetic improvement possible, or more accurate re-estima- tion of the sensor orientation -> to avoid re-estimation, estimate orientation from bundle adjustment with dense and acuurate tie points

• Reduction of misregistration errors by measuring identical points in overlap

region and computing a local tranformation between the two images. The er-

ror can be divided between both images (then both transformed), or only

one is transformed. The transformation should not be applied to the whole

image, but starting from the seam line it should fade out after a certain di-

stance from it.

Accuracy evaluation of mosaics

• If original orthoimages not available, check planimetric accuracy as with or- thoimages. Also visual control along seam lines (if visible) for radiometric and geometric differences. For latter check straight strong edges vertical to the seam lines.

• If original orthoimages available

- check radiometry by comparing statistics of grey values (e.g. histograms) from both images in the overlap region

- check geometry by measuring identical points in both orthoimages and

checking the difference of their map coordinates Xi, Yi. Measure points in a

Xo, Yo ...map coordinates of center of upper left pixel s

, s

...orthoimage pixel size in X and Y direction

X

= x

× s

+ X

Y

= – y

× s

+ Y

Orthoimage Maps

• Useful when maps do not exist or are not available, or are not accurate and/

or up-to-date, and when a fast and/or economic map production is required

• Have been often created from satellite imagery (esp. from SPOT in SE Asia, Eastern Europe etc.)

• Orthoimage maps particularly suitable for tailored specific thematic applica- tions, e.g. forestry, hiking, touristic maps etc.

• Basic elements of a topographic orthoimage map:

orthoimage mosaic, overlaid contours, map grid and tick marks, names, le- gend

• Additional useful elements: transportation lines, rivers and lake outlines

• Topographic features can be digitised from the mosaic or imported in vector form from other systems (photogrammetric, GIS)

• Orthoimage maps can be also created by combination with scanned map layers

• Orthoimage (maps) is printed in a larger scale than the original imagery (usually 2 - 5 times enlargement). The printing pixel size cannot exceed ca.

150 µ m, else it is visible by human eye. The enlargement factor of the or- thoimage map with respect to the original image (assuming that the or- thoimage grid spacing = pixel footprint) is given by the ratio

(printing pixel size) / (scan pixel size of original image)

Correction of 3-D objects (buildings etc.)

From top left clockwise: Orthoimage with no corrected build-

ings ; correct position of buildings ; buildings corrected, but

old information partially remains (occluded areas, not visible

original image

orthoimage

stereomate

Heavy lines represent corresponding distorted

and orthorectified locations.