JPS60204083A - Unequal interval interpolation system - Google Patents

Abstract

PURPOSE:To improve the speed and the accuracy by calculating and storing preliminarily weight coefficients, which are used for unequal interval interpolation, in a form of a two-dimensional table and changing the number of degree of a spline function in accordance with the width of unequal intervals. CONSTITUTION:A coefficient of distortion correction in the longitudinal direction is used to obtain a width (d) of unequal intervals, and six weight coefficients WSPi (x, d) (i=1-6) corresponding to a position (x, d) of the two-dimensional weight coefficient table in accordance with an interpolation position (x) and this width (d), and an interpolation value I(x, d) is obtained in accordance with picture intensity values I1(i=1-6) of six nearby points and these weigt coefficients by an equation. In this case, (x) and (d) are the most approximate representative points, and the value of (x) is increased by every 1/16, and that of (d) is increased by every 1/8. This processing is performed for all points of a correction picture to complete the distortion correction processing in the longitudinal direction.

Description

DETAILED DESCRIPTION OF THE INVENTION [Field of Application of the Invention] The present invention relates to an interpolation method for data arranged at irregular intervals, and particularly to an irregular interval interpolation method suitable for performing interpolation at high speed and with high precision.

[Background of the invention]

Before explaining the prior art, the geometric distortion correction processing of satellite images, which is the background of the present invention, will be explained.

Images of the earth's surface taken by artificial satellites are made up of changes in the satellite's orbit and attitude, the Earth's rotation, and imaging equipment (called sensors).
Contains geometric and strength distortions due to internal factors. When using satellite images, it is necessary to correct these distortions by digital image processing. In the following, geometric distortion correction processing will be described, and distortion simply refers to geometric distortion.

The principle of distortion correction processing called resampling, which has been carried out in satellite ground stations such as Landsat ground stations, is as shown in FIG. 1.

First, a point (u, v) on the uncorrected image 1 corresponding to an arbitrary pixel position (X, Y) on the corrected image 2 is found. This map ψ is determined from the satellite's orbit and attitude data. Point (u
, v) is between the uncorrected image data, the true image intensity value is estimated by two-dimensional interpolation from some surrounding data, and the point (X, y) on the corrected image is Let the image intensity value be .

Since this two-dimensional interpolation process can generally be performed at high speed by replacing two one-dimensional interpolations vertically, the discussion below will be limited to one-dimensional interpolation.

When the distortion changes slowly enough, each pixel on the uncorrected image can be considered to be locally arranged at equal intervals, and the interpolation is equal interval interpolation.

Regarding equally spaced interpolation, the cubic convolution method is widely known. (For example, “Digita
i ■mage l(, construction
andl(, esampling for QeOme
triCManipulation.

K.W.Simon, presented at th
e SymposjumOn Machine pro
cessing of Ren10telySense
d1) ata, June 3-5.1975 J. ) In the cubic skewer convolution method, as shown in Figure 2, the weighting coefficient W for the interpolation position
! ” (X), i = 1, 4 is made into a table in advance, and the interpolated image intensity value I (X) is
By using the intensity value I+i=1.4 of the point, high-speed interpolation processing is possible.

Here, one shift is usually taken in 1/8'' or 1/16 increments at the representative point closest to xijx.

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

By the way, one of the typical sensors installed on Landsat satellites is the multispectral scanner (IVlul).
ti 8pQCtral 5Canner, abbreviated as Mf3s) and Thematic
Mapper.

(abbreviated as TM), all of which are scanning type sensors. The imaging method of the scanning sensor will be explained according to FIG. The satellite main body 3o1 has a vibration scanning arm 302. An optical system 303 and a detector 304 are mounted. As the satellite moves forward, it images the ground surface immediately below it, but with a constant strain sensor, the vibration scanning mirror 302 repeatedly scans the ground surface in a direction perpendicular to the orbital direction 305, thereby covering the ground surface over a scanning range 306. I will continue to take images. In one scan, IjIiii element 30
Data for as many lines as there are 7 quick 07 detectors 304 can be obtained. The number of detectors is 6 for MSS and 16 for TM.

Therefore, when the ground surface is imaged with a scanning sensor, pixels aligned in the scanning direction and in the scanning perpendicular direction within one scan can be considered to be equally spaced, but they generally become unevenly spaced at the joints of the scans.

In the case of MSS, there was no need to consider these unevenly spaced portions because the scanning was done in one direction and the image interval in the scanning vertical direction was as wide as about 79 m. However, in the case of TM, unevenly spaced portions cannot be ignored because of the round-trip scanning and the pixel spacing of about 30 m.

(According to the actual TM sensor specifications, unequal intervals of about ±2 pixels at most will occur.) Therefore, TMua
1#! Regarding distortion correction, it is not possible to correct the distortion by using only equal interval interpolation methods such as the conventional cubic convolution method.

Conventional distortion correction methods for images containing unevenly spaced parts, such as the TM111li image, are the pseudo-cubic convolution method and the GE ((3eneralE1ec
The triC) company system is known. Pseudo cubic
Regarding convolution, patent application 1683-1983
No. 62 “Processing method for image data with scanning errors J”
Regarding the Gffff formula, [LANDSAT -D'
l'hematic Mapper Images
aml)ling for 5oan Geometer
y (::orrection.

B, p, Beyer, etal, MacJne
processingof l(, emotel)'
5ensed, [)ata sympos7um.

1981J and the aforementioned Japanese Patent Application No. 57-168362.

In the pseudo-cubic convolution method, unevenly spaced portions are regarded as equally spaced, and interpolation is performed using the cubic convolution method described above.

That is, as shown in FIG. 4, if the width of the irregular intervals at the scanning joint is d, and the interpolation position is X, the interpolated image intensity value I (X) is I (x) = j: W 's' (x/d)=I
+ --(2)1=1. Here, the weight WV is the same as that appearing in equation (1). Also, x/d is the representative point closest to x/d.

However, this wrinkle-like cubic convolution performs interpolation calculations by treating unequal intervals as if they were equal intervals, which has the drawback of lower interpolation accuracy.

Next, the GE method will be explained with reference to FIG.

D represents the line interval and the first position of the corrected image to be inserted into the GP.

In the method by GE, extended lines (17
The data 53 above (M, 18B) are processed by spline interpolation, and the seemingly unevenly spaced portions are eliminated and equally spaced interpolation (such as the cubic convolution method) is performed. The 3tli in the spline will be explained in detail in the example, but the feature of the GE method is that the weighting coefficient table is
For dimensioning, lines 15, 16.17B, 181:
The point is that they are spaced at equal intervals. (In other words, this corresponds to fixing the interpolation position line, and then performs evenly spaced interpolation), which takes processing time. Also, since the order of the spline function is fixed, it is not possible to use a spline function that corresponds to the width of the uneven intervals, and There was a drawback that insertion accuracy decreased.

[Purpose of the invention]

The purpose of the present invention is to improve the drawbacks of the prior art as described above,
An object of the present invention is to provide an irregular interval interpolation method that interpolates data arranged at irregular intervals at high speed and with high accuracy.

[Summary of the invention]

In order to achieve the above object A, the present invention provides a two-dimensional weighting coefficient table in the form of a two-dimensional table with two variables: the width of the waste and the interpolation position. The first feature is that it is prepared as A second feature is that when the weighting coefficients are generated according to the spline interpolation method, a spline function whose order changes depending on the width of the unequal intervals is used. Furthermore, a third feature is that the two-dimensional weighting coefficient table is divided and edited into a rectangular table in order to simplify address calculation and perform non-uniformly spaced interpolation at high speed.

[Embodiments of the invention]

Hereinafter, one embodiment of the present invention will be described by taking as an example the Sematetsukumatsuno<-(TM) image distortion correction system on board the Landput 4 satellite.

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

The CPU 101 calculates distortion correction coefficients necessary for distortion correction processing (resampling). Distortion correction processing is performed by the image correction processing device 105 (high-speed array processor AP120B)
I'll do it at The image data is recorded on a high-density magnetic tape device (HD
DT ) 102 , temporary magnetic disk 104
is stored in The corrected image data is output to the magnetic tape device 103.

Figure 7 shows what shop [! FIG. 3 is a schematic diagram of a processing flow for correcting geometric distortion of an Il image.

Low frequency attitude data 701, high frequency attitude data 702, orbit data 703, TM mirror scan center time data 704
When the uncorrected image 705 is input, a reasonable angle calculation 706, a trajectory calculation 707, and a scanning mirror correction 708 are performed based on these data, and scanning is controlled. Furthermore, a deviation amount 709 and a geometric forest correction coefficient 710 are calculated,
Based on this result, distortion correction processing 711 according to the present invention is performed, and corrected image data 712 is output.

The distortion correction process is performed in two stages, first in the scanning direction (horizontal direction) and then in the scanning direction (vertical direction). The image after distortion correction in the scanning direction is performed is called a hybrid image.

Since the pixel intervals in the scanning direction can be regarded as locally equal intervals, an equal interval interpolation method such as the conventional ab cupic convolution method is used. In the following, we will convert a hybrid image into a corrected image, in which irregularly spaced portions that occur at scan seams become a problem. The one-dimensional interpolation process in the vertical direction will be described in detail.

First, as shown in FIG. 8, a corresponding point 811 on the hybrid image 81 of a point 821 on the corrected image 82 is calculated using a vertical distortion correction coefficient. If the corresponding point, such as point 811, is in the area A shown in Fig. 8, that is, in the part where the distance between the two ends is II# in one eclipse, interpolation processing is performed using the usual four-point cubic convolution method. Let's do it.

If the corresponding point is in region B, such as point 812, non-uniform interpolation is performed. FIG. 9 shows a sectional view corresponding to AA' in FIG. 8.

First, using the vertical direction distortion correction coefficient, find the width d of the uneven intervals, and from the interpolation positions X and d, the position of the two-dimensional weighting coefficient table [
Six weighting coefficients WI (X, d) corresponding to (X, d)
)1 +=l, read 5, image enhancement of 6 nearby points 11
iI s r i=1 , 6 and interpolation 1 straight I (X,
d). However, x and d are the maximum M for x and d, respectively]: J) 9 at the second point, and X is 1/16. d is taken in 1/8 increments.

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

Below, a method of constructing a two-dimensional weight table will be described.

Interpolation is performed using a spline function. For details on spline functions, please refer to, for example, "Spline Functions and Their Applications," Kozo Yoshida et al., Kyoshi Publishing. A spline function is a piecewise polynomial that connects polynomials so as to satisfy some continuous condition.

The width d of the uneven intervals is -1 × d≦3 (according to the TM sensor specifications), and the interpolation position X is 0 (X (2 + d) (see Figure 9). When using a function, (a first-order spline function means a straight line in each interval, that is, a broken line) interpolation using a fold is as shown in Figure 10.
This is M-type interpolation. The method for generating the weighting coefficients at this time is self-evident and will be omitted here.

Next, a case will be explained in which the most commonly used cubic spline function is used.

First, consider the case when O<d.

Section kLl*几2 of Ml1 diagram. Let the cubic spline functions in R3 be respectively)'1 t)'11, yz'/I=aI"'X"+aノ"X"+aa"'x+
a4'''・=・=(4) (j=1.2.3).

12 coefficients a J"' (i = 1, 2.3, j = 1
゜2.3.4) is determined by the following conditions.

(1) Passing the end point of each section'/s (0)=L2)'t(1)=Ia yz(1)=I3 ・・・・・・・・・(5) Y 2 (1+d ) = I n Y a (1+ d) = I 4 Ys (2+ d) = Is (2) Continuous yf (1) up to the second differential coefficient at the joint of the interval
=yi(1) Yf(1) =y'! (1) ・・・・・・・・・(6) y≦(1+d)=y≦(1+d) y'l (x+d)=yf (1+d)(3)x=o
, specify the slope at 2+a. Set the conditions (1) to (3) to a
If we write it down in J and summarize it in the form of a matrix, we get D-"=F(!:'iF, then we get a as the solution to the simultaneous equations.
, lit can be expressed as follows.

j=1~4 a(?=(j-+j+4 on the right side of aj)(3
1 aj = (J→J+8 on the right side of ar)...
...Q [Therefore, #l1iI (X, d) is (quantity 1
fil fil I(x, d)=aI X3+a2 x2+a, x+a
4= Σ Wr (X, d)1+ ・・・・・・・・・
・Uυjmi1, and comparing both sides, Wl (X+ d)+ J=1
+..., you can add 6.

-1(d(0) (see FIG. 12), WJ(x, d) can be calculated using exactly the same idea.

When d=:=00 (see FIG. 13), I3 and I4 overlap, so after replacing them with an average of 11ii≦=(Is +I4)/2, the reduction WJ (X, 0) is calculated in the same way.

As described above, interpolation calculations can be performed at high speed by storing all WJ (X, d) determined in advance in the form of a two-dimensional table.

Next, the order of the spline function used for interpolation is
The order can be changed depending on the width d of the uneven intervals. According to the experiment shown in FIG. 14, interpolation using a cubic spline function is most accurate when d≦1.6, and interpolation using a linear spline function (that is, linear interpolation) is most accurate when d≦1.6. Here, the horizontal axis is the width of uneven intervals d, and wc@ is the interpolation phase change (几MS)
It is. ○ indicates the nearest neighbor interpolation method, X indicates the linear interpolation method, Δ indicates the pseudo-cubic convolution method, . indicates the cubic spline method, and ■ indicates the GE method. In this experiment, an interpolation process using the cubic convolution method was performed on a single original image, and the above-mentioned spline interpolation was performed on an image with simulated unequal intervals to compare the resulting restored image with the original image. This is due to the difference in images.

Based on the results of this experiment, weighting coefficients are generated using a spline function whose order changes according to d.

As described above, according to this embodiment, by calculating and storing the weighting coefficients used for non-uniform interpolation in the form of a two-dimensional table in advance, non-uniform interpolation can be executed at high speed. By changing the order of the spline function used for interpolation according to the width d of the unequal intervals, there is an effect that interpolation calculations can be performed with high accuracy.

Figure 15 compares the interpolation accuracy and amount of calculation for the invention method, GE method, pseudo-cubic method, and linear spline method.Here, the interpolation accuracy is
! The error (RMS) is measured experimentally, and the amount of calculation is expressed as the number of floating point operations required for one interpolation.

Next, an example will be described in which a two-dimensional static coefficient table is changed into a long sword shape.

The width d of the uneven intervals is according to the TM sensor specifications, d≦2
The interpolation position If this table representing the set of is stored in memory as it is, it will become as shown in Figure 17.At this time, the weighting coefficient for a certain (x, d) is ((D+15) (D+16)+X) relative to the beginning. ) Misery 6...
It can be seen that it is stored in the αth position. However, it is.

Since the a formula is a quadratic formula with respect to Dvc, the address calculation is complicated and takes a long processing time. This is because the table in FIG. 16 is not rectangular.

FIG. 18 shows a table configuration according to the method of the present invention. That is, in Figure 16, divide the part where X≧32 and repeat the prayer.
It has been edited into a rectangular table measuring 32 vertically by 33 horizontally. The arrangement in memory is shown in FIG.

At this time, the weighting coefficient for a certain (x, d) is the relative position from the beginning: (32D+X+512) When -6X≦31 ((-32D-X+575)) When 6X≧32...
......(2) becomes.

Since equation (2) is a linear equation regarding X and D, address MF4 is created for business luck.

〔Effect of the invention〕

According to the present invention, when interpolating data arranged at irregular intervals, the distortion coefficient used for the interpolation is set to 2' with two variables, the width d of the irregular intervals and the interpolation position X. V.

By storing it as a source table, you can perform non-uniform interpolation at any speed, and since you can change the order of the spline function used for interpolation according to d, you can perform non-uniform interpolation with high precision. It has the effect of being able to Furthermore, since the two-dimensional weighting coefficient table for non-uniformly spaced interpolation is laid out in a rectangular shape, there is an effect that the interpolation process can be executed at random.

[Brief explanation of the drawing]

Figure 1 is a diagram showing the principle of distortion correction, Figure 2 is a diagram explaining the equidistant interpolation method using the cubic convolution method, Figure 3 is a diagram showing the scanning type imaging system, and Figure 4 is a diagram showing the imaging method of the scanning type cell.
Figure 5 is a diagram explaining the pseudo-cubic convolution method, Figure 5 is a diagram explaining the GE method, Figure 6 is a diagram showing the hardware configuration for processing data on the Randonat No. 4 satellite, and Figure 7 is a diagram explaining the pseudo-cubic convolution method. The figure shows the flow of data necessary for correcting the image data of the RANDNADT 4 satellite, Fig. 8 shows an overview of the vertical distortion correction process, and Fig. 9 shows the A-A' gland cross section of Fig. 8. Figures 10 to 13 are diagrams showing the non-uniform interval interpolation method, Figure 14 is a diagram showing the interpolation accuracy experiment results,
FIG. 15 is a diagram showing comparison results between the method of the present invention and the conventional method, and FIGS. 16 to 19 are diagrams showing configuration examples in which the weighting coefficient table is rectangular. l...Uncorrected II! Pixel position in ll image, 2... Pixel position in corrected image, 301... Magic star body, 302.
... Vibrating scanning mirror, 304... Detector, 101... C
PU, 102... High density magnetic tape device, 103...
-Magnetic tape device, 104...Magnetic disk device, 1
05... Image correction processing device. Figure 1 Figure 2 Figure 3 Figure 2 Figure 'zy Figure 812 92nd 11τeVLJ IX Ll'lXlllJ% /j
mouth'fI1gmouth'fr17 figure%'t'f figure

Claims (1)

[Claims] 1. In the interpolation calculation of data arranged at irregular intervals, the weighting coefficients used for interpolation are stored in memory in the form of a two-dimensional table with the width of the irregular intervals and the interpolation position as variables. An unequal interval interpolation method characterized by storing data in . 2. The above weighting coefficients are interpolated depending on the width of the unequal intervals. , 3. 2. The nonuniform interval interpolation method according to claim 1 or 2, wherein the two-dimensional table is edited into a rectangular shape and stored in a memory.