WO1999022337A1 - Method and apparatus for enhancing cartographic images and identifying and extracting the features therein - Google Patents

Abstract

An image enhancement and feature extraction system and apparatus use a measured property of an original image (60) to select a set of operations to be performed. In a preferred method, the set of image enhancement and feature extraction operations depends on characteristics of the image histogram (62). Furthermore, parametric values (62) used in image enhancement operations may be selected according to feature type and/or image histogram. In a preferred system, a table (64) that indicates image enhancement operations is indexed according to image histogram and feature of interest (67). The table (64) may also indicate values for parameters used in image enhancement (70) or extraction (71). Additional tables (66, 68) indexed according to feature type provide additional parametric values used in feature extraction (71).

Description

METHOD AND APPARATUS FOR ENHANCING CARTOGRAPHIC IMAGES AND IDENTIFYING AND EXTRACTING THE FEATURES THEREIN

FIELD OF THE INVENTION The invention relates to the field of electronic image processing, and has particular utility in enhancing cartographic images, and identifying and extracting non- hypsographic features within such enhanced images.

BACKGROUND Systems have been proposed for making enhancements to images and for classifying and extracting the features in such images. Typically, the enhancement and extraction processes are highly specialized and suitable only for a particular image source or for a particular type of feature. As a result, commercial products for processing images and for classifying and extracting the features in such images tend to be limited in applicability, and to require significant operator skill and training. Cartographic images (geophysical imagery) exemplify certain problems in terms of image enhancement and feature extraction. There are many sources and types of images used in cartography, including images taken by satellite and aircraft. Images may be panchromatic or multispectral. The wide variation in the physical composition of the imaged geographic area often exceeds the capabilities of automatic and automated systems. A substantial need remains unsatisfied for methods and systems for enhancing cartographic images and for extracting the wide variety of features in such images in a reliable and automated fashion.

An object of the invention is to provide an improved image enhancement and feature extracting method and apparatus.

A further object of the invention is to provide an automated image enhancement and feature extraction method and apparatus. A further object of the invention is to provide a cartographic image enhancement and feature extraction method and apparatus capable of automated operation for a variety of features.

A further object of the invention is to provide a cartographic image enhancement and feature extraction method and apparatus capable of automatically selecting appropriate image processing operations from a set of multiple enhancement, classification and extraction operations.

A further object of the invention is to provide an automated cartographic image enhancement and feature extraction architecture that efficiently integrates multiple enhancement and extraction operations.

These and other objects will become apparent from the figures and written description contained herein.

BRIEF DESCRIPTION OF THE DRAWING FIGURES The present invention will be described below with reference to attached drawings in which:

FIG. 1 illustrates a suitable hardware system in which the present invention may be embodied;

FIG. 2 illustrates a flow of operations for enhancing images and detecting and extracting features in accordance with the present invention;

FIG. 3 illustrates a more detailed flow of operations for calculating image statistics, including histograms;

FIG. 4 illustrates a first table for selecting image adjustment and threshold parameters based on image statistics and feature to be extracted; FIGS. 5a and 5b illustrate a second table for selecting image enhancement and feature extraction parameters based on feature type; FIG. 6 illustrates a third table for selecting image filter parameters based on image resolution and feature type;

FIG. 7 illustrates a flow of operations for extracting a first feature based on parameters selected from the tables of FIGS. 3-5b; FIG. 8 illustrates a flow of operations for extracting a second feature based on parameters selected from the tables of FIGS. 3-5b;

FIG. 9a illustrates filter value placement for a 3 X 3 filter; FIG. 9b illustrates filter value placement for a 5 X 5 filter; and FIG. 10 illustrates a flow of operations for extracting a feature from multispectral data based on parameters selected from tables of FIGS. 4-6.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS FIG. 1 illustrates a suitable hardware system in which the present invention may be embodied. The system comprises: a Central Processing Unit 20; a Read Only Memory (ROM) 22; at least sixty-four megabytes of Random Access Memory (RAM) 24; a Local Hard Disk 26 which serves as virtual memory as well as for storing image processing routines and images; a High Resolution Display 28 displaying the input images, enhanced images, and images showing locations and extent of features of interest (as well as computer operating information); an Operator Interface element 32 to support operator commands via the mouse and keyboard 30; a suitable Graphics Card 34 with video cache to enhance display operations; an Image Data Interface element 36 which uses drivers to control importing and exporting of images to the system; and a Local Bus 40 to connect the resources within the system. A variety of peripherals 38 may be connected to the Image Data Interface element 36, such as floppy disk drive, tape drive, network hard drive, scanner, printer, etc. The Central Processing Unit 20 preferably is one or more Pentium Pro or Pentium 2 processors with a clock speed of at least 166 Mhz. The Windows NT environment is preferred to support multi-tasking, particularly if two or more processors are used. Application programs as discussed further herein may be programmed in the "C" programming language. The ROM 22 preferably contains the overall executive program, selected image processing routines and the key parameter look-up tables. The ROM 22 can be partially replaced by loading the same programs and tables on the Local Hard Disk 26 of the system. The storage space required for the executable code of this program and associated tables is contemplated to be less than 8 Mbytes. The RAM 24 should be error correction coding (ECC) enabled to reduce errors being induced into the lookup tables during processing. The High Resolution Display should be at least a 17 inch monitor with at least a 15.9 inch viewing area and 0.26 mm pixel display. A video card with 4 Mbytes of cache is essential to viewing speed, and a hard disk of a: least 4 Gbytes with LAN access to an imagery storage server is preferred.

FIG. 2 illustrates a flow of operations for enhancing images and detecting and extracting features in accordance with the present invention. In a selection step 60, an operator selects a feature of interest, and loads an image or images (in the case of, e.g., multi-spectral images) from which the feature is to be extracted. The image may, for example, be loaded from a peripheral 38 (FIG. 1) to the local hard disk 26 (FIG. 1). The operator may also be prompted to load imagery of a preferred resolution and spectral type for the feature selected. Additional processing of the image takes place by having the CPU perform operations as defined by program instructions on image pixel values, using RAM 24 as a buffer and the local hard disk 26.

In an image analysis step 62, the machine processes the pixel values of the image(s) to measure certain statistical properties. The properties may include: pixel density histogram, pixel density mean, maximum digital number (DN) in the image, minimum DN in the image, standard deviation, mode, histogram, skewness, the 0.1 percentile DN, and the 99.9 percentile DN. These properties are discussed more fully below with reference to FIG. 3.

In a first table look-up step 64, the machine selects a sequence of image processing steps to be performed, based on the histogram type and the feature to be extracted. The machine retrieves a pixel offset parameter and identifies additional adjustment, if any, to be used in an image enhancement step 70. These parameters are discussed in more detail with reference to FIGS. 5a and 5b. Such additional adjustment may include a histogram stretch, histogram equalization, histogram matching or other general manipulation of the entire image.

Histogram Stretch (also called contrast stretch) is a function whereby the range of digital intensity values is used to calculate a new value for all pixels in the image and thereby increase the contrast in the image. The new digital number in an 8 bit image for any pixel is 255 times the difference of the Old DN minus the original minimum DN in the image, and this product is divided by the intensity range of the image.

New DN = 255 * (Old DN- Min DN Original Image Range

Histogram equalization is a function whereby the image statistics are used to determine amounts of re-distribution of pixels into 'm' different intensity. Consequently to generate an equalized image, the probability of any pixel having a given gray scale level must be adjusted to be 1/m. If the number of pixels with a given original gray level F_io is represented as N_g, and the total number of pixels in the image is defined as N_t then the new pixel value for each pixel is

F_new = Cm-l) * N_s N, Histogram matching is a function whereby the image intensity distribution is adjusted so that it is the same as a specified histogram. This is done in four steps.

First, equalize the original image histogram; second, equalize the desired histogram; third, generate the inverse of the equalized desired histogram; and fourth, apply this inverse transform to the equalized original image. The resulting image is typically the result of a non-linear transformation and the final distribution of intensity values is close to the target distribution.

In a second table look-up step 66, the machine selects additional parameters to be used in enhancing the image and detecting features based on the selected feature type. The machine also retrieves physical size and sensor response parameters as well as work area expansion limits, the physical distance for line segment dilation and erosion, and the minimum acceptable segment length to be used in a feature classification step 71. These parameters are discussed more fully below with reference to FIG. 5b. In a third table look-up step 68, the machine selects the size and content parameters to be used in a filter adjustment step. The machine selects the convolution filter parameters based on the image resolution and the feature to be extracted. In a feature extraction step 71, the machine applies the selected image manipulation processes to the adjusted image, filters to the adjusted image, and generates a modified image that indicates locations and extent of classified features.

Several examples of feature extraction are discussed below with reference to FIGS. 7,

8, and 10.

In an operator acceptance step 74, the machine displays the results of the extraction step 71. The operator in this interactive step 74 views the resulting display and judges whether there is sufficient detail represented in line segments to support easy editing and completion. If the operator confirms that the detail is sufficient, then the modified image (feature image) is ready for evaluation and editing 76. If the judgment is 'no,' then the operator may provide input to the machine by clicking on an input button and invoking a second threshold generation sequence (as discussed more fully below with reference to FIG. 7). The resulting image is then merged with the first. In an output evaluation step 76, the machine displays a modified image (e.g., showing locations and extent of detected features) as an overlay to the original or adjusted image for verification and evaluation by the operator.

In a storage step 78, the machine saves the modified image in permanent or semi-permanent memory, such as the local hard disk 26 (FIG. 1) or peripheral 38 (FIG. 1).

FIG. 3 illustrates a more detailed flow of operations for calculating image statistics, including histograms. After an operator selects the feature and loads the image, the machine executes Step 1, "Get Total Pixel Count" 100, which automatically generates a total pixel count from the image. Step 2, "Initialize Histogram Array" 102, initializes an array for histogram pixel storage and manipulation. The machine then performs Step 3, "Sort Pixels into Grayscale Bins" 104, which sorts pixels by grayscale and counts them in the appropriate bins according to the intensity depth of the image. Typically this depth is 8 bits per pixel for geophysical images and results in 256 bins. In Step 4, "Compute Mean Pixel Value" 106, the machine sums the counts per bin and divides this sum by the total pixel count to determine the mean pixel value of the image intensity. In Step 5, "Determine Min and Max Value" 108, the machine determines the minimum and maximum bins containing pixels which is a measure of a range of contrast. In Step 6, "Compute Std. Dev." 110, the standard deviation (Std. Dev.) of the image about the mean is computed. The machine then computes the skewness of the distribution in Step 7, "Compute Skewness" 112. This is done by calculating the third moment of the distribution around its mean and dividing it by the cube of the Std. Dev. In other words, define 'N' as the total pixel count in the image; μ as the mean of the distribution, μ₃ as the third moment of the distribution about the mean of the histogram of the working image, and W_j as the interim calculation at each of T representing a single digital number (DN) in the pixel distribution. The machine starts at the lowest DN value in the image (not necessarily zero) and calculates ; = [(X_j - μ)³] times the count of the X; DN value. Then, the machine increments DN by one and performs the W; calculation for the next X;. Then the machine sums all the W; and divides by N to obtain μ₃. The machine calculates the Skewness and then can performs a check to determine if it is positive (dark-biased image with pixel counts trailing off to the right), or negative (light-biased image with pixel counts trailing off to the left). Skewness, designated α₃, is the third moment of the working image histogram divided by the cube of the standard deviation of the histogram. In Step 8, "Compute Mode" 114, the machine computes the mode of the distribution. In Step 9, "Compute Percentiles" 116, the machine computes the 0.1 and 99.9 percentiles of the pixel distribution. This step is performed by beginning at the lowest value bin (zero) and summing the cumulative bin pixel counts until the desired percentile is reached. After each of the calculations in steps 106, 108, 110, 112, 114 and 116, the value is stored with the actual image data in the final step 118.

FIG. 4 illustrates a first table for selecting image adjustment and threshold parameters based on image statistics and feature to be extracted. The table is multidimensional, in the sense that it provides multiple parametric values based on multiple inputs. The table is indexed according to feature type 202 and histogram type 204. Each of the histogram types is defined according to a number of parameters, including mean values 206 and skew 212. Histogram mean value 206 is calculated for the working image based on the pixel count occurring in each digital number bin (e.g., from 0 to 255 for 8 bit images). The histogram type is selected by comparing the measured mean to ranges of histogram means in the table. For instance, a working image with a mean of 120 would cause the machine to enter the table at the 64< M < 128 section 205, also called the histogram category B portion of the table. Skew 212 is the measure of the lack of symmetry of a statistical distribution. It is a secondary check for the working image statistics and (for a single mode image) the table is set up to handle positive and negative skew of the image. This information is used by the machine to determine the detailed adjustment of the image required to support the extraction of the selected feature. It is used in conjunction with the working image Mean 206 to determine the direction and size of the offset for image adjustment in subsequent image enhancement steps. The values in this field are "P" for positive skew (trailing to the right), or "N" for negative skew (trailing to the left).

For a measured histogram type, the table returns values for a number of parameters, including First (1st) Pan A/B 214, Second (2nd) Pan A/B 216, Offset 218, Adjustment 220, and Sequence 222. All PAN A/B values are predetermined independently based on typical grayscale values of the feature of choice in an 8 bit image. They are not measured from the working image in process. This PAN A/B value is used in FIG. 7, step 504, FIG. 8, step 554, and FIG. 10, step 658 in preparation for edge detection steps. For example, First PAN A/B may be used as a threshold value when First PAN A/B is larger than a value of the mean plus 1/2 standard deviation. Pixels with density values less than the threshold are replaced with a zero density value.

Second Pan A/B is used in FIG. 7, step 524 for a second threshold operation which yields an image that can be combined with the result of the first threshold operation (logical AND or logical OR) to improve continuity in primitive segments. Offset 218 is a density number which may be added to (to lighten) or subtracted from (to darken) all pixels in the image. Offset is used in steps 502 (FIG. 7), 552 (FIG. 8), and 652 (FIG. 10) to bring the image pixel values into a more efficient range for automated processing and operator verification in subsequent processing steps . Adjustment 220 indicates the type of image enhancement that may be performed in addition to simple offset. Such additional adjustment may be a histogram stretch, histogram equalization, histogram matching or other general manipulation of the entire image. If histogram matching is to be performed, the identity of the target histogram is included in the column. Sequence 222 contains identifying numbers of image processing sequences for detecting, classifying and/or extracting the related feature from an image. They may be called from within step 71 (FIG. 2) and executed as software subroutines from a subroutine library. The identifying numbers beginning with "3" are image processing for panchromatic imagery. The numbers beginning with "4" are image processing for multispectral and hyperspectral imagery. Additional number series (e.g., beginning with "5") can be used for additional image sources, such as radar and lidar imagery of the planet surface.

For example, for a histogram having a mean less than 64, and for a river/coastline feature, the table returns the following specific sequence 224: "311; 411." The first value "311 " designates a processing routine for rivers and coastlines using panchromatic imagery. The sequence 311 is shown in more detail in FIG. 7. The second alternative of sequence 224, "411 " designates a processing routine for rivers and coastlines using multispectral or hyperspectral imagery.

For this example, the table also returns an offset 226 of one ("1"). This offset is used to brighten the working image to assist with visual verification of the extracted features by the operator. This offset is multiplied by a scale factor, e.g., 18, to get an offset in terms of pixel density values. Generally, the offset is greater for images with a mean less than 64 or greater than 128. If no offset is required, a "NULL" is entered into the table.

As a second example, for a histogram having a mean of between 64 and 128, and for a lake & pond feature, the table returns a First Pan A/B value 228 of "56." This value is used in a subsequent step as an input to a threshold calculating function.

The calculating function is based on the 0.1 percentile to 99.9 percentile region of grayscale values in the image and is used to calculate a threshold in FIG. 7, steps 524 and 504.

The table also returns a Second Pan A/B value 230 of "60" which is used in an optional second pass through the image to pick up segment fragments that were not included in the first pass. After the first pass, the operator is asked if the quality of the extraction warrants an additional pass (see, e.g., FIG. 7, step 524). If the operator responds "yes," the machine processes another copy of the working image with the

Second Pan A/B value, and then merges the two resulting images. The value is generally higher than, but near, the first pan value to reduce image noise while increasing the number of primitive segments (e.g., discontinuous edge segments) to be combined during a contour tracing operation.

The "3131" sequence entry designates a process routine for identifying lakes and ponds using panchromatic images. Other than feature size parameters used in classification, the processing for lakes and ponds is similar to that performed on rivers and coastlines. The "4131" entry is for multispectral image processing routine to locate, classify and extract lakes and ponds as shown in FIG. 10.

FIG. 5 illustrates a second table for selecting line segment manipulation parameters and classification measures based on feature type. The machine enters the table of FIG. 5 with a selected feature 300. The table returns a number of parameter values. This table incorporates information about the physical characteristics of the feature of interest, such as shape and size of the features, and is independent of the ground sample distance (resolution) of the imagery processed. In subsequent processing steps, the size-dependent characteristics are scaled according to a resolution of a specific image to convert these values to pixel units. Resolution can be obtained from a header density file or other data associated with the image file. The features column 300 in this table contains entries that are generally the same entries as those in FIG. 4. One exception is "Lakes and Ponds" which is divided into two separate entries to distinguish size-related and texture-related characteristics. Ponds are smaller than lakes and may have different size-related characteristics. Major Axis 304 is the range of the largest physical dimension in meters related to the feature under consideration. In a first example, this is a "NULL" for rivers and coastlines because they typically extend beyond the region covered a single image and therefore cannot be used in processing. This may also be true for other linear features such as primary roads and railroads. In a second example, the major axis for lakes 344 is a range of 1001 to 100,000 meters. The major axis range for ponds 346 is only 10 to 100 meters.

Minor Axis 306 is the smaller cross-object dimension and is used in subsequent processing of both linear and polygonal geographic features. The Minor Axis entry 348 for rivers and coastlines is 10 to 1600 meters and is only used if both shorelines are present in an image. In a second example, the Minor Axis for lakes 350 is from 101 to 1000 meters.

Height 308 is applicable to most non-linear features such as buildings and other structures. The Height entry for the linear feature roads 352 is from 0 to 3 meters and supports analysis of the shadows cast by on-ramps/off-ramps. However, a "NULL" is used in processing rivers and coastlines, lakes and ponds.

Perimeter (Perim) 308 is applicable to polygonal features in an image and is measured in meters. During extraction of a polygonal feature, this entry is compared to a measured edge pixel count from an image. (Edge pixel count is first multiplied by the resolution of the image to match scale.) Area 310 is applicable to polygonal features and is used to discriminate among similar objects. In a first example, Area for lakes 354 is within the range 101,000 to 100 million square meters, as contrasted with the Area of ponds 356 which is only 105 to 2000 square meters.

Mean Radius (Mn Rad) 312 is the average distance from the centroid of the feature to its Perimeter and is measured in meters. This measure is used to determine the circularity of features.

Minimum Enclosing Rectangle (MER) 314 is a measure which can be determined in different ways according to the method for defining an orientation. The method used in the preferred embodiment is to make the measurement from a frame of reference oriented parallel to the bottom of the image, as opposed to a frame of reference orientated parallel to the major axis of the feature. In lakes, for example, the values for an acceptable MER range 358 is the same as the range of Area 354. In a second example, the Area 356 for ponds and the MER 360 are the same.

Percent of the MER filled by the feature, (%MER) 316 is a measure derived from the Area of the feature and the Area of the MER. This measure is not contemplated for use with lakes and ponds, but relates more appropriately to manmade feature such as buildings. Compactness (Comp) 318 is the ratio of the square of the Perimeter divided by the Area. The irregularity of lakes and ponds precludes a meaningful use of compactness for such features. However, compactness is useable for manmade features such as buildings and other structures. Compactness is used to determine the linearity and "sensed density" of an object. Expand Segment MER width (Xpnd Seg MER W) 320 is an image-resolution- independent entry. In segment joining process, the machine searches for disjoint segments in proximity to a "base" or starting segment. Xpnd Seg Mer is used to bound the expansion of the search area in a first direction for segments which can be joined to a base segment. (Breaks occur during image thresholding and initial primitive segment generation.) In an example for primary roads, the value 5 meters 362 is automatically used without operator intervention to expand the width of the MER of a base road segment to find candidate segments that can be combined. This entry is used in conjunction with Expand segment MER length (Xpnd Seg MER L) 342. It is used to bound the expansion of the search area in a second direction for segments which can be joined to a base segment. In an example for primary roads, the value 15 meters 370 is automatically used without operator intervention to expand the length of the MER of a base road segment to find candidate segments that can be combined.

Trace Length 322 is the minimum acceptable length of segments used during contour tracing to perform initial segment combinations. In one example, rivers and coastline 364, trace length is 15 meters and is the same for large water bodies such as lakes. In a second example, Ponds 366, Trace Length is only 10 meters and is the same as the initial trace length used for primary roads 368.

Dilation 324 is a physical measure expressed in meters, rather than as pixels. This value, designated "A", is divided by the image resolution in meters per pixel to determine the number of pixels to dilate a segment in order to combine segment - fragments.

Erode 326 is also a physical measure expressed in meters, rather than as pixels. This value, designated "B," is divided by the image resolution in meters per pixel to determine the number of pixels to erode a segment combination in order to retain the combination while thinning it for further processing as defined by the routine peculiar to the feature being extracted.

Red to Blue ratio "R/B" 328 is a range of the ratio of reflectance in the visible red to visible blue bands. Coarse determination of material types can be made with this ratio in conjunction with other measures (e.g., compactness or texture) during feature classification processing. Red to Green ratio "R/G" 330 is a range of the ratio of reflectance in the visible red to visible green band. Coarse determination of material types can be made with this ratio in conjunction with other measures during classification processing.

Mean First Minimum (Mean First Min) 332 is the location in nanometers of the first minimum reflectance value, which distinguishes vegetation from minerals. There is a distinctive trough related to chlorophyll. Other substances cause the minimum to occur in other regions of the spectrum. One example is class 1 & 2 roads 372 in which this spectral energy minimum is associated with concrete and asphalt falls into the range of 370 - 390 nanometers.

Mean First Following Maximum 334 is the spectral reflectance maximum which follows the first minimum in the histogram. One example is class 1 & 2 roads 374 in which this maximum is associated with concrete and asphalt and falls into the range of 700 to 750 nm.

Texture 336 is a measure of change in reflectance over a given area or length for an established pattern. To retain image resolution independence, texture is defined as the change in reflectance intensity divided by change in area (for polygonal objects) or divided by change in length (for linear objects). One example of a polygonal object is lakes 376 where the texture is a range of 0 to 100 per change in area. Area is in hundreds of square meters. A second example of a linear object is roads class 1 & 2, 378 where the texture range is 0 to 80 per change in length. Length is in tens of meters.

Period 338 is not applicable to the example features in the portion of the table reproduced in FIG. 5, but applies to features such as electrical power transmission lines and railroad tracks that have regular recurring patterns. Features of this nature have repeating characteristics, such as pole/tower placement and tie placement, which have a defined period that can be evaluated using a Fourrier transform. Period is the range of a count per 100 meters.

PAN Tolerance 340 is a measure applied to surface tracking and surface similarity processing routines. Such routines use a value from the working image and then find all pixels within a tolerance of that value (average of values) for primitive detection prior to further processing. One example is Class 1 & 2 Roads where differences in patch reflectance and other material changes may cause processing anomalies. A PAN Tolerance value 382 such as "40" reduces such anomalies in processing and can be extended to individual bands of multispectral and hyperspectral data. A second example, 380 provides a tolerance of "60" for changes in water quality and apparent depth as indicated by a single band of data for ponds. The tolerance for lakes 384 and rivers 386, which cover greater areas, is larger as seen in the table as a value of "80."

FIG. 6 illustrates a third table for selecting image filter parameters based on image resolution and feature type. This table contains the size and contents of the principal edge detection filter used for each feature at the specified image ground sample distance (GSD). Resolution 402 has six classes of GSD 0-0.5; 0.51-1.0; 1.1-5; 5.1-10; 10.1-20; and 20.1-30. This set can be extended to finer resolution in the half meter region and to coarser resolution above 30 meters. The class 1.1 - 5 meters 414 is shown as an example appropriate for satellite data. Data in the finer classes is available from commercial aerial imagery vendors (and is contemplated to be available from commercial satellite sources in the future). For other resolutions, filter values for those features would be scaled up or down from the ones shown for 1.1-5 m resolution. The machine enters the table based on the resolution of working image and the feature to be extracted. Feature 404 has the same entries used in other tables. (Lakes and Ponds both use the same entry.) Size 406 defines the filter as a 3 X 3 or 5 X 5 matrix to be used in the convolution for edge detection. FIG. 9a illustrates filter value placement for a 3x3 filter. FIG. 9b illustrates filter value placement for a 5x5 filter.

The filters in Fig 6 are varied specifically for each feature being extracted. The filter values are in the field designated "FV11" 408 to "FV55" 412 to identify the placement of the filter values (FV) in the filter mask. "FV 11 " 408 is the upper left comer of the filter mask. "FV33" 410 is the lower right comer of a 3 X 3 mask, or the center value of a 5 X 5 mask.

FIG. 7 illustrates a flow of operations for extracting rivers and coastline based on parameters selected from the tables of FIGS. 3-5b. Step 1, "Calculate Image Statistics" 500, calculates the statistics of the pixel value distribution in the image. This is the same as step 2, "Calculate Image

Statistics/Histogram" 62 illustrated in FIG. 2.

Step 2, "Adjust Image" 502, adjusts the image to optimize the segment generation and noise reduction. This is defined in look-up tables for the input parameters to the Rivers/Coastlines extraction function as well as to determine the routines to run to perform the extraction. Adjustments include an offset of "1" (i.e., an increase of magnitude of 18 in pixel density) and a histogram (contrast) stretch. Step 3, "Make Copies of Image" 506, makes three copies of the adjusted image. Step 4, "Generate SES on Copy 1" 508, applies a 3 X 3 Sobel horizontal edge detector to a first copy of the adjusted image to generate a sidereal edge segment

(SES) working image. The values for this filter are: "1, 2, 1" (top row); "0, 0, 0"

(center row); "-1, -2, -1" (bottom row). The SES results will be applied after a second pass at the shoreline definition is completed with a second set of threshold values and the two shoreline images are combined. Step 5, "Extract Contours on Copy 2" 512, applies a water contour extracting algorithm to a second copy of the adjusted images. In a threshold step, "Calculate 1st Threshold" 504, the water contour extracting routine first assigns "0" to all pixels above the Shoreline Threshold (SThld) where:

SThld =(fPan A/B. *.99.9 percentile DN- 0.1 percentile DN⁾ 1 + 6 255

(Pan A/B) is a table entry 214 in FIG. 3 for the feature being extracted (which may be First Pan A/B or Second Pan A/B, depending on when the subroutine is called). This produces a thresholded image with limited grayscale values. Then the routine applies the edge detecting filter for rivers and coastlines. The size of and values for the edge detecting filter are taken from FIG. 6. The layout of the filter mask for a 3 X 3 is shown in FIG. 9a and for a 5 X 5 filter in FIG. 9b. This process produces a binary edge image.

Step 6, "Trace Contour" 514, tracks the line segments based on the resolution of the working image. The contour tracer routine 514 looks for contiguous pixels with a DN within 18 of the threshold value. The routine assigns a density value of "255" to any pixel passing the test and "0" to any failing it. This step yields a binary file with more connected shoreline segments and noise.

Step 7, "Dilate Segments" 516, applies a dilation routine to all items in the binary shoreline file. The size of this dilation is dependent on the amount of ground area represented by each pixel and the fact that the selected feature is rivers/coastlines. The amount of this dilation is the value "A" (324 in FIG. 5b) divided by the image resolution, which is obtained from the image header file. For a resolution of 4-8 m/pixel, the result is to add from 2 to 4 pixels to all sides of the binary line segments. The machine then releases the buffer.

Step 8, "Erode Segments" 518, applies an erosion routine to the result of the dilation. The amount of erosion is the value B 326 from FIG. 5b divided by the image resolution. The erosion will remove between one and three pixels to all sides of the dilated binary image for a resolution of 4-8 m/pixel. Any line segments that were joined during dilation will remain intact, but will be thinner.

The operator is asked in Step 9, "Operator Input on Adequate Segmentation" 520, to examine the modified image, and a decision is received in an input step, "Is More Extraction Needed" 522, whether additional segment extraction is required. If the response is "yes," the routine calculates a new threshold value 524, which is the higher of: (1) the original Pan A/B 214 plus one half the absolute value of the image mean minus the image mode, and (2) the 2nd PAN A B value 216 from FIG. 4. Steps 4 through 8 (512, 514, 516, 518) are repeated by applying the new threshold value to the third copy of the image 510. This second pass is used to handle major water reflectance differences within an image such as at a confluence of a slow (dark) and a fast (light) river. The two resulting shoreline binary images are combined in Step 11 , "Merge Modified Copies 2 & 3" 526. Individual shoreline segments of the second image are used as candidates for "segment joins" with feature segments in the first feature image.

Step 12, "Generate PROI and Calculate Orientation of Segments" 528, calculates (in pixel coordinates) a zone around the binary feature (line) segments. This is a processing region of interest (PROI), and its width is defined as twice (A / image resolution), where "A" is the dilate parameter 324 of FIG. 5b. This step also calculates the orientation (in degrees from vertical in the image) of the line segments within the region.

Step 13, "Despeckle" 530, applies a 3 X 3.despeckle routine to reduce the noise away from the shoreline (outside the PROI). This despeckle routine is applied in a manner similar to soften and sharpen filters, but is applied only outside the PROI. The 3 X 3 filter values are: FV11= 1 ; FV12= 0; FV13= 1; FV21= 0; FV22= 4; FV23= 0; FV31= 1; FV32=0; FV33= 1; scale 1; offset: 0. Step 14, "Apply SES Reduction" 534, applies a Sidereal Edge Structure (SES) reduction routine to identify any horizontal edge-detected segments within the PROI that are perpendicular to a tangent of the shoreline within a small tolerance (absolute value of direction of sidereal minus direction of feature segment less than or equal to 15 degrees). The segments normal to the tangent of the shoreline are probably structures protruding from shore to water such as piers, portions of bridges, etc. and should be removed. The machine then compares with the shoreline pixel location and gives any of these suspect wharf structures a gray pixel value of 60 in the shoreline image. In an input step, "Operator Confirmation" 532, the machine then displays the gray segments, and gives the operator an option to delete them (i.e., give them a density value of zero).

Step 15, "Save Feature File" 536, prompts the operator to save the file, and then release all resources.

FIG. 8 illustrates a flow of operations for extracting a second feature, lakes and ponds, based on parameters selected from the tables of FIGS. 3-5b. The primary processing routine will call modules to perform sixteen general steps.

Step 1, "Calculate Image Statistics" 550, calculates the statistics of the pixel value distribution in the image, and is the same as the step, "Calculate Image Statistical Histogram" 62 from FIG. 2. Step 2, "Adjust Image," adjusts the image (use lighten, darken or histogram stretch) to optimize the segment generation and noise reduction as defined in look-up table FIG. 4 for the input parameters to the Lakes Ponds extraction function, as well as to determine the routines to run to perform the extraction. For lakes and ponds in an image having a "Bl" histogram (mean between 64 and 128 inclusive, positive skew), there will be no offset and no adjustment. The processing routines 3131 are described more fully below. Step 3, "Make Image Copies" 556, makes two copies of the adjusted image. A Sobel horizontal edge detector 558 is applied to a first copy of the adjusted image to generate a sidereal edge segment working image.

Step 4, "Extract Contours on Copy 2" 560, applies the water contour extracting algorithm to a copy of the adjusted image. In a threshold step, "Calculate Threshold" 554, the water contour extracting routine first assigns "0" to all pixels above the Shoreline Threshold (SThld) where:

SThld = ((Tan A/B) *.99.9 percentile DN- 0.1 percentile DN)} + 6 255

(Pan A B) is a table entry from FIG. 4 for lakes and ponds. Then, the machine applies the edge detecting algorithm. The size of and values for the edge detecting filters come from the table of FIG. 6. This process produces a binary edge image. Step 5, "Trace Contour" 562 tracks the line segments based on the resolution of the working image. The contour tracer routine looks for contiguous pixels with a DN within "18" DN of the threshold value and assigns "255" to any pixel passing the test and "0" to any failing it. This step yields a binary file with more connected shoreline segments and noise. Step 6, "Dilate Segments" 566, applies a dilation routine to all items in the binary shoreline file. The size of this dilation is dependent on the ground area represented by each pixel and the fact that the selected feature is lakes and ponds. The amount of this dilation is the value "A" in 324 in Table 5b divided by the image resolution. For a resolution of 4-8 m/pixel, the result is to add from 2 to 3 pixels to all sides of the binary line.

Step 7, "Erode Segments" 564, applies an erosion routine to the result of the dilation. The amount of erosion is the value "B" 326 from Table 5b divided by the image resolution. The erosion will remove between one and two pixels to all sides of the dilated binary image for a resolution of 4-8 m/pixel. Any line segments that were joined during dilation will remain intact but will be thinner.

Step 8, "Determine Perimeter Closure" 568, determines closure on the perimeter (edge) of the pond or lake. Step 9, "Operator Editing of Perimeters" displays perimeters for the operator. If closure does not exist, the machine draws a boundary around the perimeter break(s) in magenta to assist in operator editing of the break. The operator can accept or edit the proposed boundary. The operator is also asked if the interest is in lakes or in ponds. This allows the machine to get dimensional ranges from the table in FIG. 5. This copy of the image is then flickered as an overlay on the original working image, and a raster editing function is enabled to support editing in the highlighted work areas.

Step 10, "Generate PROI and Calculate Orientation" 570 of Segments/Tangents 572, calculates (in pixel coordinates) a zone around the binary line segments. This is a processing region of interest (PROI), and its width is defined as 2 (A / image resolution). The machine calculates the orientation (in radians from vertical in the image) of the line segments within the region. Optionally, the machine may then apply a 3 X 3 despeckle routine to reduce the noise away from the shoreline (outside the PROI). This despeckle routine is not shown in FIG. 8, but would be applied in a manner similar to the despeckle routine described with reference to FIG. 7. The filter values would be: FV11= -1; FV12= -1; FV13= -1; FV21= -1; FV22= 7; FV23= -1; FV31= -1 ; FV32= -1 ; FV33= -1. See FIG. 9a for the filter layout.

Step 11, "Generate SES on Copy 1" 558, applies a Sidereal Edge Structure (SES) reduction routine to identify any horizontal edge-detected segments within the PROI that are perpendicular to the shoreline within a small tolerance (absolute value of direction of sidereal minus direction of feature segment less than or equal to 15 degrees). The machine then compares the SES location with the shoreline pixel location and gives any of these suspect wharf structures a pixel value of 60 in the shoreline binary — it is no longer a binary image.

Step 12, "Operator Confirmation and Editing" 576, displays a message, "Gray segments are suspected structures, not shoreline — press DELETE to eliminate all gray line segments." The machine will assign all pixels with a value of "60" a value of "0," if the operator presses Delete in a manner similar to step 534 of FIG. 7.

Step 13, "Fill SES Gaps" 578, generates a yellow line to fill in the gaps created when the SES were removed. The operator is prompted to confirm the fill or draw one using a raster editing feature. Step 14, "Index Probable Water Bodies and Calculate Dimensions" 580, generates an index of all probable water bodies (closed perimeter polygons resulting from the operator editing and SES clean up) in the image and performs dimensional calculations on each. The machine selects as the desired feature, polygons that satisfy the parameter ranges defined in the tables of FIGS. 5a and 5b. Step 15, "Classify Polygons" 582, classifies the water bodies as lakes or ponds, or other. This step also outlines the perimeter(s) of ponds in yellow and lakes in magenta.

Step 16, "Save Feature File" 584, prompts the operator to save the file. If the operator indicated that ponds were of interest, the machine assigns "255" to all yellow pixels and "0" to all others. The machine then saves the image and then releases all resources. If the operator indicated that lakes were of interest, the machine assigns

"255" to all magenta pixels and "0" to all others. The machine then saves the image and then releases all resources.

A further example of the interaction of the image processing sequences and the parameter tables is the extraction of lakes and ponds from multispectral data. The primary difference between panchromatic and multispectral processing is related to the increased information available when working with ratios of reflectance in specific spectral bands. The specific feature to be illustrated is the extraction of a pond as illustrated in FIG. 10. The preferred image data (resolution and wave length) for this, in descending order of preference, are: 0-1 meter in NIR (700-1100 nm) band; 0-1 meter in SWIR (1100-3000 nm) band; 1.1-5 meter in NIR (700- 1100 nm) band; 1.1-5 meter in SWIR (1100-3000 nm) band; 5.1-10 meter in NIR (700-1100 nm) band; 5.1- 10 meter in SWIR (1100-3000 nm) band; 10.1-20 meter in NIR (700-1100 nm) band; 10.1-20 meter in SWIR (1100-3000 nm) band; 20.1-30 meter in SWIR (1100-3000 nm) band . Secondary imagery in the Blue (400-500 nm), Green (500-600 nm), and Red (600-700 nm) bands in the matching resolutions can be used to confirm classification and support final display of findings. The operations are as follows. Step 1 , "Calculate Image Statistics" 650, calculates the pixel value distribution in the NIR or SWIR image. This process is the same as "Calculate Image Statistics/Histogram" 62 (FIG. 2).

Step 2, "Adjust Image" 652, uses the offset value of FIG. 4 ("1") to lighten or darken the image and then performs a contrast stretch as also defined in FIG. 4.

Step 3, "Enter Tables 2a and 2b" 654, obtains the input parameters for the extraction function. Tables 2a and 2b are illustrated in FIGS. 5a and 5b, respectively. The Major Axis parameter 304 of FIG. 5 ( 10 to 100 meters) and the perimeter parameter 309 are used to distinguishing ponds from lakes. The ratio of Red to Blue 328 and Red to Green 330 reflectance entries of FIG. 56 are confirmatory measures, which supplement the fact that NIR and SWIR energy is absorbed by water. Also used are the trace length 364 (which is used to define the shortest line fragment allowed to be combined to generate a shoreline), .and the dilation 324 and erosion 326 parameters. The variability of pond water reflectance across the visible bands is less than such variation in lakes and may be used as another distinguishing physical phenomenon. Step 4, "Enter Table 1" 656, obtains the First PAN A B value (FIG. 4 item 214) in preparation for calculating the threshold for the NIR or SWIR image.

Step 5, "Threshold the NIR/SWIR Image" 658, calculates the threshold, and then thresholds the NIR or SWIR image by assigning the value "0" to all pixels at least 20 DN below the First PAN A/B value for Lakes and Ponds in the histogram entry in the table of FIG. 4. For instance, if First PAN A/B = 60, then the machine assigns "0" to all pixels below "40." The machine also assigns DN =255 to all values equal to or above this adjusted First PAN A/B value. Editing may be required due to occlusions such as overhanging trees, docks, bridges. Step 6, "Extract Contours" 660, applies the edge detection filter with the size and filter values in Table 3 (FIG. 6) for lakes and ponds and the associated image resolution. This processing produces an image with gray outlines for the lake and pond boundaries. The machine assigns "255" to all non-zero DN valued pixels. The resulting NIR/SWIR image now has white outlines of probable water bodies. Step 7, "Trace Contours" 662, tracks the line segments based on the resolution of the working image. The contour tracer routine looks for contiguous pixels with a DN within 18 DN of the threshold value. The machine assigns "255" to any pixel passing the test and "0" to any failing it. This step yields a binary file with more connected shoreline segments and noise. Step 8, "Dilate Segments" 668, applies a dilation routine to all items in the binary shoreline file. The size of this dilation is derived from the table of FIG. 5b, and is dependent on the ground area represented by each pixel and the selected feature of rivers/coastlines. The amount of this dilation is the value "A" 324 in FIG. 5b divided by the image resolution, which is obtained from the image header file. For a resolution of 4-8 m pixel, the result is to add from 1 to 3 pixels to all sides of the binary line segments. Step 9, "Erode Segments" 666, applies an erosion routine to the result of the dilation. The amount of erosion is the value "B" 326 from in FIG. 5b divided by the image resolution. The erosion will remove between one and two pixels to all sides of the dilated binary image for a resolution of 4-8 m/pixel. Any line segments that were joined during dilation will remain intact but will be thinner.

Step 10, "Load the Visible Image Bands (R, G, B)" 664, prompts the operator to load the associated visible band data. Typically these files are collected simultaneously with the NIR or SWIR imagery and are automatically co-registered. Step 11, "Operator Editing of Perimeters" 674, displays an editing control window for the operator and supports raster editing for closure of the perimeters of probable water bodies. This simple editing function allows only pixel inverting , i.e., white to be drawn on black, or black on black for erasure.

Step 12, "Copy Perimeter Image and Fill" 672, makes a copy of the perimeter image and then assigns an intensity value of 255 to all pixels outside the perimeter of the polygons. This image becomes a mask which, when applied to another data file, takes all values to zero for pixels outside the polygons.

Step 13, " Mask R, G, B" 670, subtracts the perimeter fill image from each of the visible band data files (red, green, and blue) and acts as a land mask to identify the suspect pond areas on the blue red and green band data sets. Step 14, "Calculate R/B; R/G" 676, determines the average ratio of red to blue and red to green reflectance of each of the water bodies viewable though the mask (non-zero).

Step 15, "Index Probable Water Bodies & Calculate Measures" 678, generates an index of probable water bodies and calculates physical dimensions of the polygons and other derived measures, such as texture.

Step 16, "Classify Polygons" 680, determines if remaining polygons fall into a specific category of polygon using values in FIGS. 5a and 5b, e.g., Area 310, Perimeter 309, R B 328, R/G 330, and Texture 336. Water bodies which have characteristics that are within the ranges displayed in the table for ponds are classified as ponds. Values such as Texture 336 are a subset of acceptable values for lakes and therefore might not be used to discriminate ponds from lakes. All parameters may be used to increase confidence in the unsupervised clustering process used for classification.

Step 17, "Save Feature File" 682, prompts the operator to save the pond boundary image as a layer to be placed under vector shorelines for revision.

Although the present invention has been described with reference to the particular preferred embodiments, various modifications and variations can be made that will be apparent to those in the art and which will fall within the scope of the present invention as defined by the appended claims.

Claims

WHAT IS CLAIMED IS:

1. A method for controlling the generation of a modified image from an original image, said modified image indicating regions representing a selected feature, said method comprising the steps of: establishing a predetermined set of a plurality of image enhancement operations; establishing a predetermined set of features; establishing a predetermined set of a plurality parametric valu s, said values varying according to an image property and according to a feature; selecting a feature of interest; measuring the image property of the original image; selecting a. subset of image enhancement operations from the predetermined set of image enhancement operations according to the measured value of the image property; selecting a subset of parametric values according to the selected image; performing the selected set of image enhancement operations to generate an enhanced image; extracting regions representing the selected feature from the enhanced image using the selected subset of parametric values; and generating a modified image indicating regions representing the selected feature.