JPH0823886B2 - Unequally spaced interpolation method - Google Patents

Description

Description: BACKGROUND OF THE INVENTION 1. Field of the Invention The present invention relates to a method of interpolating data arranged at unequal intervals, and more particularly to a suitable unequal interval interpolation method for performing interpolation with high speed and high accuracy.

[Background of the Invention]

Before describing the related art, a satellite image geometric distortion correction process, which is the background of the present invention, will be described.

The image of the surface of the earth captured by artificial satellites is the change in the orbit and attitude of the satellite, the rotation of the earth, and imaging equipment (called sensors).
Due to internal factors such as, geometrical and strength distortion is included. When using satellite images, these distortions need to be corrected by digital image processing. In the following, the geometric distortion correction processing will be extended, and simply the distortion means the geometric distortion.

The principle of distortion correction processing called re-sampling, which has been performed in satellite ground stations such as the Landsat ground station, is as shown in FIG.

First, a point (u, v) on the uncorrected image 1 corresponding to an arbitrary image position (x, y) on the corrected image 2 is obtained. This map is determined from the satellite orbit and attitude data. Since the point (u, v) lies between the uncorrected image data, the true image intensity value is estimated by two-dimensional interpolation processing from some surrounding data, and at the point (x, y) on the corrected image Image intensity value.

Since this two-dimensional interpolation processing can be generally performed at high speed in place of two vertical one-dimensional interpolations, the following description is limited to one-dimensional interpolation.

When the distortion changes sufficiently slowly, the pixels on the uncorrected image may be considered to be locally arranged at equal intervals, and the interpolation is equal interval interpolation.

For evenly spaced interpolation, the Kyubikku convolution method is widely known. (For example, "Digital Im
age Reconstruction and Resampling for Geometric Ma
nipulation, KWSimon, Presented at the Symposium on
Machine Processing of Remotely Sensed Data, June3
-5, 1975 ". ) In the queuing combo convolution method, as shown in FIG. 2, the weighting coefficient ▲ W ^cc _i ▼ (), i = 1,4 for the interpolation position x is made into a table in advance, and the interpolated image intensity value is set. I (x) is close to 4
Using the intensity values of points _i i i = 1,4 The interpolation processing can be performed at high speed by obtaining Here, the representative point closest to x usually takes a value of 1/8 or 1/16.

However, the queuing convolution method has a drawback that it can be applied only when pixels are arranged at equal intervals.

By the way, the typical sensor mounted on the Landsat satellite is
ectral Scanner, MSS) and Semathematic Mapper (abbreviated as Thematic Mapper, TM), both of which are sensors called scanning type. An imaging method of the scanning sensor will be described with reference to FIG. Vibration scanning mirror 3 on the satellite body 301
02, optical system 303 and detector 304 are mounted. The satellite captures an image of the ground surface directly below as it advances, but the scanning sensor repeatedly captures an image of the ground surface over a scanning range 306 by repeatedly scanning the vibration scanning mirror 302 in a direction perpendicular to the orbit direction 305. In one scan, the detector 3 is used for each pixel 307.
You can get data for 04 lines. The number of detectors
6 for MSS and 16 for TM.

Therefore, when imaging the ground surface with a scanning sensor,
Pixels lined up in the scanning direction and in the scanning vertical direction within one scanning can be regarded as being at equal intervals, but are generally unequal intervals at the joints of scanning. In the case of MSS, since the scanning is in one direction and the image interval in the vertical scanning direction is as wide as about 79 m, it is not necessary to consider this unequal interval portion. However, in the case of TM, since it is a reciprocating scan and the pixel interval is about 30 m, the non-uniform intervals cannot be ignored. (Actual TM
According to the sensor specifications, unequal intervals of up to ± 2 pixels will occur. Therefore, the distortion correction of the TM image cannot be corrected only by the regular interval interpolation method such as the conventional queuing convolution method.

As a distortion correction method for an image including non-equidistant portions such as a TM image, a pseudo queuing combo solution method and a GE (General Electric) method have been conventionally known. Japanese Patent Application No. 57-168362 "Processing Method of Image Data Inducing Scan Error" for the pseudo queuing convolution, and "LANDSAT-DT" for the GE method.
hematic Mapper Image Resambling for Soan Geometry
Correction, EPBeyer, etal, Machine Processing of Re
See "Motely Sensed Data Symposium, 1981" and Japanese Patent Application No. 57-168362.

In the pseudo queuing convolution method, unequal intervals are regarded as equal intervals, and interpolation is performed by the above queuing convolution method. That is, as shown in FIG. 4, when the width of the unequal interval at the scanning seam is d and the interpolation position is x, the interpolated image intensity value I
(X) Becomes Here, the weight ▲ W ^cc _i ▼ is the same as that appearing in the equation (1). Also, ▲ ▼ is the representative point closest to x / d.

However, in this pseudo-cubic convolution, since the interpolating calculation is performed by treating the unequal intervals as if they were equidistant, there is a drawback that the accuracy of the interpolation decreases.

Next, the GE method will be described with reference to FIG. D represents the line interval, and GP represents the pixel position of the corrected image to be obtained.

In the GE method, the extension line (17E,
18E) The above data 53 is obtained by spline interpolation, and apparently unequal intervals are eliminated, and even interval interpolation (Kewbit convolution method, etc.) is performed. The spline interpolation will be described in detail in Examples,
The feature of the GE method is that the weighting coefficient table is one-dimensionalized,
The points 15, 16, 17E, and 18E are evenly spaced.
(That is, it corresponds to setting the interpolation position x to a fixed value among the two variables of the unequal interval width d and the interpolation position x.) However, in the GE method, two-stage interpolation processing (first expansion It takes a lot of processing time because the line is calculated and then interpolated at regular intervals). Also, since the order of the spline function is fixed, it is not possible to use the spline function according to the width of the irregular intervals. There was a drawback that the insertion accuracy was lowered.

[Object of the Invention]

SUMMARY OF THE INVENTION An object of the present invention is to improve the above-mentioned drawbacks of the prior art and to provide an unequal interval interpolation method for interpolating data arranged at unequal intervals at high speed and with high accuracy.

[Outline of Invention]

In order to achieve the above object, the present invention provides a weighting coefficient used for interpolation as a two-dimensional weighting coefficient table in the form of a two-dimensional table in which two variables, that is, the width of unequal intervals and the interpolation position are variables. Has the first feature. The details will be described below. In the unequal interval interpolation method that images the ground surface by repeating scanning in the direction perpendicular to the orbit direction of the artificial satellite and interpolates the imaged image data arranged in unequal intervals The weighting factors used when interpolating the data are stored as a two-dimensional table in which the width of the unequal intervals and the interpolation position are variables, and the correspondence of each point on the corrected image on the input image Seek points,
When the obtained corresponding points are in the seam area of the scan of the input image, the width of the unequal interval in the area of the seam is obtained, and the obtained corresponding points and the obtained corresponding points are obtained using the two-dimensional table. It is characterized in that the weighting factors corresponding to the unequal intervals are read out, and the image data arranged at unequal intervals are interpolated using the read out weighting factors. A second feature is that when the weighting factor is generated according to the spline interpolation method, a spline function whose order changes according to the width of unequal intervals is used. Furthermore, in order to simplify the address calculation and perform unequal interval interpolation at high speed, the third feature is that the two-dimensional weighting coefficient table is divided and edited to be a rectangular table.

Example of Invention

An embodiment of the present invention will be described below by taking the Semanates Kumaper (TM) image distortion correction system mounted on the Landsat 4 satellite as an example.

FIG. 6 shows a hardware configuration for implementing the processing method of the present invention.

The CPU 101 calculates a distortion correction coefficient necessary for distortion correction processing (resampling). Distortion correction processing is an image correction processing device.
105 (High-speed array processor-AP120B). Image data is input from a high density magnetic tape device (HDDT) 102 and stored in a temporary magnetic disk 104. The corrected image data is output to the magnetic tape device 103.

FIG. 7 is a schematic diagram of a processing flow for correcting the geometric distortion of the satellite image.

When the low-frequency attitude data 701, the high-frequency attitude data 702, the orbit data 703, the TM mirror scanning center time data 704 and the uncorrected image 705 are input, the attitude angle calculation 706, the orbit calculation 707, and the scanning mirror are calculated based on these data. Correction 708 is performed and scanning is controlled. Further, the deviation amount 709 and the geometric distortion correction coefficient 710 are calculated, and the distortion correction processing 711 according to the present invention is performed on the basis of these results to obtain the corrected image data 712.
Is output.

The distortion correction processing is performed in two stages, first in the scanning direction (horizontal direction) and then in the scanning vertical direction (vertical direction). An image that has been subjected to distortion correction in the scanning direction is called a hybrid image.

Since the pixel intervals in the scanning direction can be considered locally as equal intervals, the equal interval interpolation method such as the queuing convolution method is used as usual. The vertical direction one-dimensional interpolation processing for converting a hybrid image into a corrected image, which causes a problem of non-uniform intervals occurring at the seam of scanning, will be described in detail below.

First, as shown in FIG. 8, the corresponding point 811 of the point 821 on the corrected image 82 to the hybrid image 81 is calculated by the vertical distortion correction coefficient. If the corresponding point like the point 811 is in the area A shown in FIG. 8, that is, in the portion excluding the interval of one pixel at both ends in one scan, the interpolation processing is performed by using the ordinary 4-point queuing convolution method. .

When corresponding points, such as the point 812, are in the region B, unequal interval interpolation is performed. FIG. 9 shows a sectional view corresponding to AA 'in FIG.

First, using the vertical direction distortion correction coefficient,
From the interpolated positions x and d, the six weighting factors ▲ W ^SP _i ▼ corresponding to the position (,) in the two-dimensional weighting factor table.
(,) I = 1,6 is read out, and image intensity values I of 6 points in the vicinity are read.
Interpolated value I (x, d) from _i , i = 1,6 To ask. Here, is the representative point closest to x and d, respectively, and is 1/16, and takes a value in 1/8 steps.

By performing this process for all points of the corrected image, the distortion correction process in the vertical direction is completed.

 Hereinafter, a method of constructing the two-dimensional weight table will be shown.

Interpolation is performed using a spline function. For details of the spline function, refer to, for example, “Spline function and its application, Kozo Yoshida et al., Educational Publishing”. The spline function is a piecewise polynomial in which polynomials are connected so as to satisfy some continuous condition.

The width d of unequal intervals is -1 <d≤3 (from TM sensor specifications), and the interpolation position x is 0 <x <2 + d. (Refer to FIG. 9) First, when a linear spline function is used, the linear interpolation (the linear spline function means a straight line in each section, that is, a polygonal line) is linear as shown in FIG. It becomes an interpolation. Since the method of generating the weighting coefficient at this time is self-explanatory, it is omitted.

Next, the case of using the most commonly used cubic spline function will be described.

 First, consider the case of 0 <d.

Let the cubic spline functions in the sections R ₁ , R ₂ and R ₃ in FIG. 11 be y ₁ , y ₂ and y ₃ , respectively, y _i = a ₁ ⁽ⁱ⁾ x ³ + a ₂ ⁽ⁱ⁾ x ² + a ₃ ^{(i )} x + a ₄ ⁽ⁱ⁾ (4) (i = 1,2,3).

Twelve coefficients a _j ⁽ⁱ⁾ (i = 1,2,3, j = 1,2,3,4) are obtained under the following conditions.

(1) Passing the end points of each section (2) Continuous up to quadratic differential coefficient at the joint of sections (3) Designate the slope at x = 0,2 + d When the conditions (1) to (3) are written down with ▲ a ⁽ⁱ⁾ _j ▼ and summarized in matrix form, However, If we write D ^-1 = F, then a _j ⁽ⁱ⁾ can be expressed as the solution of simultaneous equations as follows.

Therefore, the interpolated value I (x, d) Then, by comparing both sides, it is possible to obtain ▲ W ^SP _j ▼ (x, d), j = 1, ..., 6.

When -1 <d <0 (see FIG. 12), W _j (x, d) is obtained by the same idea.

When d = 0 (see FIG. 13), I ₃ and I ₄ overlap, so
After changing the average value ▲ I ^′ ₃ ▼ = (I ₃ + I ₄ ) / 2, the weight W _j (x, 0) is similarly obtained.

As described above, the previously obtained ▲ W ^SP _j ▼ (x, d) is
By storing all in the form of a two-dimensional table,
The interpolation calculation can be performed at high speed.

Next, regarding the order of the spline function used for interpolation, the order can be changed according to the width d of the unequal intervals. According to the experiment shown in FIG. 14, the cubic spline is most accurate when d ≦ 1.6, and the linear interpolation is most accurate when the d> 1.6. Here, the horizontal axis represents the width d of unequal intervals, and the vertical axis represents the interpolation accuracy (RMS). ○
Is the nearest-neighbor interpolation method, × is the linear interpolation method, △ is the pseudo queuing combo reuse method, Is the cubic spline method, Indicates the value according to the GE method. In this experiment, the original image and the restored image obtained by performing the above-mentioned spline interpolation for the image in which the unequal intervals are simulated by the interpolation processing using the Kyubik's convolution method are performed. This is the difference.

Based on the result of this experiment, a weighting coefficient is generated using a spline function of an order that changes according to d.

As described above, according to the present embodiment, the weighting factors used for unequal interval interpolation are calculated and stored in advance in the form of a two-dimensional table, so that unequal interval interpolation can be executed at high speed, and unequal By changing the order of the spline function used for interpolation according to the width d of the interval, there is an effect that the interpolation calculation can be performed accurately.

Fig. 15 shows the invention method, GE method, pseudo-qubit method,
This is a comparison of interpolation accuracy and calculation amount for linear (first-order spline) methods. Here, the interpolation accuracy is an error (RMS) measured by a simulation experiment, and the calculation amount is indicated by the number of floating point calculations required for each interpolation.

Next, an embodiment in which the two-dimensional weighting coefficient table is converted into a rectangle will be described.

The width d of the unequal interval may be −2 <d ≦ 2 according to the TM sensor specification, and the interpolation position x is 0 ≦ x <2 + d. (No. 16
Therefore, the weight coefficient table is as shown in FIG. If this table, in which one box represents a set of six weighting factors, is stored in the memory as it is, it becomes as shown in FIG. At this time, it can be seen that the weighting coefficient for a certain (,) is stored in the {(D + 15) (D + 16) + X} · 6 ... (12) th position relative to the beginning. However Is.

Since the equation (12) is a quadratic equation with respect to D, the address calculation is complicated and the processing time is long. This is because the table in Fig. 16 is not rectangular.

FIG. 18 shows a table structure according to the method of the present invention. That is, in FIG. 16, the part of X ≧ 32 is divided and folded back,
It was edited into a rectangular table of 32 x width 33. The layout in memory is shown in FIG.

At this time, the weighting coefficient for a certain (,) is a relative position from the beginning. Becomes

Since the equation (13) is a linear equation with respect to X and D, the address can be calculated at high speed.

〔The invention's effect〕

According to the present invention, when interpolating data arranged at unequal intervals, a two-dimensional table in which the weighting coefficient used for interpolation has two variables of the unequal interval width d and the interpolation position x as variables. By storing as, it is possible to perform unequal interval interpolation at high speed,
Further, since the order of the spline function used for interpolation can be changed according to d, there is an effect that unequal interval interpolation can be performed with high accuracy. Further, since the two-dimensional weighting coefficient table for unequal interval interpolation is edited into a rectangle, there is an effect that interpolation processing can be executed at high speed.

[Brief description of drawings]

FIG. 1 is a diagram showing the principle of distortion correction, FIG. 2 is a diagram explaining an equal-interval interpolation method by the queuing convolution method, and FIG. 3 is a diagram showing an imaging method of a scanning type sensor.
FIG. 7 is a diagram for explaining the pseudo queuing convolution method, FIG. 5 is a diagram for explaining the GE system, FIG. 6 is a diagram showing a hardware configuration for performing data processing of Landsat 4 satellite, FIG. 7 Is a diagram showing a flow of data necessary for image data correction of Landsat 4 satellite, FIG. 8 is a diagram showing an outline of vertical direction distortion correction processing, and FIG. 9 is A- in FIG.
FIG. 10 is a sectional view taken along line A ′, FIG. 10 to FIG. 13 are diagrams showing an unequal interval interpolation method, FIG. 14 is a diagram showing results of an interpolation accuracy experiment, and FIG. 15 is a diagram showing a method of the present invention and a conventional method. Figure showing comparison results, Figure 16 ~
FIG. 19 is a diagram showing a configuration example in which the weighting coefficient table is rectangular. 1 ... Pixel position in uncorrected image, 2 ... Pixel position in corrected image, 301 ... Satellite body, 302 ... Vibration scanning mirror, 304 ...
... Detector, 101 ... CPU, 102 ... High-density magnetic tape device, 103 ... Magnetic tape device, 104 ... Magnetic disk device, 105 ... Image correction processing device.

 ─────────────────────────────────────────────────── ─── Continuation of the front page (72) Hideo Ota 5-2-1 Omika-cho, Hitachi-shi, Ibaraki, Hitachi, Ltd. Inside the Omika Plant, Hitachi, Ltd. (56) References: 57-84058

Claims (4)

[Claims] 1. An unequal interval interpolation method for imaging the ground surface by repeating scanning in the direction perpendicular to the orbital direction of an artificial satellite and interpolating the imaged image data arranged at unequal intervals. The weighting factor used when interpolating the image data arranged at equal intervals is defined by the variables of the nonuniform spacing and the interpolation position.
It is stored as a dimension table, the corresponding points on the input image of each point on the corrected image is obtained, and if the obtained corresponding points are within the area of the seam of scanning of the input image, The widths of unequal intervals in the region of the seam are obtained, and the weighting factors corresponding to the obtained corresponding points and the obtained unequal intervals are read out from the two-dimensional table and read out. An unequal interval interpolation method characterized by interpolating image data arranged at unequal intervals using weighting factors. 2. The non-equidistant interpolation method according to claim 1, wherein the region of the seam of the scan is a portion of one pixel at both ends in one scan. Interpolation method. 3. The unequal interval interpolation method according to claim 1, wherein the weighting coefficient is generated by a splice function whose coefficient changes in accordance with the width of the unequal interval. Equal interval interpolation method. 4. The non-equidistant interpolation method according to claim 1, wherein the two-dimensional table is created by dividing and editing into a rectangular table.