Remote Sensing Image Preprocessing Guide
Remote Sensing Image Preprocessing Guide
5. Image Preprocessing
5.1 Introduction
The ideal or perfect remote sensing system has yet to be developed. Also, the terrain is amazingly
complex and does not lend itself well to being recorded by relatively simplistic remote sensing
devices that have constraints such as spatial, spectral, temporal, and radiometric resolution.
Consequently, error creeps into the data acquisition process and can degrade the quality of the
remote sensor data collected. These errors may have an impact on the accuracy of subsequent
human or machine-assisted image analysis. Therefore, it is usually necessary to preprocess the
remotely sensed data prior to analyzing it to remove some of these errors.
Image restoration is concerned with the correction of distortion, degradations, and noise
introduced during the imaging process. Image restoration produces a corrected image that is as
close as possible, both geometrically and radiometrically, to the radiant energy characteristics of
the original scene.
To correct the image data, internal and external errors must be determined. Internal errors are
created by the sensor itself. They are generally systematic (predictable) and stationary
(constant), and may be determined from prelaunch or in-flight calibration measurements.
External errors are due to platform perturbations and the modulation of scene characteristics,
which are variable in nature. Such unsystematic errors may be determined by relating points on
the ground (i.e., ground control points) to sensor system measurements. Radiometric and
geometric error are the most common types of error encountered in remotely sensed imagery.
5.2 Radiometric Correction of Remote Sensor Data
Multispectral scanning systems require calibration and data correction to provide the user with
reliable sensor data. A most important operation is to determine the relationship between the
individual detector output and the input scene radiance, and to correct the data for any detector
malfunctions or distortions.
Detector Response Error
Remote sensors use detectors that convert the sensed radiance into a voltage or brightness
value. There are several detector anomalies that can arise, including line dropouts, striping,
and/or line start problems.
If one of the detectors in the sensor fails to function during a scan, this can result in a brightness
value of zero for every pixel, j, in a particular line, i. This is often called a line dropout and may
appear as a completely black line in the band, k, of imagery. Obviously, this is a very serious
condition. There is no way to restore data that were never acquired; however, it is possible to
improve the visual interpretability of the data by introducing estimated brightness values for
each bad scan line. The first problem is to locate each bad line. A simple thresholding algorithm
1
can flag any scan line having a mean brightness value at or near zero. Once identified, it is then
possible to evaluate the output for a pixel in the preceding line (BVi-1, j, k) and succeeding line
(BVi+1, j, k ) and assign the output pixel (BVi, j, k) in the "dropout" line the average of these two
brightness values:
BVi 1, j , k BVi 1, j , k
BVijk Int
2
Sometimes, a detector does not fail completely, but simply goes out of adjustment (e.g., provides
readings perhaps twice as great as the other detectors for the same band). This is often referred
to as n-line striping or banding. The data are valid, but should be corrected (or restored) to have
the same general contrast as the other detectors per scan. First, the bad scan lines must be
identified in the scene. This is usually accomplished by computing a histogram of the values for
each of the n detectors over a homogeneous area such as a body of water. If one detector's mean
or median is significantly different from the others, it is probable that this detector is out of
adjustment. It may require a bias (additive or subtractive) correction or a more severe gain
(multiplicative) correction. Then the errors arising from the maladjusted detector are not as
noticeable.
Occasionally, scanning systems fail to collect data at the beginning of a scan line. This is called a
line start problem. Also, a detector may abruptly stop collecting data somewhere along a scan
and produce results similar to the line dropout discussed previously. Ideally, when data are not
collected the sensor system would be programmed to remember what was not collected and
place any good data in their proper location within the scan. Unfortunately, this is not always the
case. For example, the first pixel in band k, on line i (i.e., BVi, 1, k) might be improperly located at
column 50 (i.e., BVi,50, k). If the line start problem is always associated with a horizontal bias of
50 columns, it is a simple horizontal adjustment. However, if the amount of the line start
displacement is random, it is very difficult to restore the data without extensive human
interaction on a line-by-line basis. Data not recorded by the detectors can never be restored.
Atmospheric Effects
There is really no such thing as atmospheric "error" in remotely sensed data. The energy
emanating from the sky and recorded by the sensor is a true signal even when it completely
destroys our ability to measure the spectral reflectance patterns of the terrain. Nevertheless, we
often consider the deleterious effects of atmospheric scattering and absorption as being sources
of error that can minimize our ability to extract useful terrain information from remotely sensed
data.
2
Figure 5 - 1 (a) Effects of molecular (Rayleigh) scattering and sun angle on incident radiation
transmitted to the earth's surface; (b) effects of aerosol (Mie) scattering and sun angle on
radiation transmitted to the earth's surface.
Solar radiation is largely unaffected as it travels through the vacuum of space. However, when it
interacts with the earth's atmosphere it is selectively scattered and absorbed. This introduces
external (unsystematic) radiometric distortion into remotely sensed data. Atmospheric
absorption by water vapor and other gases in the atmosphere mostly affects radiation of
wavelengths longer than 0.8 µm. Thus atmospheric scattering adds brightness, whereas
atmospheric absorption subtracts brightness from landscape spectral measurements (Fig.5 - 1).
The effect of the atmosphere on spectral brightness creates problems when trying to compare in
situ measurements of spectral response and satellite sensor spectral measurements for the same
geographic area. It can also effect classification accuracy. Fortunately, it is possible to apply
atmospheric scattering adjustments to minimize some of these effects. Two techniques will be
discussed which are based primarily on the fact that infrared data (> 0.7 µm) are largely free of
atmospheric scattering effects, whereas the visible region (0.4 to 0.7 µm) is strongly influenced
by them.
3
Figure 5 - 2 Combined effects of scattering and absorption on the brightness values
Histogram adjustment
This involves evaluating the histograms of the various bands of remotely sensed data of the
desired scene. Normally, the data collected in the visible wavelengths (e.g., TM bands1 to 3) have
a higher minimum value because of the increased atmospheric scattering taking place in these
wavelengths. Conversely, atmospheric absorption subtracts brightness from the data recorded
in the longer wavelength intervals e.g., TM bands 4, 5, and 6 (band 6 in our discussion is 2.08 to
2.35 µm). This effect commonly causes data from the infrared bands to have minimums close to
zero, even when no objects in the scene truly have a reflectance of zero.
If the histograms are shifted to the left so that zero values appear in the data, the effects of
atmospheric scattering will be somewhat minimized. This simple algorithm models the first-order
effects of atmospheric scattering, or "haze." It is based on a subtractive bias established for each
spectral band. The bias is determined by evaluating the original histograms discussed previously.
The atmospheric effects correction algorithm is defined as:
Where
Input BVi, j, k = input pixel value at line i and column j of band k
Output BVi, j, k = corrected pixel value at the same location
Regression adjustment
This atmospheric (haze) correction technique requires that the analyst identify an area in an
image either in shadow or in homogeneous deep, nonturbid water. The pixel brightness values
for this area are then extracted from each of the bands of digital data under investigation. Then,
for each pixel in the area, the brightness values of a visible band (e.g., TM band 1) are plotted
against the corresponding values found at the same location in an infrared band (e.g., TM band
6).
4
On the plotted scattergram in Fig. 5 -2. A straight line is then "fit" through the distribution using
a least-squares technique. If no haze or atmospheric scattering exists, the line would be expected
to pass through the origin of the distribution. However, this is rarely the case. Instead, the line
usually intersects the x -axis. This x-intercept represents the amount of bias caused by
atmospheric scattering. It is the amount that should be subtracted from all the original data in
TM band 1 for this scene. Data from each of the other bands are then plotted against the band 6
data in a similar fashion to derive their offset values. This technique requires many more passes
through the digital data than when using the histogram adjustment method discussed previously.
There is no guarantee that it will provide superior results.
5
corrected only through the use of GCPs are sensor system attitude (roll, pitch, and yaw) and/or
altitude.
Most of the digital Landsat data acquired this days is of already systematic error removed. Unless
otherwise processed, however, the unsystematic error remains in the image, making it
nonplanimetric. This section focuses on the two geometric corrections most often used to make
the digital remote sensor data truly useful; geometric rectification and geometric registration.
Table 5-1 Sources of Image Geometry Errors in Landsat MSS and TM Data
Systematic distortions
Scan skew: Caused by the forward motion of the platform during the time required for each
mirror sweep. The ground swath is not normal to the ground track but is slightly skewed,
producing cross-scan geometric distortion.
Mirror scan velocity: The MSS mirror scanning rate is usually not constant across a given scan,
producing along-scan geometric distortion.
Panoramic distortion: The ground area imaged is proportional to the tangent of the scan angle
rather than to the angle itself. Because data are sampled at regular intervals, this produces
along scan distortion.
Platform velocity: If the speed of the platform changes, the ground track covered by successive
mirror scans changes, producing along-track scale distortion.
Earth rotation: The earth rotates as the MSS scans the terrain. This results in a shift of the
ground swath being scanned, causing along-scan distortion.
Perspective: For some applications it is desirable to have the MSS images represent the
projection of points on the earth upon a plane tangent to the earth with all projection lines
normal to the plane. This introduces along-scan distortion.
Nonsystematic distortions
Altitude: If the MSS platform departs from its normal altitude, this produces changes in scale.
Attitude: One sensor system axis is usually maintained normal to the earth's surface and the
other parallel to the spacecraft's direction of travel. If the sensor surface departs from this
attitude, geometric distortion results.
Rectification is the process by which the geometry of an image area is made planimetric. It may
not, however, remove distortion caused by topographic relief displacement in images. The
process almost always involves relating GCP pixel coordinates (row and column) with their map
coordinate counterparts (e.g., meters in northing and easting in a transverse Mercator map
projection). This is the most precise geometric correction since each pixel can be referenced not
only by its row and column in a matrix after rectification is completed, but it is also rigorously
referenced in degrees, feet, or meters in a standard map projection. Whenever accurate area,
direction, and distance measurements are required, geometric rectification is required. This is
often referred to as an image-to-map rectification.
6
Conversely, there may be times when such geometric precision is not necessary. For example,
one might want simply to compare two images obtained at different dates to see if any change
has taken place in them. Certainly, it is possible to rectify both of the images to a standard map
projection and then evaluate them (and this is often done), but this may not be necessary.
Instead, image registration might be applied, which is the translation and rotation alignment
process by which two images of like geometries and of the same set of objects are positioned
concident with respect to one another so that corresponding elements of the same ground area
appear in the same place on the registered images. This is often called image-to-image
registration.
The same general image processing principles are used in both image rectification and image
registration. The difference is that in image rectification the reference is a map in a standard map
projection, whereas in image registration the reference is another image. It should be obvious
that if an image is used as the reference base (rather than a map), any other image registered to
it will inherit the geometric errors existing in the reference image. Because of this characteristic,
most serious earth science remote sensing research is based on analysis of data that have been
rectified to a map base. Therefore, this discussion will focus on image rectification.
Geometric Rectification
Two basic operations must be performed in order to geometrically rectify a remotely sensed
image to a map coordinate system:
1. The geometric relationship between the input pixel location (row and column) and the
associated map coordinate of this same point (x, y) must be identified. This will establish
the nature of the geometric coordinate transformation that must be applied to rectify or
relocate every pixel in the original input image (x' , y') to its proper position in the rectified
output image (x, y). This process is called spatial interpolation.
2. Pixel brightness values must be determined. Unfortunately, there is no direct one-to-one
relationship between the movements of input pixel values to output pixel locations. It will
be shown that often a pixel in the rectified output image requires a value from the input
pixel grid that does not fall neatly on a row-and-column coordinate. When this occurs
there must be some mechanism for determining the brightness value (BV) to be assigned
to the new rectified pixel. This process is called intensity interpolation.
Spatial interpolation using coordinate transformations
As discussed earlier, some distortions in MSS data may be removed or mitigated using techniques
that model systematic orbital and sensor characteristics. Unfortunately, this does not remove
error produced by changes in attitude (roll, pitch, and yaw) or altitude. Such errors are generally
unsystematic and are best removed by identifying ground control points (GCPs) in the original
imagery and on the reference map and then mathematically modeling the geometric distortion
present. This procedure generally requires that polynomial equations be fit to the control point
data using least-squares criteria to model the corrections directly in the image domain without
7
explicitly identifying the source of the distortion. For moderate distortions in a relatively small
area of an image (e.g., a quarter of a Landsat scene), a six parameter affine transformation is
sufficient to rectify the imagery to a geographic frame of reference.
This type of transformation can model six kinds of distortion in the remote sensor data, including:
translation in x and y, scale changes in x and y, skew, and rotation. When all six operations are
combined into a single expression, it becomes
x' ao a1 x a2 y
y ' bo b1 x b2 y
Where x and y are positions in the output-rectified image or map, and x' and y' represent
corresponding positions in the original input image. By the use of these six coordinate transform
coefficients that model distortions in the original scene, it is then possible to transfer (relocate)
pixel values from the original distorted image x', y' to the grid of the rectified output image, x, y.
However, before applying the rectification to the entire set of data, it is important to determine
how well the six coefficients derived from the least-squares regression of the initial GCPs
accounts for the geometric distortion in the input image. The method used most often involves
computation of the root mean square error (RMS error) for each of the ground control points.
Where xorig and yorig are the original row and column coordinates of the GCP in the image, and x'
and y' are the computed or estimated coordinates in the original image. The square root of the
squared deviations represents a measure of the accuracy of this GCP in the image. By computing
RMS error for all GCPs, it is possible:
(1) To see which GCPs exhibit the greatest error, and
(2) To sum all the RMS error.
Normally, the user specifies a certain amount (a threshold) of acceptable total RMS error. If an
evaluation of the total RMS error reveals that a given set of control points exceeds this threshold,
it is common practice to:
(1) Delete from the analysis the GCP that has the greatest amount of individual error,
(2) Recompute the six coefficients, and
(3) Recompute the RMS error for all points.
This process continues until one of the following occurs: the total RMS error is less than the
threshold specified, or too few points remain to perform a least-squares regression to compute
the coefficients. Once the acceptable RMS error is reached, one can proceed to the intensity
interpolation phase of geometric rectification, which attempts to fill an output grid (x, y) with
brightness values found within the original input grid (x', y').
8
Figure 5 - 4 Residuals and RMS Error Per Point
Total RMS Error
From the residuals, the following calculations are made to determine the total RMS error, the X
RMS error, and the Y RMS error:
1 n 1 n
2 2
RX X Ri RY YRi
n i 1 n i1
1 n
RXi RYi
2 2 2 2
T RX RY
n i1
Where
RX = X RMS e rror
RY = Y RMS e rror
T = total RMS error
n = the nu mber of GCPs
i = GCP number
X Ri = X residual for GCPi
YRi = Y resid ual for GCPi
9
Minimum Number of GCPs
Higher orders of transformation can be used to correct more complicated types of distortion.
However, to use a higher order of transformation, more GCPs are needed. For instance, three
points define a plane. Therefore, to perform a 1st order polynomial transformation, which is
expressed by the equation of a plane, at least three GCPs are needed. Similarly, the equation
used in a 2nd order polynomial transformation is the equation of a paraboloid. Six points are
required to define a paraboloid. Therefore, at least six GCPs are required to perform a 2nd-order
transformation. The minimum number of points required to perform a transformation of order t
equals:
((t 1)(t 2))
Minimum [Link] GCPs
2
Use more than the minimum number of GCPs whenever possible. Although it is possible to get a
perfect fit, it is rare, no matter how many GCPs are used.
Intensity Interpolation
This process involves the extraction of a brightness value from an x', y' location in the original
(distorted) input image and its relocation to the appropriate x, y coordinate location in the
rectified output image. This "pixel filling" logic may be used to produce the output image line by
line, column by column. Most of the time, the x' and y' coordinates to be sampled in the input
image are real numbers (i.e., they are not integers). For example, in Fig. 5-4 we see that pixel 5,
4 (x, y) in the output image is to be "filled" with the value from coordinates 2.4, 2.7 (x', y') in the
input image. When this occurs, there are several methods of brightness value (BV) interpolation
that can be applied, including nearest-neighbor, bilinear interpolation, and cubic convolution.
The practice of brightness value interpolation is commonly referred to as resampling.
Figure 5 - 5 illustration of how the rectified output grid (X, Y) is filled with brightness values
from the geometrically distorted input grid (X', Y').
10
In zero-order, or nearest-neighbor interpolation, the brightness value of the pixel closest to the
x', y' coordinate specified is assigned to the output x, y coordinate. For example, in Fig. 5-4, the
output pixel 5,4 (x, y) requests the brightness value in the original input image at location 2.4, 2.7
(x', y'). There is no value at this location; however, there are nearby values at the integer grid
intersections. A nearest neighbor rule would assign the output pixel (x, y) the value of 15, which
is the value found at the nearest input pixel.
This is a computationally efficient procedure. It is especially liked by earth scientists because it
does not alter the pixel brightness values during resampling. It is often the very subtle changes
in brightness value that make all the difference when discriminating between one type of
vegetation versus another, an edge associated with a geologic lineament, or different levels of
turbidity or temperature in a lake. Other interpolation techniques to be discussed use averages
to compute the output intensity value, often removing valuable spectral information.
First-order, or bilinear interpolation assigns output pixel values by interpolating brightness
values in two orthogonal directions in the input image. It basically fits a plane to the four pixel
values nearest the desired position (x', y') in the input image and then computes a new brightness
value based on the weighted distances to these points. For example, the distances from the
requested x', y' position at 2.4, 2.7 in the input image in Fig. 5-4 to the closest four input pixel
coordinates (2,2; 3,2; 2,3; 3,3) are computed in Table 5-2. The closer a pixel is to the desired x',
y' location, the more weight it will have in the final computation of the average. The weighted
average of the new brightness value (BVwt) is computed according to the equation
4
(Z
2
k / Dk )
BVwt k 1
4
(1 / D
2
k )
k 1
Where Zk are the surrounding four data point values, and Dk are the distances squared from the
point in question (x', y') to these data points. In our example, the weighted average of BV wt is
13.54 (truncated to 13), as shown in Table 5-2. The average without weighting is 12. In many
respects this method acts as a spatial moving filter that subdues extremes in brightness value
throughout the output image. The method is also more computationally demanding than the
nearest-neighbor method.
Cubic convolution resampling assigns values to output pixels in much the same manner as
bilinear interpolation, except that the weighted values of 16 input pixels surrounding the location
of the desired x', y' pixel are used to determine the value of the output pixel.
Table 5-2 Bilinear Interpolation of a Weighted Brightness Value (BVwt) at Location x', y' Based on
the Analysis of Four Sample Points in Fig. 6-7
11
Value at Distance from x', y'
Sample point
sample point, to the sample point
location column row
z D Dk2 Z/ Dk2 1/ Dk2
2,2 9 0.806 0.65 13.85 1.53
3,2 6 0.922 0.85 7.05 1.176
2,3 15 0.5 0.25 60 4
3,3 18 0.67 0.45 40 2.222
∑ 120.00 8.928
BVwt = 120.00/8.928 =
13.54
A threshold of 0.5 was not satisfied until was not satisfied until 13 GCPs were deleted
from the analysis (the order is summarized in Table 5-3).
The seven GCPs finally selected that produced an acceptable RMS error are shown in
Table 5-4. The six coefficients derived from the seven suitable GCPs are found in Table 5-
5.
Table 5-3 Characteristics of Twenty Ground Control Points Used in the computation
Total RMS
Order of
Point Easting on Northing on error after
Points X’ Pixel Y’ Pixel
number map, X1 map, X2 this point
deleted
deleted
1 12 597,120 3,627,050 150 185 0.501
2 9 597,680 3,627,800 166 165 0.663
3 Kept 598,285 3,627,280 191 180 -
4 Kept 595,650 3,627,730 98 179 -
5 2 596,750 3,625,600 123 252 6.569
6 13 597,830 3,624,820 192 294 0.435
7 Kept 596,250 3,624,380 137 293 -
8 Kept 602,200 3,628,530 318 115 -
9 Kept 600,350 3,629,730 248 83 -
10 5 600,680 3,629,340 259 93 1.291
11 Kept 600,440 3,628,860 255 113 -
12 10 599,150 3,626,990 221 186 0.601
13 8 600,300 3,626,030 266 211 0.742
12
14 6 598,840 3,626,460 211 205 1.113
15 3 598,940 3,623,430 214 295 4.773
16 Kept 600,540 3,626,450 272 196 -
17 4 596,985 3,629,350 134 123 1.95
18 7 596,035 3,627,880 109 174 0.881
19 11 600,995 3,630,000 269 71 0.566
20 1 601,700 3,632,580 283 12 8.542
Table 5-4 Information concerning the Final Seven Ground Control Points
Solution
Let
Y = either the x' or y' location in the image, depending on which is being evaluated; in this
example it will represent the x' values
X1 = easting coordinate (x) of the map GCP
X2 = northing coordinate (y) of the map GCP
However, for purposes of analysis it is convenient to express the general equation in matrix
form. Let
13
Y1 1 X 11 . . . X k1 a0 1
. . . . . . . .
Y . a . .
, X . . . . .
, ,
. . . . . . . .
Y 1 X
n 1n . . . X kn a
k
k
In our matrix notation, we can rewrite the general equation as Y Xa . Solving for the
Constants (a0, a1 and a2) and Residuals (ɛ 0, ɛ 1 and ɛ 2) based on the least square principle we
will have the following matrix notation.
Y Xa and
a ( X T X )1 X T Y
Now, let us develop the matrix first for the equation x' ao a1 x a2 y as shown below.
1 1 1 1 1 1 1
X 598285 595650 596250 602200 600350 600440 600540
T
7 4193715 25392960
X T X 4193715 2512499373925 15212991317200
25392960 15212991317200 92114649871600
14
900704.68 0.090562077 0.263251046
1
( X X ) 0.090562077
T
0.000000041 0.000000032
0.263251046 0.000000032 0.000000078
1519
X Y 911177865
T
5510667360
382.063954220
1
a ( X X ) X Y 0.034187679
T T
0.005481044
a0 382.06395422, a1 0.034187679 and a 2 0.005481044
180
179
293
Y 115
83
113
196
Since the remaining quantities are similar, let us only solve for
1519
X T Y 693651510
4203611330
130162716072864
1
a ( X X ) X Y 0.005576180
T T
0.034914996
a0 130162716072864, a1 0.005576180 and a 2 0.034914996
15