mSSBBHSSSSHm 
 
 nBRgmraHH 
 
 H 
 
 DUMB 
 
 ""ar 
 
 BnSBnsKmt 
 
 BEnHMQ 
 
 ■ 
 
 AVV" 
 
 
 ■HAM KAaiiVw^ig' 
 
 BBS fto8 
 
 Hon 
 
LIBRARY OF THE 
 
 UNIVERSITY OF ILLINOIS 
 
 AT URBANA-CHAMPAIGN 
 
 5IO. 84 
 y>o. 800 - 8o*3 
 
 Cop. 2^ 
 
Digitized by the Internet Archive 
 in 2013 
 
 http://archive.org/details/proceduresforana801haid 
 
fU r 
 
 nc ft/ UIUCDCS-R-76-801 
 
 l v [<iik- 
 
 7 
 
 T 
 
 
 PROCEDURES FOR ANALYZING DATA OBTAINED FROM 
 "LANDSAT" MULTISPECTRAL IMAGERY 
 
 by 
 
 MAHMOUD HAIDAR 
 
 May 1976 
 
UIUCDCS-R-76-801 
 
 PROCEDURES FOR ANALYZING DATA OBTAINED FROM 
 "LANDSAT" MULTISPECTRAL IMAGERY 
 
 BY 
 
 MAHMOUD HAIDAR 
 
 May 1976 
 
 Department of Computer Science 
 University of Illinois 
 Urbana, Illinois 61081 
 
 This work was supported in part by Contract No. USNASANGR14-005-202 
 at the Center for Advanced Computation and by the Department of 
 Computer Science and was submitted in partial fulfillment of the 
 requirements for the degree of Master of Science in Computer Science 
 at the University of Illinois. 
 
Ill 
 
 ACKNOWLEDGMENT 
 
 I wish to thank my thesis advisor, Professor D. Slotnick, 
 and Professor R. Ray for their guidance and encouragement in the 
 preparation of this thesis. 
 
IV 
 
 TABLE OF CONTENTS 
 
 Chapter Page 
 
 1. INTRODUCTION 1 
 
 1.1 Remote Sensing of Earth's Resources 1 
 
 2. THE LANDSAT PROJECT 4 
 
 2.1 LANDSAT Mission 5 
 
 3. COMPILATION AND GEOMETRIC REGISTRATION OF LANDSAT 
 MULTITEMPORAL IMAGERY 14 
 
 3.1 A Basic LANDSAT Data Preprocessing Problem: 
 
 Image Rectification 15 
 
 3.2 Oblique Transformation of LANDSAT Images to 
 Approximate North-South Orientation 21 
 
 4. PATTERN RECOGNITION: A BASIS FOR REMOTE SENSING 
 
 DATA ANALYSIS 48 
 
 4.1 Classification Procedure 48 
 
 4.2 Clustering Procedure 60 
 
 LIST OF REFERENCES 66 
 
LIST OF FIGURES 
 
 Figure Page 
 
 1 Overall LANDSAT System ... 6 
 
 2 RBV Camera Head Orientation 9 
 
 3 MSS Ground Scan Pattern 11 
 
 4 Buffer Requirement for Core-Oriented Rectification 
 
 of LANDSAT Images 19 
 
 5 The Oblique Transformation Process 22 
 
 6 A Simplified Model of a Pattern Recognition System . 49 
 
 7 Decision Regions and Surfaces 50 
 
 8 A Pattern Classifier Defined in Terms of 
 
 Discriminant Function 50 
 
 9 Clustering Algorithm 61 
 
 10 Separability of Clusters 64 
 
1. INTRODUCTION 
 
 Human beings are impelled to develop means for automatic image 
 processing for at least two major reasons: 
 
 1. As man expands into portions of the universe formerly 
 denied him, the need for expanding his senses likewise 
 increases. 
 
 2. Our need and our ability to collect and transmit image data 
 are increasing at a rapid rate c 
 
 In fact, the data collection rate is currently held at a level well 
 below capacity simply because of the limits imposed by processing. Thus 
 the extremely large volume of data generated- by a multi spectral source 
 will overburden the available computer resources. For example, it has 
 
 been estimated that data rates of the order of 10 to 10 bits per 
 
 14 
 year can be expected in the near future. A data rate of 10 bits per 
 
 year is equivalent to covering the surface of the earth every 18 days in 
 six spectral regions with a 200 foot digital sample step. The digital 
 processing of this volume of data will generate a large computer require- 
 ment. 
 
 The main objective of this thesis is the study of the specific 
 problems and procedures relating to the LANDSAT multi spectral imagery-- 
 such as geometric correction of the image, classification and clustering 
 procedures. 
 
 1.1 Remote Sensing of Earth's Resources 
 
 Remote sensing of the environment is rapidly becoming a major 
 
area of research for engineers and scientists throughout the world. 
 Satellites and high altitude aircraft provide ideal platforms from 
 which earth resource data may be gathered. The data is gathered in the 
 form of several spectral images of a particular area of the earth under 
 observation. Each image represents the spatial distribution of electro- 
 magnetic energy as seen through a given spectral "window." Information 
 about a particular area is obtained through the study of the spatial 
 and spectral variations of the data for that area. Thus an important 
 subject before the engineering and scientific community at the present 
 time is the processing of senses which represent tracts of the earth's 
 surface as viewed from above. 
 
 A typical scene consists primarily of regular and/or irregular 
 regions arranged in a patchwork manner and each containing one class 
 of surface cover type. These homogeneous regions are the "objects" 
 in the scene. A basic processing goal is to locate and classify the 
 objects and produce a description of the scene in terms of tabulated 
 results and/or a "type map." 
 
 As in the image processing applications, the locations and 
 spatial features (e.g., size, shape, orientation) of objects are re- 
 vealed by changes in average spectral properties that occur at boundaries 
 But unlike most other applications, the spatial features of an object 
 often have only a weak relationship to its class. Research has shown, 
 however, that many classes can be distinguished reasonably well on the 
 basis of their spectral features, using statistical pattern classifica- 
 tion techniques. 
 
Remote sensing combined with a variety of surface, airborne, 
 and space platforms can now give us an integrated overview previously 
 difficult, if not impossible, to obtain. In turn we are learning to 
 extract from these data a new class of information about the earth-- 
 information necessary for making decision in what we might term "planetary 
 engineering"; we foresee the development of interrelated physical methods 
 of the earth and its environment that will predict both the natural 
 course of resources and environmental conditions as well as the effects 
 of human actions. 
 
2. THE LANDSAT PROJECT 
 
 The primary objective of the project [1] is to conduct experi- 
 ments in the interpretation and utilization of earth-resources survey 
 data from space. Doing this means collecting substantial quantities of 
 data about the earth's surface from orbital altitudes, and then trans- 
 formitting it in suitable format to the scientists who will analyze it 
 and determine if it can yield useful information. Thus the knowledge 
 gained from the application of data acquired by the LANDSAT will point 
 the way toward development of fully operational and more effective sys- 
 tems for earth resources management. 
 
 In late July 1972, the National Aeronautics and Space Admin- 
 istration (NASA) launched a one-ton satellite into a sun-synchronous 
 orbit some 560 miles up with a period of 103 minutes. Every period the 
 satellite crosses the equator at about 9:30 a.m. local time on the day- 
 light side. The spin of the earth during each satellite period causes 
 the ground under the orbit to slip westward about 1,550 nautical miles 
 each period. After a full day the orbit precision combined with the 
 fact that the orbital period is not an exact submultiple of the day, 
 insures that the satellite traces out a track on the ground in which each 
 orbital sweep lies about 85 nautical miles west of the corresponding 
 sweep of the day before. Thus every 18 days, the satellite will revisit 
 nearly the same ground points, always by morning. The meteorological 
 cloud-cover photographs made by the satellite are higher-grade versions 
 of the ones that are now familiar on the daily television weather pro- 
 grams; they have a resolution of about one kilometer on the surface by 
 
day and five to ten kilometers even at night with the de!ep infrared. 
 2.1 LANDSAT Mission 
 
 To achieve its broad objectives, the mission of the LANDSAT 
 provides for the repetitive acquisition of high resolution data of the 
 earth's surface on a global basis. LANDSAT carries three 5+inch lenses 
 that produce television images of the same scene in three distinct color 
 bands. This is exceptionally good television, with 4,000 lines, about 
 10 times the detail of home or weather-satellite television. It is 
 the RCA developed system, the Return Beam Vidicon (RBV) system. The 
 satellite also has a second image system--!' t is the Hughes Aircraft de- 
 velopment called Multispectral Scanner (MSS). Here the field is scanned 
 not electronically but mechanically by an oscillating mirror. The 
 scheme allows the use of many detectors, one for each of four wave bands, 
 without any problem of registration since eyery detector samples one and 
 the same scanning spot. In addition, the LANDSAT observatory will be 
 utilized as a relay system to gather data from remote, widely distri- 
 buted, earth-based sensor platforms equipped by individual investigators. 
 The data acquired by the total LANDSAT systems will thus permit quanti- 
 tative measurements to be made of earth-surface characteristics on a 
 spectral, spatial, and temporal basis. 
 
 The overall system is illustrated in Figure 1. The Operations 
 Control Center (OCC) is the focal point of the mission orbital operations, 
 There the overall system is scheduled, spacecraft commands are originated 
 and orbital operations are monitored and evaluated. Data Collection 
 
Commands 
 Tracking 
 
 ^^^^ DCS data, 
 
 ^f^ housekeeping RBV 
 
 ■ I .Toiomptry. tracking MSS 
 
 / j jy data, payload video images 
 
 Earth based 
 
 DCS 
 sensing plat- 
 forms 
 
 data 
 
 Remote Ground Receiving Site 
 
 Goldstone 
 
 NTTF 
 
 ALASKA 
 
 ORBIT 
 Determination 
 
 NASA 
 
 LANDSAT 
 
 PROJECT 
 
 OFFICE 
 
 Command 
 
 DCS, 
 
 TLM 
 
 NASCOM 
 and 
 DTS 
 
 Command 
 
 Payload video tapes 
 
 (mailed from Alaska and Goldstone) 
 
 User 
 
 Evalua- 
 tion 
 
 Ground data handling system 
 
 Operations 
 
 Control 
 
 Center 
 
 NASA Data 
 
 Processing 
 
 Facility 
 
 j 
 
 Figure 1 Overall LANDSAT System 
 
System (DCS), telemetry and command data transfer between the OCC and 
 remote sites is accomplished by NASA communications. The NASA Data Pro- 
 cessing Facility (NDPF) accepts payload video data in the form of mag- 
 netic tapes received in real time at the NTTF station via the OCC or 
 by mail from Alaska and Goldstone. The NDPF then performs the video-to- 
 film conversion and correction, producing black and white images from 
 individual spectral bands and color composites from several spectral 
 bands. 
 
 2.1.1 Observatory System 
 
 The elements of the observatory system include the payload 
 systems and the various support subsystems comprising the spacecraft 
 vehicle. Control of observatory attitude to the local, vertical and 
 orbit velocity vectors within 0.7 degree of each axis is achieved by 
 a three-axis active attitude control subsystem. It uses horizon scan- 
 ners for pitch [2] (pitch variations cause a change in the vertical 
 scale by changing the vertical size of the frame. The magnitude is 
 Ay = he while 9 is the pitch variation from the top to bottom of the 
 frame) and rol 1 (roll causes a skew in the horizontal direction of magni- 
 tude Ax = he where h is the nominal altitude, 6 is the roll variation 
 r r 
 
 over the frame, and Ax is the amount of horizontal skew) control, and 
 a gyro-compassing for yaw (yaw variations are due to errors in the alti- 
 tude control system) orientation. An independent passive Attitude Mea- 
 suring Subsystem (AMS), operating over a narrow range of about 2 degrees 
 provides pitch and roll attitude data accurate to within 0.07 degree to 
 
8 
 
 aid in image location. Orbit-adjustment capability is furnished by a 
 monopropellant hydrazine subsystem employing one pound force thrusters. 
 This system is used to remove launch vehicle injection errors, and to 
 provide periodic trim to maintain a precise orbit. 
 
 Payload video data are transmitted to ground stations over 
 two wideband S-band data links. 
 
 2.1.2 Payloads 
 
 2 . 1 . 2 o 1 Return Beam Vidicon Camera 
 
 The RBV camera system operates by shuttering three independent 
 cameras simultaneously, each sensing a different spectral band in the 
 range of 0.48 to 0.83 micrometers. Since these are visible wavelengths, 
 the RBV is operated only in daylight. The viewed ground scene, 100 by 
 100 nautical miles in area, is stored on the photosensitive surface of 
 the camera tube and, after shuttering, the image is scanned by an elec- 
 tron beam to produce a video signal output. Each camera is read out 
 sequentially, requiring about 3.5 seconds for each of the three spectral 
 images. To produce overlapping images along the direction of spacecraft 
 motion, the cameras are reshuttered eyery 25 seconds. The video band- 
 width during readout is 3.5 MHz. Orientation of the three camera heads 
 is shown in Figure 2. 
 
 2.1.2.2 Multi spectral Scanner 
 
 The Multispectral Scanner (MSS) is a line scanning device which 
 uses an oscillating mirror to continuously scan perpendicular to the 
 
Three RBV 
 
 Camera Mounted in 
 Spacecraft 
 
 185 km x 185 km 
 (100 urn x 100 urn) 
 
 Direction of 
 Flight 
 
 Figure 2 RBV Camera Head Orientation 
 
10 
 
 spacecraft velocity as shown in Figure 3. Six lines, with the same 
 bandpass, are scanned simultaneously in each of the four spectral bands 
 for each mirror sweep. Spacecraft motion provides the along-track pro- 
 gression of the six scanning lines. Optical energy is sensed simultan- 
 eously by an array of detectors in four visible spectral bands from 0.5 
 to 1.1 micrometer for daylight operations of LANDSAT. A fifth band in 
 the near (thermal) infrared from 10.4 to 12.6 micrometers is installed. 
 The detector outputs are sampled, encoded to six bits and formatted into 
 a continuous data stream of 15 megabits per second. During image data 
 processing in the ground data handling system facility, the continuous 
 strip imagery is transformed to framed images with a 10 percent overlap 
 of consecutive frames and an area coverage approximately equal to that 
 of the RBV images. 
 
 2.1.3 Wideband Video Tape Recorders 
 
 The use of data from the RBV and MSS sensors are complementary 
 in several respects, and both sensors are generally operated simultan- 
 eously over the same terrain during daylight hours. When operated over 
 a ground receiving station, their data are transmitted in real time to 
 the ground receiving site and recorded there on magnetic tape. 
 
 2.1.4 Data Collection System 
 
 The Data Collection System (DCS) obtains data from remote, 
 automatic data collection platforms, which are equipped by specific in- 
 vestigators, and relays the data to ground stations whenever the LANDSAT 
 
optics 
 
 11 
 
 -scan mirror 
 (oscillates 
 nominally ± 2.88°) 
 
 6 Lines/Scan/Band 
 
 'it. ■' 
 
 Path of spacecraft 
 
 Figure 3 MSS Ground Scan Pattern 
 
12 
 
 spacecraft can mutually view any platform and any one of the ground sta- 
 tions. Each DCS platform collects data from as many as eight sensors, 
 supplied by the cognizant investigator, sampling such local environ- 
 mental conditions as temperature, stream flow, snow depth, or soil mois- 
 ture. Data from any platform is available to investigators within 24 
 hours from the time the sensor measurements are relayed by the spacecraft, 
 
 2.1.5 LANDSAT Multi spectral Scanner Data 
 
 LANDSAT multi spectral scanner data is received from the satel- 
 lite by NASA, processed, and delivered to users recorded on computer 
 compatible tape and in photographic form. The computer tape form of 
 the data is calibrated and line length adjusted by NASA but no geometric 
 corrections are applied. The system-corrected photographic products 
 are corrected for many geometric distortions including earth rotation 
 effects in addition to the above two corrections. Also, these images 
 are rescaled so that the horizontal and vertical scales are the same. 
 Thus, the digital form of the MSS data contains many geometric distor- 
 tions and users of this data are faced with the problem of compensating 
 for these errors. 
 
 When digital MSS data is reproduced in image form on a stan- 
 dard IBM computer line printer the resulting scale factor is approxi- 
 mately 1 in. = 22,400 in. in the horizontal direction and 1 in. = 
 24,200 in. in the vertical direction. This scale differential exists 
 in addition to the skew due to earth's rotation and all other geometric 
 errors. Similarly when this data is reproduced on a video display 
 
13 
 
 device the horizontal scale is 56 meters per point and the vertical 
 scale is 79 meters per point. 
 
 2.1.6 MSS Digital Data Geometric Characteristics 
 
 The LANDSAT MSS System produces four spectral band digitized 
 imagery of approximately 100 nautical mile (185 km) wide strips beneath 
 the satellite path. The scanner has an instantaneous resolution of 79 
 meters and scan lines are sampled at a rate such that samples are spaced 
 approximately 56 meters apart and successive scan lines are spaced ap- 
 proximately 79 meters apart as determined by the forward motion of the 
 satellite. The image data are edited so that the along size is approxi- 
 mately 96.3 nautical miles (155 km). The resulting data set consists 
 of nominally 3240 samples horizontally (E-W) and 2340 samples vertically 
 (N-S). The geometric distortions are due to sensor, satellite, and 
 earth effects. 
 
14 
 
 3. COMPILATION AND GEOMETRIC REGISTRATION 
 OF LANDSAT MULTI TEMPORAL IMAGERY 
 
 Within the context of computer-assisted earth-resources data 
 analysis, image registration objectives fall generally into two separate 
 (but procedurally related) categories: 
 
 1. A first set of image registration objectives has to do 
 with the need for accurate spatial registration of imagery with respect 
 to conventional geographic coordinate systems. Generally, the goal of 
 procedures here is to relate earth-resources imagery of interest to map 
 bases in such a manner that they can be associated with each data ele- 
 ment of a multispectral scanner (MSS) image (multivariate observation, 
 resolution element, pixel) an accurate ground position in terms of a 
 conventional map coordinate system, typically latitude-longitude. Such 
 image-to-map registration is necessary, for example, to enable conven- 
 ient cross referencing of image data with other geographically-referenced 
 ground-collected information within automated-classifier training opera- 
 tions and other image analysis procedures. Image-to-map registration is 
 also required where image data are to be output in graphical format 
 suitable for photo-reproduction as topo map overlays. In this instance, 
 not only must the mathematical function relating image coordinates to 
 map coordinates be determined, but also the data of the image must be 
 physically reformatted such that a graphical display of the output data 
 via some photo-recording device conforms in area, orientation, and pro- 
 jection to a specific map. 
 
 2. A second category of objectives within image registration 
 operations relates to the automatic overlay of multiple images (usually 
 
15 
 
 two) taken at different times over the same ground area in creating 
 multitemporal multispectral image data sets. In the case of imagery 
 collected by LANDSAT, conventional mathematical correlation techniques 
 may be employed directly to determine points of optimal match between 
 temporally separate images. The goal of digital processing here of 
 course is to determine automatically an element-by-element matching 
 of the pixels of one image with those of the other image(s) such that 
 a composite multitemporal data set might be produced for a specified 
 ground area common to both (all) images. 
 
 Usually, it is also desirable that the resulting mutli temporal 
 imagery be accurately registered to some specified geographic coordinate 
 system. Where one map-registered image already exists for a specific 
 ground area, subsequent images can be registered automatically to the 
 same map base by simply correlating later images with the initial map 
 registered image and then applying identical image-to-map spatial trans- 
 formation. Such combined uses of image- to- image and image-to-map regis- 
 tration procedures suggest the practicality of automated systems for 
 monitoring natural resource and land use boundary changes over time. 
 
 The remainder of this chapter describes the approach taken at 
 the Center of Advanced Computation, University of Illinois [33 in 
 adapting image registration procedures previously researched at LARS 
 [2,4,5] Purdue University, for more economical processing using the 
 ILLIAC IV and other computers of the ARPA network. 
 
 3.1 A Basic LANDSAT Data Prepr ocessing Problem: Image Rectification 
 
 For a number of reasons, it is generally most convenient to 
 
16 
 
 work with the digital data of LANDSAT MSS imagery in units of map- 
 oriented rectangular image areas. For example, the image analysis area 
 of interest is often the geographic region of a particular map bounded 
 by lines of latitude and longitude. However, the characteristics of the 
 near-polar orbit of LANDSAT are such that the grid of data elements 
 scanned for each MSS image is always rotational ly displaced to a consid- 
 erable degree from alignment with geodetic meridians and parallels. Also, 
 the rotation of the earth under LANDSAT during the scanning of each MSS 
 image introduces an oblique skew to the grid geometry of each image 
 digitalized. Hence, to enable convenient retrieval, analysis, and dis- 
 play of selected MSS imagery in terms of map-oriented rectangular areas, 
 cumbersome digital preprocessing is usually required to transform geo- 
 metrically and reformat the computer tape data of each LANDSAT image to 
 north-vertical image grid orientation. 
 
 This computer preprocessing of LANDSAT MSS data from initial 
 scanner format to map-oriented image analysis format is known generally 
 as image rectification or geometric correction . Again, the two principle 
 objectives of this data preprocessing step are: 
 
 1. To remove the oblique skew within the MSS image grid 
 caused by the earth's rotation during image scanning, 
 
 2. To transform the MSS image grid to a vertical north, map 
 orientation. 
 
 A third objective to be served by geometric correction is sometimes 
 the differential rescaling of the MSS image grid along its transformed 
 axes to achieve equality of horizontal and vertical image scales. 
 
17 
 
 In any case via image-rectification data preprocessing, a new 
 image data tape is created from which map-oriented rectangular areas of 
 data can be retrieved and displayed directly for all subsequent image 
 interpretation purposes. 
 
 To understand the nature of the data processing problem as- 
 sociated with geometric rectification of LANDSAT MSS imagery, consider 
 the following typical image rectification task. We wish to extract from 
 a single LANDSAT MSS image those data elements corresponding to ground 
 points lying within a specific USGS 1° x 2° map quadrangle, and to trans- 
 form that image subarea such that ordered rectangular blocks of data of 
 the rectified image tape conform closely in geometric area and orienta- 
 tion for the USGS 7-1/2 minutes, 15 minutes, and 1° x 2° topo maps exist- 
 ing for that quadrangle. 
 
 Assume that a mathematical function is given that, for each 
 element of the rectified image designated by output image row and column 
 coordinates (r , c ), determines that single element of the original 
 scanner image (r, c) that is most nearly geographically equivalent. We 
 assume here an integer function, say (r,c) = f(r ,c ), that establishes 
 an integer-pair mapping from rectified pixel indices back to "best- 
 corresponding" scanner image indices. (Note that such a mapping between 
 rectified and original image pixels will not in general be one-to-one 
 [next paragraph describes an oblique transformation method to achieve 
 a one-to-one mapping of LANDSAT input picture grid onto a LANDSAT output 
 picture grid] and thus some data elements are usually lost or duplicated 
 in tran format ion.) 
 
18 
 
 Having a particular rule for locating pixels of the original 
 image to be copied into specific pixel positions within the rectified 
 image format, we then wish to apply this reformatting rule for all pixels 
 lying within the given 1° x 2° map quadrangle. Figure 4 illustrates the 
 impractical ity of core-contained computational strategies for this 
 problem. 
 
 We would like to input from magnetic tape successive portions 
 of the oribinal scanner image, and following transformation, output se- 
 quentially onto another tape individual complete lines of the rectified 
 image. The main difficulty arises with the size of the real -core memory 
 required to hold that portion of the input image containing all data ele- 
 ments needed to compose the next line of the output image. Rectifica- 
 tion of the 1° x 2° quadrangle of LANDSAT data as depicted in Figure 4 
 would require typically a real core input data buffer on the order of 
 7.10 bytes . 
 
 Confronted with the sheer quantity of data comprising each 
 LANDSAT MSS image and given such computational problems, most LANDSAT 
 data analysts reexamine the necessity of preprocessing LANDSAT image 
 tapes to achieve map-oriented data format. At this point, it is usually 
 decided that the goals of image-to-map registration should be limited to 
 only the determination of those numerical transformation necessary for 
 referencing individual data elements by map coordinates leaving completely 
 undone the actual reformatting of LANDSAT image data tapes to map orien- 
 tations. In many instances where geo-referencing and graphical output is 
 only identical to data analysis objectives, such a strategy seems entirely 
 appropriate. 
 
19 
 
 7M bytes 
 (7.106) 
 of data 
 
 A 
 
 Original LANDSAT MSS Image 
 
 Rectified 1° x 2° Quad Image 
 
 Figure 4 Buffer Requirement for Core-Oriented Rectification of 
 LANDSAT Images 
 
20 
 
 However, a fundamental problem remains In order to determine 
 adequately (to one-pixel average accuracy) the numerical parameters 
 needed for accurate geo-referencing, some initial set of control points 
 must be identified simultaneously in terms of both LANDSAT image pixel 
 position and map position. Such geographic control points most conven- 
 iently obtained for LANDSAT imagery over the U.S. by manually correlat- 
 ing small area line-printer displays of rectified LANDSAT image data 
 with USGS 7-1/2 minute topo maps and digitizing some small set of points 
 denoting overlay positions of maximum geographic correspondence. 
 
 Here a reasonable approach might seem simply to retrieve from 
 tape and rectify only those small areas of image data needed for manual 
 correlation with quadrangle maps. However, the rectangular latitude- 
 longitude grid of the USGS map series is such that a uniformly distri- 
 buted set of image sub-areas, selected to fall within 7-1/2 minute map 
 quadrangles, is such more efficiently retrieved from an image data tape 
 of recti fied-image format. Since also there will be needed some sys- 
 tematic grid of image subareas for automated digital correlation within 
 any subsequent mul ti temporal image registration procedures, it would 
 seem that a number of image registration objectives would be served by 
 some efficient computational strategy for initial transformation of all 
 LANDSAT imagery to vertical-north orientation. Such an efficient LANDSAT 
 data preprocessing operation is provided by the oblique (skew) trans- 
 formation. 
 
21 
 
 3.2 Oblique Transformation of LANDSAT Images to Approximate North- 
 South Orientation 
 
 The oblique transform [6] achieves a one-to-one mapping of a 
 LANDSAT input picture grid onto a LANDSAT output picture grid that 
 closely approximates a skewed rotation of the input picture with scale 
 factored out. A one-to-one mapping is desirable since a further cor- 
 rection step, dependent on a particular map projection and some set of 
 specific image-to-map control points, will be applied to the obliquely- 
 transformed output picture. The major purpose of the oblique transform 
 is to reduce the size of the random-access memory needed for application 
 of that final correction step. As a bonus, it turns out that a line 
 printed display (assuming 10 characters per inch horizontally and 8 
 characters per inch vertically) of the obliquely transformed LANDSAT 
 image has a scale and geometry quite close to that of the USGS seven- 
 and-a-half-minute quadrange map series. The oblique transformation is 
 equivalent to a vertical skew of the input picture followed by a hori- 
 zontal skew. See Figure 5. Each input column is shifted an integral 
 number of pixels, the shift distance being the rounded off values of 
 - 3 x for that column. The resulting rows are now shifted horizontally 
 by the rounded off value of ay for that row. It is clear that with this 
 process, no input points are thrown away or duplicated. The question is 
 now the relationship between the oblique transform and a general linear 
 transform. 
 
 The two parameters of the oblique transform are the skew fac- 
 tors a and 3. Since all input and output pixels must fall within integral 
 
22 
 
 • • ■ » • m 
 
 (a) 
 
 (ay,y) 
 
 (c) 
 
 Figure 5 The Oblique Transformation Process 
 

 23 
 
 grid location, the integral or "digitized" oblique transform is: 
 
 x = x - ay + 1/2 
 
 ' I ' I (1) 
 
 y = [Bx + 1/2J+ y 
 
 with x and y output pixel coordinates, x and y input pixel coordinates, 
 and where |_aj denotes the greatest integer ^ a. Thus, [a_ + 1/2J is theo- 
 retically equivalent to rounding a^ to the nearest integer. The trans- 
 formation (1) is the inverse of the process displayed in Figure 5 to 
 reflect the way picture transformations are normally obtained: for each 
 output pixel, the input pixel that transforms to it is found and assessed. 
 The continuous equivalent of (1) is the linear transform: 
 
 1 
 x = x - ay 
 
 (2) 
 
 y = Bx + (1 - a3) y 
 
 Note that rounding off x and y values in (2) does not necessarily re- 
 sult 1n an information-lossless transform. For example, with a = 1/2 and 
 
 1 1 
 B = - 2, both coordinates (1,1) and 2,2) are transformed into (x ,y ) = 
 
 (l,o). 
 
 The general linear transform may be factored into an oblique 
 transform (2) followed by a scale transformation. One factorization is: 
 
 (3a) 
 
 If there is no solution to (3a), i.e., a = then the alternate factor- 
 ization is: 
 
24 
 
 (3b) 
 
 Given any non-singular linear transform, i.e., ad-bc f o, a, 3> m and n 
 ncessarily exist and can be easily determined. For more details see 
 reference [2], 
 
 3.2.1 Error Bounds 
 
 The maximum errors incurred by using (1) instead of (2) are: 
 
 6x = ay - [ay + 1/2J (4a) 
 
 6y = 3(x - ay) - [g(x - [ay + 1/2_J ) + 1/2] (4b) 
 
 Any real number (u) can be expressed as m + e, where m is an integer, 
 and -1/2 ^ e < 1/2. Thus, we can set ay = m + e, and obtain from (4a): 
 
 <5x = m+£-[m+£+l/2j =m+e-m = e 
 Thus, J 6x J 1 1/2. For (4b), we set x - ay = m + e-| 
 where - 1/2 < £] - 1/2 instead. Then 
 
 x - [ay + 1/2] = x - [x-m-E 1 + 1/2J = m, and (5a) 
 
 we get 6y = 3m - 3m + 1/2 + 3£i = e + $e, 
 
 1 — — ' ' ' (5b) 
 
 |5y| < 1/2 (|3|+D 
 
 From theoretical consideration of the LANDSAT orbit, the actual 3 re- 
 quired to bring LANDSAT images over the continental United States into 
 north-south alignment does not exceed .17 in absolute value, giving a 
 maximum of error of 0.59 of a pixel's vertical dimension. 
 
25 
 
 3.2.2 Implementation Consideration 
 
 Magnetic tape used at 1600 bps is currently the only reasonable 
 way to store LANDSAT images. Small disk packs (IBM 2314 equivalent) will 
 not hold a complete LANDSAT picture and large ones (IBM- 3330 equivalent) 
 are relatively scarce and expensive. Since typically images are mailed 
 via tape anyway, the present implementation of the oblique transformation 
 procedure assumes magnetic tape storage of both the output and the input 
 images. The main restriction, therefore, is that all access to the data 
 of both input and output pictures is necessarily sequential , i.e., a 
 row or line at a time. 
 
 Whenever the 3 parameter of the oblique transformation (1) is 
 non-zero, the information of more than one input line must be more or 
 less randomly accessible for the creation of one output line of data. 
 That means input lines must be buffered. The size of that buffer is 
 clearly the maximum number of input lines needed for any output line 
 multiplied by the memory needed per input line, i.e., (3. width) (width) 
 (Kbytes/pixel). For a typical U. S. LANDSAT frame, 3 = 0.17, width = 
 3240 pixels per row, and K = 4, 8-bit bytes per pixel, yielding a 
 required buffer size of slightly more than 7.10 bytes. Note if we 
 perform the transformation or quarter-frames rather than the entire 
 image, all operations are done four times, but memory requirements will 
 be reduced by a factor of approximately 16. Furthermore, computation 
 time using such multiple passes will not be significantly higher since 
 the same number of inner-loop operations are done in either case. 
 
26 
 
 Specific implementation of the oblique transformation pro- 
 cedure should be designed with respect to "real core" availability, since 
 assessing the full 7M bytes buffer as "virtual core" on most virtual 
 machines would create a large number of page transfers or faults. The 
 great increase in system overhead that causes would reduce enormously 
 the efficiency of the procedure. 
 
 It is interesting to compute the approximate number of page 
 faults that would occur; each output line causes at least 3- width-w 
 faults, when w is the working set size or the number of pages storable 
 in real core. If we assume 2340 output lines, use the same width and 
 3 as before, and use w = 100 (equivalent to a 400K byte working set in 
 real core on an IBM 360/67), about 1,300,000 faults would occur. 
 
 Finally, even if the input width is cut down as above, we 
 are still faced with the condition that tape input and output is more 
 efficient as a real-core machine, especially when large tape input and 
 output buffers are used. 
 
 3.2.3 Details of the Real -Core Machine Implementation 
 
 The process of oblique transformation is broken up into several 
 processing steps, each step transforming a successive part of a LANDSAT 
 input picture. Thus, each pass has one output (the freshly-updated out- 
 put picture) and two inputs (the appropriate LANDSAT MSS quarter-frame 
 tape as received from Goddard and the output picture from the immediately 
 preceding pass). 
 
27 
 
 The part of a LANDSAT tape that is transformed in each pass is 
 a vertical strip of it, e.g., the entire Goddard tape if sufficient real 
 core is available. Within each step then, for each output line, the 
 following things are done. First, the output line is "clipped" into the 
 largest segment that, transformed, lies wholly within the current input 
 strip. 
 
 Second, if that segment is non-empty, then the required lines 
 of the input picture are read in and queued into the "back" part of the in- 
 put buffer. Note that due to the linear nature of the oblique transfor- 
 mation, lines previously queued which are not needed for the current output 
 line will necessarily be at the "front" part of the input buffer and can 
 be discarded systematically with the knowledge that their information 
 will not be needed again in the current pass. That fact allows us to 
 efficiently keep the size of the queue (the input lines buffer as dis- 
 cussed above) just that needed to hold all input line data for a single 
 output segment. 
 
 Finally, the accumulated data for the current output line is 
 read from the tape output from the last pass, and for each point in the 
 output segment, the inner loop of the oblique transformation updates 
 that part of the output line by moving the corresponding image data from 
 the input line. The new output line is then written out onto another 
 tape. 
 
 3.2.4 Precision Geographic Registration of LANDSAT Imagery 
 
 Initial rectification of all LANDSAT imagery to map-oriented 
 
28 
 
 data formats greatly facilitates all subsequent precision image-to-map 
 registration efforts. Once approximate registration has been determined 
 for an image and data tapes re-formulated via skew transformation, pre- 
 cision image-to-map registration may be conveniently accomplished inthe 
 following manner: 
 
 The central problem confronting all precision image-to-map 
 registration efforts is to identify accurately, in terms of both digital 
 image pixel coordinates and map coordinates, a geographic network of 
 contral points visible both on the image and on maps. For LANDSAT image- 
 ry over the U.S., experience suggests that such a system of control 
 points can be established by overlaying on a light table small line- 
 printer displays of rectified imagery with USGS 7 1/2 quadrangle maps 
 at 1:24,000 to determine manually image-to-map overlay positions of 
 maximum geographic correspondence. 
 
 The fact that, the horizontal and vertical scales of standard 
 line-printer displays of rectified LANDSAT imagery (8 pixels to the 
 inch vertically and 10 pixels per inch horizontally) overlay USGS 7 1/2 
 minute quadrangles to approximately one-pixel accuracy contributes much 
 toward the convenience of such a strategy. Following quick visual cor- 
 relation of line printer display and map quadrangle, a single point of 
 geographic correspondence between image and map can be established simply 
 by recording the row and column coordinates of a central pixel of the 
 display and measuring the position of the center of this pixel in terms 
 of map coordinates. The availability of a data-tablet digitizer for 
 such a procedure greatly facilitates the establishment of a network of 
 
29 
 
 control points by this method. In selecting image areas for determina- 
 tion of control points, essentially two strategies are available. 
 
 We have noted that the skew-transformed image tape now stores 
 approximate map quadrangles of imagery as ordered rectangular blocks of 
 data since we will need some grid of image blocks (64 lines by 64 col- 
 umns) for any subsequent image-to-image registration procedures that may 
 be required. It is convenient in several respects to establish for each 
 rectified image a grid of image correlation blocks centered approximately 
 on USGS 7 1/2 minute quadrangles to be used throughout all image-to-map 
 and image-to-image registration procedures. Once this quad-centered 
 grid of image blocks has been determined and retrieved from tape, line- 
 printer displays (MSS channel two) can be produced for only a small addi- 
 tional cost. Quick visual inspection of these displays with reference 
 to available quadrangle maps, should suggest some subset of image blocks 
 to be used in establishing precision geographic control. However, such 
 a strategy is appropriate only for those geographic regions for which 
 USCS 7 1/2 mapping is relatively complete. 
 
 In those geographic regions for which 7 1/2 quadrangle maps 
 are unavailable, image-to-map registration must proceed in more conven- 
 tional fashion. Some sets of candidate control points must be selected 
 manually by the data analyst with reference to suitable terrain features 
 identifiable on both film image and available maps. A set of small 
 image blocks containing these candidate control points is retrieved from 
 tape and displayed in line-printer image format. Then, whenever it is 
 meaningful to do so given available materials, individual image pixels 
 
30 
 
 are recorded in terms of both yow and column positions and geographic 
 coordinates. 
 
 Our initial experience suggests that to achieve one-pixel 
 average registration accuracy over a full LANDSAT image, approximately 
 two dozen widely distributed control points should be established by 
 either of those methods. 
 
 3.2.5 Compilation of LANDSAT Multitemporal Imagery 
 
 In overlaying multiple LANDSAT images to achieve multitemporal 
 data sets, points of best match between the two images to be overlaid 
 must be determined. These match points may be determined digitally by 
 computing points of maximum image correlation between small blocks of 
 data paired between the two images. Following the LARS image- to- image 
 registration methodology [5] such a procedure may be implemented as fol- 
 lows. 
 
 Given two temporally separate LANDSAT images over the same 
 ground area, initial geographic control allows pairs of image correla- 
 tion blocks to be selected from the two images such that spatial dis- 
 placement between paired blocks are only on the order of 10 pixels. 
 One "floater" block of data taken from one image is translated hori- 
 zontally and vertically through all possible positions of pixel-to- 
 pixel overlay with the corresponding (and slightly larger) "stationary" 
 block taken from the other image. Numerical values of image correla- 
 tion are computed for all overlay positions, and thus a two-dimensional 
 correlation surface is determined. A single point of spatial correspond- 
 ence between the two image blocks is then established automatically at 
 
31 
 
 the peak of this surface. Repeated for a complete grid of blocks paired 
 between two images such a procedure outputs a computer-determined system- 
 atic set of control points which can then be employed directly in cali- 
 brating the spatial transformation required for overlaying the two images 
 as a single multi temporal image. 
 
 For typical applications, it is sufficient to choose floater 
 blocks from one image as 32-line by 32-column data rectangles and 
 stationary blocks from the other images as 64 by 64 rectangles. Allow- 
 ing all possible translations of floaters with respect to stationary 
 blocks, 1089 (33 x 33) individual image correlations are computed for 
 each block pair. 
 
 Also, it is sufficient to space correlation blocks approximate- 
 ly 180 pixel positions apart, both horizontally and vertically, i.e., 
 the approximate dimensions of a 7 1/2 minute map quadrangle. Therefore, 
 for two LANDSAT images perfectly coincident over the same geographic 
 area, correlation would be computed for approximately 230 block pairs. 
 This corresponds to more than 250,000 individual image correlation cal- 
 culations. Since there image correlation calculations represent the 
 majority of all computations needed for registration of two images, and 
 since these calculations can be performed in parallel computation for- 
 mat, it is convenient to isolate these correlation calculations for 
 execution on the ILL (AC IV. This can be done by retrieving from both 
 image tapes of UCLA two quad-centered grids of image correlation blocks 
 paired between the two images and transmitting these two files via the 
 ARPA network to ILLIAC IV at NASA/ Ames. Following the ILLIAC IV calcu- 
 lations of all block correlation peaks, only minor TENEX processing is 
 
32 
 
 required to calculate the numerical parameters of the two-dimensional 
 least squares quadratic polynomial needed for pixel-to-pixel overlay of 
 the images. Given the parameters of this polynomial, the UCLA IBM 360/91 
 may be again employed efficiently for re-formatting the two separate 
 images as a single 8-channel multi temporal image. 
 
 3.2.6 The EDITOR . . .ERTS Data interpreter and TENEX Operation Recorder 
 
 The EDITOR system represents a component within a more exten- 
 sive computer software development effort undertaken at CAC to facili- 
 tate large-scale ILLIAC IV - ARPA Network analysis of multispectral 
 reconnaissance imagery such as that collected by the Earth Resources 
 Technology Satellite--ERTS (currently called LANDSAT) of NASA and USGS/ 
 DI. [9,10] This software effort has been undertaken by CAC, in collab- 
 oration with the Laboratory for Applications of Remote Sending (LARS) 
 of Purdue University, in support of the LANDSAT/EROS programs of NASA 
 and USGS/DI, in support of the ILLIAC IV and ARPA Network applications 
 programs of NASA and ARPA/DOD, and in support of U. S. agricultural 
 crop-acreage monitoring objectives of SRS/USDA. Image analysis algorithms 
 initially implemented within these efforts follow closely software pre- 
 viously researched at LARS as multispectral image interpretation tech- 
 niques. [11,12] 
 
 EDITOR has been conceived as an interactive job-editing and 
 file-management interface to ERTS data analysis systems being implemented 
 for the ILLIAC IV hardware complex at NASA's Ames Research Center (Ames) 
 at Moffett Field, California. Nationally accessible via the ARPA Network, 
 
33 
 
 the EDITOR system allows portable, dial-up terminal command of ILLIAC 
 IV image processing facilities at Ames. Multispectral image data manage- 
 ment services within EDITOR have been designed to optimize the efficiency 
 of ERTS data file transfers between the computers of the ILLIAC IV com- 
 plex at Ames (the ILLIAC IV parallel processor, its peripheral TENEX 
 PROCESSORS, AND THE UNICON Data Computer) and other computers on the ARPA 
 Network. 
 
 Since the present EDITOR system, using tape inputs, has been 
 
 developed for stimulation of data management systems to be implemented 
 
 12 
 using the 10 -bit laser memory of the UNICON Data Computer at Ames, 
 
 specific data tape formatting is required for all multispectral images 
 
 to be processed by EDITOR. Required tape formats are described in CAC 
 
 Technical Memorandum No. 19. [13] 
 
 While conceived initially as a peripheral TENEX file-management 
 and job-editing interface to ILLIAC IV batch image interpretations, the 
 scope of EDITOR functions has been widened subsequently to include pro- 
 cedures for direct TENEX interpretation of small-scale ERTS data samples 
 in preparation for ILLIAC IV batch analyses. Thus, where only small- 
 scale image interpretations are needed and where portable terminal output 
 is sufficient, EDITOR may be used alone as a self-sufficient tape- 
 operating multispectral image analysis system. 
 
 Since EDITOR has been developed as an interactive interface to 
 ILLIAC IV image processing facilities, the system should also prove 
 advantageous to a wider community of ARPA Network picture processing 
 researchers. Where image disk files conform to EDITOR data file formats, 
 
34 
 
 the system may be used for image file management, editing, and display 
 in conjunction with more general ILLIAC IV image processing activities. 
 Specific file formats assumed by EDITOR are described in CAC Technical 
 Memorandum No. 19. [13] 
 
 The EDITOR system described here may be executed on any ARPA 
 Network TENEX system having 7-track or 9-track magnetic tape drives. 
 Currently, a copy of the EDITOR system is maintained by CAC at Bolt, 
 Beranek and Newman (BBN) in Cambridge, Massachusetts, in directory 
 <RAY>. The BBN EDITOR can be made available to any user having access 
 to the ARPA Network. The file <RAY>TAPES. SHARED is maintained there by 
 CAC describing all tape imagery available at BBN for experimental system 
 usage. 
 
 3.2.7 System Overview 
 
 The actual use of EDITOR is best described by way of examples. 
 The following verbal description of EDITOR simply sketches major system 
 functions. (More detailed descriptions of all presently available 
 system procedures are given in Section 3.2.9.) A familiarity with basic 
 multispectral image interpretation procedures is assumed. 
 
 EDITOR is "implemented in modular format with each system 
 module addressed by a single command. For each command, only that num- 
 ber of characters sufficient to distinguish the command from all others 
 need be typed by the user. 
 
 The RETRIEVE command is used to create TENEX disk files con- 
 taining multispectral image data from single rectangular areas of an 
 LANDSAT image, or from sets of rectangular areas within an image. Use 
 
35 
 
 of RETRIEVE assumes that, prior to entering EDITOR, the user has request- 
 ed the TENEX operator to mount a specific image tape. Analysis areas to 
 be retrieved from tape may then be specified within EDITOR by the user 
 in terms of image line and column coordinates and line and column sampl- 
 ing increments. 
 
 If, on the other hand, the user has an image file already 
 resident on the TENEX disk and wishes to use only some portion of it for 
 analysis, he may use the CLIP A WINDOW command. This command allows the 
 user to construct new windows by specifying new image line and column 
 coordinates for the subwindows desired. Additional image sampling (by 
 additional line and column sampling increments) is also possible within 
 CLIP. 
 
 Once image analysis areas have been selected and stored on the 
 TENEX disk, the PRINT command within EDITOR allows a user to view the 
 image disk files created to verify that the desired data windows have 
 indeed been retrieved. Gray-scale displays of any of the individual 
 spectral bands of data within an image file may be achieved by overprint- 
 ing characters on the user's terminal. A data histogramming routine is 
 available for displaying the proportional distribution of pixel intensi- 
 ties within any spectral band of an image file. This feature of the 
 system can be used for terminal display enhancement purposes. With 
 reference to image data histograms, a user can select manually a partic- 
 ular range of image point intensities to be assigned to each gray-scale 
 shade available within the character-overprinting terminal display facil- 
 ity. 
 
36 
 
 Once a user has verified that the desired imagery has been 
 selected, several options are available for data analysis. 
 
 The CLUSTER command within EDITOR evokes a multivariate cluster 
 analysis of all multispectral image data points within a particular file 
 to determine a logical partitioning or grouping of all data points into 
 a user-specified small number of spectral clusters. [7,14] The cluster 
 analysis procedure may be regarded as an unsupervised image interpre- 
 tation procedure categorizing a set of image points into most-separable 
 spectral classes. Also, following any cluster analysis procedure, a set 
 of summary statistics characterizing the multivariate distributions of 
 data within all clusters (i.e., mean vectors and variance-covariance 
 matrices) is computed and, at the request of the user, displayed on the 
 terminal device. 
 
 The CLASSIFY command within EDITOR allows statistical classifi- 
 cation of all data points within a specified disk file into a particular 
 set of spectral ly-distinct categories determined by a previous cluster 
 analysis. 
 
 This methodology of classification assumes that "ground-truth" 
 samples will always be cluster-analyzed prior to class signature compu- 
 tation in order to verify the existence of spectral homogeneity within 
 nominal terrain classes and spectral separability between classes. Where 
 spectral properties of ground- truth samples are not homogeneous within 
 nominal classes (i.e., multivariate distributions are non-spherical), the 
 multiple spectral classes determined by cluster analysis for such terrain 
 classes may be used for analysis purposes as spectrally-distinct, but 
 
37 
 
 nominally-synonymous, classes to be re-grouped following classification 
 into single nominal terrain classes. 
 
 In accordance with the LARS multispectral image interpretation 
 methodology, the statistical classification algorithm evoked by the com- 
 mand CLASSIFY employs the Gaussian maximum-likelihood decision rule . 
 [11,7,8] For each data point of a file of data points to be classified, 
 discriminant functions for all classes are computed and the point is 
 assigned to that class for which the discriminant is largest. The param- 
 eters of the discriminant function for each class are the estimates of 
 the means, variances and covariances of the spectral properties of all 
 data points belonging to that class. These are the spectral class sta- 
 tistics computed and saved by EDITOR immediately following any cluster 
 analysis for the set of spectral classes determined by that particular 
 CLUSTER procedure. 
 
 The ENTER STATISTICS FILE EDITOR command provides a set of 
 file-management procedures convenient for interactive display and edit- 
 ing of the spectral class statistics accumulated from previous cluster 
 analyses. Within this mode, class signatures of different statistics 
 files may be re-combined in any manner to produce new composite sets of 
 signatures to be used for classification. Also within this mode, meas- 
 ures of spectral separability or distinctness may be computed and dis- 
 played for all pairs of signatures within any statistics file. A pro- 
 cedure is also implemented (again following LARS methods) for "pooling" 
 the statistics of several different signatures into a single new compos- 
 ite signature. 
 
38 
 
 The EXECUTE ANALYSIS ON ILLIAC IV command instructs EDITOR to 
 use the ILLIAC IV at Ames for a particular cluster analysis or classifi- 
 cation task specified by a user. Given this command, EDITOR immediately 
 requires the user's I4-TENEX user code and password. If such a code can 
 be supplied, EDITOR will automatically set up the task for ILLIAC IV 
 processing, transmit the job via the ARPA Network to Ames, and enter the 
 job into the ILLIAC IV batch stream. Procedures are also provided within 
 this mode for inquiring on the status of ILLIAC IV jobs and for conven- 
 ient retrieval by EDITOR of cluster analysis and classification output 
 files resulting from previous ILLIAC IV executions. Specific ILLIAC IV 
 mul tispectral image cluster analysis and classification procedures acces- 
 sible to EDITOR users have been documented elsewhere. [15,16,17,18] 
 
 In addition to the procedures for displaying raw mul tispectral 
 imagery, the PRINT command within EDITOR also provides procedures for 
 displaying the interpreted imagery that results from either cluster 
 analysis or classification operations. Interpreted image data files may 
 be output as terminal printer maps where the symbols 1 - 9 and A - Z are 
 used to represent terrain classes. Display of only selected subsets of 
 classes is permitted. Following again the LARS methodology, provisions 
 are included for reliability thresholding of classification results with- 
 in PRINT mode. Where a particular confidence level is specified by the 
 user, EDITOR displays only those image data points whose classification 
 reliability meets the specified level. Also, procedures are included 
 within PRINT mode to allow re-grouping of spectral classes into composite 
 
39 
 
 nominal classes, and to allow manual re-assignment of printer symbols 
 among these nominal classes for more meaningful display of interpreta- 
 tions. 
 
 3.2.8 Getting Started with EDITOR and TENEX 
 
 In this section a brief description of how to run EDITOR 
 under TENEX will be presented. TENEX is a time-sharing system running 
 on a modified DEC SYSTEM-! as developed by Bolt, Beranek and Newman 
 (BBN). [19] When connecting to a TENEX system, the user should estab- 
 lish a character-at-a-time, full-duplex connection. Commands to TENEX 
 may be abbreviated and if followed by ESC (ESCAPE or ALT MODE on some 
 terminals) the remainder of the command will be typed out followed by 
 a prompt for any parameters required. All commands should be terminated 
 by a carriage return. 
 
 To begin an EDITOR session, one must first LOGIN. The se- 
 quence is LOG <user code> <password> <account number> followed by a 
 carriage return. (The password will not be echoed to preserve secur- 
 ity.) If a space i c . substituted for <account number>, that account num- 
 ber associated with the user code given will be assumed by default. 
 
 To get an ERTS image tape mounted, a request must be made to 
 the operator who should be logged in as OPERATOR. To do this, the LINK 
 command is used. The LINK command links two terminals together so that 
 whatever is typed on either may be viewed on both. The syntax of the 
 command is LINK OPERATOR. To communicate with the operator, a message 
 should then be typed. Each line of this message must begin with a semi- 
 colon (";") to indicate to TENEX that the line is not a command. The 
 
40 
 
 request should be to mount a specific tape, without write ring (to in- 
 sure that the data on the tape cannot be accidentally over-written), on 
 a 9-track or 7-track drive. Once the operator has mounted the tape, he 
 will inform the user of the drive. The tape drives are MTA0:, MTA1 : , 
 etc. The drive used should be remembered, since it will be needed later 
 in using EDITOR to read the tape. When discussion with the operator has 
 been completed, the BREAK command should be used to break the link. 
 
 In order to secure use of the tape, two more TENEX commands 
 must be used. The ASSIGN <tape drive> (e.g., ASSIGN MTA0:) command 
 assigns the tape drive to the user's job, preventing other users from 
 accessing it. The MOUNT <tape drive> (e.g., MOUNT MTA0:) positions the 
 tape for reading. 
 
 To evoke EDITOR, the user simply types <RAY> EDITOR followed 
 by one or two carriage returns. The two-carriage-returns form provides 
 a shorter introductory herdld and is recommended for regular users. The 
 user is then at the top command level of EDITOR and may enter commands 
 as detailed in the following section. A list of legitimate commands is 
 available by entering question mark ("?") followed by a carriage return 
 (CR). At the top level, any command may be abbreviated—as short as 
 required to distinguish it from all other commands—and followed by a 
 carriage return or ESC. If followed by ESC, the remainder of the abbre- 
 viated command will be completed by EDITOR. In certain of the analysis 
 modules within EDITOR, a somewhat different command structure is used in 
 which the user is prompted for short inputs by a series of questions. 
 To return to the top level of EDITOR from any lower-level module, an 
 
41 
 
 exclamation point ("!" and CR) is used instead of a command. (One is 
 at the top level of EDITOR when one is prompted by a single exclama- 
 tion point. ) 
 
 Most of the EDITOR modules require TENEX disk files for input 
 and/or output of data. The simplest form of a file designator is FILE- 
 NAME. EXTENSION where FILENAME and EXTENSION are identifiers and the 
 period is part of the file designator. Note that the period (".") is 
 always required between FILENAME and EXTENSION. 
 
 In entering commands through a terminal, errors will occa- 
 sionally be made. To delete a single character, CONTROL-A is used and 
 to delete an entire line, CONTROL-X is used. Similar conventions are 
 also available at the level of TENEX commands. For example, TENEX pro- 
 vides CONTROL-C for terminating any program and returning to TENEX 
 EXEC. (In some cases, two CONTROL-C characters must be entered.) It 
 should be noted that, wherever one might be in EDITOR, the CONTROL-C 
 completely exits EDITOR. However, this is the only way to terminate a 
 long program, such as a large cluster or classify, which was somehow 
 initiated erroneously. 
 
 Finally, when the EDITOR session is complete, EDITOR is exited 
 using the QUIT EDITOR command. If a tape was used, the user should 
 relinquish control of the tape drive. First the DEASSIGN <tape drive> 
 (e.g., DEASSIGN MTA0:) command is used to release exclusive control of 
 the tape drive. Next the UNLOAD <tape drive> (e.g., UNLOAD MTA0:) 
 command is used to rewind the tape and position it for dismounting, hence 
 signaling the operator that he may remove the tape and put it away. 
 
42 
 
 Thus, it is not necessary to LINK to the operator again at the end of an 
 EDITOR session. 
 
 Before logging out, the user should check to see if there are 
 any files on TENEX disk that will be no longer needed. A listing of all 
 files in the user's directory may be obtained via the TENEX command 
 DIRECTORY. It is prudent to delete files since all disk space used costs 
 money, and since once a directory is filled to its allocation (the 
 DSKSTAT command will indicate the space remaining) no new files may be 
 created. 
 
 To delete files, the command DELETE <file name> is used. 
 Where it is desirable to delete all files of a common FILENAME (or prefix 
 within the FILENAME. EXTENSION convention), this may be done conveniently 
 by using the DELETE FILENAME.* form of the DELETE command. This option 
 for TENEX file deletion should be noted by EDITOR users since, if care 
 is given to naming files within EDITOR sessions, file directory mainte- 
 nance can be greatly simplified. In any case, if obsolete TENEX files 
 are being deleted in order to regain disk space needed for immediate 
 use, the TENEX EXPUNGE command should be issued by the user following 
 all DELETES. 
 
 When the TENEX session is complete, the LOGOUT command is used 
 to exit the system before closing the connection. At this point, the 
 TENEX system will report to the user the total amount of processing (CPU) 
 time and terminal-connect time used during the complete session. 
 
43 
 
 3.2.9 Summary of EDITOR Commands 
 
 In this section, a brief description of the various EDITOR 
 commands and the modules they evoke will be given. The use of most of 
 these commands is illustrated in the example EDITOR sessions given in 
 the Appendix. Several commands have already been described partially 
 in Section 3.2.7 above. The commands are discussed here in alphabetical 
 order to assist subsequent referencing. 
 
 3.2.9.1 CLASSIFY 
 
 The CLASSIFY command is used to perform statistical classifica- 
 tion of an image data file using a TENEX routine. The number of channels 
 will be requested. At present, only 4- or 8-channel data may be handled 
 so "4" or "8" should be entered. The user will be prompted for other 
 required inputs. For large-scale classifications, the ILLIAC IV may be 
 used via the EXECUTE ANALYSIS ON ILLIAC IV command described below. 
 
 3.2.9.2 CLIP A WINDO W 
 
 The CLIP A WINDOW command is used to create subwindows of 
 windows in an image file already on TENEX disk. A question mark ("?") 
 upon entering CLIP lists the commands. CLIP may be used for editing 4- 
 or 8-channel raw or categorized data files. The subwindows to be created 
 may be specified in terms of new lines and columns (contained within the 
 input window) and/or additional image sampling. Commands within the 
 CLIP A WINDOW module allow selection of the image window to be edited, 
 and the new coordinates and sampling increments to be used in creating 
 
44 
 
 the new subwindows. After all window editing parameters have been 
 specified, a command is available to write the output file and return 
 to the top level of EDITOR. 
 
 3.2.9.3 CLUSTER 
 
 The CLUSTER command is used to perform cluster analysis on 
 image data using TENEX processing. The number of channels will be 
 requested. At present, only 4- or 8-channel data may be handled so 
 "4" or "8" must be entered. The user will be prompted for other required 
 inputs. For large cluster analyses, the ILLIAC IV may be used via the 
 EXECUTE ANALYSIS ON ILLIAC IV command (again, described below). 
 
 3.2.9.4 ENTER STATISTICS FILE EDITOR 
 
 The ENTER STATISTICS FILE EDITOR command enters the statistics 
 editor module. Subcommands are of the form one letter followed by a 
 carriage return. A list of commands may be displayed by typing question 
 mark ("?") followed by carriage return. Within this module, the summary 
 statistics for sets of spectral categories generated by cluster analyses 
 may be displayed. Also, the statistics of individual categories of 
 various files may be grouped in new combinations to form new statistics 
 files. In addition, the statistics of several categories may be pooled 
 together to create a composite spectral "signature." Finally, where a 
 classification is to be done using the ILLIAC IV, the variance-covariance 
 matrices of the appropriate statistics file may be inverted. 
 
45 
 
 3.2.9.5 EXECUTE ANALYSIS ON ILLIAC IV 
 
 The EXECUTE ANALYSIS ON ILLIAC IV command enables automatic 
 submission and retrieval of ILLIAC IV batch jobs from a remote TENEX 
 site. The ARPA Network is used to submit the jobs and retrieve the 
 results from the ILLIAC IV at NASA-Ames. To use the ILLIAC IV, the user 
 must have a valid user code and password at I4-TENEX (the TENEX front- 
 end subsystem of the ILLIAC IV complex). EDITOR will ask for these be- 
 fore allowing any commands to be entered. At present, both cluster 
 analysis and statistical classification procedures may be performed on 
 the ILLIAC IV remotely through EDITOR. 
 
 3.2.9.6 IDENTIFY A WINDOW 
 
 The IDENTIFY A WINDOW command displays information extracted 
 from the header of an image file. This information includes the image 
 file type, the line and column coordinates of all windows in the file, 
 the line and column sampling increments used in creating all windows in 
 the file, the number of channels of the image data, the cluster analysis 
 parameters (DELTA or separability threshold and the NUMBER OF CLUSTERS), 
 and certain descriptive information supplied by NASA for the ERTS image 
 from which the data have been extracted. 
 
 3.2.9.7 INDICATE A PRIORI PROBABILITIES 
 
 The INDICATE A PRIORI PROBABILITIES command allows creation of 
 a file of a priori probabilities of various categories for classification. 
 
46 
 
 3.2.9.8 MODIFY WINDOW HEADER 
 
 The MODIFY WINDOW HEADER command allows insertion of neces- 
 sary parameters (NUMBER OF CLUSTERS and DELTA or separability threshold) 
 to be used for a particular cluster analysis. Also there the alpha- 
 numeric information recorded in the header of an image file may be modi- 
 fied. Once all desired changes are made, the CLOSE FILE command should 
 be issued to return to the top level of EDITOR. 
 
 3.2.9.9 PRINT A WINDOW 
 
 The PRINT A WINDOW command allows display of image windows. 
 The display may be directed either to the user's terminal or to a back-up 
 disk file for later printing using a conventional line printer. Raw 
 data windows are displayed with an overprinting routine allowing up to 
 eleven (11) levels of gray-scale. A histogram of the distribution of 
 spectral intensities for any spectral channel may also be displayed. For 
 categorized windows, the character a-ssigned to each class may be selected 
 manually as any printing character to achieve more meaningful displays. 
 This feature facilitates the re-grouping of spectral classes by nominal 
 categories for tabulation and display purposes. 
 
 NOTE: When PRINT A WINDOW output is directed to a terminal, 
 the line width used for displays of windows is the line width assigned 
 to that terminal by TENEX. The default value is 72 characters, but this 
 may be changed by the TENEX WIDTH <number> command, where <number> is 
 the width desired. Since WIDTH is a TENEX command it must be performed 
 before entering EDITOR. 
 
47 
 
 3.2.9.10 QUIT EDITOR 
 
 The QUIT EDITOR command provides for a proper exit from the 
 top EDITOR command level back to TENEX EXEC. The QUIT command should 
 always be used whenever possible (instead of CONTROL-C) to exit EDITOR. 
 
 3.2.9.11 RETRIEVE WINDOW FROM TAPE 
 
 The RETRIEVE WINDOW FROM TAPE command is used to copy windows 
 from a specific tape (that has already been mounted) into a disk file. 
 Where desired, line and column sampling increments should be supplied 
 via the SAMPLE command. Then the coordinates of the windows to be 
 retrieved are entered using the INSERT COORDINATES command. Additional 
 facilities are available to allow the user to delete window coordinates 
 entered in error. 
 
 Prior to any EDITOR analysis of a specific ERTS image, specific 
 geographic areas of interest are known to the user, typically, only in 
 terms of latitude and longitude boundaries. On the other hand, all image 
 windows to be retrieved from tape by EDITOR must be specified in terms 
 of image line and column coordinates. Therefore as a convenience to the 
 user, there is included within the RETRIEVE module a subprocedure for 
 transformation of window corners specified in terms of latitude and 
 longitude coordinates into "nearest-approximation" image pixels identi- 
 fied by image line and column numbers. (The inverse transformation is 
 also possible; that is, for image pixels identified by line and column 
 coordinates, approximate latitude and longitude coordinates may be com- 
 puted. ) 
 
48 
 
 4. PATTERN RECOGNITION: A BASIS FOR 
 REMOTE SENSING DATA ANALYSIS 
 
 4.1 Classification Procedure 
 
 Pattern recognition plays a central role in numerically 
 oriented remote sensory systems. It provides an automatic procedure 
 for deciding to which class any given ground resolution elements 
 should be assigned. 
 
 Figure 6 shows a model of a general pattern recognition 
 system. In the LARS [7] context the receptor or sensor is usually a 
 multispectral scanner. For each ground resolution element the sensor 
 produces n_ numbers or measurements corresponding to the n channels of 
 the scanner. It is convenient to think of the n measurements as defin- 
 ing a point in n-dimensional Eucl idean Space which is referred to as 
 the Measurement Space . Any particular measurement can be represented by 
 the vector: 
 
 ~ x l 
 
 X = 
 
 The feature extractor transforms the n-dimensional measurement vector 
 
 i 
 into an n -dimensional feature vector. The decision maker in Figure 7 
 
 performs calculations on the feature vectors presented to it and, based 
 
49 
 
 sensor 
 
 x 2 
 
 xn 
 
 decision maker 
 
 result 
 
 Figure 6 A Simplified Model of a Pattern Recognition Syst 
 
 em 
 
class 1 
 
 50 
 
 class 2 
 
 decision 
 surface 
 
 ^— x- 
 
 Figure 7 Decision Regions and Surfaces 
 
 Result 
 
 Figure 8 A Pattern Classifier Defined in Terms of Discriminant Functi 
 
 on 
 

 51 
 upon a decision rule, assigns the "unknown" data point to a particular 
 class . 
 
 4.1.1 Discriminant Functions: Quantifying the Decision Procedure 
 
 Patterns arising in remote sending problems exhibit some ran- 
 domness due to the randomness of nature. Vectors from the same class 
 tend to form a "cloud of points" as shown in Figure 7. The job of the 
 pattern classifier is to divide the feature space into decision regions, 
 each region corresponding to a specific class. Any data point falling 
 in a particular region is assigned to the class associated with that 
 region. The surface separating the decision regions are known as deci - 
 sion surface . Designing a pattern recognizer really boils down to form- 
 ing a procedure for determining the decision surface so as to optimize 
 some performance criterion, such as maximizing the frequency of correct 
 classification. 
 
 These concepts can be put on a quantitative basis by intro- 
 ducing Discriminant Functions . Assume there are m pattern classes. Let: 
 
 g 1 (x), g 2 (x),..., g m (x) 
 
 be scaler single valued functions of x such that 
 
 g i (x) > g, (x) for all x 
 
 in the corresponding region corresponding to the i class (J f i). If 
 the discriminant functions are continuous across the decision boundaries, 
 the decision surface are given by equation of the form. 
 
 9^ (x) - g. (x) = 
 
52 
 
 A pattern classifier can then be represented by the block diagram of 
 Figure 8. By taking this approach the pattern classifier design 
 problem is reduced to the problem of how to reflect the discriminant 
 functions in an optimal fashion. 
 
 4.1.2 "Training" the Classifier 
 
 In some cases it is possible to select discriminant functions 
 on the basis of theoretical considerations or experience. More com- 
 monly the discriminant functions are based upon a set of training pat - 
 terns . Training patterns which are typical of those to be classified 
 are "shown to the classifier together with the identity of each pattern 
 and based on this information the classifier establishes its discriminant 
 functions. 
 
 g i (x), i = 1,2, .... n. 
 
 Example : Consider a two dimensional, two class problem in which the 
 discriminant functions are assumed to have the form 
 
 g 1 (x) = a n x ] + 9 ]? x 2 + b ] 
 
 g 2 (x) = a 21 x 1 + g 22 x 2 + b 2 
 
 Then g, (x) - g 2 (x) = is the equation of a straight line dividing x, , 
 x« plane. Given a set of training patterns, how should the constants 
 a,-., a,2» b-i > etc. be chosen? It can be proven that if the training pat- 
 tern are indeed separable by a straight ine, then the following procedure 
 will converge. [8] 
 
53 
 
 Initially select a's and b's arbitrarily. For example, let: 
 a ll = a 12 = b l = ] 
 
 Oiry-t ~ ®00 ~ 9 _ ~ 
 
 Then take the first training pattern (say it is from w, , i.e., from 
 class 1) and calculate g, (x) and g 2 (x). If g, (x) > g 2 (x) the deci- 
 sion is correct; go on to the next training sample. If g-, (x) < g 2 (x) 
 a wrong decision would be made. In this case, alter the coefficients so 
 as to increase the discriminant function associated with the correct 
 class and decrease the discriminant function associated with the incor- 
 rect class. If x is from to-, but co 2 was decided, let 
 
 i i 
 
 a -in - aii cxXi ^9i — ^9i ™ o'X i 
 
 i i 
 
 3i2 = 3-j n ■ CXX o ^99 = ^99 ~ CXX^ 
 
 i i 
 
 b-, = b, + a bp = b - a 
 
 where a is a convenient positive constant. If x is from w^ but w-, was 
 decided, change the signs in the above equation so as to increase g 2 and 
 decrease g-,. Then go to the next training pattern. Cycle through the 
 training patterns until all are correctly classified. 
 
 4.1.3 The Statistical Approach 
 
 Remote sensing is typical of many practical applications of 
 pattern recognition for which statistical methods are appropriate in the 
 following respects: 
 
 1. The data exhibit many incidental variations (noise) which 
 tend to obscure differences between the pattern classes. 
 
54 
 
 2. The pattern classes of interest may actually overlap in 
 the measurement space (may not always be discriminable), suggesting the 
 use of an approach which leads to decisions which are "most likely" 
 correct. 
 
 Statistical pattern recognition techniques often make use of 
 the probability density functions associated with the pattern classes 
 (including the approach to be described here). However, the density 
 functions are usually unknown and must be estimated from a set of train- 
 ing patterns. In some cases, the form of the density functions is 
 assumed and only certain parameters associated with the functions are 
 estimated. Such methods are called " parametri c. " Methods for which 
 not even the form of the density functions is assumed are called " non - 
 parametric . " The parametric case requires more a priori knowledge or 
 some basic assumptions regarding the nature of the patterns. The non- 
 parametric case requires less initial knowledge and fewer assumptions 
 but is generally more difficult to implement. Let there be m classes 
 characterized by the conditional probability density functions. 
 
 p (X^.) i = 1,2,. . ., m (1) 
 
 The function p (XJoo.) gives the probability of occurrence of pattern X, 
 given that X is in fact from class i. 
 
 An important assumption in the LARSYS (Purdue University) 
 algorithms is that the p (X|gj. ) are each multivanate Gaussian (or normal) 
 distributions. This is a parametric assumption which leads to a form of 
 classifier which is relatively simple to implement. Under this assumption, 
 
55 
 
 a mean vector and covariance matrix are sufficient to characterize the 
 probability distribution of any pattern class. 
 
 Returning to the problem of how to specify the discriminant 
 functions, an approach based on statistical decision theory is taken. 
 A set of less functions is defined: 
 
 X (i|j) i,j = 1,2,. . ., m (2) 
 
 where A (i|j) is the loss (or cost) incurred when a pattern is classified 
 into class j_ when it is actually from class j_. If the pattern classifier 
 is designed so as to minimize the average (expected) loss, then the 
 classifier is said to be Bayes optimal . For a given pattern X, the 
 expected loss resulting from the decision Xeo). is given by 
 
 m 
 
 L x (D = I A (i|j) pU |x) (3) 
 
 j = l J 
 
 where p(w.|x) is the probability that a pattern X is from class j. 
 Applying Bayes' rule, i.e., 
 
 p (x,Wj) = p (x|w.) P (wj) = p(o).jx) p (x) (4) 
 
 the expected loss can be: 
 
 L x (D = I A Ci id) P(x|a).) p(u).)/p(x) (5) 
 
 j = l J J 
 
 where p (go.) is the priori probability of to., 
 j j 
 
 Note the minimizing L x (i) with respect to i_ is the same as 
 maximizing - L x (i). Thus a suitable set of discriminant functions is 
 
 9,- (x) = -L x (1) i = 1,2,. . .,m (6) 
 
56 
 
 A simple (a reasonable) loss function is 
 X (i|j) = i = j 
 
 A (i|j) = 1 i t J 
 
 m 
 Then g.(x) = - \ p(x|u>.) p(ui.)/p(x) (7) 
 
 J7i 
 Note that from any set of discriminant functions, another set of dis- 
 criminant functions can be formed by taking the same monotonic functions 
 of each of the original discriminant functions. For example, if: 
 
 g i (x), 1 = 1,2, • . -,m 
 
 is a set of discriminant functions, then so are the sets 
 
 9,- U) = 9,- ( x ) + constant i = 1,2,. . . ,m 
 
 and 
 
 9-| (x) = log [g^x)] i = 1,2,. . . ,m 
 
 Examing (7) note that p (x) is not a function of i so it is just as well 
 
 to maximize. 
 
 m 
 
 9,- (x) = - I P (xLo) p (o>.) = - [p(x) - p(x to. ) p (u>.)] 
 i j =1 I J J « l i 
 
 But this is maximum if 
 
 q (x) = p (x|oo i ) p (w.) 
 
 is maximum. Thus, the decision rule is: 
 decide x ecu. if and only if 
 
 p (x|uj i ) p (x^ > p (x|oa.) p (go.) for all j 
 
57 
 
 that is deciding x eu. if g. (x) = g. (x) and i > j. This is commonly 
 referred to as the Maximum Likelihood Decision Rule . 
 
 Consider the maximum likelihood discriminant function as it 
 applies to remote sensing. The p (w. ) represent the priori probability 
 
 t h 
 
 of the i class. This can often be estimated. Taking agricultural 
 crop types as an example, the p (to. ) may be estimated from previous year 
 yields, seed sales records or statistical reporting service information. 
 The densities p (xju).), on the other hand, generally have to be esti- 
 mated from training samples. 
 
 The assumptions upon which the classification algorithms are 
 based is that p (x|w. ) is a multivariate Gaussian probability density 
 function. Examining the maximum likelihood decision criterion in 
 one-dimensional Gaussian case will serve both as a review of Gaussian 
 density functions and as a means of illustrating the principles of pat- 
 
 tern classification. In this case (e.g., one spectral channel) 
 
 ,2 
 
 exp 
 
 P (x|oa i ) = __ 
 
 1/2 
 (2tt) a f 
 
 1/2 (x - y i V 
 
 2 2 
 
 when u. = E [x] and a. = E [(x-u^) ] are the mean and variance for 
 
 2 
 class i. In practice u. and a. are unknown and must be estimated from 
 
 training samples. From statistical theory, we have 
 
 n, 
 
 /\ -I L 
 
 rt -i — I <J 
 
 "t 
 
 o, = Si = 
 
 1 n t 
 
 V 1 
 
 I (x.- m .)' 
 
58 
 
 (n. = number of training patterns in class 1) are unbiased estimators of 
 the mean and variance. Thus, the estimated density function is: 
 
 2 
 
 P (xjoj^ = 1 
 
 1/2 
 (2tt) S i 
 
 exp 
 
 -1/2 < X - m i>' 
 
 1 _J 
 
 Following the decision theory approach the discriminant function is 
 
 9i (x) 
 
 p(u).j) 
 
 exp 
 
 (2ir) 
 
 1/2 
 
 S i 
 
 - 1/2 (x - m i }< 
 
 S 2 
 b i 
 
 and since a monotonic function of a discriminant function may also be 
 used as a discriminant function, we shall take the logarithm of the 
 previous function to obtain: 
 
 ( \ 2 
 
 g. (x) = log p (u>.) - 1/2 log 2 tt - logs.- 1/2 u " V 
 
 g.' (x) = log p (a).) - log S i - 1/2 ^" m 1 ) ' 
 
 Thus the decision rule becomes: 
 
 Decide Xcw. if and only if 
 
 2 2 
 
 (x-m.) (x-m.) 
 
 log p(o) i ) -log S-i/2 — ^ — > log p(u.) - log S-. 1/2 —^ — 
 
 s i S. 
 
 In the two dimensional case 
 
X = 
 
 59 
 
 and 
 
 exp 
 
 P (x|a3 i ) = 1 
 
 (x l - Vj} y 
 a ill 
 
 2 tt ( 0ill a i22 - a il2 ) 
 
 2a il2 (x l " y il } (x 2 " y i2 ) + (x 2 " Hl ){ 
 
 (a ill • a i22 } 
 
 -1/2 
 
 1 - 
 
 a il2 2 
 
 a ill a i22 
 
 122 
 
 where 
 
 y il = E l- x l | 0J i-'' y i2 = E ^ x 2 | a) i-' 
 Q ijk = E[(x i " »*1J> (x k - y ik>h ] 
 
 by defining a mean vector and covariance matrix 
 
 j,k = 1,2 
 
 i = 1,2, . . .R 
 
 Ui = 
 
 II- 
 
 M il 
 y i2 
 
 °ill a il2 
 
 a i21 a i22 
 
 The density can be rewritten in the simple form: 
 
 P (xjo)^ = 1_ 
 
 2ir|I.| 
 
 1/2 
 
 exp 
 
 1/2 (x 
 
 U,) T I- 1 (x -Uj) 
 
60 
 
 Wh 
 
 ere |J.| is the determinant of J. and (x-U.) is the transpose of 
 
 (x - U.). Thus this matrix formulation holds for n_ dimensions as well 
 as for two dimensions. 
 
 4.2 Clustering Procedure 
 
 Clustering is a data analysis technique by which one attempts 
 to determine the "natural" relationships in a set of observations or 
 data points. It is sometimes referred to as unsupervised classifica- 
 tion because the end product is generally a classification of each 
 observation into a "class" which has been established by the analysis 
 procedure, based on the data, rather than by the person interested in 
 the analysis. Clustering has been applied as a means of data compres- 
 sion and for the purpose of determining differentiating characteristics 
 in complex data sets. An increasingly important application is unsuper- 
 vised classification, in which the clustering algorithm determines the 
 classes based on the clustering tendencies in the data. Essentially, the 
 definition of a clustering algorithm depends on the specification of two 
 distance measures: a measure of distance between data points or indi- 
 vidual observations, and a measure of a distance between groups of obser- 
 vations. Figure 9 is a block diagram for a typical clustering algorithm. 
 The point-to-point distance measure is used in the step labeled "assign 
 each vector to nearest cluster center." The distance between groups of 
 points (clusters, in this case) is calculated in the step "compute 
 separability information." Euclidean distance, the most familiar point- 
 to-point distance measure, is defined for two-n-dimensional points or 
 
61 
 
 Initialize 
 
 cluster 
 
 centers 
 
 Assign each vector 
 to nearest 
 cluster center 
 
 Calculate means of 
 
 new "clusters" 
 
 (new cluster centers) 
 
 Yes 
 
 Compute 
 
 separability 
 
 information 
 
 No 
 
 Figure 9 Clustering Algorithm 
 
M ? 
 
 Euclidean distance : D = [£ (x.-y.) ] 
 
 i=l 1 n 
 
 62 
 1/2 
 
 Several alternatives are available as candidate measures of distance 
 between clusters. One possibility is the divergence or transformed 
 divergence used for feature selection. In LARSYS, a measure called 
 "Swain-Fu distance" has been implemented, which compares the separation 
 of cluster centers to the dispersion of the data in the clusters. The 
 dispersion of the data in a cluster is measured in terms of the 
 "ellipsoid of concentration" associated with the cluster. 
 
 Ellipsoid of concentration : Let the random vector X have a 
 distribution with mean vector U and covariance matrix I = [<?■,•■;]• If 
 Z is another random vector uniformly distributed over the volume of the 
 ellipsoid given by 
 
 Q (Z) = I I \hA (Z.-U.) (Z.-U.) = n + 2 
 1-1 J-l m J 3 
 
 where n is the number of components in 1, and |T. . | is the cofactor of 
 ff. ., then Z also has zero mean and covariance matrix Z. The ellipsoid 
 Q is called the ellipsoid of concentration of the distribution of X. 
 
 Consider two clusters and their respective ellipsoids of con- 
 centration as shown in Figure 10. D-.^ is the distance between the 
 cluster centers. D-. is the distance from the center of cluster 1 to the 
 surface of its ellipsoid of concentration along the line connecting the 
 cluster centers similarly to D«. In terms of these distances, D, , Dp, 
 D-|2» the Swain in distance is proven by 
 
 A= D 12 
 
 D l +D 2 
 
63 
 
 In terms of the cluster centered (cluster means) and the covariance 
 matrices associated with the clusters, the Swain-Fu distance can be 
 expressed as 
 
 A = TzTcZ 
 
 /c7 + /c\. 
 
 where C k = tr {£~ ] (Uj-u 2 ) (u r u 2 ) T } 
 
 tr {A} = trace of matrix A 
 
 J, = covariance matrix for cluster K 
 
 U. = mean vector for cluster K. 
 
 An illustration of how the clustering algorithm works refers 
 to Figure 10. The first step is to select initial cluster centers. The 
 analyst must specify how many clusters are to be isolated; the algorithm 
 determines (arbitrarily) where the initial centers are to be located (the 
 final results are relatively insensitive to the initial selection. Each 
 data point is then labeled as "belonging" to the nearest cluster center 
 (using Euclidean distance), effectively creating a cluster of data points 
 associated with each center. The boundaries between clusters are formed 
 by the lines (planes in n-dimensional space) which are the perpendicular 
 bisectors of the lines connecting the centers. Next, new cluster cen- 
 ters are calculated. The new center for each cluster is the mean vector 
 of all points just assigned to that cluster. A check to see whether the 
 algorithm has achieved the final results, which is the care when the new 
 
64 
 
 -*- x- 
 
 Figure 10 Separability of Clusters 
 
65 
 
 cluster centers are identical with the previous centers. When no 
 further change is detected, the pair wise distances (Swain-Fu distance) 
 between the resulting clusters are computed and all results are printed 
 for evaluation by the analyst. 
 
 
66 
 
 LIST OF REFERENCES 
 
 1. NASA. Earth Resources Technology Satellite—Data Users Handbook. 
 
 "Goddard Space, NASA Flight Center." 
 
 2. Paul E. Anuta. Geometric Correction of ERTS-1 Digital Multispectral 
 
 Scanner Data. LARS, Purdue University 1973. 
 
 3. W. E. Donovan, M. Ozga, R. M. Ray. CAC Technical Memorandum No. 52. 
 
 4. P. E. Anuta. "Spatial Registration of Multispectral and Multi- 
 
 tempora- Digital Imagery Using Fast Fourier Transform Tech- 
 niques. IEEE Transaction on Geoscience Electronics, Vol. 
 GE-8, pp. 353-368, Oct. 1970. 
 
 5. P. E. Anuta, M. Bauer. An Analysis of Temporal Data for Crop 
 
 Species Classification and Urban Change Detection. LARS 
 Information Note 110873. Purdue University, West Lafayette, 
 Indiana, 1973. 
 
 6. W. Donovan. Oblique Transformation, CAC Technical Memorandum No. 
 
 38. 
 
 7. Philip H. Swain. Pattern Recognition. LARS Information Note 11572. 
 
 Purdue University, Indiana. 
 
 8. Nilsson, N. J. Learning Machines. McGraw-Hill, 1965. 
 
 9. R. Ray, J. D. Thomas, W. E. Donovan, P. H. Swain. "Implementation 
 
 of ILLIAC IV Algorithms for Multispectral Image Interpretation. 
 CAC Document No. 112, University of Illinois 1974. 
 
 10. P. H. Swain. "Implementation and Evaluation of ILLIAC IV Algorithms 
 
 for Multispectral Image Processing." LARS Technical Memorandum 
 T-16 111273, Purdue University, Indiana, July 1974. 
 
 11. K. S. Fu, D. A. Landgrebe, and T. L. Phillips. "Information Proces- 
 
 sing of Remotely Sensed Agricultural Data," Proceedings IEEE, 
 Vol. 57, No. 4, pp. 639-653, April 1969. 
 
 12. "Remote Multispectral Sensing in Agriculture," Research Bulletin 
 
 873, Laboratory for Agricultural Remote Sensing, Purdue Uni- 
 versity, West Lafayette, Indiana, December 1970. 
 
 13. Martin Ozga, Walter E. Donovan, Robert M. Ray, John D. Thomas, and 
 
 Marvin L. Graham. "Data File Formats for Processing of Multi- 
 spectral Image Data," CAC Technical Memorandum No. 19 (revised), 
 Center for Advanced Computation, University of Illinois at 
 Urbana-Champaign, Urbana, Illinois, July 1975. 
 
67 
 
 14. G. H. Ball and D. J. Hall. "ISODATA, A Novel Technique for Data 
 
 Analysis and Pattern Classification," Technical Report, 
 Stanford Research Institute, Menlo Park, California, May 1965. 
 
 15. John Thomas. "An ILLIAC IV Algorithm for Cluster Analysis of 
 
 ERTS-1 Data," CAC Technical Memorandum No. 17, Center for 
 Advanced Computation, University of Illinois at Urbana- 
 Champaign, Urbana, Illinois, May 1974. 
 
 16. John Thomas. "An ILLIAC IV Algorithm for Statistical Classifica- 
 
 tion of ERTS-1 Imagery," CAC Technical Memorandum No. 18, 
 Center for Advanced Computation, University of Illinois at 
 Urbana-Champaign, Urbana, Illinois, May 1974. 
 
 17. Martin Ozga. "An ILLIAC IV Algorithm for Cluster Analysis of 
 
 8-Channel Multi spectral Image Data," CAC Technical Memorandum 
 No. 53, Center for Advanced Computation, University of 
 Illinois at Urbana-Champaign, Urbana, Illinois, May 1975. 
 
 18. Martin Ozga and John Thomas, "An ILLIAC IV Algorithm for Statistical 
 
 Classification of 8-Channel Multispectral Image Data," CAC 
 Technical Memorandum No. 54, Center for Advanced Computation, 
 University of Illinois at Urbana-Champaign, Urbana, Illinois, 
 June 1975. 
 
 19. Theodore H. Myer and John R. Barnaby, TENEX Executive Language - 
 
 Manual for Use rs, Bolt, Beranek and Newman, Inc., Cambridge, 
 Massachusetts, April 1973. 
 
BIBLIOGRAPHIC DATA 
 SHEET 
 
 1. Report No 
 
 UI 
 
 -R-76-801 
 
 3. Recipient's Accession No. 
 
 4. Title and Subtitle 
 
 PROCEDURES FOR ANALYZING DATA OBTAINED FROM 
 "LANDSAT" MULTISPECTRAL IMAGERY 
 
 5. Report Date 
 
 May 1976 
 
 7. Author(s) 
 
 MAHMOUD HAIDAR 
 
 8. Performing Organization Rept. 
 
 No - UIUCDCS-R-76-801 
 
 9. Performing Organization Name and Address 
 
 Department of Computer Science 
 
 University of Illinois at Urbana-Champaign 
 
 Urbana, Illinois 61801 
 
 10. Project/Task/Work Unit No. 
 
 11. Contract /Grant No. 
 
 USNASANGR-14-005-202 
 
 12. Sponsoring Organization Name and Address 
 
 13. Type of Report & Period 
 Covered 
 
 Master's Thesis 
 
 14. 
 
 15. Supplementary Notes 
 
 16. Abstracts 
 
 Procedures for analyzing data obtained from the "LANDSAT" Satellite Multispectral 
 Imagery: 
 
 1) Image Correction and rectification procedure. 
 
 2) Clustering procedure. 
 
 3) Classification procedure. 
 
 17. Key Words and Document Analysis. 17o. Descriptors 
 
 LANDSAT, Multispectral, pixel, EDITOR 
 
 17b. Identifiers /Open-Ended Terms 
 
 17c. COSATI Field/Group 
 
 18. Availability Statement 
 
 Release Unlimited 
 
 FORM NTIS-35 (10-70) 
 
 19. Security Class (This 
 Report) 
 
 UNCLASSIFIED 
 
 20. Security Class (This 
 
 Page 
 UNCLASSIFIED 
 
 21. No. of Pages 
 
 73 
 
 22. Price 
 
 USCOMM-DC 40329-P71 
 
<fc 
 
 v^ 
 
JUN 2 1 REC'O