WO2007088759A1 - 変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、特徴点マッチング処理方法、特徴点マッチングプログラム - Google Patents

変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、特徴点マッチング処理方法、特徴点マッチングプログラム Download PDF

Info

Publication number
WO2007088759A1
WO2007088759A1 PCT/JP2007/051095 JP2007051095W WO2007088759A1 WO 2007088759 A1 WO2007088759 A1 WO 2007088759A1 JP 2007051095 W JP2007051095 W JP 2007051095W WO 2007088759 A1 WO2007088759 A1 WO 2007088759A1
Authority
WO
WIPO (PCT)
Prior art keywords
displacement
phase
procedure
phase singularity
singularity
Prior art date
Application number
PCT/JP2007/051095
Other languages
English (en)
French (fr)
Inventor
Mitsuo Takeda
Wei Wang
Tomoaki Yokozeki
Reika Ishijima
Yu Qiao
Original Assignee
National University Corporation The University Of Electro-Communications
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National University Corporation The University Of Electro-Communications filed Critical National University Corporation The University Of Electro-Communications
Priority to US12/162,370 priority Critical patent/US7813900B2/en
Priority to JP2007556824A priority patent/JP4385139B2/ja
Publication of WO2007088759A1 publication Critical patent/WO2007088759A1/ja

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B11/00—Measuring arrangements characterised by the use of optical techniques
    • G01B11/002—Measuring arrangements characterised by the use of optical techniques for measuring two or more coordinates

Definitions

  • Displacement detection method displacement detection apparatus, displacement detection program, feature point matching processing method, feature point matching program
  • the present invention relates to a displacement detection method, a displacement detection device, a displacement detection program, a feature point matching processing method, and a feature point matching program, and in particular, displacement detection that detects displacement from phase singularities before and after displacement.
  • the present invention relates to a method, a displacement detection device, a displacement detection program, a feature point matching processing method, and a feature point matching program.
  • non-destructive inspection and material strength have been used in non-destructive inspection and material strength to measure minute displacement of an object using a pattern or texture having a spatially random structure such as a laser speckle pattern or random dot pattern as an index. It occupies an important position in industrial applications such as testing.
  • the correlation method is known as a signal processing technique for detecting minute displacements using noturns and textures.
  • FIG. 38 and FIG. 39 are diagrams for explaining the operation of the correlation method.
  • Fig. 38 (A) shows the field displacement
  • Fig. 38 (B) shows a conceptual diagram of the correlation method
  • Fig. 38 (C) shows a conceptual diagram of this measurement method.
  • FIG. 39 (A) shows the displacement of the object A in one direction
  • FIG. 39 (B) shows a diagram for explaining the rotational displacement of the object A.
  • phase-only correlation method is known as a method of calculating a correlation function using phase information.
  • This correlation method suppresses the amplitude of the composite spectrum, which is a complex function, by a constant key or a logarithmic function, etc., creates an amplitude-limited complex composite spectrum that includes only phase information, and performs inverse Fourier transform on this. A cross-correlation function is obtained.
  • Patent Document 1 Patent No. 3035654
  • the present invention has been made in view of the above points, and is an element that characterizes the structure of the phase singularity.
  • the phase singularity before and after the displacement can be reliably identified, and the displacement
  • a displacement detection method, a displacement detection device, a displacement detection program, a feature point matching processing method, and a feature point matching program that can detect various displacements including rotational displacement easily and reliably. The purpose is to provide.
  • the present invention is a displacement detection method for detecting a phase singularity force displacement before and after displacement, acquiring a predetermined element of phase structure force of a phase singularity, and specifying a phase singularity based on the element. The position of the phase singularity after the displacement is detected.
  • the acquired element interpolates the real part and imaginary part values near the phase singularity in a plane, and as a result, the line segment in which the real part becomes zero and the line segment in which the imaginary part becomes zero. It is characterized by the angle at which and intersect, the eccentricity of the shape of the intensity distribution near the phase characteristic point, and the intensity gradient near the phase singularity.
  • the position of the specified phase singularity before and after the displacement is specified based on the element, and the displacement of the phase singularity is measured based on the position of the phase singularity before and after the specified displacement. And At this time, the square root of the difference in coordinates of the position of the phase singularity before and after the displacement is calculated as a displacement amount of the phase singularity.
  • the position of the phase singularity before and after the displacement is obtained by obtaining the point where the sum of the distances of the plurality of perpendicular bisectors of the line segment connecting the positions of the phase singularities before and after the plural sets of displacements as the rotation center.
  • the angle of the side that sandwiches the rotation center of the triangle consisting of the rotation center and the rotation center is obtained as the rotation angle. To do.
  • the present invention performs a Fourier transform on the acquired intensity signal of the image and performs a filtering process that breaks the symmetry of the complex conjugate of the Fourier-transformed signal, thereby acquiring and acquiring the pseudo-phase information. It is characterized by acquiring phase singularities based on pseudo-phase information
  • filtering is performed by a Laguerre Gaussian beam filter characteristic having a helical phase characteristic having a phase singular point at the center and a donut-shaped amplitude characteristic with respect to a Fourier-transformed signal.
  • a Laguerre Gaussian beam filter characteristic having a helical phase characteristic having a phase singular point at the center and a donut-shaped amplitude characteristic with respect to a Fourier-transformed signal.
  • an image obtained by subtracting the average value of the intensity signals of the image from the intensity signals of the image is used as the original image.
  • the present invention is characterized in that a pseudo phase information is obtained by performing filtering by a filter characteristic having a helical phase characteristic having a phase singular point at the center of a Fourier transformed signal or by Hilbert transform. To do.
  • the present invention is a point matching processing method for matching a feature point specified before displacement with a feature point specified after displacement, wherein the feature parameters of the feature point before displacement are Based on the difference from the feature parameter of the feature point after the displacement, the feature point after the displacement corresponding to the feature point before the displacement is found, and the similarity force of the geographical position relation of the feature point is found The feature point before the displacement The feature points after the displacement corresponding to are matched.
  • the present invention is also characterized in that the difference between the total distance of feature parameters between a plurality of feature points before displacement and the total distance of feature parameters between the plurality of feature points after displacement is minimized. It is characterized by specifying a set of characteristic points as characteristic points that match each other before and after displacement.
  • the feature parameter of the predetermined feature point before the displacement is calculated after the Tth place.
  • a feature point set that minimizes is extracted as a set of feature points that match before and after displacement.
  • the feature point is a phase singularity acquired based on the intensity information of the image.
  • the phase singularity is characterized by positional information, eccentricity, vorticity, and alignment as feature parameters. And having an angle at which the real axis and the imaginary axis intersect.
  • each phase singularity before and after displacement is uniquely identified by acquiring a predetermined element of the phase structural force of the phase singularity and identifying the phase singularity based on the element.
  • FIG. 1 is a system configuration diagram of an embodiment of the present invention.
  • FIG. 2 is a processing flowchart of a displacement detection program.
  • FIG. 3 is a process flowchart of phase singularity acquisition processing.
  • FIG. 4 is an operation explanatory diagram of phase singularity acquisition processing.
  • FIG. 5 is an operation explanatory diagram of phase singularity acquisition processing.
  • FIG. 6 is a diagram for explaining the operation of determining the position of the phase singularity with subpixels.
  • FIG. 7 is a diagram for explaining an operation of determining the position of a phase singularity by subpixels.
  • FIG. 8 is a diagram for explaining an operation of determining the position of a phase singularity by subpixels.
  • FIG. 9 is a process flowchart of a phase singularity specifying process.
  • FIG. 10 is a diagram for explaining an operation for obtaining an eccentricity e and a vorticity ⁇ .
  • FIG. 11 is a diagram for explaining an operation for obtaining an eccentricity e and a vorticity ⁇ .
  • FIG. 12 is a diagram for explaining an operation for obtaining an eccentricity e and a vorticity ⁇ .
  • FIG. 13 This figure shows the difference between the positions of all phase singularities and shows the displacement in the X and y directions in a histogram.
  • FIG. 14 is a diagram for explaining a rotational displacement detection operation.
  • FIG. 15 is a diagram for explaining a simulation experiment.
  • FIG. 16 is a histogram showing the number of phase singularities before and after rotation.
  • FIG. 17 is a histogram showing rotation angles ⁇ i of phase singularities before and after rotation.
  • FIG. 18 is a histogram showing the difference in angle ⁇ before and after rotation.
  • FIG. 19 is a histogram showing the eccentricity e of the corresponding phase singularities before and after rotation.
  • FIG. 20 A histogram showing the difference in vorticity ⁇ of corresponding phase singularities before and after rotation.
  • FIG. 21 is an operation explanatory diagram of point matching processing.
  • FIG. 22 is a flowchart of point matching processing.
  • FIG. 23 is a diagram for explaining the processing operation of the matching algorithm.
  • FIG. 24 is a diagram for explaining the processing operation of the matching algorithm.
  • FIG. 25 is a diagram for explaining the processing operation of the matching algorithm.
  • FIG. 26 is an amplitude characteristic diagram of the LG filter used for the experiment 'verification.
  • FIG. 27 is a diagram for explaining the operation of matching processing.
  • FIG. 28 is a diagram for explaining the operation of matching processing.
  • FIG. 29 is a diagram for explaining the results of the experiment and verification.
  • FIG. 30 is a diagram for explaining the results of the experiment and verification.
  • FIG. 31 is a diagram for explaining the results of the experiment / verification.
  • FIG. 32 is a diagram for explaining the results of the experiment and verification.
  • FIG. 33 is a diagram for explaining the results of the experiment and verification.
  • FIG. 34 is a diagram for explaining the results of the experiment / verification.
  • FIG. 35 is a diagram for explaining the results of the experiment and verification.
  • FIG. 36 is a diagram for explaining an operation of performing matching based on a difference between the total distances between points.
  • FIG. 37 is a diagram for explaining an operation of performing matching based on a difference between the total distances between points. [38] FIG. 38 is a diagram for explaining the operation of the correlation method.
  • FIG. 39 is a diagram for explaining the operation of the correlation method.
  • FIG. 1 shows a system configuration diagram of an embodiment of the present invention.
  • the displacement detection device 101 of this embodiment includes an image detection unit 111 and a signal processing unit 112.
  • the image detection unit 111 includes a light source 121 and an imaging device 122.
  • the image detection unit 111 irradiates the measurement target 102 with light from the light source 121, and images the surface image of the measurement target 102 with the imaging device 122.
  • the image captured by the image detection unit 111 is supplied to the signal processing unit 112.
  • the signal processing unit 112 is, for example, a computer system, and includes a frame memory 131, a processing unit 132, a Fourier conversion unit 133, a ROM 134, a hard disk drive 135, a display 136, and an input device 137.
  • the frame memory 131 stores the image from the image detection unit 111 for each frame.
  • the processing unit 132 executes processing for detecting a minute displacement of the measurement object 102 based on a displacement detection program installed in the hard disk drive 135 in advance.
  • the Fourier transform unit 133 performs Fourier transform on the image captured by the image detection unit 111.
  • the memory 134 is used as a working storage area for the processing unit 132.
  • the hard disk drive 135 stores a displacement detection program, image data detected by the image detection unit 111, or data such as an intermediate result in the displacement detection program and a displacement detection result.
  • the display 136 is configured with a force such as an LCD and a CRT, and displays processing results in the processing unit 132 such as images, converted images, intermediate results, and displacement detection results.
  • the input device 137 includes a keyboard, a mouse, and the like, and is used to input data and various commands.
  • the displacement detection apparatus 101 of the present embodiment analyzes an image acquired from the measurement object 102 based on the installed displacement detection program, and detects the displacement. [Displacement detection program]
  • FIG. 2 shows a processing flowchart of the displacement detection program.
  • the processing unit 132 first reads out the target images before and after displacement from the node disk drive 135 to the memory 134 or the like in step S1-1 according to the displacement detection program.
  • step S1-2 the processing unit 132 executes phase singularity acquisition processing for acquiring phase singularities for the read target images before and after displacement.
  • the processing unit 132 performs phase singularity based on factors such as the zero cross angle between the real part and the imaginary part, the eccentricity, and the vorticity of the phase singularity acquired from the target image before and after displacement in step S1-3.
  • a phase singularity specifying process for specifying a point is executed.
  • step S1-4 the processing unit 132 determines the phase singularities corresponding to before and after the displacement by factors unique to the phase singularities such as the zero cross angle between the real part and the imaginary part, the eccentricity, and the vorticity. Identify
  • the displacement of the target image can be accurately identified.
  • step S 1-2 Next, the phase singularity acquisition process in step S 1-2 will be described.
  • FIG. 3 is a process flowchart of the phase singularity acquisition process
  • FIGS. 4 and 5 are operation explanatory diagrams of the phase singularity acquisition process.
  • step S12 the phase singularity acquisition process performs a Fourier transform on the intensity signal of the acquired image and performs a filtering process that breaks the symmetry of the complex conjugate of the Fourier-transformed signal. Is obtained, and the phase singularity is obtained based on the obtained pseudo-phase information.
  • the processing unit 132 first performs preprocessing in step S2-1.
  • the pre-processing in step S2-1 is a process of obtaining an average value of the intensity information of the entire target image and subtracting the obtained average value of the intensity information of the entire image from the intensity information of the target image. As a result, the DC component in the image data intensity information is removed. By removing the DC component, phase singularities are emphasized.
  • the processing unit 132 controls the Fourier transform unit 133 in step S2-2 to perform a two-dimensional fast Fourier transform on the intensity information to obtain a spatial frequency spectrum.
  • the function 201 of the texture intensity information 211 of the object obtained by the image detection unit 111 is obtained.
  • 212 indicates the amplitude distribution of the spatial frequency spectrum
  • 213 indicates the phase distribution. Since the spatial frequency spectrum is a Fourier transform of the intensity information with only the real part, it has complex conjugate symmetry with respect to the frequency origin, that is, Hermitian symmetry. If the FFT is a function F, for example, it can be expressed as Equation 202.
  • the processing unit 132 generates a imaginary part based on the real part and obtains a complex analysis signal in step S2-3, and has a frequency domain having a distribution as indicated by 221 to 225 in FIG.
  • the Hilbert transform, spiral phase filter, and LG (Laguerre Gauss) filter are used in the filtering process.
  • 221 indicates the amplitude distribution of the Hilbert transform
  • 222 indicates the phase distribution of the spiral phase filter
  • 224 indicates the amplitude distribution of the LG filter
  • 225 indicates the phase distribution of the LG filter.
  • An arithmetic expression of the Hilbert transform can be expressed by, for example, Figs.
  • the arithmetic expression of the spiral phase filter can be expressed by, for example, FIGS.
  • the calculation formula of the amplitude transmission distribution 224 of the LG filter can be expressed by FIG. 4, 205 (A).
  • the LG filter obtains a donut-shaped amplitude transmittance distribution by shifting the complex conjugate by ⁇ on the phase symmetrically with the center.
  • the arithmetic expression of the phase distribution 225 of the LG filter is expressed as shown in FIGS. 4 and 205 ( ⁇ ) by introducing the angle ⁇ into FIGS. 4 and 205 ( ⁇ ), for example.
  • an analysis signal required in this embodiment can be obtained.
  • high-frequency components that cause errors can be removed, and DC components can be zeroed automatically and continuously.
  • the number of phase singularities can be increased by removing low-frequency components from the donut-shaped intensity distribution.
  • analysis is performed using intensity information of a random pattern, and there are many patterns such as speckles depending on the structure of the measurement object 102 and the like.
  • This speckle has a corresponding relationship with the phase singularity used as an index in this embodiment.
  • the frequency including a lot of phase singularity information in the frequency domain when Fourier transformed is changed.
  • the frequency component including phase singularity information is high, and conversely if the size is large, the frequency component including phase singularity information is low.
  • an LG filter with a doughnut-shaped intensity distribution is optimal for obtaining phase singularity information with high accuracy, and the phase of the phase can be increased by adjusting the size of the donut-shaped width according to the speckle size. Singularity information can be obtained with high accuracy.
  • the low-pass filter emphasizes low-frequency components and cannot extract phase singularity information well.
  • FIG. 4 214 indicates the amplitude distribution after the filtering process
  • 215 indicates the phase distribution after the filtering process
  • the arithmetic expression thereof can be expressed by, for example, FIGS. [0053] It should be noted that the same kind of analysis signal can be obtained even if a Hilbert transform or spiral phase filter as shown in 221 and 223 is applied. In order to use the Hilbert transform or spiral phase filter, the DC component is removed in advance. I have to keep it.
  • the processing unit 132 performs Fourier inverse transform on the spatial frequency spectrum filtered in step S2-4 to perform a complex analysis signal.
  • 216 indicates the amplitude distribution after the inverse Fourier transform
  • 217 indicates the phase characteristic after the inverse Fourier transform
  • the arithmetic expression thereof can be expressed by, for example, FIGS.
  • FIG. 5A is an enlarged amplitude distribution after inverse Fourier transform
  • FIG. 5B is an enlarged view of a predetermined portion of FIG. 5A.
  • a phase singular point P as shown in FIG. 5 (B) can be obtained from the pseudo phase information of this complex analysis signal.
  • step S2-5 the processing unit 132 performs a process of determining the position of the phase singularity with the resolution of the sub-pixel having a resolution smaller than one pixel.
  • FIG. 6, FIG. 7, and FIG. 8 are diagrams for explaining the operation of determining the position of the phase singularity by the subpixel.
  • FIG. 6A shows pseudo-phase information
  • FIG. 6B shows an enlarged view of pixels around the phase singular point
  • FIG. 6C shows the real part
  • FIG. 6D shows the imaginary part
  • Figure 7 (A) shows the quasi-phase information near the phase singularity
  • Figure 7 (B) shows the imaginary part of the quasi-phase information near the phase singularity
  • Figure 7 (C) shows the quasi-phase information near the phase singularity
  • Figure 7 (D) is a linear interpolation of the pseudo phase information near the phase singularity
  • Figure 7 (E) is a linear interpolation of the imaginary part of the pseudo phase information near the phase singularity
  • (F) shows a linear interpolation of the real part of the pseudo-phase information near the phase singularity.
  • the pseudo-phase information takes values up to (- ⁇ , ⁇ ), and increases in value from black to white in Figs.
  • phase singularity ⁇ has the property that the integral of the phase gradient on the closed path surrounding the point is an integral multiple of ⁇ 2 ⁇ , and the real and imaginary parts of the complex number are both zero. .
  • the integral of the phase gradient is an integral multiple of ⁇ 2 ⁇ for a square closed path connecting four adjacent pixels as shown in Fig. 6 ( ⁇ ) and Fig. 7 (C).
  • the position where the phase singularity exists was determined. For this reason, the position of the phase singularity cannot be determined with a resolution of 1 pixel or more.
  • the position of the phase singularity is determined with subpixel resolution smaller than 1 pixel resolution by using the property that both the real part and imaginary part of the complex number are zero.
  • Fig. 7 (E) and Fig. 7 (F) are performed by linearly interpolating the discretized real part and imaginary part in the vicinity of the singular points shown in Fig. 7 (B) and Fig. 7 (C).
  • the real part and imaginary part are both smooth, and the line segment Lre0 where the real part becomes zero and the line segment LimO where the imaginary part becomes zero are obtained.
  • the position of the phase singularity is determined by setting the point where the line segment LreO where the real part is zero and the line segment LimO where the imaginary part is zero as the phase singular point P0.
  • the phase singularity can be determined with resolution.
  • FIG. 9 shows a process flowchart of the phase singularity specifying process.
  • the phase singularity specifying process of step S1-3 is a process of acquiring a predetermined element from the phase structure of the phase singularity and specifying the phase singularity based on the element.
  • the processing unit 132 acquires a zero cross angle as an element for specifying a phase singular point in step S3-1.
  • the zero-crossing angle is the line segment LreO where the real part shown in Fig. 7 (D) and Fig. 8 (B) is zero and imaginary.
  • the angle ⁇ ⁇ intersects with the line segment LimO where the part becomes zero.
  • the processing unit 132 acquires the eccentricity e as an element for specifying the phase singularity in step S3-2, and acquires the vorticity ⁇ in step S3-3.
  • FIG. 10, FIG. 11, and FIG. 12 are diagrams for explaining the operation for obtaining the eccentricity e and the vorticity ⁇ .
  • Fig. 10 (A) shows the pseudo-phase information
  • Fig. 10 (B) shows the amplitude distribution near the phase singularity of the pseudo-phase information
  • Fig. 10 (C) shows the phase distribution!
  • FIG. 11 shows a structural diagram of the amplitude distribution near the singular point.
  • Figure 11 ( ⁇ ) shows the amplitude displayed with contour lines
  • Fig. 11 ( ⁇ ) shows the three-dimensional display.
  • the amplitude distribution in the vicinity of the phase singularity is substantially elliptical, and its shape varies depending on the phase singularity, and can be used as an element for specifying the phase singularity.
  • a phase singularity can be identified by indicating this elliptical feature by the eccentricity, which is a value representing the characteristic of the quadratic curve.
  • the eccentricity e of the ellipse shown in Fig. 12 is a: major axis radius, b: minor axis radius, f: focal length.
  • the vorticity in the vicinity of the phase singularity is obtained as one of the elements that specify each phase singularity.
  • the vorticity is a quantity that represents the strength of local rotation and the direction of rotation in a small part.
  • the vorticity increases as the interval between the contour lines of the amplitude becomes narrower.
  • the vorticity increases as the gradient (degree of change) near the phase singularity displayed in three dimensions becomes steeper. . In other words, the vorticity depends on the amplitude near the phase singularity.
  • the phase singularity has the elements of position (x, y), zero cross angle ⁇ , eccentricity e, and vorticity ⁇ , and each phase singularity is uniquely identified using these elements. be able to . At this time, the phase singularity can be reliably identified by using a plurality of elements.
  • step S1-4 the phase singularity displacement detection process in step S1-4 will be described.
  • phase displacement detection process the phase singularities before and after the displacement are identified by the above identification method. Thus, the displacement amount is detected.
  • the phase singularity is determined using the zero cross angle ⁇ , the eccentricity e, and the vorticity ⁇ obtained as described above as parameters for uniquely identifying the phase singularity before and after the displacement. Identify. The method assumes that the parameters before displacement are 0 1, el, ⁇ 1 and the parameters after displacement are 0 2, e2, ⁇ 2.
  • FIG. 13 is a graph showing the difference in the positions of all phase singularities and the displacement in the X and y directions in a histogram.
  • Fig. 13 (A) shows the displacement in the X direction
  • Fig. 13 (B) shows the displacement in the y direction.
  • each phase singularity can be specified before and after rotational displacement, so that rotation can be detected.
  • FIG. 14 is a view for explaining the rotational displacement detection operation.
  • the center of rotation can be determined as described above.
  • the strength information obtained in the experiment was rotated 90 degrees inside the computer, and a simulation experiment was performed using the strength information before and after the rotation.
  • FIG. 15 is a diagram for explaining a simulation experiment.
  • Figure 15 shows the number of pixels
  • the intensity information power before and after rotation at 1024 X 1024 is also the amplitude distribution and phase distribution when two analytical signals are obtained using the LG filter.
  • Fig. 15 (A) shows the pseudo-phase information before rotational displacement
  • Fig. 15 (B) shows the amplitude distribution of the pseudo-phase information before rotational displacement
  • Fig. 15 (C) shows the phase distribution of the pseudo-phase information before rotational displacement
  • Fig. 15 (D) shows the pseudo phase information after rotational displacement
  • FIG. 15 (E) shows the amplitude distribution of the pseudo phase information after rotational displacement
  • FIG. 15 (F) shows the phase distribution of the pseudo phase information after rotational displacement.
  • the intensity information power before and after rotation at 1024 X 1024 is also the amplitude distribution and phase distribution when two analytical signals are obtained using the LG filter.
  • Fig. 15 (A) shows the pseudo-phase information before rotational displacement
  • FIG. 16 is a histogram showing the number of phase singularities after rotation corresponding to the phase singularities before rotation.
  • FIG. 17 shows a histogram representing the rotation angle ⁇ i obtained for the phase singularity after rotation corresponding to the phase singularity before rotation and the central force of rotation.
  • the histogram of the rotation angle which is the simulation result, shows that all phase singularities are rotated 90 degrees with an accuracy of ⁇ 0.001 degrees.
  • the measurement object 102 is rotated 90 degrees, and the measurement result shown in FIG. 17 also indicates that the measurement result is rotated 90 degrees. It can be confirmed that it can be measured.
  • FIG. 18 is a histogram showing the difference in angle ⁇ 0 between corresponding phase singularities before and after rotation in a histogram
  • FIG. 19 is a histogram showing the eccentricity e of the corresponding phase singularity before and after rotation in a histogram
  • 20 shows a histogram representing the difference in vorticity ⁇ at the corresponding phase singularity before and after rotation.
  • the point matching process is a process for detecting a pattern force before displacement using a positional relationship between arbitrary points and the similarity of the spatial structure of the points.
  • a point having a similar phase structure is set as a candidate point, and a point having the smallest distance between points in the pattern including the candidate point is determined as a pattern after displacement.
  • phase singularity is extracted by performing the filter operation for extracting the above-described phase singularity on the two intensity distribution images before and after the displacement.
  • the position of all phase singularities is specified with subpixel accuracy.
  • the eccentricity e, the angle ⁇ between the real axis and the imaginary axis, and the vorticity ⁇ are obtained and labeled for each phase singularity whose position is specified.
  • FIG. 21 shows an operation explanatory diagram of the point matching process.
  • Fig. 21 ( ⁇ ) shows the state of the phase singularity before displacement
  • Fig. 21 ( ⁇ ) shows the state of the phase singularity after displacement.
  • phase singularity indicated by a black circle is selected as the phase singularity force shown in FIG.
  • the phase singularities to be selected should be those whose characteristic parameters are the eccentricity e, the angle between the real and imaginary axes ⁇ ⁇ , and the vorticity ⁇ that are clearly different from many other phase singularities. . this This makes it easier to specify the phase singularity after displacement.
  • phase singularity the characteristic parameters of the phase singularity, the eccentricity e, the vorticity ⁇ , and the angle ⁇ between the real axis and the imaginary axis are the same before and after the displacement, or hardly change. Therefore, it is possible to accurately measure movement, rotation, scale, etc. by accurately identifying the phase singularities before and after the displacement.
  • FIG. 22 shows a flowchart of the point matching process.
  • step S4-1 the processing unit 132 acquires and labels phase singularities before and after displacement by the phase singularity acquisition process described with reference to FIG.
  • the processing unit 132 replaces the phase singularities having the characteristic parameters that are substantially the same as the phase singularities before the displacement in step S4-2 or that have characteristic parameters approximated within a predetermined range. Search from special points.
  • the processing unit 132 selects two predetermined points as the phase singularity force before displacement in step S4-3, and combines them to obtain a combination ⁇ . Next, the processing unit 132 assigns “3” to the variable k in step S4-4.
  • processing section 132 determines whether or not variable k is smaller than constant n in step S4-5.
  • n is the number of phase singularities to be searched.
  • the processing unit 132 extracts a phase singularity as a next search target from the plurality of phase singularities selected in step S4-6.
  • the processing unit 132 searches for the phase singularity after the displacement approximated by the feature meter with respect to the phase singularity newly extracted in step S4-7, and adds it as the combination M.
  • step S4-8 the processing unit 132 calculates the characteristic parameter of the phase singularity of the combination M by the following equation (1), and calculates a value less than the threshold value as the phase singularity before displacement. This is left as a candidate for a powerful combination of a point and a phase singularity after displacement. [0123] [Equation 10]
  • t represents the points before and after displacement. Note that t is w (i) w (i) corresponding to.
  • ⁇ ⁇ , ⁇ , ⁇ , ⁇ ' ⁇ , ⁇
  • I t t I is the distance between point t and point t
  • ⁇ is the scale (degree) of scale il i2 il i2
  • variable w represents the correspondence between ⁇ and ⁇ .
  • the processing unit 132 substitutes (k + 1) for the variable k in step S4-9, and returns to step S4-5.
  • the processing unit 132 repeats the above steps S4-6 to S4-9 until the variable k becomes larger than the constant n in step S4-5.
  • the processing unit 132 extracts, as a matching pattern, the one having the smallest distance between the feature parameters of the phase singularities of the combination M remaining in step S4-10.
  • the processing unit 132 converts the phase singularity conversion parameter that minimizes the distance between the feature parameters extracted in step S4-11, for example, translation (shift), rotation (rotation), scale (s cale).
  • the target can find the best match before and after the displacement in the positional relationship of a plurality of points in the pattern.
  • Figs. 23, 24, and 25 are diagrams for explaining the processing operation of the matching algorithm.
  • phase singularities are selected from the phase singularities before displacement.
  • Eccentricity e which is a characteristic parameter of each phase singularity t in the pattern before displacement, zero cross angle ⁇ , vorticity ⁇ , and eccentricity e, which is a characteristic parameter of phase singularity pj after displacement, zero cross angle
  • ⁇ and vorticity ⁇ are compared with a threshold value. For example, it can be expressed by the following equation.
  • E is the eccentricity threshold
  • ⁇ is the zero-crossing angle threshold
  • ⁇ is the vorticity threshold
  • a circle indicated by a solid line indicates a positive phase singularity
  • a circle indicated by a broken line indicates a negative phase singularity by a dotted line.
  • the positive phase singularity and the negative phase singularity are positive phase singularities when the characteristic parameter vorticity value is positive, and negative ones. Negative phase singularity.
  • phase singularity t and t in space T are selected in step S4-3. This phase singularity
  • Points t and t have at least one matching candidate point.
  • the phase singularity t has n matching points il i2 il 1 and the phase singularity t has n matching candidate points. Therefore,
  • the temporal matching is repeated until all matching points of the phase singularities are finally determined while acquiring new matching points.
  • equation (1) the similarity of temporal matching is examined, and at the same time, the operation of excluding matching points detected by mistake in the distance threshold is also performed.
  • FIG. 26 and FIG. 27 are diagrams for explaining the operation of the matching process.
  • Figures 26 (A) to (C) and 27 (A) to (C) show the distribution of phase singularities before displacement, that is, the distribution of phase singularities in T-space
  • Figure 26 (D) to (F) show the distribution of phase singularities after displacement, that is, the distribution of phase singularities in P space.
  • FIGS. 26 (A) and 26 (D) feature parameters of the phase singularities of T space and P space are calculated, and the point of interest is selected.
  • Fig. 26 (A) the combination of tl, t2, and t4 was selected.
  • a set of phase singularities smaller than the threshold is extracted from the calculation results. For example, here, it is assumed that the combinations Ml (pl, p2, p4) and M2 (pl, p2, p3) are extracted. The extracted phase singularity pairs M1 and M2 are extracted as matching combinations.
  • the lowest cost that is, the set Ml having the lowest value of Equation (1) is extracted as a matched set. To do.
  • the movement of the phase singularity can be expressed by three parameters: parallel movement ⁇ , Ay, rotation 0, and scale ⁇ . These three parameters are called conversion parameters and are calculated in steps S4-11.
  • FIG. 28 shows the experiment.
  • the amplitude characteristics of the LG filter used for the verification, and Figs. 29 to 35 show the experiment.
  • FIG. 29 shows an experimental image
  • FIG. 29A shows an image before displacement
  • FIG. 29B shows an image after displacement
  • Fig. 30 is an image showing the phase singularities from which the image forces shown in Fig. 29 are also extracted.
  • Fig. 29 (A) shows the distribution of phase singularities before displacement
  • Fig. 29 (B) shows the distribution of phase singularities after displacement.
  • Fig. 31 is an enlarged view of the selected phase singularity.
  • Fig. 31 (A) shows an enlarged view before displacement
  • Fig. 31 (B) shows an enlarged view after displacement.
  • LG filter operation is performed on the captured image to obtain a phase singularity.
  • the LG filter is a spiral phase filter with the characteristics of a bandpass filter, and has an amplitude characteristic as shown in FIG. If the amplitude width ⁇ of the LG filter is increased, the number of phase singularities that appear will increase, and conversely, it will decrease as the amplitude narrows.
  • phase singularity force before displacement shown in Fig. 30 (A) Six phase singularities indicated by black circles were selected and matched.
  • the phase singularity indicated by the black circle corresponds to the phase singularity indicated by the black circle before displacement, and is a pattern actually detected by the matching algorithm.
  • FIG. 31 is an enlarged view of the combination of the phase singularities selected before the displacement and the combination of the phase singularities after the displacement detected by the point matching algorithm. In this way, when the patterns before and after displacement are compared, deformation that combines movement, rotation, and scale has occurred, but this was verified by a matching algorithm.
  • FIG. 32 is a diagram showing the distribution of phase singularities before and after displacement when the point matching process of the present embodiment is applied focusing on different patterns for two different images.
  • 32 (A) shows the distribution before displacement
  • FIG. 32 (B) shows the distribution after displacement.
  • Figure 33 is shown in Figure 32
  • FIG. 33 (B) is an enlarged view of the selected phase singularity.
  • FIG. 34 is a diagram showing an image when an experiment “verification example” is performed on continuous images.
  • FIG. 34A shows the first image
  • FIG. 34B shows the last image
  • Figure 35 focuses on the six phase singularities from the image shown in Figure 34 and tracks their trajectories.
  • phase singularities can be continuously tracked by using the point matching processing of the present embodiment.
  • phase singularity pair is detected based on the phase structure of the phase singularity and the distance between the other phase singularities, and the phase singularities are detected. It can be seen that by tracking the pair, random movement can be tracked. It has become possible to measure the movement and displacement of an object with a large displacement that was impossible to measure because the characteristic parameters changed.
  • matching was performed by detecting the combination of the phase singularities that minimizes Equation (1), but the difference between the point-to-point total distances of the phase singularities was calculated. Then, matching may be performed by detecting a phase singularity that minimizes the difference in the total distance between points between the pattern before displacement and the pattern after displacement.
  • Figs. 36 and 37 are diagrams for explaining the operation of performing the matching based on the difference between the total points.
  • Figs. 36 (A) and 37 (A) show the distribution of phase singularities before displacement
  • Figs. 36 (B) and 37 (A) show the distribution of phase singularities after displacement.
  • phase singularities p to p are extracted from the image before displacement, and from the image after displacement.
  • phase singularities p ', p', p ', p' are different from the image after displacement.
  • the point matching process of the present embodiment uses the positional relationship of a plurality of phase singularities, the point matching is performed even if the positional relationship does not completely match before and after displacement due to scale, movement, and rotation.
  • the closest set can be selected, so that the present invention can measure any displacement.
  • the phase singularity having, as the feature parameter, the position information, the eccentricity, the vorticity, and the zero cross angle that is the crossing angle between the real axis and the imaginary axis is described as the feature point.
  • the feature point is not limited to this, and may be only position information, for example.
  • spatial random pattern intensity information is associated with phase information without an interferometer, and before and after displacement based on the characteristics of all phase singularities existing in the phase information.

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Length Measuring Devices By Optical Means (AREA)
  • Image Analysis (AREA)

Abstract

　本発明は、変位前後の位相特異点から変位を検出する変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、特徴点マッチング処理方法、特徴点マッチングプログラムに関し、変位前後の位相特異点から変位を検出する変位検出方法であって、位相特異点の位相構造から所定の要素を取得し、要素に基づいて位相特異点を特定し、変位後の位相特異点の位置を検出することにより、変位前後の位相特異点を確実に特定し、その変位を容易、かつ、確実に、さらには、回転変位を含む種々の変位を検出可能とすることを特徴とする。

Description

変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、特 徴点マッチング処理方法、特徴点マッチングプログラム
技術分野
[0001] 本発明は変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、特徴 点マッチング処理方法、特徴点マッチングプログラムに係り、特に、変位前後の位相 特異点から変位を検出する変位検出方法、及び、変位検出装置、変位検出プロダラ ム、並びに、特徴点マッチング処理方法、特徴点マッチングプログラムに関する。 背景技術
[0002] 近年、レーザースペックルパターンやランダムドットパターンなどのように空間的にラ ンダムな構造をもつパターンやテクスチャを指標として物体の微小変位を非接触で 計測する技術は非破壊検査や材料強度試験など産業応用上の重要な位置を占め て 、る。ノターンやテクスチャを用いて微小変位を検出する際の信号処理技術として 相関法が知られている。
[0003] 図 38、図 39は相関法の動作を説明するため図を示す。図 38 (A)は場の変位、図 38 (B)相関法の概念図、図 38 (C)は本測定法の概念図を示している。また、図 39 ( A)は物体 Aの一方向への変位、図 39 (B)は物体 Aの回転変位を説明するための図 を示している。
[0004] 相関法を用いて変位を検出する場合、図 38 (A)に矢印で示すように場の変位が位 置により異なるときには図 38 (B)に示すように各領域に分けて計算する必要があり、 各点での相関値を決定するごとに 2次元相関積分計算または 2次元フーリエ変換'逆 変換を行う必要があつたため、 1ピクセル以下での変位検出の分解能を上げようとし てピクセル内部での相関値の計算点数を増やすと計算時間が膨大になっていた。
[0005] また、従来の相関法では、図 39 (A)に示すように一方向の変位に対しては相関を とることにより変位量を求めることができる力 図 39 (B)に示すように回転に対しては 適用できなかった。
[0006] なお、位相情報を利用する相関関数の計算方法としては、位相限定相関法が知ら れている（特許文献 1参照)。この相関法は、複素関数である合成スペクトルの振幅を 一定ィ匕又は対数関数等により抑制し、位相情報のみカゝらなる振幅限定複素合成スぺ タトルを作成してこれをフーリエ逆変換して相互相関関数を求めるものである。
特許文献 1：特許第 3035654号
発明の開示
発明が解決しょうとする課題
[0007] し力るに、従来の変位検出方法では、計算時間が膨大になるとともに、回転変位の 検出が容易に行えな!/、などの問題点があった。
[0008] 本発明は上記の点に鑑みてなされたもので、位相特異点の構造を特徴付ける要素 力 所望の位相特異点を特定することにより変位前後の位相特異点を確実に特定し 、その変位を容易、かつ、確実に、さらには、回転変位を含む種々の変位を検出でき る変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、特徴点マッチ ング処理方法、特徴点マッチングプログラムを提供することを目的とする。
課題を解決するための手段
[0009] 本発明は、変位前後の位相特異点力 変位を検出する変位検出方法であって、位 相特異点の位相構造力 所定の要素を取得し、要素に基づいて位相特異点を特定 し、変位後の位相特異点の位置を検出することを特徴とする。
[0010] このとき、取得される要素は、位相特異点近傍の実部と虚部の値を平面で補間し、 その結果、実部がゼロとなる線分と虚部がゼロとなる線分とが交差する角度、位相特 異点近傍の強度分布の形状の離心率、位相特異点近傍の強度の勾配であることを 特徴とする。
[0011] また、特定された位相特異点の変位前後の位置を要素に基づいて特定し、特定さ れた変位前後の位相特異点の位置に基づいて位相特異点の変位を測定することを 特徴とする。このとき、変位前後の位相特異点の位置の座標の差分の二乗平方根を 位相特異点の変位量として算出することを特徴とする。
[0012] 複数組の変位前後の位相特異点の位置を結ぶ線分の複数の垂直二等分線の距 離の和が最小となる点を回転中心として求め、変位前後の位相特異点の位置と回転 中心とからなる三角形の回転中心を挟む辺の角度を回転角として求めることを特徴と する。
[0013] さらに、本発明は、取得した画像の強度信号をフーリエ変換し、フーリエ変換された 信号の複素共役の対称性を崩すようなフィルタリング処理を行 ヽ、擬似位相情報を 取得し、取得した擬似位相情報に基づ ヽて位相特異点を取得することを特徴とする
[0014] このとき、フーリエ変換された信号に対して中心に位相特異点を持つらせん状の位 相特性を持ち、かつ、ドーナツ状の振幅特性を持つラゲールガウスビームフィルタ特 性によりフィルタリングすることを特徴とする。また、画像の強度信号の平均値を画像 の各強度信号から減算した画像を元画像とすることを特徴とする。
[0015] さらに、フーリエ変換された信号に対して中心に位相特異点を持つらせん状の位相 特性を持つフィルタ特性、又は、ヒルベルト変換によりフィルタリングを行い、擬似位 相情報を取得することを特徴とする。
[0016] また、本発明は、変位前に特定された特徴点と変位後に特定された特徴点とをマツ チングを行うポイントマッチング処理方法であって、変位前の特徴点の特徴パラメ一 タと変位後の特徴点の特徴パラメータとの差異に基づいて変位前の特徴点に対応す る変位後の特徴点を探しだし、その特徴点の地理的な位置関係の類似性力 変位 前の特徴点に対応する変位後の特徴点のマッチングを行うことを特徴とする。
[0017] また、このとき、本発明は、変位前の複数の特徴点間の特徴パラメータの総和距離 と変位後の複数の特徴点間の特徴パラメータの総和距離との差分とが最小となる特 徴点の組を変位前と変位後とで互いにマッチングする特徴点として特定することを特 徴とする。
[0018] さらに、本発明は、変位前の所定の特徴点の特徴パラメータを T 位後の所
il、 T、変
i2
定の特徴点の特徴パラメータを P P α
w(il)、 w(i2)、 を縮尺としたとき、
[0019] [数 1]
Figure imgf000005_0001
\≤i\,il≤n が最小となる特徴点の組を変位前後でマッチングする特徴点の組として抽出すること を特徴とする。
[0020] なお、特徴点は、画像の強度情報に基づ 、て取得される位相特異点であることを 特徴とし、位相特異点は特徴パラメータとして位置情報、離心率、及び、渦度、並び に、実軸と虚軸とが交わる角度を有することを特徴とする。
発明の効果
[0021] 本発明によれば、位相特異点の位相構造力 所定の要素を取得し、要素に基づい て位相特異点を特定することにより、変位前後の各位相特異点を一意に特定するこ とができ、これによつて、変位後の位相特異点を追跡できるなどの効果を得られる。 図面の簡単な説明
[0022] [図 1]本発明の一実施例のシステム構成図である。
[図 2]変位検出プログラムの処理フローチャートである。
[図 3]位相特異点取得処理の処理フローチャートである。
[図 4]位相特異点取得処理の動作説明図である。
[図 5]位相特異点取得処理の動作説明図である。
[図 6]位相特異点の位置をサブピクセルで決定する動作を説明するための図である。
[図 7]位相特異点の位置をサブピクセルで決定する動作を説明するための図である。
[図 8]位相特異点の位置をサブピクセルで決定する動作を説明するための図である。
[図 9]位相特異点特定処理の処理フローチャートである。
[図 10]離心率 e、渦度 ωを取得する動作を説明するための図である。
[図 11]離心率 e、渦度 ωを取得する動作を説明するための図である。
[図 12]離心率 e、渦度 ωを取得する動作を説明するための図である。
[図 13]全位相特異点の位置の差を求め X方向， y方向の変位をヒストグラムに表した 図である。
[図 14]回転変位の検出動作を説明するための図である。
[図 15]シミュレーション実験を説明するための図である。
[図 16]回転前後の位相特異点の数をヒストグラムで表した図である。
[図 17]回転前後の位相特異点の回転角 Θ iをヒストグラムで表した図である。 [図 18]回転前後の角度 Δ Θの差をヒストグラムで表した図である。
[図 19]回転前後の対応する位相特異点の離心率 eをヒストグラムで表した図である。
[図 20]回転前後の対応する位相特異点の渦度 ωの差をヒストグラムで表した図であ る。
[図 21]ポイントマッチング処理の動作説明図である。
[図 22]ポイントマッチング処理のフローチャートである。
[図 23]マッチングアルゴリズムの処理動作を説明するための図である。
[図 24]マッチングアルゴリズムの処理動作を説明するための図である。
[図 25]マッチングアルゴリズムの処理動作を説明するための図である。
[図 26]実験'検証に用いた LGフィルタの振幅特性図である。
[図 27]マッチング処理の動作説明図である。
[図 28]マッチング処理の動作説明図である。
[図 29]実験 ·検証結果を説明するための図である。
[図 30]実験 ·検証結果を説明するための図である。
[図 31]実験 ·検証結果を説明するための図である。
[図 32]実験 ·検証結果を説明するための図である。
[図 33]実験 ·検証結果を説明するための図である。
[図 34]実験 ·検証結果を説明するための図である。
[図 35]実験 ·検証結果を説明するための図である。
[図 36]点間総和距離の差によりマッチングを行う動作を説明するための図である。
[図 37]点間総和距離の差によりマッチングを行う動作を説明するための図である。 圆 38]相関法の動作を説明するため図である。
圆 39]相関法の動作を説明するため図である。
符号の説明
101 変位検出装置、 102 測定対象
111 画像検出部、 112 信号処理部
121 光源、 122 撮像装置
131 フレームメモ!;、 132 ^Μ , 133 フーリエ変 咅 134 ROM
135 ハードディスクドライブ、 136 ディスプレイ、 137 入力装置
発明を実施するための最良の形態
[0024] 〔システム構成〕
図 1は本発明の一実施例のシステム構成図を示す。
[0025] 本実施例の変位検出装置 101は、画像検出部 111及び信号処理部 112から構成 されている。
[0026] 画像検出部 111は、光源 121及び撮像装置 122から構成されている。画像検出部 111は、光源 121から測定対象 102に光を照射し、撮像装置 122により測定対象 10 2の表面画像を撮像する。
[0027] 画像検出部 111で撮像された画像は、信号処理部 112に供給される。信号処理部 112は、例えば、コンピュータシステムであり、フレームメモリ 131、処理部 132、フー リエ変換部 133、 ROM134、ハードディスクドライブ 135、ディスプレイ 136,入力装 置 137から構成されている。
[0028] フレームメモリ 131は、画像検出部 111からの画像をフレーム毎に記憶する。処理 部 132は、ハードディスクドライブ 135に予めインストールされた変位検出プログラム に基づいて測定対象 102の微小変位を検出するための処理を実行する。
[0029] フーリエ変換部 133は、画像検出部 111で撮像された画像に対してフーリエ変換を 行う。メモリ 134は、処理部 132の作業用記憶領域として用いられる。ハードディスク ドライブ 135は、変位検出プログラムや画像検出部 111で検出された画像データ、あ るいは、変位検出プログラムでの途中結果、変位検出結果などのデータが記憶され る。
[0030] ディスプレイ 136は、 LCD、 CRTなど力も構成されており、画像、変換画像、途中 結果、変位検出結果などの処理部 132での処理結果を表示する。入力装置 137は、 キーボード、マウスなど力 構成されており、データ入力や各種コマンドの入力するた めに用いられる。
[0031] 本実施例の変位検出装置 101は、インストールされた変位検出プログラムに基づい て測定対象 102から取得した画像を解析して、その変位を検出する。 [0032] 〔変位検出プログラム〕
図 2は変位検出プログラムの処理フローチャートを示す。
[0033] 処理部 132は、まず、変位検出プログラムに従って、ステップ S1— 1で変位前後の 対象画像をノヽードディスクドライブ 135からメモリ 134などに読み出す。
[0034] 次に処理部 132は、ステップ S1— 2で、読み出された変位前後の対象画像に対し て各々位相特異点を取得するための位相特異点取得処理を実行する。
[0035] 次に処理部 132は、ステップ S1— 3で変位前後の対象画像から取得した位相特異 点が持つ、実部と虚部とのゼロクロス角度、離心率、渦度などの要素により位相特異 点を特定する位相特異点特定処理を実行する。
[0036] 次に処理部 132は、ステップ S1— 4で、実部と虚部とのゼロクロス角度、離心率、渦 度などの位相特異点に特有の要素により変位前後で対応する位相特異点を特定し
、変位前後の位相特異点の位置力 対象画像の変位を検出する変位検出処理を実 行する。
[0037] 位相特異点によって、対象画像の変位を検出することで対象画像の変位を正確に 特定することが可能となる。
[0038] 〔複素解析信号取得処理〕
次にステップ S 1 - 2の位相特異点取得処理にっ 、て説明する。
[0039] 図 3は位相特異点取得処理の処理フローチャート、図 4、図 5は位相特異点取得処 理の動作説明図を示す。
[0040] ステップ S 1 2の位相特異点取得処理は、取得した画像の強度信号をフーリエ変 換し、フーリエ変換された信号の複素共役の対称性を崩すようなフィルタリング処理 を行い、擬似位相情報を取得し、取得した擬似位相情報に基づいて位相特異点を 取得する処理である。
[0041] 処理部 132は、まず、ステップ S2— 1で前処理を行う。ステップ S2— 1の前処理は 、対象画像の全体の強度情報の平均値を求め、求めた全体画像の強度情報の平均 値を対象画像の強度情報から減算する処理を行う処理であり、この処理によって、画 像データ強度情報中の直流成分が除去される。直流成分が除去されることにより、位 相特異点が強調される。 [0042] 次に、処理部 132は、ステップ S2— 2でフーリエ変換部 133を制御して、強度情報 に対して 2次元高速フーリエ変換を行 、、空間周波数スペクトルを得る。
[0043] 例えば、図 4にお 、て、画像検出部 111で得られた物体のテクスチャの強度情報 2 11の関数 201を
[0044] [数 2]
g0 ( ，
とする。この強度情報を 2次元フーリエ変換することにより、空間周波数スペクトル [0045] [数 3]
F {g0O， }
を求めることができる。なお、図 4において 212は空間周波数スペクトルの振幅分布、 213は位相分布を示している。空間周波数スペクトルは、実部のみをもつ強度情報を フーリエ変換したものなので、周波数原点に対して複素共役の対称性、すなわちェ ルミート対称性をもつ。 FFTを関数 Fとすると、例えば、式 202のように表せる。
[0046] そこで、処理部 132は、ステップ S2— 3で実部を基に虚部を作り出して複素解析信 号を得るために、図 4に 221〜225で示すような分布を持つ、周波数領域上の複素 共役の対称性を崩すようなフィルタリング処理を行う。フィルタリング処理では、例え ば、ヒルベルト変換、スパイラル位相フィルタ、 LG (ラゲールガウス）フィルタなどが用 いられる。なお、図 4において 221はヒルベルト変換の振幅分布、 222はスパイラル位 相フィルタの位相分布、 224は LGフィルタの振幅分布、 225は LGフィルタの位相分 布を示している。 [0047] ヒルベルト変換の演算式は、例えば、図 4、 203で表すことができる。また、スパイラ ル位相フィルタの演算式は、例えば、図 4、 204で表すことができる。 LGフィルタの振 幅透過分布 224の演算式は、図 4、 205 (A)で表すことができる。 LGフィルタは、中 心を対称に複素共役を位相上で πだけずらすことにより，ドーナツ状の振幅透過率 分布が得られる。なお、 LGフィルタの位相分布 225の演算式は、例えば、図 4、 205 (Α)に角度 φを導入することにより図 4、 205 (Β)に示すように表される。
[0048] 以上のように、演算式 202に演算式 203、 204、 205などを乗算することにより、本 実施例で必要となる解析信号を得ることができる。また、エラーの原因となる高周波 数成分を除去でき、直流成分を自動的、かつ、連続的にゼロにできる。さらに、ドー ナツ状の強度分布により低周波数成分を取り除くことにより、位相特異点の数を増や すことができる。
[0049] さらに、本実施例ではランダムなパターンの強度情報を用いて解析を行っており、 測定対象 102の構造などによってスペックという丸い斑点のようなパターンが多数存 在することになる。このスペックルは、本実施例で指標としている位相特異点に対応 関係がある。
[0050] このスペックルのサイズによりフーリエ変換したときの周波数領域に位相特異点の 情報を多く含む周波数が変わる。例えば、スペックルのサイズが小さいと位相特異点 の情報を含む周波数成分は高くなり、逆にサイズが大きいと位相特異点の情報を含 む周波数成分は低くなる。つまり、位相特異点の情報を高精度に求めるにはドーナ ッ状の強度分布をもつ LGフィルタが最適であり、ドーナツ状の幅のサイズをスペック ルのサイズに合わせて調節することにより、より位相特異点の情報を高精度に求める ことができる。一方、ローパスフィルタでは、低周波数成分を強調し、位相特異点の情 報をうまく取り出すことができな 、。
[0051] 以上のように、本実施例における解析信号の生成過程では LGフィルタを用いること が最適となる。
[0052] なお、図 4において 214はフィルタリング処理後の振幅分布、 215はフィルタリング 処理後の位相分布を示しており、その演算式は、例えば、図 4、 206で表すことがで きる。 [0053] なお、 221、 223に示すようなヒルベルト変換またはスパイラル位相フィルタをかけ ても同種の解析信号を得ることはできる力 ヒルベルト変換またはスパイラル位相フィ ルタを使うためには予め直流成分を除去しておかなければならない。
[0054] 処理部 132は、ステップ S2— 4でフィルタをかけた空間周波数スペクトルをフーリエ 逆変換して複素解析信号
[0055] [数 4]
g0O,
を得る。なお、図 4において 216はフーリエ逆変換後の振幅分布、 217はフーリエ逆 変換後の位相特性を示しており、その演算式は、例えば、図 4、 207で表すことがで きる。
[0056] 図 5 (A)はフーリエ逆変換後の振幅分布を拡大したものであり、図 5 (B)は図 5 (A) の所定の部分を更に拡大したものである。
[0057] この複素解析信号の擬似位相情報から図 5 (B)に示すような位相特異点 Pを取得 することができる。
[0058] 次に処理部 132は、ステップ S2— 5で位相特異点の位置を 1ピクセルより小さい分 解能であるサブピクセルの分解能で決定する処理を行う。
[0059] 図 6、図 7、図 8は位相特異点の位置をサブピクセルで決定する動作を説明するた めの図を示す。
[0060] 図 6 (A)は擬似位相情報、図 6 (B)は位相特異点周辺画素の拡大図、図 6 (C)は 実部、図 6 (D)は虚部を示している。また、図 7 (A)は位相特異点近傍の擬似位相情 報、図 7 (B)は位相特異点近傍の擬似位相情報の虚部、図 7 (C)は位相特異点近傍 の擬似位相情報の実部、図 7 (D)は位相特異点近傍の擬似位相情報を直線補間し たもの、図 7 (E)は位相特異点近傍の擬似位相情報の虚部を直線補間したもの、図 7 (F)は位相特異点近傍の擬似位相情報の実部を直線補間したものを示して ヽる。 [0061] 擬似位相情報は (- π , π )までの値をとり、図 6、図 7において黒から白になるにした がって値が大きくなる。
[0062] 位相特異点 Ρは、その点を囲む閉経路上の位相勾配の積分が ± 2 πの整数倍に なり、また、複素数の実部と虚部が共にゼロとなるという性質を持っている。
[0063] 図 6 (Β)、図 7 (C)に示すような互いに隣接しあう 4画素を結ぶ正方形閉経路に対し て位相勾配の積分が ± 2 πの整数倍になるという性質を用いて位相特異点の存在し ている位置を決定していた。このため、 1ピクセル以上の分解能をもって位相特異点 の位置を決定することができな 、。
[0064] 本実施例では、複素数の実部と虚部が共にゼロとなるという性質を用いて 1ピクセ ルの分解能より小さ 、サブピクセルの分解能で位相特異点の位置を決定するように している。本実施例では、図 7 (B)、図 7 (C)に示す特異点近傍の離散化された実部 と虚部を平面で線形補間することにより図 7 (E)、図 7 (F)に示すように実部及び虚部 を共に滑らかな構造として、実部がゼロとなる線分 Lre0、虚部がゼロとなる線分 LimO を求める。図 7 (D)に示すように位相特異点の位置は実部がゼロとなる線分 LreOと虚 部がゼロとなる線分 LimOの交わる点を位相特異点 P0として設定することによりサブピ クセルの分解能をもって位相特異点を決定することができる.
図 8 (A)に示すように位相勾配の積分が ± 2 πの整数倍になる性質を用いる方法 では、 1ピクセルの面積を持つ正方形の閉経路内のどこ力に位相特異点が存在する ことはわかってもその位置を特定することができな力つた。これに対して図 8 (Β)に示 すように複素数の実部と虚部が共にゼロとなると 、う性質を用いて実部と虚部の内挿 と実部と虚部のゼロの交差点から位相特異点を求めることにより 1ピクセルより小さい サブピクセルの分解能をもって位相特異点の位置を決定できている。
[0065] 〔位相特異点特定処理〕
図 9は位相特異点特定処理の処理フローチャートを示す。
[0066] ステップ S1— 3の位相特異点特定処理は、位相特異点の位相構造から所定の要 素を取得し、要素に基づ 、て位相特異点を特定する処理である。
[0067] 処理部 132は、ステップ S3— 1で位相特異点を特定する要素としてゼロクロス角を 取得する。ゼロクロス角は、図 7 (D)、図 8 (B)に示す実部がゼロとなる線分 LreOと虚 部がゼロとなる線分 LimOとの交差する角度 Δ Θである。
[0068] また、処理部 132は、ステップ S3— 2で位相特異点を特定する要素として離心率 e を取得し、ステップ S3— 3で渦度 ωを取得する。
[0069] 図 10、図 11、図 12は離心率 e、渦度 ωを取得する動作を説明するための図を示す
。図 10 (A)は擬似位相情報、図 10 (B)は擬似位相情報の位相特異点近傍の振幅 分布、図 10 (C)は位相分布を示して!/ヽる。
[0070] 図 10 (B)に示すように振幅分布は位相特異点近傍に近づくにつれて暗くなつてい て、これは値がゼロに近づ 、て 、ることを示して 、る。
[0071] また、図 11は特異点近傍の振幅分布の構造図を示して!/、る。図 11 (Α)は振幅を 等高線で表示したものを示しており、図 11 (Β)は三次元表示したものを示す。
[0072] 図 11 (Α)に示すように位相特異点近傍の等振幅の分布は略楕円形状をしている。
一般に位相特異点近傍の振幅分布は略楕円形をしており、その形状は位相特異点 に応じて異なっており、位相特異点を特定する要素として用いることができる。この楕 円形状の特徴を二次曲線の特徴を表す値である離心率によって示すことによって、 位相特異点を特定することが可能となる。
[0073] 例えば、図 12に示すような楕円の場合、一般に数式で、
[0074] [数 5]
X ヌ ，
a b
で表せる。図 12に示す楕円の離心率 eは a :長軸半径， b :短軸半径， f :焦点距離と すると、離心率 eは
[0075] [数 6] a a
ただし、 0 <離心率 eく 1
で表させる。
[0076] 図 11 (B)に示される振幅を三次元で表示した図を見ると位相特異点近傍の形状は 楕円形の円錐のような構造をしていることが分かる。したがって、本実施例では、各位 相特異点を特定する要素の一つとして位相特異点近傍の渦度を求める。渦度とは微 小部分の局所的な回転の強さと回転の方向を表す量である。図 11 (A)では振幅の 等高線の間隔が狭いほど渦度は大きくなり、図 11 (B)では三次元表示した位相特異 点近傍の勾配 (変化の度合い）が急になるほど渦度は大きくなる。すなわち、渦度は 位相特異点近傍の振幅の大きさに依存する。
[0077] 渦度 ωは次式で表される。
[0078] [数 7]
Figure imgf000015_0001
(r:実部、 i :虚部である。 )
以上のように位相特異点は位置 (x、 y) ,ゼロクロス角 Δ θ ,離心率 e,渦度 ωの要 素を持っていて、これらの要素を用いて各位相特異点を一意に特定することができる 。このとき、複数の要素を用いることにより確実に位相特異点を特定することができる
[0079] 〔直線変位検出〕
次にステップ S 1—4の位相特異点変位検出処理について説明する。
[0080] 位相変位検出処理では、変位前後の位相特異点を上記の特定方法により特定す ることにより変位量を検出する。
[0081] このとき、まず、変位前後で位相特異点をそれぞれ一意に特定するためのパラメ一 タとして前述のように求めたゼロクロス角 θ ,離心率 e,渦度 ωを用いて位相特異点を 特定する。その方法は変位前の各パラメータを 0 1、 el、 ω 1、変位後の各パラメータ を 0 2、 e2、 ω 2とする。
[0082] 上記パラメータ（Δ 0 1、 el、 ω 1)、（Δ 0 2、 e2、 ω 2)のそれぞれの差（Δ θ 2 - A θ 1、 e2— el、 ω 2— ω 1)の絶対値がある値より小さい時に同一の位相特異点である とみなす。各位相特異点をそれぞれ一意に特定できるため、変位前後で各位相特異 点がどこに移動したのか分かる。このとき、変位前の位相特異点の位置 (xl、 yl)と変 位後の位相特異点の位置 (x2、 y2)の差力 変位量を求めることができる。
[0083] 実際の変位量 A Lは、 Δ χ= (χ2— xl)、 A y= (y2— yl)とすると、
[0084] [数 8]
AL = ^Ax2 + Ay2
で表すことができる。
[0085] なお、図 13は、全位相特異点の位置の差を求め X方向， y方向の変位をヒストグラ ムに表した図を示す。図 13 (A)は X方向の変位、図 13 (B)は y方向の変位を示して いる。
[0086] 図 13に示すように全位相特異点の位置が X方向に略平行に変位したことがわかる
[0087] 〔回転変位検出〕
また、本実施例では、各位相特異点を回転変位前後で特定できるため、回転の検 出も可能である。
[0088] 図 14は回転変位の検出動作を説明するための図を示す。
[0089] 図 14 (A)に示すように回転前の 2組の位相特異点 Pl l、 P12の xy座標をそれぞれ (xl、yl)、 (x2、y2)
とし，回転後の対応する位相特異点 P21、 P22の xy座標をそれぞれ
(x'l、y'l)、（x'2、y'2)
とすると、回転の中心
(xc、 ycノ
は図 14 (A)に示すように回転前後の位相特異点 Pl l、 P12と位相特異点 P21、 P22の 垂直二等分線
[0090] [数 9]
> 'c 2
の交わる点となる。
[0091] 一般的に回転の中心を求めるには 2組の回転前後の位相特異点が分かれば十分 である力 実際は多数の位相特異点が存在するので回転前後の各位相特異点から 得られる垂直二等分線が図 14 (B)に示すに必ずしも 1点では交わらない。このような 場合には、回転の中心 (xc、 yc)と各垂直二等分線との距離の和が最小となるような 回転の中心 (xc、 yc)を最小二乗法により求めるようにする。
[0092] 以上により回転中心を決定できる。
[0093] 次に図 14 (C)に示すように回転前後の各位相特異点 Pl l、 P12と回転中心 (xc、 yc )力 なる三角形に対して余弦定理を用いることによって、回転角 Θを求めることがで きる。
[0094] 以上により、本実施例によれば、回転変位の検出を行うことが可能となる。
[0095] 〔シミュレーション実験結果〕
実験で得た強度情報をコンピュータ内部で 90度だけ回転させて，回転させる前後 の強度情報を用いてシミュレーション実験を行った。
[0096] 図 15はシミュレーション実験を説明するための図を示す。なお、図 15は、画素数が 1024 X 1024で、回転させる前後の強度情報力も LGフィルタを用いて 2つの解析信 号を得たときの振幅分布と位相分布を表示したものである。図 15 (A)は回転変位前 の擬似位相情報、図 15 (B)は回転変位前の擬似位相情報の振幅分布、図 15 (C) は回転変位前の擬似位相情報の位相分布、図 15 (D)は回転変位後の擬似位相情 報、図 15 (E)は回転変位後の擬似位相情報の振幅分布、図 15 (F)は回転変位後 の擬似位相情報の位相分布を示して 、る。
[0097] この回転前後の位相特異点の数は両方とも正の位相特異点が 322個，負の位相区 異点 321個の合計 643個で同数の位相特異点が存在していた。図 15 (B)及び図 15 ( E)に示される振幅分布からも明らかなように回転前の三角で囲った 2つのリング状の 振幅が回転後には回転の中心 (全画素の中心）に対して 90° 回転している。
[0098] 次に位相特異点を特定するための要素である Re=0と lm=0の直線の交わる角度： Δ
θ ,離心率 (eccentricity)： e,渦度 (vorticity)： ωを回転前後の全位相特異点に対し て求める。回転前後で対応する位相特異点を探す条件は以下の条件を用いた。
[0099] 条件 l : Re=0と lm=0の直線の交わる角度 Δ 0の差が ±0.5度以下
条件 2 :回転前後の離心率 eの差が ±0.005以下
条件 3：回転前後の渦度 ωの差が ±0.05以下
回転前の位相特異点に対応する回転後の位相特異点が 1対 1の対応になっている のかを調べた。
[0100] 図 16は回転前の位相特異点に対応する回転後の位相特異点の数をヒストグラムで 表した図を示す。
[0101] 回転前後でほぼ全ての位相特異点が一意に特定できた。
シミュレーションで得られた測定結果の回転の中心は
(xcゝ yc) = (511.5、 511.5)
となった。これは設定した回転の中心
(xcゝ yc) = (511.5、 511.5)
と一致する。
[0102] 図 17は回転前の位相特異点に対応する回転後の位相特異点と回転の中心力も求 めたそれぞれの回転角 Θ iをヒストグラムで表した図を示す。 [0103] シミュレーション結果である回転角のヒストグラムは ± 0.001度の精度で全位相特異 点が 90度回転していることを示している。シミュレーション実験では、測定対象 102を 90度回転させており、図 17に示す測定結果も 90度回転していることを示しているの で、この測定結果力も本発明が回転に対しても有効に計測できることを確認できる。
[0104] 図 18は回転前後の対応する位相特異点の角度 Δ 0の差をヒストグラムで表した図 、図 19は回転前後の対応する位相特異点の離心率 eをヒストグラムで表した図、図 2 0は回転前後の対応する位相特異点の渦度 ωの差をヒストグラムで表した図を示す。
[0105] 図 18、図 19、図 20に示すようにゼロ付近に固まっており、誤差が少なぐ位相特異 点を確実に特定できて 、ることがわ力る。
[0106] なお、上記実施例では、位相特異点を用いて回転変位を検出する処理につ！、て 説明したが縮尺、移動、回転を簡単な処理により検出することも可能である。
[0107] 〔ポイントマッチング処理〕
ポイントマッチング処理は、任意の複数点の位置関係と、その点の空間構造の類似 性を用いて変位前のパターン力 変位後のパターンを検出する処理である。
[0108] すなわち、位相構造が類似する点を候補点として、その候補点を含んだパターン 内の点間の距離が最小となる点を変位後のパターンとして決定するというものである
[0109] 以下にポイントマッチング処理の手順を説明する。
[0110] まず、変位前後の 2つの強度分布画像に対して各々前述の位相特異点を抽出する ためのフィルタ演算を行うことにより位相特異点を抽出する。抽出した位相特異点近 傍にある実部及び虚部を線形補間することによって、サブピクセルの精度で全位相 特異点の位置を特定する。更に、位置が特定された各位相特異点に対して離心率 e 、実軸と虚軸との角 Δ Θ、渦度 ωを求め、ラベリングを行う。
[0111] 図 21は、ポイントマッチング処理の動作説明図を示す。図 21 (Α)は変位前の位相 特異点の状態、図 21 (Β)は変位後の位相特異点の状態を示している。
[0112] まず、図 21 (Α)に示す位相特異点力も黒丸で示す位相特異点を選択する。選択 する位相特異点はその特徴パラメータである離心率 e、実軸と虚軸との角 Δ Θ、渦度 ωが他の多くの位相特異点とは明らかにことなるものを抽出することが望ましい。これ によって、変位後の位相特異点を特定しやすくなる。
[0113] 図 21 (A)に示す変位前の画像力も選択した位相特異点の特徴パラメータと位置と 図 21 (B)に示す変位後の画像力も選択した位相特異点の特徴パラメータと位置とか ら図 21 (B)に黒丸で示すように図 21 (A)に黒丸で示す位相特異点に対応する位相 特異点を検出する。
[0114] 位相特異点の特徴パラメータ、離心率 e、渦度 ω、実軸と虚軸との角 Δ Θは、変位 前後で同じであるか、または、ほとんど変化しない。このため、変位前後の位相特異 点を正確に特定することによって移動、回転、縮尺などを正確に測定することができ る。
[0115] 図 22はポイントマッチング処理のフローチャートを示す。
[0116] まず、処理部 132は、ステップ S4— 1で、図 3とともに説明した位相特異点取得処 理によって変位前後の位相特異点を取得し、ラベリングする。
[0117] 次に、処理部 132は、ステップ S4— 2で変位前の各位相特異点にほぼ同じ、あるい は、所定の範囲内に近似した特徴パラメータを持つ位相特異点を変位後の位相特 異点から探索する。
[0118] 次に処理部 132は、ステップ S4— 3で変位前の位相特異点力も所定の 2点を選択 し、それを組み合わせて組み合わせ Μとする。次に、処理部 132は、ステップ S4—4 で変数 kに「3」が代入される。
[0119] 次に、処理部 132は、ステップ S4— 5で変数 kが定数 nより小さいか否かを判定する
。なお、定数 nは、探索する位相特異点の数である。
[0120] 次に処理部 132は、ステップ S4— 6で選択された複数の位相特異点から次の探索 対象となる位相特異点を抽出する。
[0121] 次に処理部 132は、ステップ S4— 7で新たに抽出された位相特異点に対して特徴 ノ メータが近似する変位後の位相特異点を探索し、組み合わせ Mとして追加する
[0122] 次に処理部 132は、ステップ S4— 8で組み合わせ Mの位相特異点の特徴パラメ一 タに対して下記の式（1)によって計算を行い、閾値以下のものを変位前の位相特異 点と変位後の位相特異点との有力な組み合わせの候補として残す。 [0123] [数 10]
min = ^ (I Tn一 Ta \ -a | Pw( ) - Pw{n) |)2 * · · (丄 )
l≤il,;2≤n
式（1)にお 、て tと は変位前と変位後の点を表して 、る。なお、 tは と対応して w(i) w(i) いる。
T= {t , t , t , - - - .t }
1 2 3 n
は変位前の位相特異点の組み合わせを表して 、る。そして、
Ρ= {ρ , ρ , ρ , · ' · ,ρ }
1 2 3 η
は変位後の位相特異点の組み合わせを表している。
[0124] また、 I t t I は、点 tと点 tとの間の距離， αは縮尺のスケール (度合い)を表して il i2 il i2
いる。さら〖こ、変数 wは Τ→Ρの対応関係を表している。
[0125] なお、上記式（1)の関数をコスト関数（cost fonction)と呼ぶ。
[0126] 処理部 132は、ステップ S4— 9で変数 kに（k+ 1)を代入し、ステップ S4— 5に戻る
。処理部 132は、上記ステップ S4— 6〜S4— 9をステップ S4— 5で変数 kが定数 nよ り大きくなるまで繰り返す。
[0127] 次に処理部 132は、ステップ S4— 10で残った組み合わせ Mの位相特異点の特徴 ノ ラメータ間の距離が最小となるものをマッチングパターンとして抽出する。
[0128] 次に、処理部 132は、ステップ S4— 11で抽出した特徴パラメータ間の距離が最小 となる位相特異点の変換パラメータ、例えば、平行移動 (shift),回転 (rotation),縮尺 (s cale)を計算する。
[0129] 上記ポイントマッチング処理により、 目標はパターン内の複数の点の位置関係にお いて変位前後で最も良く一致するものを見つけ出すことができる。
[0130] 次に、上記マッチングアルゴリズムの処理動作を、図面を用いてさらに詳細に説明 する。
[0131] 図 23、図 24、図 25はマッチングアルゴリズムの処理動作を説明するための図を示 す。
[0132] 図 23 (A)は変位前の位相特異点の中力も複数の位相特異点の組を選んだもので あり、これを T空間とし、 ti(i=l、 ...,η),Τί(Χί,Υί,Θί,ωί,Δ θί)で表す。図 23(B)は 変位後の位相特異点で Ρ空間とし、 pj(j =1、 ...,k)、 Pj(Xj,Yj,ej, ω j, Δ θ j)で表され る。する。
[0133] まず、図 23 (A)に黒丸で示すように、変位前の位相特異点から任意の個数の位相 特異点を選ぶ。
[0134] 変位前のパターン内のそれぞれの位相特異点 tの特徴パラメータである離心率 e、 ゼロクロス角 Θ、渦度 ωと変位後の位相特異点 pjの特徴パラメータである離心率 e、 ゼロクロス角 Θ、渦度 ωとの差を閾値と比較する。それは、例えば、次式で表せる。
[0135] [数 11]
： {\ eTj― ePj \< etk,\\ ωτ> - Pj \i\ ωτ< + Pj |< ω,,,\ ¾ \≤0lh) . . . (2)
なお、式（2)で 2番目の分母の I ω +ω |は正規ィ匕するために付けカ卩えたもの
Ti Pj
である。実際のプログラム内ではこの形で使われている。また、 e は離心率の閾値、
th
Θ はゼロクロス角の閾値、 ω は渦度の閾値を示している。
th th
[0136] 次に T空間の位相特異点の組と P空間の位相特異点の組とのマッチングを行う。
[0137] 図 24において、実線で示す丸印は正の位相特異点を示しており、破線で示した丸 印は負の位相特異点を点線で示して 、る。
[0138] まず、図 24 (A)に示す T空間の 2個の位相特異点の極性、正か負力符号と、その 位置関係とに基づいて図 24(B)〖こ破線の直線で示すように P空間に存在しうる可能 性のある 2個の位相特異点の組を探索する。
[0139] このマッチングでは移動 (shift),回転 (rotation),だけでなく縮尺 (scale)も考慮してい るので、 T空間での初めの 2つの組を P空間の中から 2つを選ぶ組み合わせは、この 例では正の位相特異点と負の位相特異点の組み合わせ全てであり、正が 7個、負が 6個なので 42通り存在する。
[0140] なお、ここで、正の位相特異点と負の位相特異点というのは特徴パラメータの渦度（ vorticity)の値が正のものを正の位相特異点、また、負のものを"負の位相特異点とし ている。
[0141] この得られた位相特異点の組み合わせに対して，位相特異点の特徴パラメータの 条件に合わない組み合わせを切り捨てる。このとき、各特徴パラメータは実験環境に より変化しやすいことを考え、条件を緩くとることが好ましい。ここで言う条件とは、式（ 2)の閾値 e 、 0 、 ω のことである。
th th th
[0142] 次に、ステップ S4— 3で空間 T内の 2つの位相特異点 t、tを選ぶ。この位相特異
il i2
点 t、tは少なくとも 1つのマッチング候補点を持っている。位相特異点 tは n個のマ il i2 il 1 ツチング候補の点を、位相特異点 tは n個のマッチング候補点を持つ。したがって、
Ϊ2 2
変位後の空間 P内の候補点の組み合わせは n n通りあることがわ力る。それをテンポ
1 2
ラノレマッチング (temporal matching)と呼ぶこと〖こする。
[0143] 次に、図 25のように T空間の位相特異点を 1つ増やし、 P空間の中から 3点目として 可能性のある点を探す。これにより 3点目を含む組み合わせが得られる。なお、特徴 ノ メータを用いて明らかに間違った点は除外される。このような操作を T空間で選 択される位相特異点が n個になるまで繰り返す。
[0144] つまり、新しいマッチングポイントをカ卩えつつ、上記テンポラルマッチングを最終的 にすベての位相特異点のマッチング点が決定されるまで繰り返す。
[0145] この操作を行うためには、新しく追加されていく点の地理的な位置関係が変位前と 変位後で類似性があるかどうか調べる必要がある。
[0146] 式（1)では、テンポラルマッチングの類似性を調べると同時に距離の閾値は間違つ て検出されたマッチングポイントを除外する操作も行っていることになる。
[0147] 最終的にで得られた位相特異点の組み合わせ Mの中力 式（1)の評価式により T 空間の全位相特異点間の距離と対応する P空間の位相特異点間の空間的な距離及 び特徴点の距離の差の総和が最小になるものを選び出してマッチングを行っている
[0148] さらに、図面を用いて説明を続ける。 [0149] 図 26、図 27はマッチング処理の動作説明図を示す。図 26 (A)〜（C)、図 27 (A) 〜（C)は変位前の位相特異点の分布、すなわち、 T空間における位相特異点の分 布、図 26 (D)〜(F)、図 27 (D)〜(F)は変位後の位相特異点の分布、すなわち、 P 空間における位相特異点の分布を示す。
[0150] まず、図 26 (A)、 (D)に示すように T空間及び P空間のそれぞれの位相特異点の 特徴パラメータを計算し、注目ポイントを選択する。図 26 (A)では、 tl、 t2、 t4の組み 合わせを選択した。
[0151] 次に、図 26 (B)、（E)に示すように T空間内のそれぞれの点の特徴パラメータと似 た P空間の候補点 P1〜P4を探索する。
[0152] 次に、図 26 (C)、 (E)に示すように T空間内の 2つの位相特異点 tl、 t2を選択して、
P空間内の対応する 2つの位相特異点 pl、 p2を選択し、組み合わせ Mとする。
[0153] 次に、図 27 (A)、（D)に示すように既に選択した位相特異点以外で候補となる位 相特異点 t3、 t4及びそれに対応する P空間内の位相特異点 p3、 p4を探索する。
[0154] 次に、図 27 (B)、 (F)に示すように T空間の位相特異点 tl、 t2、 t4と P空間の位相 特異点 pl、 p2、 p4の特徴パラメータを式（1)に代入して、式（1)の値を算出する。同 様に、 T空間の注目位相特異点 tl、 t2、 t4と P空間の位相特異点 pl、 p2、 p3の特徴 ノラメータを式（1)に代入して、式（1)の値を算出する。
[0155] 算出結果のうち閾値より小さい位相特異点の組を抽出する。例えば、ここでは、組 み合わせ Ml (pl、 p2、 p4)及び M2 (pl、 p2、 p3)が抽出されたとする。抽出された位 相特異点の組 M1、M2をマッチングした組み合わせとして抽出する。
[0156] 次に図 27 (C)、 (G)に示すようにマッチングした組 Ml、 M2のうち最もコスト、すな わち、式（1)の値が低い組 Mlをマッチングした組として抽出する。
[0157] 以上は、説明を簡単にするために注目する位相特異点の数を 3点とした場合につ いて説明したが、 3点以上の位相特異点に注目して同様な処理を行うことも可能であ る。
[0158] 次に、ステップ S4— 11で算出される変換パラメータについて説明する。すなわち、 位相特異点の移動が平行移動、回転、縮尺の 3つのパラメータに分解できることにつ いて説明する。 [0159] 一般に、流体力学の観点力 どんな複雑な動きもそれを平行移動 (transition),回 転 (rotation),縮尺 (scale)で表すことができることが知られており、変位前後座標の回 転関係式は次式で表すことができる。
[0160] [数 12]
X,— xt = cos0(xi -x) +な sin ( ,.ー ,·)，
,. ~y. =-a sin 9{xi ~x) + cos 9(yi - ,),
ここで、 （x;、 y.)は変位前の位相特異点の座標、 (x\ y.')は変位後の位相特異点の 座標、 αは縮尺、 Θは回転角である。
[0161] [数 13]
x =∑xi !N,y =
i i
は変位前の位相特異点の重心、
[0162] [数 14]
Figure imgf000025_0001
は変位後の位相特異点の重心を表して 、る。
まず、はじめに、マッチングを表すのに適した関数について説明する。
[0163] この関数はそれぞれ変位前後で対応している位相特異点の差を足し合わせた形で 書ける。 [0164] [数 15]
E
Figure imgf000026_0001
- y)f }
ここで、 C= acos 0、 S = asin 0を表している。 Eの値の最小値を見つけるために はそれぞれのパラメータをまず計算しなければならな 、。例えば、
[0165] [数 16]
dE/dC = ,dE/dS = 0
などである。
[0166] ここで、
[0167] [数 17]
Figure imgf000026_0002
とお <と、
[0168] [数 18] 8E/dC =
力 次式が得られる。
[0169] [数 19]
(X; - CXi - ; )( ) +∑(y1+Sxi -Cyl.)(yl) = 0. i /
これを解くと次のようになる。
[0170] [数 20]
Figure imgf000027_0001
同様にして
[0171] [数 21]
dE/dS二 Q
力 次式が得られる。
[0172] [数 22] ∑ (x; - x,― 5^ ) (- y i ) +∑ (y； + Sxt - Cyt )(x, ) = 0 i i
これを解くと次式が得られる。
[0173] [数 23]
Figure imgf000028_0001
Sと Cの定義より以下のパラメータが得られる。
[0174] [数 24]
a = s2+C2,
<9 = arctan(5/C);
Ax = x - ,
y = ~y-
このように、位相特異点の移動は、平行移動 Δχ、 Ay、回転 0、縮尺 αの 3つのパ ラメータで表せることがわかる。この 3つのパラメータを変換パラメータと呼び、ステツ プ S4— 11で算出される。
[0175] 〔実験'検証〕
次に、上記マッチング処理方法の実験'検証結果について説明する。
[0176] 図 28は実験.検証に用いられた LGフィルタの振幅特性図、図 29〜図 35は実験. 検証結果を説明するための図を示す。図 29は実験画像を示す図であり、図 29 (A) は変位前の画像、図 29 (B)は変位後の画像を示している。図 30は図 29に示す画像 力も抽出された位相特異点を示す画像であり、図 29 (A)は変位前の位相特異点の 分布、図 29 (B)は変位後の位相特異点の分布を示している。図 31は選択した位相 特異点を拡大した図であり、図 31 (A)は変位前、図 31 (B)は変位後の拡大図を示し ている。
[0177] まず、今回の実験では、対象としてみどりふぐを使用し、高速 CCDカメラを用いて 秒間 30枚、連続撮影した画像を用いてマッチング処理を検証して 、る。
[0178] 撮影した画像に対して LGフィルタ演算を行 ヽ、位相特異点を得る。このとき、 LGフ ィルタのフィルタ幅は、 ω 2 = 30としている。なお、 LGフィルタは、 spiral phase filterに バンドパスフィルタの特性をカ卩えたものであり、図 28に示すような振幅特性とした。な お、 LGフィルタの振幅幅 ωを広くとれば、出てくる位相特異点の数は多くなり、逆に 狭くすれば少なくなる。
[0179] 図 29に示す変位前後の画像に対して今回提案したマッチングアルゴリズムを適用 した例にっ 、て説明する。図 29に示す画像の変位前後の位相特異点の分布は図 3 0に示すようになる。図 30 (A)に示す変位前の位相特異点力 黒丸で示される 6つ の位相特異点を選択してマッチングを行った。
[0180] 図 30 (Β)にお 、て、黒丸で示す位相特異点は変位前の黒丸で示す位相特異点に 対応しており、実際にマッチングアルゴリズムにより検出されたパターンである。
[0181] 図 31は変位前に選んだ位相特異点の組み合わせと、ポイントマッチングァルゴリズ ムによって検出された変位後の位相特異点の組み合わせを拡大図示したものである 。このように変位前，変位後のパターンを見比べると移動，回転，縮尺を組み合わせ た変形が起こっているがマッチングアルゴリズムによって検証された。
[0182] 異なる 2枚の画像に対して異なったパターンに着目して、本実施例のポイントマッチ ング処理を採用した場合の結果にっ 、ても検証した。
[0183] 図 32は異なる 2枚の画像に対して異なったパターンに着目して、本実施例のポイン トマッチング処理を採用した場合の変位前後の位相特異点の分布を示す図であり、 図 32 (A)は変位前の分布、図 32 (B)は変位後の分布を示している。図 33は図 32に 示した位相特異点の変位を示す図であり、図 33 (B)は選択された位相特異点を拡 大した図である。
[0184] この検証例にぉ 、ても図 32 (A)に示すように変位前の 4つの位相特異点が変位後 、図 32 (B)に示す位置に移動し、図 33に示すように注目した 4つの位相特異点が変 位して 、ることが検証できた。
[0185] さらに同じ作業をふぐの連続撮影画像 100枚に対して行った。
[0186] 図 34は連続した画像について実験'検証例を行った場合の画像を示す図である。
図 34 (A)は最初の画像、図 34 (B)は最後の画像を示している。図 35は図 34に示す 画像から 6つの位相特異点に注目し、その軌跡を追跡したものである。
[0187] このように、本実施例のポイントマッチング処理を用いることにより、位相特異点を連 続して追跡できる。
[0188] さらに、任意の位相特異点の組み合わせに対して位相特異点が持つ位相構造と他 の位相特異点との距離とによって対応する位相特異点の組を検出し、それらの位相 特異点の組を追跡することによってランダムな動きを追跡できることがわかる。特徴パ ラメータが変化してしまう理由で計測不可能であった大きな変位を伴った物体の運動 変位計測が可能となった。
[0189] また、本実施例では、位相特異点の組のうち式（1)が最小となる組み合わせを検出 することにより、マッチングを行ったが、位相特異点の点間総和距離の差を算出し、 変位前のパターンと変位後のパターンとで点間総和距離の差が最小となる位相特異 点を検出することによりマッチングを行うようにしてもよい。
[0190] 図 36、図 37は点間総和距離の差によりマッチングを行う動作を説明するための図 を示す。
[0191] 図 36 (A)、図 37 (A)は変位前の位相特異点の分布、図 36 (B)、図 37 (A)は変位 後の位相特異点の分布を示す。
[0192] まず、変位前の画像から位相特異点 p〜pを抽出するとともに、変位後の画像から
1 4
位相特異点 P '、 p '、 p '、 p
1 2 4 5 'を抽出し、変位前の位相特異点 p〜p
1 4tと変位後の位相 特異点 P '、 p '、 p '、 p 'との点間総和距離の差 Sを求める。点間総和距離の差 Sは S = ( I P P I I P 'P ' I )2+ ( I P P I I P 'P ' I )2
1 1 3 1 5 1 2 1 2
+ ( I P P I I P 'P ' I )2+ ( I P P I I P 'P ' I )2
2 3 2 5 4 3 4 5
+ ( I P P I I P 'P ' I )2+ ( I P P I I P 'P ' I )2
1 4 1 4 2 4 2 4
で求められる。
[0193] また、変位後の画像から位相特異点 p '、 p '、 p '、 p 'とは異なる位相特異点 p "、 p
1 2 4 5 1 2
"、 P "、 P
4 5〃を抽出し、変位前の位相特異点 P〜P
1 4tと変位後の位相特異点 P "、 P "、 P
1 2
"、 p "との点間総和距離の差 Sを求める。点間総和距離の差 Sは、
4 5 2 2
S = ( I P P I Ι Ρ "Ρ )2+ ( Ι Ρ Ρ I Ι Ρ 'Ρ )2
2 1 4 1 4 1 2 1 2
+ ( Ι Ρ Ρ I Ι Ρ "Ρ )2+ ( Ι Ρ Ρ I Ι Ρ 'Ρ )2
2 3 2 3 3 4 3 4
+ ( Ι Ρ Ρ I Ι Ρ "Ρ )2+ ( Ι Ρ Ρ I Ι Ρ 'Ρ )2
1 3 1 3 2 4 2 4
で求められる。
[0194] このとき、例えば、 Sく Sであれば、変位前の画像力 位相特異点 ρ〜ρは、位相
2 1 1 4 特異点 Ρ "、 Ρ "、 Ρ "、 Ρ
1 2 4 5 "に対応すると判定する。
[0195] なお、ここでは、説明を簡単にするために 4点について説明した力 2点以上であれ ば、点間総和距離の差によりマッチングを行うことが可能であることは言うまでもな 、。
[0196] 本実施例のポイントマッチング処理は、複数の位相特異点の位置関係を用いるの で縮尺、移動そして、回転により変位前後で完全にその位置関係が一致しなくてもポ イントマッチングを行うことで一番近い組を選び出すことができるので、本発明はあら ゆる変位を計測することができる。
[0197] なお、本実施例では、特徴点として位置情報、離心率、渦度、実軸と虚軸との交叉 角度であるゼロクロス角を特徴パラメータとして持つ位相特異点を例に説明したが、 特徴点はこれに限定されるものではなぐ例えば、位置情報だけであってもよい。
[0198] 以上、本実施例によれば、干渉計なしで空間的なランダムパターンの強度情報を 位相情報に対応づけ、位相情報の中に存在する全位相特異点の特徴を基に変位 前後で各位相特異点がどこに移動したか特定することにより、計測物体の全体の変 位計測に対しても、物体の局所領域における変位計測のいずれに対しても適用する ことができ、さらに場の回転の計測をも実現した変位検出装置を提供することできる。
[0199] 本国際出願は、 2006年 2月 1日に出願した日本国特許出願 2006— 024846号、 及び、 2006年 11月 8日に出願した日本国特許出願 2006— 303238号に基づく優 先権を主張するものであり、 日本国特許出願 2006— 024846号、及び、 日本国特 許出願 2006— 303238号の全内容を本国際出願に援用する。

Claims

請求の範囲
[1] 変位前後の位相特異点力 変位を検出する変位検出方法であって、
前記位相特異点の位相構造力 所定の要素を取得し、前記要素に基づいて前記 位相特異点を特定する位相特異点特定手順と、
前記位相特異点特定手順で取得された要素に基づいて変位後の位相特異点の位 置を検出する変位検出手順とを有することを特徴とする変位検出方法。
[2] 前記位相特異点特定手順で取得される要素は、前記位相特異点近傍の実部と虚部 の値を平面で補間し、その結果、実部がゼロとなる線分と虚部がゼロとなる線分とが 交差する角度であることを特徴とする請求項 1記載の変位検出方法。
[3] 前記位相特異点特定手順で取得される要素は、前記位相特異点近傍の強度分布 の形状の離心率であることを特徴とする請求項 1記載の変位検出方法。
[4] 前記位相特異点特定手順で取得される要素は、前記位相特異点近傍の強度の勾 配であることを特徴とする請求項 1記載の変位検出方法。
[5] 前記位相特異点特定手順は、前記位相特異点近傍の実部と虚部の値を平面で補 間し、その結果、実部がゼロとなる線分と虚部がゼロとなる線分とが交差する角度と、 前記位相特異点近傍の強度分布の形状の離心率と、前記位相特異点近傍の強度 の勾配である渦度とを含むことを特徴とする請求項 1記載の変位検出方法。
[6] 前記変位検出手順は、前記位相特異点特定手順で特定された位相特異点の変位 前後の位置を前記要素に基づいて特定する変位前後位置検出手順と、
前記変位前後位置検出手順で検出された変位前後の位相特異点の位置に基づ いて前記位相特異点の変位を測定する変位測定手順とを有することを特徴とする請 求項 1記載の変位検出方法。
[7] 前記変位測定手順は、前記変位前後位置検出手順で検出された変位前後の位相 特異点の位置の座標の差分の二乗平方根を位相特異点の変位量として算出するこ とを特徴とする請求項 6記載の変位検出方法。
[8] 前記変位測定手順は、前記変位前後位置検出手順で検出された複数組の変位前 後の位相特異点の位置を結ぶ線分の複数の垂直二等分線の距離の和が最小となる 点を回転中心として求める回転中心取得手順と、 前記変位前後位置検出手順で検出された変位前後の位相特異点の位置と前記回 転中心取得手順で取得した前記回転中心とからなる三角形の前記回転中心を挟む 辺の角度を回転角として求める回転角度取得手順とを有することを特徴とする請求 項 6記載の変位検出方法。
[9] 取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、
前記フーリエ変換手順でフーリエ変換された信号の複素共役の対称性を崩すよう なフィルタリング処理を行!、、擬似位相情報を取得するフィルタリング処理手順と、 前記フィルタリング処理手順で取得した擬似位相情報に基づいて位相特異点を取 得する位相特異点取得手順を有することを特徴とする請求項 1記載の変位検出方法
[10] 前記フィルタリング処理手順は、前記フーリエ変換手順でフーリエ変換された信号に 対して中心に位相特異点を持つらせん状の位相特性を持ち、かつ、ドーナツ状の振 幅特性を持つラゲールガウスビームフィルタ特性により前記擬似位相情報を取得す ることを特徴とする請求項 9記載の変位検出方法。
[11] 前記画像の強度信号の平均値を前記画像の各強度信号から減算した画像を前記画 像とする前処理手順を有することを特徴とする請求項 9記載の変位検出方法。
[12] 前記フィルタリング処理手順は、前記フーリエ変換手順でフーリエ変換された信号に 対して中心に位相特異点を持つらせん状の位相特性を持つフィルタ特性、又は、ヒ ルペルト変換により前記擬似位相情報を取得することを特徴とする請求項 9記載の変 位検出方法。
[13] 変位前の特徴点の特徴パラメータと変位後の特徴点の特徴パラメータとの差異に基 づいて変位前の特徴点に対応する変位後の特徴点を探しだし、その特徴点の地理 的な位置関係の類似性力 変位前の特徴点に対応する変位後の特徴点のマツチン グを行うマッチング手順を、更に、有することを特徴とする変位検出方法。
[14] 前記マッチング手順は、変位前の複数の位相特異点間の特徴パラメータの総和距 離と変位後の複数の位相特異点間の特徴パラメータの総和距離との差分とが最小と なる特徴点の組を変位前と変位後とで互いにマッチングする特徴点として特定するこ とを特徴とする請求項 13記載の変位検出方法。
[15] 前記マッチング手順は、変位前の所定の位相特異点の特徴パラメータを T、 T、変 位後の所定の位相特異点の特徴パラメータを P 、 P 、 αを縮尺としたとき、
[数 1]
Figure imgf000035_0001
l≤il,i2≤n
が最小となる位相特異点の組を変位前後でマッチングする位相特異点の組として抽 出することを特徴とする請求項 13記載の変位検出方法。
[16] 変位前後の位相特異点力 変位を検出する変位検出装置であって、
前記位相特異点の位相構造力 所定の要素を取得し、前記要素に基づいて前記 位相特異点を特定する位相特異点特定手段と、
前記位相特異点特定手段で取得された要素に基づいて変位後の位相特異点の位 置を検出する変位検出手段とを有することを特徴とする変位検出装置。
[17] 前記位相特異点特定手段で取得される要素は、前記位相特異点近傍の実部と虚部 の値を平面で補間し、その結果、実部がゼロとなる線分と虚部がゼロとなる線分とが 交差する角度であることを特徴とする請求項 16記載の変位検出装置。
[18] 前記位相特異点特定手段で取得される要素は、前記位相特異点近傍の強度分布 の形状の離心率であることを特徴とする請求項 16記載の変位検出装置。
[19] 前記位相特異点特定手段で取得される要素は、前記位相特異点近傍の強度の勾 配であることを特徴とする請求項 16記載の変位検出装置。
[20] 前記位相特異点特定手段は、前記位相特異点近傍の実部と虚部の値を平面で補 間し、その結果、実部がゼロとなる線分と虚部がゼロとなる線分とが交差する角度と、 前記位相特異点近傍の強度分布の形状の離心率と、前記位相特異点近傍の強度 の勾配である渦度とを含むことを特徴とする請求項 16記載の変位検出装置。
[21] 前記変位検出手段は、前記位相特異点特定手段で特定された位相特異点の変位 前後の位置を前記要素に基づいて特定する変位前後位置検出手段と、 前記変位前後位置検出手段で検出された変位前後の位相特異点の位置に基づ いて前記位相特異点の変位を測定する変位測定手段とを有することを特徴とする請 求項 16記載の変位検出装置。
[22] 前記変位測定手段は、前記変位前後位置検出手段で検出された変位前後の位相 特異点の位置の座標の差の二乗平方根を位相特異点の変位量として算出すること を特徴とする請求項 21記載の変位検出装置。
[23] 前記変位測定手段は、前記変位前後位置検出手段で検出された複数組の変位前 後の位相特異点の位置を結ぶ線分の複数の垂直二等分線の距離の和が最小となる 点を回転中心として求める回転中心取得手段と、
前記変位前後位置検出手段で検出された変位前後の位相特異点の位置と前記回 転中心取得手段で取得した前記回転中心とからなる三角形の前記回転中心を挟む 辺の角度を回転角として求める回転角度取得手段とを有することを特徴とする請求 項 21記載の変位検出装置。
[24] 取得した画像の強度信号をフーリエ変換するフーリエ変換手段と、
前記フーリエ変換手段でフーリエ変換された信号の複素共役の対称性を崩すよう なフィルタリング処理を行!ヽ、擬似位相情報を取得するフィルタリング処理手段と、 前記フィルタリング処理手段で取得した擬似位相情報に基づいて位相特異点を取 得する位相特異点取得手段を有することを特徴とする請求項 16記載の変位検出装 置。
[25] 前記フィルタリング処理手段は、前記フーリエ変換手段でフーリエ変換された信号に 対して中心に位相特異点を持つらせん状の位相特性を持ち、かつ、ドーナツ状の振 幅特性を持つラゲールガウスビームフィルタ特性により前記擬似位相情報を取得す ることを特徴とする請求項 24記載の変位検出装置。
[26] 前記画像の強度信号の平均値を前記画像の各強度信号から減算した画像を前記画 像とする前処理手段を有することを特徴とする請求項 24記載の変位検出装置。
[27] 前記フィルタリング処理手段は、前記フーリエ変換手段でフーリエ変換された信号に 対して中心に位相特異点を持つらせん状の位相特性を持つフィルタ特性、又は、ヒ ルペルト変換により前記擬似位相情報を取得することを特徴とする請求項 24記載の 変位検出装置。
[28] 変位前の特徴点の特徴パラメータと変位後の特徴点の特徴パラメータとの差異に基 づいて変位前の特徴点に対応する変位後の特徴点を探しだし、その特徴点の地理 的な位置関係の類似性力 変位前の特徴点に対応する変位後の特徴点のマツチン グを行うマッチング手段を、更に、有することを特徴とする変位検出装置。
[29] 前記マッチング手段は、変位前の複数の位相特異点間の特徴パラメータの総和距 離と変位後の複数の位相特異点間の特徴パラメータの総和距離との差分とが最小と なる特徴点の組を変位前と変位後とで互いにマッチングする特徴点として特定するこ とを特徴とする請求項 28記載の変位検出装置。
[30] 前記マッチング手段は、変位前の所定の位相特異点の特徴パラメータを T、 T、変 位後の所定の位相特異点の特徴パラメータを P 、 P 、 αを縮尺としたとき、
[数 2]
min = 〉 (I Tn - Τη \ -a \
Figure imgf000037_0001
\<i\,i2≤n
が最小となる位相特異点の組を変位前後でマッチングする位相特異点の組として抽 出することを特徴とする請求項 28記載の変位検出装置。
[31] コンピュータに、
位相特異点の位相構造力も所定の要素を取得し、前記要素に基づ 、て前記位相 特異点を特定する位相特異点特定手順と、
前記位相特異点特定手順で取得された要素に基づいて変位後の位相特異点の位 置を検出する変位検出手順とを実行させる変位検出プログラム。
[32] 前記位相特異点特定手順で取得される要素は、前記位相特異点近傍の実部と虚部 の値を平面で補間し、その結果、実部がゼロとなる線分と虚部がゼロとなる線分とが 交差する角度であることを特徴とする請求項 31記載の変位検出プログラム。
[33] 前記位相特異点特定手順で取得される要素は、前記位相特異点近傍の強度分布 の形状の離心率であることを特徴とする請求項 31記載の変位検出プログラム。
[34] 前記位相特異点特定手順で取得される要素は、前記位相特異点近傍の強度の勾 配であることを特徴とする請求項 31記載の変位検出プログラム。
[35] 前記位相特異点特定手順は、前記位相特異点近傍の実部と虚部の値を平面で補 間し、その結果、実部がゼロとなる線分と虚部がゼロとなる線分とが交差する角度と、 前記位相特異点近傍の強度分布の形状の離心率と、前記位相特異点近傍の強度 の勾配である渦度とを含むことを特徴とする請求項 31記載の変位検出プログラム。
[36] 前記変位検出手順は、前記位相特異点特定手順で特定された位相特異点の変位 前後の位置を前記要素に基づいて特定する変位前後位置検出手順と、
前記変位前後位置検出手順で検出された変位前後の位相特異点の位置に基づ いて前記位相特異点の変位を測定する変位測定手順とを有することを特徴とする請 求項 31記載の変位検出プログラム。
[37] 前記変位測定手順は、前記変位前後位置検出手順で検出された変位前後の位相 特異点の位置座標の差分の二乗平方根を位相特異点の変位量として算出すること を特徴とする請求項 35記載の変位検出プログラム。
[38] 前記変位測定手順は、前記変位前後位置検出手順で検出された複数組の変位前 後の位相特異点の位置を結ぶ線分の複数の垂直二等分線の距離の和が最小となる 点を回転中心として求める回転中心取得手順と、
前記変位前後位置検出手順で検出された変位前後の位相特異点の位置と前記回 転中心取得手順で取得した前記回転中心とからなる三角形の前記回転中心を挟む 辺の角度を回転角として求める回転角度取得手順とを有することを特徴とする請求 項 36記載の変位検出プログラム。
[39] 取得した画像の強度信号をフーリエ変換するフーリエ変換手順と、
前記フーリエ変換手順でフーリエ変換された信号の複素共役の対称性を崩すよう なフィルタリング処理を行!、、擬似位相情報を取得するフィルタリング処理手順と、 前記フィルタリング処理手順で取得した擬似位相情報に基づいて位相特異点を取 得する位相特異点取得手順を有することを特徴とする請求項 31記載の変位検出プ ログラム。
[40] 前記フィルタリング処理手順は、前記フーリエ変換手順でフーリエ変換された信号に 対して中心に位相特異点を持つらせん状の位相特性を持ち、かつ、ドーナツ状の振 幅特性を持つラゲールガウスビームフィルタ特性により前記擬似位相情報を取得す ることを特徴とする請求項 39記載の変位検出プログラム。
[41] 前記画像の強度信号の平均値を前記画像の各強度信号から減算した画像を前記画 像とする前処理手順を有することを特徴とする請求項 39記載の変位検出プログラム
[42] 前記フィルタリング処理手順は、前記フーリエ変換手順でフーリエ変換された信号に 対して中心に位相特異点を持つらせん状の位相特性を持つフィルタ特性、又は、ヒ ルペルト変換により前記擬似位相情報を取得することを特徴とする請求項 39記載の 変位検出プログラム。
[43] 変位前の位相特異点の特徴パラメータと変位後の位相特異点の特徴パラメータとの 差に基づいて前記変位前に特定された位相特異点と変位後に特定された位相特異 点とをマッチングを行うマッチング手順を、更に、有することを特徴とする請求項 31記 載の変位検出プログラム。
[44] 前記マッチング手順は、変位前の複数の位相特異点間の特徴パラメータの総和距 離と変位後の複数の位相特異点間の特徴パラメータの総和距離との差分とが最小と なる特徴点の組を変位前と変位後とで互いにマッチングする特徴点として特定するこ とを特徴とする請求項 43記載の変位検出プログラム。
[45] 前記マッチング手順は、変位前の所定の位相特異点の特徴パラメータを T、 T、変 位後の所定の位相特異点の特徴パラメータを P 、 P 、 αを縮尺としたとき、
[数 3]
min = > , (I T - Τ \ -a \ Pw( - Ρ^α) \f
が最小となる位相特異点の組を変位前後でマッチングする位相特異点の組として抽 出することを特徴とする請求項 43記載の変位検出プログラム。
[46] 変位前に特定された特徴点と変位後に特定された特徴点とをマッチングを行うポイン トマッチング処理方法であって、
変位前の特徴点の特徴パラメータと変位後の特徴点の特徴パラメータとの差異に 基づいて変位前の特徴点に対応する変位後の特徴点を探しだし、その特徴点の地 理的な位置関係の類似性力 変位前の特徴点に対応する変位後の特徴点のマッチ ングを行うことを特徴とすることを特徴とする特徴点マッチング処理方法。
[47] 変位前の複数の特徴点間の特徴パラメータの総和距離と変位後の複数の特徴点間 の特徴パラメータの総和距離との差分とが最小となる特徴点の組を変位前と変位後 とで互いにマッチングする特徴点として特定することを特徴とする請求項 37記載の特 徴点マッチング処理方法。
[48] 変位前の所定の特徴点の特徴パラメータを T、 T、変位後の所定の特徴点の特徴
il i2
ノ ラメータを P 、 P 、 ひを縮尺としたとき、
w(il) w(i2)
[数 4]
min = 〉 （I T - Ta \ -a |
Figure imgf000040_0001
|):
l≤il,i2≤"
が最小となる特徴点の組を変位前後でマッチングする特徴点の組として抽出すること を特徴とする請求項 38記載の特徴点マッチング処理方法。
[49] 前記特徴点は、画像の強度情報に基づいて取得される位相特異点であることを特徴 とする請求項 37記載の特徴点マッチング処理方法。
[50] 前記位相特異点は、前記特徴パラメータとして位置情報、離心率、及び、渦度、並び に、実軸と虚軸とが交わる角度を有することを特徴とする請求項 40記載の特徴点マ ツチング処理方法。
[51] コンピュータに、
変位前の特徴点の特徴パラメータと変位後の特徴点の特徴パラメータとの差異に 基づいて変位前の特徴点に対応する変位後の特徴点を探しだし、その特徴点の地 理的な位置関係の類似性力 変位前の特徴点に対応する変位後の特徴点のマッチ ングを行うことを特徴とするコンピュータにより読み取り可能な特徴点マッチング処理 プログラム。
[52] 変位前の複数の特徴点間の特徴パラメータの総和距離と変位後の複数の特徴点間 の特徴パラメータの総和距離との差分とが最小となる特徴点の組を変位前と変位後 とで互いにマッチングする特徴点として特定することを特徴とする請求項 51記載の特 徴点マッチング処理プログラム。
[53] 変位前の所定の特徴点の特徴パラメータを T、 T、変位後の所定の特徴点の特徴
il i2
パラメータを P 、 P 、縮尺をひとしたとき、
w(il) w(i2)
[数 5]
min = 2¾卜 " | " I)2
1≤1, 2<«
が最小となる特徴点の組を変位前後でマッチングする特徴点の組として抽出すること を特徴とする請求項 52記載の特徴点マッチング処理プログラム。
[54] 前記特徴点は、画像の強度情報に基づいて取得される位相特異点であることを特徴 とする請求項 51記載の特徴点マッチング処理プログラム。
[55] 前記位相特異点は、前記特徴パラメータとして位置情報、離心率、及び、渦度、並び に、実軸と虚軸とが交わる角度を有することを特徴とする請求項 54記載の特徴点マ ツチング処理プログラム。
PCT/JP2007/051095 2006-02-01 2007-01-24 変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、特徴点マッチング処理方法、特徴点マッチングプログラム WO2007088759A1 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
US12/162,370 US7813900B2 (en) 2006-02-01 2007-01-24 Displacement detection method, displacement detection device, displacement detection program, phase singularity matching method and phase singularity matching program
JP2007556824A JP4385139B2 (ja) 2006-02-01 2007-01-24 変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、位相特異点マッチング処理方法、位相特異点マッチング処理プログラム

Applications Claiming Priority (4)

Application Number Priority Date Filing Date Title
JP2006024846 2006-02-01
JP2006-024846 2006-02-01
JP2006-303238 2006-11-08
JP2006303238 2006-11-08

Publications (1)

Publication Number Publication Date
WO2007088759A1 true WO2007088759A1 (ja) 2007-08-09

Family

ID=38327338

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2007/051095 WO2007088759A1 (ja) 2006-02-01 2007-01-24 変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、特徴点マッチング処理方法、特徴点マッチングプログラム

Country Status (3)

Country Link
US (1) US7813900B2 (ja)
JP (1) JP4385139B2 (ja)
WO (1) WO2007088759A1 (ja)

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2009145208A (ja) * 2007-12-14 2009-07-02 National Institute Of Advanced Industrial & Technology 高さを測定する方法及び高さ測定装置
JP2010210571A (ja) * 2009-03-12 2010-09-24 Mitsutoyo Corp 画像相関変位計、及び変位測定方法
JP2017015742A (ja) * 2013-06-29 2017-01-19 堀　健治 位相変換作用を持つフィルター、レンズ、結像光学系及び撮像システム
CN109883333A (zh) * 2019-03-14 2019-06-14 武汉理工大学 一种基于图像特征识别技术的非接触位移应变测量方法
KR20190124026A (ko) * 2018-04-25 2019-11-04 경북대학교 산학협력단 구조물 변형 측정 장치 및 방법, 기록 매체
CN112926356A (zh) * 2019-12-05 2021-06-08 北京沃东天骏信息技术有限公司 一种目标跟踪方法和装置

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8116553B2 (en) * 2007-10-03 2012-02-14 Siemens Product Lifecycle Management Software Inc. Rotation invariant 2D sketch descriptor
US8730332B2 (en) 2010-09-29 2014-05-20 Digitaloptics Corporation Systems and methods for ergonomic measurement
US8917950B2 (en) * 2011-01-18 2014-12-23 Sony Corporation Simplifying parametric loop filters
US8587666B2 (en) * 2011-02-15 2013-11-19 DigitalOptics Corporation Europe Limited Object detection from image profiles within sequences of acquired digital images
US8705894B2 (en) 2011-02-15 2014-04-22 Digital Optics Corporation Europe Limited Image rotation from local motion estimates
US9285206B1 (en) * 2012-02-07 2016-03-15 Pile Dynamics, Inc. Measurement device for pile displacement and method for use of the same
JP2016009043A (ja) * 2014-06-24 2016-01-18 ソニー株式会社 イメージセンサ、演算方法、および電子装置
CN112614180B (zh) * 2020-12-31 2022-05-31 中国科学院长春光学精密机械与物理研究所 一种位移测量方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2002171395A (ja) * 2000-11-29 2002-06-14 Canon Inc 画像処理装置及びその方法並びに記憶媒体
JP2002190020A (ja) * 2000-12-20 2002-07-05 Monolith Co Ltd 映像効果方法および装置

Family Cites Families (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5118935A (en) * 1991-02-28 1992-06-02 Eastman Kodak Company Composite imaging mask
DE69631845T2 (de) 1996-11-01 2005-01-05 Yamatake Corp. Musterextrahierungsgerät
US6577966B2 (en) * 2000-06-21 2003-06-10 Siemens Corporate Research, Inc. Optimal ratio estimator for multisensor systems
US20020076116A1 (en) * 2000-12-15 2002-06-20 Xerox Corporation Fast implementation of homomorphic filters for image enhancement
JP3871309B2 (ja) * 2001-01-31 2007-01-24 フジノン株式会社 位相シフト縞解析方法およびこれを用いた装置
FR2863080B1 (fr) * 2003-11-27 2006-02-24 Advestigo Procede d'indexation et d'identification de documents multimedias
KR100981401B1 (ko) * 2004-04-22 2010-09-10 고쿠리쓰다이가쿠호진 덴키쓰신다이가쿠 미소 변위 계측법 및 장치
JP2006038779A (ja) * 2004-07-30 2006-02-09 Hitachi High-Technologies Corp パターン形状評価方法、評価装置、及び半導体装置の製造方法
US7539354B2 (en) * 2004-08-25 2009-05-26 Canon Kabushiki Kaisha Image database key generation method

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2002171395A (ja) * 2000-11-29 2002-06-14 Canon Inc 画像処理装置及びその方法並びに記憶媒体
JP2002190020A (ja) * 2000-12-20 2002-07-05 Monolith Co Ltd 映像効果方法および装置

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2009145208A (ja) * 2007-12-14 2009-07-02 National Institute Of Advanced Industrial & Technology 高さを測定する方法及び高さ測定装置
JP2010210571A (ja) * 2009-03-12 2010-09-24 Mitsutoyo Corp 画像相関変位計、及び変位測定方法
JP2017015742A (ja) * 2013-06-29 2017-01-19 堀　健治 位相変換作用を持つフィルター、レンズ、結像光学系及び撮像システム
KR20190124026A (ko) * 2018-04-25 2019-11-04 경북대학교 산학협력단 구조물 변형 측정 장치 및 방법, 기록 매체
KR102107145B1 (ko) * 2018-04-25 2020-05-06 경북대학교 산학협력단 구조물 변형 측정 장치 및 방법, 기록 매체
CN109883333A (zh) * 2019-03-14 2019-06-14 武汉理工大学 一种基于图像特征识别技术的非接触位移应变测量方法
CN112926356A (zh) * 2019-12-05 2021-06-08 北京沃东天骏信息技术有限公司 一种目标跟踪方法和装置

Also Published As

Publication number Publication date
US20090037137A1 (en) 2009-02-05
US7813900B2 (en) 2010-10-12
JPWO2007088759A1 (ja) 2009-06-25
JP4385139B2 (ja) 2009-12-16

Similar Documents

Publication Publication Date Title
WO2007088759A1 (ja) 変位検出方法、及び、変位検出装置、変位検出プログラム、並びに、特徴点マッチング処理方法、特徴点マッチングプログラム
Bhowmick et al. Identification of full-field dynamic modes using continuous displacement response estimated from vibrating edge video
Radkowski Object tracking with a range camera for augmented reality assembly assistance
de Oliveira et al. Improved impact damage characterisation in CFRP samples using the fusion of optical lock-in thermography and optical square-pulse shearography images
Guo et al. A novel multi-scale edge detection technique based on wavelet analysis with application in multiphase flows
Smith et al. CoLux: Multi-object 3d micro-motion analysis using speckle imaging
Javed et al. Vibration measurement of a rotating cylindrical structure using subpixel-based edge detection and edge tracking
Zhang et al. Accuracy improvement in laser stripe extraction for large-scale triangulation scanning measurement system
Gomit et al. Large-scale free surface measurement for the analysis of ship waves in a towing tank
Zheng et al. Image processing and edge detection techniques to quantify shock wave dynamics experiments
Wang et al. Target-less approach of vibration measurement with virtual points constructed with cross ratios
Abdallah et al. Three-dimensional point cloud analysis for automatic inspection of complex aeronautical mechanical assemblies
Guo et al. Visual pattern recognition supporting defect reporting and condition assessment of wastewater collection systems
Luo et al. Modeling and detection of heat haze in computer vision based displacement measurement
Tehsin et al. Improved maximum average correlation height filter with adaptive log base selection for object recognition
Wang et al. Similarity evaluation of 3D surface topography measurements
Prasad et al. A real-time feature-based clustering approach for vibration-based SHM of large structures
Kong et al. Automated fatigue crack identification through motion tracking in a video stream
Tsukamoto et al. Application of feature matching trajectory detection algorithm for particle streak velocimetry
Ramasamy et al. Data fusion strategy for multiscale surface measurements
Zhang et al. Flow visualization based on a derived rotation field
Boucheron et al. Dynamic trajectory measurement of a strained cable excited by mean flow using an improved optical stereovision system
JP4820998B2 (ja) 位相特異点検出方法、及び、位相特異点検出装置、並びに、プログラム
Al-Temeemy et al. Laser radar invariant spatial chromatic image descriptor
Zhang et al. Fast automatic multi-defects recognition based on the spatial characteristics of Freeman chain code

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application
DPE1 Request for preliminary examination filed after expiration of 19th month from priority date (pct application filed from 20040101)
ENP Entry into the national phase

Ref document number: 2007556824

Country of ref document: JP

Kind code of ref document: A

WWE Wipo information: entry into national phase

Ref document number: 12162370

Country of ref document: US

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 07707342

Country of ref document: EP

Kind code of ref document: A1

DPE1 Request for preliminary examination filed after expiration of 19th month from priority date (pct application filed from 20040101)