WO2016042779A1 - 三角測量装置、三角測量方法およびそのプログラムを記録した記録媒体 - Google Patents

三角測量装置、三角測量方法およびそのプログラムを記録した記録媒体 Download PDF

Info

Publication number
WO2016042779A1
WO2016042779A1 PCT/JP2015/004798 JP2015004798W WO2016042779A1 WO 2016042779 A1 WO2016042779 A1 WO 2016042779A1 JP 2015004798 W JP2015004798 W JP 2015004798W WO 2016042779 A1 WO2016042779 A1 WO 2016042779A1
Authority
WO
WIPO (PCT)
Prior art keywords
coordinates
correction amount
triangulation
correction
constraint condition
Prior art date
Application number
PCT/JP2015/004798
Other languages
English (en)
French (fr)
Inventor
中野　学
Original Assignee
日本電気株式会社
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 日本電気株式会社 filed Critical 日本電気株式会社
Priority to US15/510,737 priority Critical patent/US10126115B2/en
Priority to JP2016548565A priority patent/JP6547754B2/ja
Publication of WO2016042779A1 publication Critical patent/WO2016042779A1/ja

Links

Images

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
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C3/00—Measuring distances in line of sight; Optical rangefinders
    • G01C3/10—Measuring distances in line of sight; Optical rangefinders using a parallactic triangle with variable angles and a base of fixed length in the observation station, e.g. in the instrument
    • G01C3/14—Measuring distances in line of sight; Optical rangefinders using a parallactic triangle with variable angles and a base of fixed length in the observation station, e.g. in the instrument with binocular observation at a single point, e.g. stereoscopic type
    • 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
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C3/00—Measuring distances in line of sight; Optical rangefinders
    • G01C3/02—Details
    • G01C3/06—Use of electric means to obtain final indication
    • G01C3/08—Use of electric radiation detectors
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/50—Depth or shape recovery
    • G06T7/55—Depth or shape recovery from multiple images
    • G06T7/593—Depth or shape recovery from multiple images from stereo images
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/70—Determining position or orientation of objects or cameras
    • G06T7/73—Determining position or orientation of objects or cameras using feature-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2200/00—Indexing scheme for image data processing or generation, in general
    • G06T2200/04—Indexing scheme for image data processing or generation, in general involving 3D image data
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10004—Still image; Photographic image
    • G06T2207/10012—Stereo images

Definitions

  • the present invention relates to a triangulation apparatus or the like that calculates a three-dimensional position of a measurement target point using a stereo method.
  • Triangulation is the measurement of a measurement target point from two different calibrated cameras from different positions, and the point corresponding to the measurement target point in each image (hereinafter simply referred to as the corresponding point). This refers to estimating the optimal three-dimensional coordinates. “Calibrated” means that internal parameters such as the focal length or optical center of each camera and external parameters such as the position or orientation of each camera are known.
  • the optimum three-dimensional coordinates are the sum of squares of the distances between the two-dimensional image coordinates obtained by projecting the estimated three-dimensional coordinates onto each camera and the image coordinates of the observed corresponding points (hereinafter simply referred to as reprojection errors). It is the minimum point. Finding the point with the smallest reprojection error is equivalent to finding one point where a straight line passing through the center of each camera and the corresponding point intersects in the three-dimensional space. At this time, it is known that the relative geometric relationship between the corresponding point and each camera satisfies a constraint condition called an epipolar equation. When the corresponding point satisfies the epipolar equation, the optimal three-dimensional coordinate can be obtained by the linear least square method. However, since the observed corresponding points include errors due to the effects of various noises such as sensors, the epipolar equation is generally not satisfied, and the three-dimensional coordinates estimated by the linear least square method are not optimal.
  • a method called a bundle adjustment method is widely known as a general method for estimating three-dimensional coordinates.
  • the bundle adjustment method uses non-optimal 3D coordinates obtained by the linear least-squares method as an initial value, and estimates the optimal 3D coordinates by minimizing reprojection errors using a non-constrained nonlinear optimization method such as the Levenberg-Market method. is doing.
  • a method for estimating an initial value by a linear least square method is described in, for example, Non-Patent Document 1 and Non-Patent Document 2.
  • Non-Patent Documents 3, 4, and 5 propose a method of correcting the coordinates of corresponding points instead of directly optimizing the three-dimensional coordinates. That is, first, the coordinates of the optimum corresponding point satisfying the epipolar equation are obtained, and then the optimum three-dimensional coordinates are estimated by the linear least square method. In Non-Patent Document 3, the calculation of the globally optimum correction amount is reduced to a six-variable polynomial with one variable. It is known that the optimal correction amount is equivalent to the reprojection error.
  • Non-Patent Document 4 and Non-Patent Document 5 by linearly approximating an epipolar equation that is originally a quadratic equation, the correction amount of each corresponding point is calculated by repeating linear operations without solving the sixth-order equation. A scheme has been proposed.
  • the triangulation method by the bundle adjustment method may converge to a local solution that is not globally optimal because the search range of the three-dimensional coordinates by the unconstrained nonlinear optimization method is the entire three-dimensional space. That is, undesired three-dimensional coordinates may be estimated.
  • the triangulation method described in Non-Patent Document 3 can always estimate globally optimal three-dimensional coordinates.
  • the sixth-order equation used in the triangulation method cannot generally represent a solution by a finite number of algebraic operations (hereinafter simply referred to as a closed form) like a quadratic equation. Therefore, it cannot be solved without using an iterative method such as DKA (Durand-Kerner-Aberth) method or eigenvalue decomposition method.
  • the triangulation method of Non-Patent Document 3 has a large amount of calculation, and particularly when the number of corresponding points is large, an increase in the amount of calculation becomes a problem.
  • Non-Patent Document 4 and Non-Patent Document 5 have a very small amount of calculation for each iteration compared to the method used in Non-Patent Document 3. Therefore, even when the amount of calculation is small and the number of corresponding points is large, processing can be performed with a calculation time shorter than that of Non-Patent Document 3.
  • the triangulation methods of Non-Patent Document 4 and Non-Patent Document 5 are approximate methods of Non-Patent Document 3, there is no guarantee that globally optimal three-dimensional coordinates can be obtained.
  • the objective of this invention is providing the triangulation apparatus etc. which solve the said subject.
  • the triangulation device is: A triangulation device that calculates the three-dimensional position of a measurement point using a stereo method, The coordinates of the corresponding point corresponding to the measurement target point in the two images having the image of the measurement target point, the internal parameter and the external parameter of the optical apparatus that generated the two images are input and corrected.
  • a correction vector satisfying an epipolar equation in which the coordinates of the corresponding point are composed of the internal parameter and the external parameter is calculated, and the correction is calculated based on an eigen polynomial using the correction amount of the corresponding point coordinate or the reciprocal of the correction amount as a variable.
  • a three-dimensional coordinate calculation unit configured to calculate a three-dimensional coordinate of the measurement target point based on the corrected coordinates of the corresponding point, the internal parameter, and the external parameter;
  • the triangulation method is: A triangulation method executed by a triangulation device that has an optimal image coordinate estimation unit and a three-dimensional coordinate calculation unit and calculates a three-dimensional position of a measurement target point using a stereo method
  • the optimum image coordinate estimation unit includes coordinates of corresponding points corresponding to the measurement target points in two images having images of the measurement target points, internal parameters and external parameters of the optical apparatus that generated the two images.
  • the correction vector satisfying the epipolar equation composed of the internal parameter and the external parameter is used as a correction vector, and the correction amount of the coordinate of the corresponding point or the reciprocal of the correction amount is a variable.
  • the three-dimensional coordinate calculation unit calculates the three-dimensional coordinates of the measurement target point based on the corrected coordinates of the corresponding point, the internal parameter, and the external parameter.
  • the program according to the third aspect of the present invention is: A computer that calculates the three-dimensional position of the measurement point using the stereo method, The coordinates of the corresponding point corresponding to the measurement target point in the two images having the image of the measurement target point, the internal parameter and the external parameter of the optical apparatus that generated the two images are input and corrected.
  • the present invention can output the three-dimensional position without increasing the amount of calculation when calculating the three-dimensional position of the measurement target point using the stereo method.
  • FIG. 1 is a block diagram showing the configuration of the triangulation apparatus according to the first embodiment.
  • the direction of the arrow in the drawing shows an example, and does not limit the direction of the signal between the blocks.
  • the triangulation apparatus shown in FIG. 1 includes an optimum corresponding point coordinate estimation unit 1 and a three-dimensional coordinate calculation unit 2.
  • the optimum corresponding point coordinate estimation unit 1 includes an eigen polynomial calculation unit 11 and a correction vector calculation unit 12.
  • the optimal corresponding point coordinate estimation unit 1 inputs image coordinates of a set of corresponding points in two images, and internal parameters and external parameters of the camera that captured the corresponding points. Then, the optimum corresponding point coordinate estimation unit 1 calculates a correction vector that satisfies the epipolar equation in which the corrected corresponding point is composed of the internal parameter and the external parameter, using the correction amount of the corresponding point or the reciprocal of the correction amount as a variable. .
  • the eigenpolynomial calculation unit 11 receives a corresponding point, an internal parameter, and an external parameter, and calculates a predetermined constraint condition matrix that represents a constraint condition to be satisfied by a correction vector having a correction value of the corresponding point or a reciprocal of the correction amount as an eigenvalue. To do. Subsequently, the eigenpolynomial calculation unit 11 calculates the minimum correction amount by solving the eigen equation of the constraint condition matrix, and outputs the constraint condition matrix and the minimum correction amount to the correction vector calculation unit 12.
  • the correction vector calculation unit 12 receives the image coordinates of the corresponding points, the constraint condition matrix output from the eigenpolynomial calculation unit 11 and the minimum correction amount, and calculates the eigenvector of the constraint condition matrix corresponding to the minimum correction amount. . Subsequently, the correction vector calculation unit 12 performs a predetermined scaling on the eigenvector to calculate a correction vector, calculates the image coordinates of the optimum corresponding point from the correction vector, and outputs them to the three-dimensional coordinate calculation unit 2.
  • the three-dimensional coordinate calculation unit 2 receives the corrected image coordinates of the optimum corresponding point, the internal parameter, and the external parameter, calculates the three-dimensional coordinate of the corresponding point by the linear least square method, and outputs it.
  • the optimum corresponding point coordinate estimation unit 1 (more specifically, the eigen polynomial calculation unit 11 and the correction vector calculation unit 12) and the three-dimensional coordinate calculation unit 2 perform, for example, specific calculation processing and the like.
  • This is realized by designed hardware or an information processing apparatus such as a CPU (Central Processing Unit) that operates according to a program.
  • CPU Central Processing Unit
  • FIG. 2 is a flowchart showing an example of the operation of the first embodiment of the triangulation apparatus according to the present invention.
  • the image coordinates of a set of corresponding points in the two images and the internal parameters and external parameters of the camera that captured the corresponding points are input to the eigen polynomial calculation unit 11.
  • the eigenpolynomial calculation unit 11 calculates a predetermined constraint condition matrix that represents a constraint condition to be satisfied by a correction vector having a correction value of a corresponding point or a reciprocal of the correction amount as an eigenvalue.
  • the eigenpolynomial calculation unit 11 calculates the minimum correction amount by solving the eigen equation of the constraint condition matrix, and outputs the constraint condition matrix and the minimum correction amount to the correction vector calculation unit 12 (step S11).
  • the correction vector calculation unit 12 sets the minimum correction amount to the minimum correction amount. Compute the eigenvector of the corresponding constraint matrix. Subsequently, the correction vector calculation unit 12 calculates a correction vector by performing predetermined scaling on the eigenvector, calculates the image coordinates of the optimum corresponding point from the correction vector, and outputs the image coordinates to the three-dimensional coordinate calculation unit 2 (step S12). ).
  • the three-dimensional coordinate calculating unit 2 calculates the three-dimensional coordinate of the corresponding point by the linear least square method and outputs it. (Step S13).
  • the image coordinates of the corresponding points observed by the respective cameras are x, x ′
  • the image coordinates of the optimum corresponding points are x 1 , x ′
  • the correction vectors that optimize the observed corresponding points are ⁇ x, ⁇ x ′.
  • the internal parameters of the camera that captured the corresponding points are represented as K and K ′
  • the external parameters are represented as [Rt] and [R′t ′]
  • the optimum three-dimensional coordinates are represented as X.
  • K and K ′ are 3 ⁇ 3 (3 rows ⁇ 3 columns) matrices representing focal lengths and optical centers.
  • R and R ′ are 3 ⁇ 3 (3 rows ⁇ 3 columns) rotation matrices, and t and t ′ are three-dimensional translation vectors. Since the generality is not lost even if the first camera is the center of the three-dimensional coordinate system, R is a 3 ⁇ 3 (3 ⁇ 3) unit vector, and t is a zero vector in the following.
  • 0 mxn is zero matrix or the zero vector of m ⁇ n (m rows and n columns)
  • I n x n is a unit matrix of n ⁇ n (n rows and n columns)
  • ) is the L 2 norm of the vector
  • [] x is the matrix representation of the cross product of the three-dimensional vector
  • det is the determinant
  • superscript T is the transpose of the matrix or vector
  • the superscript -T is the transpose of the inverse matrix
  • ⁇ is indefinite.
  • x and x ′ can be detected by using a known matching point search technique such as template matching or SIFT (Scale-Invariant Feature Transform). Since x and x ′ include an error due to camera noise or the like, as shown in FIG. 3, the line of sight (solid line) has no intersecting point in the three-dimensional space, and the optimum three-dimensional coordinate X cannot be obtained. Therefore, as described in the triangulation method described in Non-Patent Document 3, the optimal x 1 , x ′, or correction from observation coordinates to optimum coordinates so that the line of sight (broken line) intersects at one point Vectors ⁇ x and ⁇ x ′ may be obtained.
  • ⁇ x and ⁇ x ′ are solutions to the constrained nonlinear optimization problem expressed by Equation 1 in FIG.
  • Equation 1 Next, a method for solving the constrained nonlinear optimization problem of Equation 1 will be described.
  • Equation 1 uses a vector b and a matrix A expressed by using F ij that is an (i, j) component of the matrix F obtained by multiplying the basic matrix F by T and T ′ T , as shown in FIG. It can be expressed as Equation 2. Since Expression 2 does not change even if both sides are multiplied by a constant because the right side becomes 0 due to the constraint condition, Expression 1 can be expressed by using a four-dimensional vector v T ⁇ [ ⁇ x T , ⁇ x ′ T ] and a certain constant k. This is equivalent to the constrained nonlinear optimization problem shown in Equation 3 of FIG.
  • Equation 4 When Expression 3 is represented by a Lagrangian function L (v, k, ⁇ ), Expression 4 in FIG. 6 is obtained.
  • the optimum solution of Equation 4 satisfies the following Equations 5 and 6, which are first-order optimality conditions (Karush-Kuhn-Tucker Condition: KKT conditions).
  • Equations 5 and 6 which are first-order optimality conditions (Karush-Kuhn-Tucker Condition: KKT conditions).
  • Equation 8 A ⁇ bb T / c Equation 9
  • D is shown as a constraint condition matrix for obtaining the optimum correction vector v.
  • Equation 10 in FIG. 6 is obtained.
  • Expression 10 does not change even if v is multiplied by a constant. That is, ⁇ can be interpreted as representing the relative magnitude of the correction amount of the corresponding point in Equation 3.
  • D is a real symmetric matrix
  • three non-zero eigenvalues are real numbers.
  • the eigenvalue having the minimum absolute value among the three eigenvalues may be selected, and in the case of ⁇ , the eigenvalue having the maximum absolute value among the three eigenvalues may be selected.
  • the eigenpolynominal calculation unit 11 may use, for example, a Cardano formula or a Lagrange formula expressed in a closed form.
  • Newton's method or power method may be used to obtain only the absolute value minimum or absolute value maximum solution.
  • the correction vector v can be calculated as an eigenvector corresponding to the correction vector v from Equation 8 or Equation 11.
  • the three-dimensional coordinate calculation unit 2 may use a method described in Non-Patent Document 1 or Non-Patent Document 2.
  • step S11 the eigenpolynomial calculator 11 calculates the image coordinates x, x ′ of the observed corresponding points, the internal parameters K, K ′ of the camera that captured the corresponding points, and the external parameters [Rt], [R′t ′]. ]
  • a constraint condition matrix D expressed by Equation 9 is calculated.
  • the eigen polynomial represented by Expression 12 or Expression 13 is solved, the correction amount ⁇ or the inverse number ⁇ of the correction amount is calculated, and output to the correction vector calculation unit 12 together with the constraint condition matrix D.
  • the internal parameters and external parameters of the camera can be acquired by, for example, calibration using a planar grid image in which the three-dimensional coordinates of the grid points are known, which are reflected in two cameras. Further, for example, a plurality of corresponding points obtained by template matching, SIFT, or the like may be used to obtain a method called self-calibration.
  • step S ⁇ b> 12 the correction vector calculation unit 12 first uses the constraint condition matrix D output from the eigenpolynomial calculation unit 11 and the correction amount ⁇ or the reciprocal ⁇ of the correction amount to first perform the correction expressed by Equation 8 or Equation 11.
  • the eigenvector v corresponding to the quantity or the inverse of the correction quantity is calculated.
  • k is obtained from Equation 7, and a correction vector is calculated with ⁇ x as the first and second components of the eigenvector multiplied by 1 / k of v ⁇ (1 / k) v and ⁇ x ′ as the third to fourth components.
  • the optimum corresponding points x 1 and x ′ are calculated from the correction vector and the image coordinates of the observed corresponding points, and output to the three-dimensional coordinate calculation unit 2.
  • step S13 the three-dimensional coordinate calculation unit 2 inputs the optimum corresponding points output from the correction vector calculation unit 12, and the internal parameters and external parameters of the camera that captured the corresponding points.
  • the three-dimensional coordinate calculation unit 2 calculates the three-dimensional coordinates X of the corresponding points by the linear least square method using the methods described in Non-Patent Document 1 and Non-Patent Document 2.
  • the optimum corresponding points x 1 and x ′ are calculated.
  • the optimum corresponding points x 1 and x ′ may be calculated without calculating the correction vectors ⁇ x and ⁇ x ′.
  • the [Delta] x in the above described x- x, also replacing each 'the x'x' [Delta] x in, may be the same calculation in each section.
  • a threshold such as 1e- 10 may be determined in advance, or computer epsilon may be used as the threshold.
  • the eigen polynomial calculation unit 11 receives a threshold value as an input.
  • the optimal corresponding point coordinate estimation unit 1 can obtain eigenvalues of the constraint matrix D by solving, for example, a cubic eigen equation, and then obtain eigenvectors corresponding to the eigenvalues.
  • the present invention is not limited to this.
  • the Jacobian method or the house holder method for obtaining eigenvalues and eigenvectors simultaneously may be used.
  • the triangulation device can reduce the order of the eigen polynomial in the correction vector calculation required for the image coordinates of the observed corresponding points to satisfy the epipolar equation.
  • the solution can be obtained by a cubic equation. That is, the triangulation device of this embodiment does not need to solve the sixth-order equation using a complicated iterative method as in Non-Patent Document 3.
  • the cubic equation can use a known closed-form solution, it is not necessary to use an iterative method such as the DKA method or the eigenvalue decomposition method that increases the amount of calculation.
  • the global optimality is guaranteed for the correction vector obtained by the triangulation device of the present embodiment. This is because, unlike Non-Patent Document 4 and Non-Patent Document 5, the epipolar equation is not approximated, and a solution that minimizes the reprojection error globally can be selected from the solution of the cubic equation. In other words, it is possible to select the eigenvalue of the constraint condition matrix that minimizes the correction amount.
  • the triangulation device of the present embodiment it is possible to realize a globally optimal triangulation with a low calculation amount.
  • FIG. 8 is a block diagram showing the configuration of the second embodiment of the triangulation apparatus according to the present invention.
  • the triangulation apparatus of this embodiment differs from the configuration of the first embodiment shown in FIG. 1 in that the observation corresponding point coordinates (2) are two-dimensional coordinates projected by the projector. Further, the optimum corresponding point coordinate estimation unit 1 and the three-dimensional coordinate calculation unit 2 are different in that they accept the internal parameters and external parameters of the projector. Since other configurations are the same as those in the first embodiment, the description thereof is omitted.
  • an image is projected by a calibrated projector, and the projected image is photographed by a calibrated camera to generate corresponding points.
  • the triangulation device detects the corresponding points from the photographed image by the camera using a known corresponding point search technique such as template matching or SIFT.
  • the triangulation apparatus may use the image coordinates as the observation corresponding point coordinates (1) and the direction projected by the projector as the observation corresponding point coordinates (2).
  • the direction projected by the projector is a three-dimensional vector, but its length is indefinite, so it is independent of only two dimensions and can be treated as corresponding point coordinates.
  • the operation of the triangulation apparatus of this embodiment is a two-dimensional coordinate on which the observation corresponding point coordinate (2) is projected by the projector. Since the operations are the same as those in the first embodiment except that the projector's internal parameters and external parameters are input to the optimum corresponding point coordinate estimation unit 1 and the three-dimensional coordinate calculation unit 2, description thereof will be omitted.
  • the triangulation device can realize globally optimum triangulation with a low amount of computation in a so-called active stereo composed of a camera and a projector.
  • the reason for this is that the same method as in the first embodiment can be applied to projecting light and photographing because they can be treated mathematically equivalently only with the difference in direction.
  • FIG. 9 is a block diagram when the triangulation apparatus according to the present invention is mounted on an information processing system.
  • the information processing system shown in FIG. 9 is a general information processing system including a processor 31, a program memory 32, and a storage medium 33.
  • the storage medium 33 may be a storage area composed of separate storage media, or may be a storage area composed of the same storage medium.
  • a magnetic storage medium such as a RAM (Random Access Memory) or a hard disk can be used.
  • the processing of each part of the optimum corresponding point coordinate estimation unit 1 (more specifically, the proper polynomial calculation unit 11 and the correction vector calculation unit 12) and the three-dimensional coordinate calculation unit 2 described above is performed.
  • a program to be executed by the computer 31 is stored.
  • the processor 31 operates according to this program.
  • the processor 31 may be a processor that operates according to a program such as a CPU, for example.
  • this embodiment can also be realized by a computer program. Note that it is not necessary to operate all the parts that can be operated by the program (such as the optimum corresponding point coordinate estimation unit 1 and the three-dimensional coordinate calculation unit 2) by the program, and some of them may be configured by hardware. Moreover, you may implement
  • the triangulation device 4 uses the passive stereo method to calculate the three-dimensional position of the measurement target point 71 on the structure 7 such as a bridge. Device.
  • the triangulation device 4 is connected to two cameras 5 and 6 that capture images of the structure 7 including the measurement target point 71 from different positions by wire or wirelessly.
  • the two cameras 5 and 6 are constituted by a CCD (Charge Coupled Device) camera or the like, and both have been calibrated.
  • CCD Charge Coupled Device
  • the triangulation device 4 includes a communication I / F unit (communication interface unit) 41, an operation input unit 42, a screen display unit 43, a storage unit 44, and an arithmetic processing unit 45.
  • the communication I / F unit 41 includes a dedicated data communication circuit, and has a function of performing data communication with various devices such as cameras 5 and 6 connected via a wired or wireless communication line. .
  • the operation input unit 42 includes an operation input device such as a keyboard and a mouse, and has a function of detecting an operator operation and outputting it to the arithmetic processing unit 45.
  • the screen display unit 43 includes a screen display device such as an LCD (Liquid Crystal Display) or a PDP (Plasma Display Panel), and displays the three-dimensional position of the measurement target point 71 on the screen according to an instruction from the arithmetic processing unit 45. It has a function to do.
  • a screen display device such as an LCD (Liquid Crystal Display) or a PDP (Plasma Display Panel)
  • LCD Liquid Crystal Display
  • PDP Plasma Display Panel
  • the storage unit 44 includes a storage device such as a hard disk or a memory, and has a function of storing processing information and programs 441 necessary for various processes in the arithmetic processing unit 45.
  • the program 441 is a program that realizes various processing units by being read into the arithmetic processing unit 45 and executed, and an external device (not shown) or the like via a data input / output function such as the communication I / F unit 41. It is read in advance from a storage medium (not shown) and stored in the storage unit 44.
  • Main processing information stored in the storage unit 44 includes images 442 and 443, internal parameters 444, and external parameters 445.
  • the image 442 is an image of the structure 7 including the measurement target point 71 photographed by the camera 5.
  • the image 443 is an image of the structure 7 including the measurement target point 71 photographed by the camera 5.
  • the internal parameter 444 is a parameter indicating the imaging characteristics of each camera, such as the focal length, the image center, and the pixel size of each of the cameras 5 and 6.
  • the external parameter 445 is a parameter indicating the positions, postures, and the like of the two cameras 5 and 6.
  • the arithmetic processing unit 45 includes a microprocessor and its peripheral circuits, and reads and executes the program 441 from the storage unit 44, thereby causing the hardware and the program 441 to cooperate to realize various processing units. Have. As main processing units realized by the arithmetic processing unit 45, there are a corresponding point detection unit 451, an optimum image coordinate estimation unit 452, and a three-dimensional coordinate calculation unit 453.
  • Corresponding point detection unit 451 reads images 442 and 443 photographed by cameras 5 and 6 from storage unit 44. Then, the corresponding point detection unit 451 uses a known corresponding point search method such as template matching or SIFT, the image coordinates corresponding to the measurement target point 71 in the image 442 of the camera 5, and the measurement target in the image 443 of the camera 6. It has a function of detecting image coordinates corresponding to the point 71.
  • a known corresponding point search method such as template matching or SIFT
  • the optimal image coordinate estimation unit 452 inputs the coordinates of the corresponding point corresponding to the measurement target point 71 in the images 442 and 443 from the corresponding point detection unit 451, and inputs the internal parameter 444 and the external parameter 445 from the storage unit 44.
  • the optimum image coordinate estimation unit 452 has a function of outputting the coordinates of the corrected corresponding point corresponding to the measurement target point 71 in the images 442 and 443.
  • the optimum image coordinate estimation unit 452 calculates a correction vector that satisfies the epipolar equation in which the coordinates of the corrected corresponding points are composed of the internal parameters 444 and the external parameters 445, and outputs the corrected coordinates of the corresponding points.
  • the correction vector is calculated based on a predetermined polynomial having a variable of the correction amount of the coordinates of the corresponding point or the reciprocal of the correction amount.
  • the optimal image coordinate estimation unit 452 includes an eigen polynomial calculation unit 4521 and a correction vector calculation unit 4522.
  • the eigenpolynomial calculation unit 4521 calculates a predetermined constraint condition matrix that represents a constraint condition to be satisfied by a correction vector having a coordinate value of the corresponding point or a reciprocal of the correction value as an eigenvalue.
  • the predetermined constraint condition matrix is calculated based on the coordinates of the corresponding point input from the corresponding point detection unit 451, the internal parameter 444 and the external parameter 445 input from the storage unit 44. It has a function of calculating the minimum correction amount by solving the eigen equation of the constraint condition matrix.
  • the correction vector calculation unit 4522 corresponds to the minimum correction amount based on the coordinates of the corresponding point input from the corresponding point detection unit 451, the constraint condition matrix calculated by the eigenpolynomial calculation unit 4521, and the minimum correction amount.
  • the eigenvector of the constraint matrix is calculated.
  • the correction vector calculation unit 4522 has a function of performing a predetermined scaling on the eigenvector to calculate a correction vector, and calculating and outputting the coordinates of the corresponding point corrected from the correction vector.
  • the three-dimensional coordinate calculation unit 453 receives the corrected coordinates of the corresponding points from the optimum image coordinate estimation unit 452, inputs the internal parameters 444 and the external parameters 445 from the storage unit 44, and the three-dimensional coordinates of the measurement target point 71. It has a function to calculate The three-dimensional coordinate calculation unit 453 has a function of outputting the calculated three-dimensional coordinates of the measurement target point 71 to the screen display unit 43 and transmitting it to another terminal (not shown) through the communication I / F unit 41.
  • FIG. 11 is a flowchart showing an example of the operation of the triangulation device 4 according to the present embodiment.
  • the operation of the triangulation apparatus 4 according to the present embodiment will be described with reference to FIG.
  • the corresponding point detection unit 451 of the triangulation device 4 captures an image including the measurement target point 71 with the cameras 5 and 6, and receives the image obtained by capturing from the cameras 5 and 6 with the communication I / F unit 41. Then, the images are stored as images 442 and 443 in the storage unit 44 (step S21). Next, the corresponding point detection unit 451 detects the image coordinates of the corresponding point corresponding to the measurement target point 71 in the images 442 and 443 captured by the cameras 5 and 6 and transmits the detected image coordinates to the optimum image coordinate estimation unit 452 ( Step S22).
  • the eigen polynomial calculation unit 4521 of the optimum image coordinate estimation unit 452 calculates a predetermined constraint condition matrix (corresponding to D in Expression 9 described above) from the internal parameters 444, the external parameters 445, and the image coordinates of the detected corresponding points. ) Is calculated (step S23).
  • the eigenpolynomial calculation unit 4521 calculates the minimum correction amount (corresponding to the above-described correction amount ⁇ or its inverse ⁇ ) by solving the eigen equation of this constraint condition matrix (step S24).
  • the correction vector calculation unit 4522 calculates an eigenvector (corresponding to the eigenvector v represented by the above-described equation 8 or 11) corresponding to the minimum correction amount (step S25).
  • the correction vector calculation unit 4522 performs predetermined scaling on the eigenvector to calculate a correction vector (corresponding to the correction vector v (or ⁇ x, ⁇ x ′) described above) (S26).
  • the correction vector calculation unit 4522 calculates the coordinates of the corresponding point corrected from the coordinates of the corresponding point and the correction vector, and transmits them to the three-dimensional coordinate calculation unit 453 (step S27).
  • the three-dimensional coordinate calculation unit 453 calculates and outputs the optimum three-dimensional coordinates of the corresponding points from the corrected corresponding point coordinates, internal parameters, and external parameters using the linear least square method. (Step S28).
  • the globally optimal three-dimensional position can be calculated with a small amount of calculation.
  • the reason is the same as described in the effect of the first embodiment. For example, it is sufficient to solve a cubic equation instead of a hexagonal equation, and a solution that minimizes the reprojection error globally can be selected from the solution of the cubic equation.
  • the triangulation device 8 uses the active stereo method to calculate the three-dimensional position of the measurement target point 71 on the structure 7 such as a bridge.
  • the triangulation device 8 is connected to a projector 9 that projects a predetermined image on the structure 7 and a camera 6 that captures an image projected on the structure 7 by the projector 9 by wire or wirelessly.
  • the camera 6 is the same as the camera 6 shown in FIG. 10, and has been calibrated.
  • the projector 9 can use a spot laser projector or a slit laser projector, but is not limited thereto.
  • a spot laser projector When a spot laser projector is used, a spot image with high brightness is projected on the surface of the structure 7, and the spot image location becomes the measurement target point 71.
  • a slit laser projector When a slit laser projector is used, a slit image with high brightness is projected on the surface of the structure 7, and the position of the slit image becomes the measurement target point 71.
  • the projector 9 has been calibrated.
  • the triangulation device 8 includes a communication I / F unit (communication interface unit) 81, an operation input unit 82, a screen display unit 83, a storage unit 84, and an arithmetic processing unit 85.
  • the communication I / F unit 81, the operation input unit 82, and the screen display unit 83 are the communication I / F unit 41, the operation input unit 42, and the screen display unit 43 in the third embodiment shown in FIG. Has the same function.
  • the storage unit 84 includes a storage device such as a hard disk or a memory, and has a function of storing processing information and programs 841 necessary for various processes in the arithmetic processing unit 85.
  • the program 841 is a program that implements various processing units by being read and executed by the arithmetic processing unit 85, and an external device (not shown) or the like via a data input / output function such as the communication I / F unit 81. It is read in advance from a storage medium (not shown) and stored in the storage unit 84.
  • Main processing information stored in the storage unit 84 includes images 842 and 843, internal parameters 844, and external parameters 845.
  • the image 842 is an image projected onto the structure 7 including the measurement target point 71 by the projector 9.
  • An image 843 is an image obtained by photographing the image projected on the structure 7 including the measurement target point 71 by the projector 9 with the camera 6.
  • the internal parameter 844 is a parameter indicating the imaging characteristics of the camera and the projector, such as the focal length, the image center, and the pixel size of the camera 6 and the projector 9.
  • the external parameter 845 is a parameter indicating the position, posture, etc. of the camera 6 and the projector 9.
  • the arithmetic processing unit 85 includes a microprocessor and its peripheral circuits.
  • the arithmetic processing unit 85 reads the program 841 from the storage unit 84 and executes it, thereby causing the hardware and the program 841 to cooperate with each other to realize various processing units. Have.
  • main processing units realized by the arithmetic processing unit 85 there are a corresponding point detection unit 851, an optimum image coordinate estimation unit 852, and a three-dimensional coordinate calculation unit 853.
  • the corresponding point detection unit 851 has a function of reading the image 842 projected on the structure 7 by the projector 9 from the storage unit 44 and detecting the image coordinates corresponding to the measurement target point 71 in the image 842 as the coordinates of the corresponding points.
  • the corresponding point detection unit 851 has a function of reading an image 843 captured by the camera 6 from the storage unit 84 and detecting the coordinates of the highest luminance location in the image as the coordinates of the corresponding point.
  • the image projected by the projector 9 may be known, and the corresponding point may be detected from the image captured by the camera 6 using a known corresponding point search technique such as template matching or SIFT.
  • the optimal image coordinate estimation unit 852 has the following functions. That is, the optimum image coordinate estimation unit 852 inputs the coordinates of the corresponding points corresponding to the measurement target point 71 in the images 842 and 843 from the corresponding point detection unit 851, and inputs the internal parameter 844 and the external parameter 845 from the storage unit 84. To do. Furthermore, the optimum image coordinate estimation unit 852 outputs the coordinates of the corrected corresponding point corresponding to the measurement target point 71 in the images 842 and 843.
  • the optimum image coordinate estimation unit 852 is based on a predetermined polynomial whose variable is the correction amount of the coordinate of the corresponding point or the inverse of the correction amount, and the coordinate of the corrected corresponding point is configured by the internal parameter 844 and the external parameter 845.
  • a correction vector that satisfies the epipolar equation is calculated, and the coordinates of the corrected corresponding point are output.
  • the optimal image coordinate estimation unit 852 includes an eigen polynomial calculation unit 8521 and a correction vector calculation unit 8522.
  • the eigenpolynomial calculation unit 8521 calculates a predetermined constraint condition matrix representing a constraint condition to be satisfied by the correction vector having the correction amount of the corresponding point coordinate or the reciprocal of the correction amount as the eigenvalue, and the eigen equation of the constraint condition matrix is calculated. It has a function to solve and calculate the minimum correction amount.
  • the predetermined constraint condition matrix is calculated based on the coordinates of the corresponding points input from the corresponding point detection unit 851 and the internal parameters 844 and external parameters 845 input from the storage unit 84.
  • the correction vector calculation unit 8522 calculates an eigenvector of the constraint condition matrix corresponding to the minimum correction amount, calculates a correction vector by applying a predetermined scaling to the eigenvector, and the corresponding points corrected from the correction vector It has a function to calculate and output the coordinates.
  • the eigenvector of the constraint condition matrix is calculated based on the coordinates of the corresponding point input from the corresponding point detection unit 851, the constraint condition matrix calculated by the eigenpolynomial calculation unit 8521, and the minimum correction amount.
  • the three-dimensional coordinate calculation unit 853 receives the coordinates of the corrected corresponding points from the optimum image coordinate estimation unit 852, inputs the internal parameters 844 and the external parameters 845 from the storage unit 84, and the three-dimensional coordinates of the measurement target point 71. It has a function to calculate The three-dimensional coordinate calculation unit 853 has a function of outputting the calculated three-dimensional coordinates of the measurement target point 71 to the screen display unit 83 and transmitting it to another terminal (not shown) through the communication I / F unit 81.
  • FIG. 13 is a flowchart showing an example of the operation of the triangulation device 8 according to the present embodiment.
  • the operation of the triangulation device 8 according to the present embodiment will be described with reference to FIG.
  • the corresponding point detection unit 851 of the triangulation device 4 transmits the image 842 read from the storage unit 84 to the projector 9 through the communication I / F unit 81, and the projector 9 projects the image 842 onto the structure 7 (step S31). ).
  • the corresponding point detection unit 851 captures an image including the image 842 projected on the structure 7 with the camera 6, receives the image obtained by capturing from the camera 6 with the communication I / F unit 81, and stores it.
  • the image is stored in the unit 84 as an image 843 (step S32).
  • the corresponding point detection unit 851 detects the image coordinates of the corresponding point corresponding to the measurement target point 71 in the images 842 and 843, and transmits it to the optimum image coordinate estimation unit 852 (step S33). Thereafter, processing similar to steps S23 to S28 shown in FIG. 11 is performed (steps S34 to S39).
  • the globally optimal three-dimensional position can be calculated with a small amount of calculation.
  • the reason is the same as described in the effect of the first embodiment. For example, it is sufficient to solve a cubic equation instead of a hexagonal equation, and a solution that minimizes the reprojection error globally can be selected from the solution of the cubic equation.
  • a triangulation device that calculates a three-dimensional position of a measurement target point using a stereo method
  • a storage unit that stores two images having an image of the measurement target point, an internal parameter and an external parameter of the optical device that generated the image, and a processor connected to the storage unit;
  • the processor is The coordinates of the corresponding point corresponding to the measurement target point in the two images, the internal parameter, and the external parameter are input, and the corrected coordinate of the corresponding point is configured by the internal parameter and the external parameter.
  • a correction vector satisfying the epipolar equation is calculated based on a correction coefficient of the coordinates of the corresponding point or a cubic polynomial having a variable of the reciprocal of the correction amount, and the corrected coordinates of the corresponding point are calculated,
  • a triangulation device programmed to receive the coordinates of the corrected corresponding point, the internal parameter, and the external parameter as input, and to calculate a three-dimensional coordinate of the measurement target point.
  • the processor is A predetermined constraint condition matrix representing a constraint condition to be satisfied by a correction vector having the coordinates of the corresponding point, the internal parameter, and the external parameter as input and the correction amount of the coordinate of the corresponding point or the reciprocal of the correction amount as an eigenvalue , Calculate the minimum correction amount by solving the eigen equation of the constraint matrix, The coordinates of the corresponding points, the constraint condition matrix, and the minimum correction amount are input, the eigenvector of the constraint condition matrix corresponding to the minimum correction amount is calculated, and the eigenvector is subjected to predetermined scaling to obtain a correction vector.
  • the triangulation device according to claim 1, wherein the triangulation device is programmed to calculate the coordinates of the corresponding points corrected from the correction vector.
  • the correction amount of the corresponding point is described in the supplementary note 1 or 2, which represents a Lagrange undetermined multiplier in a constrained nonlinear optimization problem in which an epipolar equation having an indefiniteness of a constant multiple is a constraint condition with a correction vector size as an objective function Triangulation instrument.
  • the stereo method is a passive stereo method in which two images of the measurement target point are captured by a camera from different positions.
  • the triangulation device which is calculated by [Appendix 7] A storage unit that stores two images having an image of the measurement target point, an internal parameter and an external parameter of the optical device that generated the image, and a processor connected to the storage unit, the measurement target A triangulation method executed by a triangulation device that calculates a three-dimensional position of a point using a stereo method,
  • the processor is The coordinates of the corresponding point corresponding to the measurement target point in the two images, the internal parameter, and the external parameter are input, and the corrected coordinate of the corresponding point is configured by the internal parameter and the external parameter.
  • a correction vector satisfying the epipolar equation is calculated based on a correction coefficient of the coordinates of the corresponding point or a cubic polynomial having a variable of the reciprocal of the correction amount, and the corrected coordinates of the corresponding point are calculated,
  • a triangulation method that receives the corrected coordinates of corresponding points, the internal parameters, and the external parameters as input, and calculates three-dimensional coordinates of the measurement target points.
  • the processor is A predetermined constraint condition matrix representing a constraint condition to be satisfied by a correction vector having the coordinates of the corresponding point, the internal parameter, and the external parameter as input and the correction amount of the coordinate of the corresponding point or the reciprocal of the correction amount as an eigenvalue , Calculate the minimum correction amount by solving the eigen equation of the constraint matrix, The coordinates of the corresponding points, the constraint condition matrix, and the minimum correction amount are input, the eigenvector of the constraint condition matrix corresponding to the minimum correction amount is calculated, and the eigenvector is subjected to predetermined scaling to obtain a correction vector.
  • the correction amount of the corresponding point is described in appendix 7 or 8, which represents a Lagrangian undetermined multiplier in a constrained nonlinear optimization problem in which an epipolar equation having a constant multiple indefiniteness is a constraint condition with a correction vector size as an objective function.
  • Triangulation method. [Appendix 10] The triangulation method according to any one of appendixes 7 to 9, wherein the stereo method is a passive stereo method in which two images of the measurement target point are captured by a camera from different positions.
  • the triangulation method which is calculated by: [Appendix 13] A storage unit that stores two images having an image of the measurement target point, an internal parameter and an external parameter of the optical device that generated the image, and a processor connected to the storage unit; To the processor of the triangulation device that calculates the three-dimensional position of using the stereo method, The coordinates of the corresponding point corresponding to the measurement target point in the two images, the internal parameter, and the external parameter are input, and the corrected coordinate of the corresponding point is configured by the internal parameter and the external parameter.
  • a correction vector that satisfies the epipolar equation is calculated based on a third-order polynomial having the correction amount of the coordinates of the corresponding point or the inverse of the correction amount as a variable, and calculates the coordinates of the corrected corresponding point,
  • the processor is A predetermined constraint condition matrix representing a constraint condition to be satisfied by a correction vector having the coordinates of the corresponding point, the internal parameter, and the external parameter as input and the correction amount of the coordinate of the corresponding point or the reciprocal of the correction amount as an eigenvalue , Calculate the minimum correction amount by solving the eigen equation of the constraint matrix, The coordinates of the corresponding points, the constraint condition matrix, and the minimum correction amount are input, the eigenvector of the constraint condition matrix corresponding to the minimum correction amount is calculated, and the eigenvector is subjected to predetermined scaling to obtain a correction vector.
  • the program according to attachment 13 wherein the coordinates of the corresponding point corrected from the correction vector are calculated.
  • Appendix 15 The correction amount of the corresponding point is described in appendix 13 or 14, which represents a Lagrange undetermined multiplier in a constrained nonlinear optimization problem in which an epipolar equation having an indefiniteness of a constant multiple is used as a constraint condition with a correction vector size as an objective function.
  • Program. [Appendix 16] The program according to any one of appendices 13 to 15, wherein the stereo method is a passive stereo method in which two images of the measurement target point are captured by a camera from different positions.

Abstract

　ステレオ法を使用して計測対象点の三次元位置を決定する際に、演算量を増加させることなく三次元位置を決定することができる三角測量装置を提供する。計測対象点の画像を有する２枚の画像における前記計測対象点に対応する対応点の座標と前記２枚の画像を生成した光学機器の内部パラメータと外部パラメータとに基づき、補正された前記対応点の座標が前記内部パラメータと前記外部パラメータとで構成されるエピポーラ方程式を満たす補正ベクトルを計算し、対応点の座標の補正量または該補正量の逆数を変数とする固有多項式を用いて、補正された前記対応点の座標を算出する最適画像座標推定手段と、前記補正された対応点の座標と前記内部パラメータと前記外部パラメータとに基づいて、前記計測対象点の三次元座標を計算する三次元座標計算手段と、を有する。

Description

三角測量装置、三角測量方法およびそのプログラムを記録した記録媒体
　本発明は、ステレオ法を使用して計測対象点の三次元位置を算出する三角測量装置等に関する。
　三角測量とは、キャリブレーション済みの２つのカメラで計測対象点を異なる位置から撮影し、それぞれの画像中の計測対象点に対応する点（以下、単に対応点と示す）から、計測対象点の最適な三次元座標を推定することをいう。キャリブレーション済みとは、各カメラの焦点距離又は光学中心等の内部パラメータと、各カメラの位置又は姿勢等の外部パラメータとが既知であることをいう。
　最適な三次元座標とは、推定した三次元座標を各カメラへ射影した２次元の画像座標と、観測した対応点の画像座標との距離の二乗和（以下、単に再投影誤差と示す）が最小となる点のことである。再投影誤差が最小となる点を求めることは、各カメラの中心と対応点を通る直線が、三次元空間中で交わる一点を求めることと等価である。このとき、対応点と各カメラの相対的な幾何関係は、エピポーラ方程式と呼ばれる拘束条件を満たすことが知られている。対応点がエピポーラ方程式を満たすとき、最適な三次元座標は線形最小二乗法により求めることができる。しかしながら、観測された対応点はセンサなどの様々なノイズの影響による誤差を含むため、一般にエピポーラ方程式は満たされず、線形最小二乗法により推定される三次元座標は最適ではない。
　バンドル調整法と呼ばれる方法が、三次元座標を推定するための一般的な方法として広く知られている。バンドル調整法は、線形最小二乗法によって得られる最適でない三次元座標を初期値として、レーベンバーグマーカート法などの制約なし非線形最適化法により再投影誤差を最小化して最適な三次元座標を推定している。線形最小二乗法により初期値を推定する方法は、例えば非特許文献１や非特許文献２に記載されている。
　非特許文献３、４、５では、直接に三次元座標を最適化するのではなく、各対応点の座標を補正する方式が提案されている。つまり、まずエピポーラ方程式を満たす最適な対応点の座標を求め、次に、線形最小二乗法により最適な三次元座標を推定する。そして、非特許文献３では、大域的に最適な補正量の計算を１変数の６次多項式に帰着させている。
最適な補正量は再投影誤差と等価なことが知られている。他方、非特許文献４や非特許文献５では、本来２次式であるエピポーラ方程式を１次近似することで、６次方程式を解かずに線形演算の反復で各対応点の補正量を計算する方式が提案されている。
徐剛、「写真から作る３次元ＣＧ」、近代科学社、２００１年、ｐ６９－７３ 佐藤淳、「コンピュータビジョン　視覚の幾何学」、コロナ社、１９９９年、ｐ１２２－１２３ Ｒ．Ｈａｒｔｌｅｙ　ａｎｄ　Ｐ．Ｓｔｕｒｍ，"Ｔｒｉａｎｇｕｌａｔｉｏｎ"、Ｉｎｔｅｒｎａｔｉｏｎａｌ　Ｊｏｕｒｎａｌ　ｏｆ　Ｃｏｍｐｕｔｅｒ　Ｖｉｓｉｏｎ　ａｎｄ　Ｉｍａｇｅ　Ｕｎｄｅｒｓｔａｎｄｉｎｇ．ｖｏｌ．６８，Ｉｓｓｕｅ２，ｐｐ．１４６－１５７，１９９７ 金谷、菅谷、新妻、"２画像からの三角測量 : Ｈａｒｔｌｅｙ－Ｓｔｕｒｍ　ｖｓ．最適補正"、電子情報通信学会技術研究報告 ＩＥ，画像工学、ｖｏｌ．１０７、ｎｏ５３８、ｐｐ．２２５－３４２、２００８年 Ｐ．Ｌｉｎｄｓｔｒｏｍ，"Ｔｒｉａｎｇｕｌａｔｉｏｎ　Ｍａｄｅ　Ｅａｓｙ"，ＩＥＥＥ　Ｃｏｍｐｕｔｅｒ　Ｖｉｓｉｏｎ　ａｎｄ　Ｐａｔｔｅｒｎ　Ｒｅｃｏｇｎｉｔｉｏｎ　２０１１，ｐｐ１５５４－１５６１，２０１０
　バンドル調整法による三角測量方法は、制約なし非線形最適化法による三次元座標の探索範囲が無限大の三次元空間全体であるため、大域的に最適でない局所解に収束する可能性がある。つまり、望ましくない三次元座標が推定されることがある。
　非特許文献３に記載されている三角測量方法は、常に大域的に最適な三次元座標を推定できる。しかしながら、その三角測量方法で使われている６次方程式は、２次方程式のように有限回の代数的な演算（以下、単に閉形式と示す）で解を表すことが一般にできない。このため、ＤＫＡ（Ｄｕｒａｎｄ－Ｋｅｒｎｅｒ－Ａｂｅｒｔｈ）法又は固有値分解法といった反復法を用いなければ解けない。非特許文献３の三角測量方法は、演算量が多くなり、特に対応点の数が多い場合、演算量の増加が問題になる。
　非特許文献４および非特許文献５に記載の三角測量方法は、反復毎の演算量が非特許文献３で用いられる方法に比べて非常に少ない。そのため、演算量は少なく、対応点数が多い場合でも、非特許文献３よりも少ない計算時間で処理ができる。しかしながら、非特許文献４および非特許文献５の三角測量方法は、非特許文献３の近似的な方法であるため、大域的に最適な三次元座標が得られる保証はない。
　このように、ステレオ法を使用した計測対象点の三次元位置を算出する場合、演算量を増加させることなく三次元位置を出力することは困難である。
［発明の目的］
　本発明の目的は、上記課題を解決する三角測量装置等を提供することにある。
　本発明の第１の観点に係る三角測量装置は、
　ステレオ法を使用して計測対象点の三次元位置を算出する三角測量装置であって、
　前記計測対象点の画像を有する２枚の画像における前記計測対象点に対応する対応点の座標と前記２枚の画像を生成した光学機器の内部パラメータと外部パラメータとを入力し、補正された前記対応点の座標が前記内部パラメータと前記外部パラメータとで構成されるエピポーラ方程式を満たす補正ベクトルを計算し、前記対応点の座標の補正量または該補正量の逆数を変数とする固有多項式に基づき計算し、補正された前記対応点の座標を算出する最適画像座標推定部と、
　前記補正された対応点の座標と前記内部パラメータと前記外部パラメータとに基づいて、前記計測対象点の三次元座標を計算する三次元座標計算部と
を有する。
　本発明の第２の観点に係る三角測量方法は、
　最適画像座標推定部と三次元座標計算部とを有し、ステレオ法を使用して計測対象点の三次元位置を算出する三角測量装置が実行する三角測量方法であって、
　前記最適画像座標推定部が、前記計測対象点の画像を有する２枚の画像における前記計測対象点に対応する対応点の座標と前記２枚の画像を生成した光学機器の内部パラメータと外部パラメータとを入力し、補正された前記対応点の座標が前記内部パラメータと前記外部パラメータとで構成されるエピポーラ方程式を満たす補正ベクトルを、前記対応点の座標の補正量または該補正量の逆数を変数とする固有多項式に基づき計算し、補正された前記対応点の座標を算出し、
　前記三次元座標計算部が、前記補正された対応点の座標と前記内部パラメータと前記外部パラメータとに基づいて、前記計測対象点の三次元座標を計算する。
　本発明の第３の観点に係るプログラムは、
　ステレオ法を使用して計測対象点の三次元位置を算出するコンピュータを、
　前記計測対象点の画像を有する２枚の画像における前記計測対象点に対応する対応点の座標と前記２枚の画像を生成した光学機器の内部パラメータと外部パラメータとを入力し、補正された前記対応点の座標が前記内部パラメータと前記外部パラメータとで構成されるエピポーラ方程式を満たす補正ベクトルを、前記対応点の座標の補正量または該補正量の逆数を変数とする固有多項式に基づき計算し、補正された前記対応点の座標を算出する最適画像座標推定部と、
　前記補正された対応点の座標と前記内部パラメータと前記外部パラメータとに基づいて、前記計測対象点の三次元座標を計算し三次元座標計算部と
して機能させる。
　本発明は、ステレオ法を使用した計測対象点の三次元位置を算出する際に、演算量を増加させることなく三次元位置を出力することができる。
第１の実施形態の三角測量装置の構成を示すブロック図である。 第１の実施形態の三角測量装置の動作の一例を示すフローチャートである。 本発明における三角測量の模式図である。 第１の実施形態の三角測量装置で使用する数式を示す図である。 第１の実施形態の三角測量装置で使用する数式を示す図である。 第１の実施形態の三角測量装置で使用する数式を示す図である。 第１の実施形態の三角測量装置で使用する数式を示す図である。 第２の実施形態による三角測量装置の構成を示すブロック図である。 本発明による三角測量装置を情報処理システムに実装した場合のブロック図である。 第３の実施形態による三角測量装置の構成を示すブロック図である。 第３の実施形態による三角測量装置の動作の一例を示すフローチャートである。 第４の実施形態による三角測量装置の構成を示すブロック図である。 第４の実施形態による三角測量装置の動作の一例を示すフローチャートである。
 [第１の実施形態]
　次に、本発明の第１の実施形態について図面を参照して詳細に説明する。
[構成の説明]
　図１は、第１の実施形態による三角測量装置の構成を示すブロック図である。図面中の矢印の方向は、一例を示すものであり、ブロック間の信号の向きを限定するものではない。なお、第２、第３及び第４の実施形態におけるブロック図でも同様である。図１に示す三角測量装置は、最適対応点座標推定部１と、三次元座標計算部２を備える。最適対応点座標推定部１は、固有多項式計算部１１と、補正ベクトル計算部１２とを含む。
　最適対応点座標推定部１は、２枚の画像における１組の対応点の画像座標と、対応点を撮影したカメラの内部パラメータと外部パラメータとを入力する。そして、最適対応点座標推定部１は、対応点の補正量もしくは補正量の逆数を変数として、補正された対応点が内部パラメータと外部パラメータとで構成されるエピポーラ方程式を満たす補正ベクトルを計算する。
　固有多項式計算部１１は、対応点と内部パラメータと外部パラメータとを入力とし、対応点の補正量もしくは補正量の逆数を固有値とする補正ベクトルが満たすべき拘束条件を表す所定の拘束条件行列を計算する。続いて、固有多項式計算部１１は、拘束条件行列の固有方程式を解いて最小の補正量を算出し、拘束条件行列と最小の補正量を補正ベクトル計算部１２に出力する。
　補正ベクトル計算部１２は、対応点の画像座標と、固有多項式計算部１１が出力する拘束条件行列と最小の補正量とを入力とし、最小の補正量に対応する拘束条件行列の固有ベクトルを計算する。続いて、補正ベクトル計算部１２は、固有ベクトルに所定のスケーリングを施して補正ベクトルを計算し、補正ベクトルから最適な対応点の画像座標を計算し、三次元座標計算部２へ出力する。
　三次元座標計算部２は、補正された最適な対応点の画像座標と、内部パラメータと外部パラメータとを入力とし、線形最小二乗法により対応点の三次元座標を計算し、出力する。
　本実施形態において、最適対応点座標推定部１（より具体的には、固有多項式計算部１１、補正ベクトル計算部１２）、三次元座標計算部２は、例えば、特定の演算処理等を行うよう設計されたハードウェア、またはプログラムに従って動作するＣＰＵ（Ｃｅｎｔｒａｌ　Ｐｒｏｃｅｓｓｉｎｇ　Ｕｎｉｔ）等の情報処理装置によって実現される。
[動作の説明]
　次に、本実施形態の三角測量装置の動作を説明する。図２は、本発明による三角測量装置の第１の実施形態の動作の一例を示すフローチャートである。
　まず、２枚の画像における１組の対応点の画像座標と、対応点を撮影したカメラの内部パラメータと外部パラメータとが固有多項式計算部１１に入力される。固有多項式計算部１１は、対応点の補正量もしくは補正量の逆数を固有値とする補正ベクトルが満たすべき拘束条件を表す所定の拘束条件行列を計算する。続いて、固有多項式計算部１１は、拘束条件行列の固有方程式を解いて最小の補正量を算出し、拘束条件行列と最小の補正量を補正ベクトル計算部１２に出力する（ステップＳ１１）。
　次いで、対応点の画像座標と、固有多項式計算部１１が出力する拘束条件行列と最小の補正量とが補正ベクトル計算部１２に入力されると、補正ベクトル計算部１２は、最小の補正量に対応する拘束条件行列の固有ベクトルを計算する。続いて、補正ベクトル計算部１２は、固有ベクトルに所定のスケーリングを施して補正ベクトルを計算し、補正ベクトルから最適な対応点の画像座標を計算し、三次元座標計算部２へ出力する（ステップＳ１２）。
　次いで、三次元座標計算部２は、補正された最適な対応点の画像座標と、内部パラメータと外部パラメータとが入力されると、線形最小二乗法により対応点の三次元座標を計算し、出力する（ステップＳ１３）。
　以下、具体的な例を用いて第１の実施形態を更に詳細に説明する。
　まず、最適な対応点の画像座標が満たすべきエピポーラ方程式を、図３を用いながら説明する。以下では、それぞれのカメラにより観測された対応点の画像座標をｘ、ｘ’、最適な対応点の画像座標をｘ、ｘ’、観測された対応点を最適にする補正ベクトルをΔｘ、Δｘ’と表す。また、対応点を撮影したカメラの内部パラメータをＫ、Ｋ’、外部パラメータを［Ｒｔ］、［Ｒ’ｔ’］、最適な三次元座標をＸと表す。ＫとＫ’は、焦点距離や光学中心を表す３×３（３行３列）の行列である。ＲとＲ’は、３×３（３行３列）の回転行列、ｔとｔ’は、三次元の平行移動ベクトルである。一番目のカメラが三次元座標系の中心としても一般性は失わないため、以下では、Ｒは３×３（３行３列）の単位ベクトル、ｔはゼロベクトルとする。また、０mｘnは、ｍ×ｎ（ｍ行ｎ列）のゼロ行列またはゼロベクトル、Ｉn x nは、ｎ×ｎ（ｎ行ｎ列）の単位行列、前後に二重線のある記号（||　||）はベクトルのＬ２ノルム、［　］xは、三次元ベクトルのクロス積の行列表現、ｄｅｔは行列式、上付き添字のＴは行列もしくはベクトルの転置、上付き添字の－１は逆行列、上付き添字の－Ｔは逆行列の転置、∝は定数倍に不定なことを表す。
　ｘ、ｘ’は、テンプレートマッチングやＳＩＦＴ（Ｓｃａｌｅ－Ｉｎｖａｒｉａｎｔ　Ｆｅａｔｕｒｅ　Ｔｒａｎｓｆｏｒｍ）といった、公知の対応点探索技術を用いて検出できる。ｘ、ｘ’は、カメラのノイズなどにより誤差を含むため、図３に示すように、視線（実線）は三次元空間で交わる点がなく、最適な三次元座標Ｘを求めることができない。そこで、非特許文献３に記載されている三角測量方法で述べられているように、視線（破線）が１点で交わるように最適なｘ、ｘ’、もしくは、観測座標から最適座標への補正ベクトルΔｘ、Δｘ’を求めればよい。ΔｘとΔｘ’は、図４の式１で表される制約付き非線形最適化問題の解である。
　式１中のＦは、基礎行列やｆｕｎｄａｍｅｎｔａｌ　ｍａｔｒｉｘなどと呼ばれ、その行列式ｄｅｔ（Ｆ）＝０、階数は２という特徴を持つことが知られている。Ｆは、基礎行列ＦにＴ、Ｔ’Ｔを乗じた行列なので、基礎行列Ｆと同様にその行列式はｄｅｔ（Ｆ）＝０、階数は２である。式１の拘束条件は、エピポーラ方程式と呼ばれ、これが満たされているとき、ｘ＝ｘ－Δｘ、ｘ’＝ｘ’－Δｘ’は、三次元空間中で一致することを示している。
　次に、式１の制約付き非線形最適化問題を解く方法について説明する。
　式１の拘束条件は、基礎行列ＦにＴ、Ｔ’Ｔを乗じた行列Ｆの（ｉ，ｊ）成分であるＦ ｉｊを用いて表されるベクトルｂと行列Ａを用いて、図５の式２のように表すことが出来る。
式２は、拘束条件により右辺が０となるため両辺を定数倍しても変化しないから、４次元ベクトルｖT ∝ ［ΔｘT，Δｘ’T］と、ある定数ｋを用いると、式１は、図６の式３に示す制約付き非線形最適化問題と等価である。
　式３をラグランジュ関数Ｌ（ｖ，ｋ，λ）で表すと、図６の式４を得る。式４の最適解は、１次最適性条件（Ｋａｒｕｓｈ－Ｋｕｈｎ－Ｔｕｃｋｅｒ Ｃｏｎｄｉｔｉｏｎ：ＫＫＴ条件）である以下の式５、式６を満たす。
∂Ｌ／∂ｖ＝ｖ－λ（Ａｖ＋ｋｂ）＝０4 x 1　　　・・・式５
∂Ｌ／∂ｋ＝λ（ｂTｖ＋ｋｃ）＝０　　　　　　・・・式６
　λは非ゼロであるから、式６より式７を得る。
ｋ＝－ｂTｖ／ｃ　　　・・・式７
　式７を式５に代入すると、式８を得る。
λＤｖ＝ｖ　　　・・・式８
ここで、Ｄ＝Ａ－ｂｂT／ｃ　　　・・・式９
以下では、Ｄを最適な補正ベクトルｖを求めるための拘束条件行列と示す。
　式８の両辺の内積を取ると、図６の式１０を得る。式１０は、ｖを定数倍しても値は変わらない。つまり、λは式３における対応点の補正量の相対的な大きさを表していると解釈できる。
　また、式８において、λの逆数をμ＝１／λとおくと、式１１となる。
Ｄｖ＝μｖ　　　・・・式１１
従って、μはＤの固有値である。
　Ｄは４×４の拘束条件行列であり、図７に示すように、その行列式ｄｅｔ（Ｄ）＝０となり、階数は３以下である。
　以下の式１２、式１３は、行列式ｄｅｔ（Ｄ）＝０により解の１つがそれぞれ∞、０と既知となる。つまり、補正量λまたは補正量の逆数μは、いずれかの固有方程式の解として次数を減らして計算することができる。
ｄｅｔ（λＤ－Ｉ4 x 4）＝０　　　・・・式１２
ｄｅｔ（Ｄ－μＩ4 x 4）＝０　　　・・・式１３
　Ｄは実対称行列であるから、０ではない３つの固有値は実数である。λの場合は、３つの固有値の中で絶対値最小となる固有値を、μの場合は、３つの固有値の中で絶対値最大となる固有値を選べばよい。
　式１２または式１３のすべての解を求めるために、固有多項式計算部１１は、例えば、閉形式で表されるカルダノの公式やラグランジュの公式を用いてもよい。また、絶対値最小もしくは絶対値最大の解のみを求めるために、ニュートン法やべき乗法を用いてもよい。
　λまたはμが求まれば、式８または式１１から補正ベクトルｖはそれに対応する固有ベクトルとして計算できる。固有値が既知であるときに固有ベクトルを計算する方法として、補正ベクトル計算部１２は、例えば、（λＤ－Ｉ4 x 4）ｖ＝０4 x 1、もしくは（Ｄ－μＩ4 x 4）ｖ＝０4 x 1を連立１次方程式として解いてもよいし、逆べき乗法を用いて解いてもよい。求めた固有ベクトルはｋ倍の不定性を持つので、式７より１／ｋ倍するとΔｘ、Δｘ’が得られる。
　補正された最適な対応点ｘ、ｘ’は、ｘ＝ｘ－Δｘ、ｘ’＝ｘ'－Δｘ’で算出される。補正された最適な対応点ｘ、ｘ’は、エピポーラ方程式を満たすため、視線の先は三次元空間中で一点に交わる。線形最小二乗法により三次元座標Ｘを推定するために、例えば、三次元座標計算部２は、非特許文献１や非特許文献２に記載されている方法を利用してもよい。
　次に、本実施形態における各部の動作を具体的に説明する。
　ステップＳ１１において、固有多項式計算部１１は、観測された対応点の画像座標ｘ、ｘ’と、対応点を撮影したカメラの内部パラメータＫ、Ｋ’と外部パラメータ［Ｒｔ］、［Ｒ’ｔ’］とを入力とし、まず、式９で表される拘束条件行列Ｄを計算する。次に、式１２または式１３で表される固有多項式を解き、補正量λまたは補正量の逆数μを計算し、拘束条件行列Ｄとともに補正ベクトル計算部１２に出力する。
　カメラの内部パラメータと外部パラメータは、例えば、２つのカメラに映る、格子点の三次元座標が既知の平面格子画像を用いるキャリブレーションによって取得できる。また、例えば、テンプレートマッチングやＳＩＦＴなどで取得した複数の対応点を用いて、セルフキャリブレーションと呼ばれる方法を用いて取得してもよい。
　ステップＳ１２において、補正ベクトル計算部１２は、固有多項式計算部１１から出力された拘束条件行列Ｄと補正量λまたは補正量の逆数μを用いて、まず、式８または式１１で表される補正量または補正量の逆数に対応する固有ベクトルｖを計算する。次に、式７よりｋを求め、ｖ←（１／ｋ）ｖと１／ｋ倍した固有ベクトルの第１～２成分をΔｘ、その第３～４成分をΔｘ’として補正ベクトルを計算する。次に、補正ベクトルと観測された対応点の画像座標から最適な対応点ｘ、ｘ’を計算し、三次元座標計算部２へ出力する。
　ステップＳ１３において、三次元座標計算部２は、補正ベクトル計算部１２から出力された最適な対応点と、対応点を撮影したカメラの内部パラメータと外部パラメータとを入力する。例えば、三次元座標計算部２は、非特許文献１や非特許文献２に記載されている方法を利用して、線形最小二乗法により対応点の三次元座標Ｘを計算する。
　本実施形態では、補正ベクトルΔｘ、Δｘ’を計算してから最適な対応点ｘ、ｘ’を計算した。しかし、本発明では、補正ベクトルΔｘ、Δｘ’を計算せずに、最適な対応点ｘ、ｘ’を計算してもよい。この場合、上記の説明中のΔｘをｘ－ｘに、またΔｘ’をｘ'－ｘ’にそれぞれ置き換えて、各部で同様の計算をすればよい。
　本実施形態では、固有多項式計算部１１は、拘束条件行列Ｄを式９にもとづいて計算したが、ｃ＝２Ｆ 3,3の絶対値が小さい場合に数値計算が不安定になる可能性があるため、式９の右辺をｃ倍したものを拘束条件行列Ｄとしてもよい。この場合、固有値である補正量λまたは補正量の逆数μもｃ倍されるが、固有ベクトルは変化しないため、補正ベクトル計算部１２の動作は変わらない。
　また、固有多項式計算部１１は、ｃの絶対値が０に非常に近い場合、以降の処理を省略し、ｖ＝０4 x 1を補正ベクトル計算部へ出力してもよい。０に近いかどうかの判定は、例えば、１ｅ-10のような閾値を事前に決めてもよいし、計算機イプシロンを閾値としてもよい。この場合、固有多項式計算部１１は、閾値を入力として受け付ける。
　本実施形態では、最適対応点座標推定部１は、例えば、３次固有方程式を解いて拘束条件行列Ｄの固有値を求め、次に固有値に対応する固有ベクトルを求めることができる。但し、これに限られるものではない。例えば、固有値と固有ベクトルを同時に求めるヤコビ法又はハウスホルダー法を用いても良い。
[第１の実施形態の効果]
　以上のように、本実施形態の三角測量装置は、まず、観測された対応点の画像座標がエピポーラ方程式を満たすのに必要な補正ベクトルの計算において、固有多項式の次数を減らすことができ、例えば、３次方程式によってその解を求めることができる。すなわち、本実施形態の三角測量装置は、非特許文献３のように複雑な反復法を用いて６次方程式を解く必要がない。また、３次方程式は公知の閉形式で表される解法が利用できるため、ＤＫＡ法や固有値分解法といった演算量が増加するような反復法を用いる必要がない。
　さらに、本実施形態の三角測量装置で求めた補正ベクトルは、大域的な最適性が保証されている。これは、非特許文献４や非特許文献５のようにエピポーラ方程式を近似せず、３次方程式の解から再投影誤差が大域的に最小となる解を選択できるためである。換言すれば、補正量が最小となる拘束条件行列の固有値を選択できるためである。
　よって、本実施形態の三角測量装置によれば、低演算量で、大域的に最適な三角測量を実現できる。
[第２の実施形態]
　図８は、本発明による三角測量装置の第２の実施形態の構成を示すブロック図である。本実施形態の三角測量装置は、図１に示した第１の実施形態の構成と比べて、観測対応点座標（２）が、プロジェクタにより投影される二次元座標であることが異なる。また、最適対応点座標推定部１と三次元座標計算部２が、プロジェクタの内部パラメータと外部パラメータを受け付ける点が異なる。それ以外の構成に関しては、第１の実施形態と同様であるため、説明を省略する。
　本実施形態では、キャリブレーション済みのプロジェクタで映像が投影され、投影された映像をキャリブレーション済みのカメラで撮影し、対応点が生成される。例えば、投影される映像を既知として、三角測量装置は、カメラによる撮影画像の中からテンプレートマッチングやＳＩＦＴといった公知の対応点探索技術を用いて対応点を検出する。三角測量装置は、その画像座標を観測対応点座標（１）、プロジェクタで投影した方向を観測対応点座標（２）とすればよい。プロジェクタで投影した方向は三次元ベクトルであるが、その長さは不定であるため、独立なのは二次元しかないため、対応点座標として扱える。
　本実施形態の三角測量装置の動作は、観測対応点座標（２）がプロジェクタにより投影される二次元座標である。最適対応点座標推定部１と三次元座標計算部２にプロジェクタの内部パラメータと外部パラメータが入力される点以外は、第１の実施形態の動作と同様であるため、説明を省略する。
　本実施形態の三角測量装置は、カメラとプロジェクタで構成された、いわゆるアクティブステレオにおいて低演算量で大域的に最適な三角測量を実現できる。その理由は、光を投影することと撮影することは方向の正負が異なるだけで数学的には等価に扱うことができるため、第１の実施形態と同様の方法が適用できるためである。
　なお、上記の各実施形態の三角測量装置は、例えば、各部に対応したハードウェア等により実現可能であり、情報処理システムによっても実現可能である。図９は、本発明による三角測量装置を情報処理システムに実装した場合のブロック図である。図９に示す情報処理システムは、プロセッサ３１と、プログラムメモリ３２と、記憶媒体３３とを備える一般的な情報処理システムである。記憶媒体３３は、別個の記憶媒体からなる記憶領域であってもよいし、同一の記憶媒体からなる記憶領域であってもよい。記憶媒体として、ＲＡＭ（Ｒａｎｄｏｍ　Ａｃｃｅｓｓ　Ｍｅｍｏｒｙ）や、ハードディスク等の磁気記憶媒体を用いることができる。
　プログラムメモリ３２には、上述した、最適対応点座標推定部１（より具体的には、固有多項式計算部１１、補正ベクトル計算部１２）と、三次元座標計算部２の各部の処理を、プロセッサ３１に行わせるためのプログラムが格納されている。このプログラムに従ってプロセッサ３１は動作する。プロセッサ３１は、例えば、ＣＰＵ等のプログラムに従って動作するプロセッサであればよい。このように、本実施形態は、コンピュータプログラムにより実現することも可能である。なお、プログラムによる動作が可能な部分（最適対応点座標推定部１や三次元座標計算部２など）のすべてをプログラムで動作させる必要はなく、一部をハードウェアで構成してもよい。また、それぞれ別々のユニットとして実現されていてもよい。
[第３の実施形態]
　図１０を参照すると、本発明の第３の実施形態に係る三角測量装置４は、パッシブステレオ法を使用して橋梁などの構造物７上の計測対象点７１の三次元位置を算出する情報処理装置である。三角測量装置４は、計測対象点７１を含む構造物７の画像をお互いに異なる位置から撮影する２台のカメラ５、６に有線または無線により接続されている。２台のカメラ５、６は、ＣＣＤ（Ｃｈａｒｇｅ　Ｃｏｕｐｌｅｄ　Ｄｅｖｉｃｅ）カメラ等で構成され、双方ともキャリブレーション済みである。
　三角測量装置４は、通信Ｉ／Ｆ部（通信インターフェイス部）４１、操作入力部４２、画面表示部４３、記憶部４４、および演算処理部４５を有する。
　通信Ｉ／Ｆ部４１は、専用のデータ通信回路からなり、有線または無線による通信回線を介して接続されたカメラ５、６などの各種装置との間でデータ通信を行う機能を有している。
　操作入力部４２は、キーボードやマウスなどの操作入力装置からなり、オペレータの操作を検出して演算処理部４５に出力する機能を有している。
　画面表示部４３は、ＬＣＤ（Ｌｉｑｕｉｄ　Ｃｒｙｓｔａｌ　Ｄｉｓｐｌａｙ）やＰＤＰ（Ｐｌａｓｍａ　Ｄｉｓｐｌａｙ　Ｐａｎｅｌ）などの画面表示装置からなり、演算処理部４５からの指示に応じて、計測対象点７１の三次元位置などを画面表示する機能を有している。
　記憶部４４は、ハードディスクやメモリなどの記憶装置からなり、演算処理部４５における各種処理に必要な処理情報やプログラム４４１を記憶する機能を有している。プログラム４４１は、演算処理部４５に読み込まれて実行されることにより各種処理部を実現するプログラムであり、通信Ｉ／Ｆ部４１などのデータ入出力機能を介して外部装置（図示せず）や記憶媒体（図示せず）から予め読み込まれて記憶部４４に保存される。記憶部４４で記憶される主な処理情報として、画像４４２、４４３、内部パラメータ４４４、外部パラメータ４４５がある。
　画像４４２は、カメラ５により撮影された計測対象点７１を含む構造物７の画像である。画像４４３は、カメラ５により撮影された計測対象点７１を含む構造物７の画像である。内部パラメータ４４４は、カメラ５、６それぞれの焦点距離、画像中心、画素サイズなど、カメラそれぞれの結像の特性を示すパラメータである。外部パラメータ４４５は、２つのカメラ５、６の位置、姿勢などを示すパラメータである。
　演算処理部４５は、マイクロプロセッサとその周辺回路を有し、記憶部４４からプログラム４４１を読み込んで実行することにより、上記ハードウェアとプログラム４４１とを協働させて各種処理部を実現する機能を有している。演算処理部４５で実現される主な処理部として、対応点検出部４５１と最適画像座標推定部４５２と三次元座標計算部４５３とがある。
　対応点検出部４５１は、カメラ５、６で撮影された画像４４２、４４３を記憶部４４から読み出す。そして、対応点検出部４５１は、テンプレートマッチング又はＳＩＦＴといった公知の対応点探索手法を用いて、カメラ５の画像４４２中の計測対象点７１に対応する画像座標、カメラ６の画像４４３中の計測対象点７１に対応する画像座標を検出する機能を有する。
　最適画像座標推定部４５２は、画像４４２、４４３における計測対象点７１に対応する対応点の座標を対応点検出部４５１から入力し、内部パラメータ４４４と外部パラメータ４４５とを記憶部４４から入力する。最適画像座標推定部４５２は、画像４４２、４４３における計測対象点７１に対応する補正された対応点の座標を出力する機能を有する。最適画像座標推定部４５２は、補正された対応点の座標が内部パラメータ４４４と外部パラメータ４４５とで構成されるエピポーラ方程式を満たす補正ベクトルを計算し、補正された対応点の座標を出力する。なお、当該補正ベクトルは、対応点の座標の補正量または該補正量の逆数を変数とする所定の多項式に基づいて計算される。
　最適画像座標推定部４５２は、固有多項式計算部４５２１と補正ベクトル計算部４５２２とを有する。固有多項式計算部４５２１は、対応点の座標の補正量または該補正量の逆数を固有値とする補正ベクトルが満たすべき拘束条件を表す所定の拘束条件行列を計算する。所定の拘束条件行列は、対応点検出部４５１から入力した対応点の座標と記憶部４４から入力した内部パラメータ４４４と外部パラメータ４４５に基づき計算される。その拘束条件行列の固有方程式を解いて最小の補正量を算出する機能を有する。また補正ベクトル計算部４５２２は、対応点検出部４５１から入力した対応点の座標と固有多項式計算部４５２１で計算された拘束条件行列と最小の補正量とに基づいて、上記最小の補正量に対応する上記拘束条件行列の固有ベクトルを計算する。さらに、補正ベクトル計算部４５２２は、上記固有ベクトルに所定のスケーリングを施して補正ベクトルを計算し、上記補正ベクトルから補正された対応点の座標を計算し出力する機能を有する。
　三次元座標計算部４５３は、補正された対応点の座標を最適画像座標推定部４５２から入力し、内部パラメータ４４４と外部パラメータ４４５とを記憶部４４から入力し、計測対象点７１の三次元座標を計算する機能を有する。また、三次元座標計算部４５３は、計算した計測対象点７１の三次元座標を画面表示部４３に出力し、また通信Ｉ／Ｆ部４１を通じて図示しない他の端末へ送信する機能を有する。
　図１１は、本実施形態に係る三角測量装置４の動作の一例を示すフローチャートである。以下、図１１を参照して、本実施形態に係る三角測量装置４の動作について説明する。
　三角測量装置４の対応点検出部４５１は、計測対象点７１を含む画像をカメラ５、６により撮影し、撮影して得られた画像を通信Ｉ／Ｆ部４１でカメラ５、６から受信し、記憶部４４に画像４４２、４４３として記憶する（ステップＳ２１）。次に、対応点検出部４５１は、カメラ５、６で撮影された画像４４２、４４３中の計測対象点７１に対応する対応点の画像座標を検出し、最適画像座標推定部４５２に伝達する（ステップＳ２２）。
　次に、最適画像座標推定部４５２の固有多項式計算部４５２１は、内部パラメータ４４４と外部パラメータ４４５と検出された対応点の画像座標とから、所定の拘束条件行列（前述した式９のＤに相当）を計算する（ステップＳ２３）。次に固有多項式計算部４５２１は、この拘束条件行列の固有方程式を解いて、最小の補正量（前述した補正量λあるいはその逆数μに相当）を計算する（ステップＳ２４）。なお、拘束条件行列の固有多項式となる式１２、式１３は、行列式ｄｅｔ（Ｄ）＝０により解の１つがそれぞれ∞、０と既知となる。つまり、補正量λまたは補正量の逆数μは、いずれかの固有多項式に対して次数を減らして解を求めることができる。
　次に、補正ベクトル計算部４５２２は、最小の補正量に対応する拘束条件行列の固有ベクトル（前述した式８または式１１で表される固有ベクトルｖに相当）を計算する（ステップＳ２５）。次に、補正ベクトル計算部４５２２は、その固有ベクトルに所定のスケーリングを施して補正ベクトル（前述した補正ベクトルｖ（あるいはΔｘ、Δｘ’）に相当）を計算する（Ｓ２６）。次に補正ベクトル計算部４５２２は、対応点の座標と補正ベクトルとから補正された対応点の座標を計算し、三次元座標計算部４５３へ伝達する（ステップＳ２７）。
　次に、三次元座標計算部４５３は、補正された対応点の座標と内部パラメータと外部パラメータとから、線形最小二乗法を使用して、対応点の最適な三次元座標を計算し、出力する（ステップＳ２８）。
　このように本実施形態によれば、パッシブステレオ法を使用して計測対象点の三次元位置を算出する三角測量装置において、大域的に最適な三次元位置を少ない演算量で計算することができる。その理由は、第１の実施形態の効果における説明と同様である。例えば、６次方程式ではなく３次方程式を解けばよく、また３次方程式の解から再投影誤差が大域的に最小となる解を選択できる。
[第４の実施形態]
　図１２を参照すると、本発明の第４の実施形態に係る三角測量装置８は、アクティブステレオ法を使用して橋梁などの構造物７上の計測対象点７１の三次元位置を算出する情報処理装置である。三角測量装置８は、構造物７に所定の画像を投影するプロジェクタ９と、プロジェクタ９によって構造物７に投影された画像を撮影するカメラ６とに有線または無線により接続されている。
　カメラ６は、図１０に示したカメラ６と同じであり、キャリブレーション済みである。
プロジェクタ９は、スポットレーザ投光器やスリットレーザ投光器などを使用できるが、それらに限定されない。スポットレーザ投光器を使用した場合、構造物７の面上に高輝度のスポット像が投影され、そのスポット像の箇所が計測対象点７１となる。またスリットレーザ投光器を使用した場合、構造物７の面上に高輝度なスリット像が投影され、そのスリット像の箇所が計測対象点７１となる。プロジェクタ９は、キャリブレーション済みである。
　三角測量装置８は、通信Ｉ／Ｆ部（通信インターフェイス部）８１、操作入力部８２、画面表示部８３、記憶部８４、および演算処理部８５を有する。このうち、通信Ｉ／Ｆ部８１、操作入力部８２、および画面表示部８３は、図１０に示した第３の実施形態における通信Ｉ／Ｆ部４１、操作入力部４２、および画面表示部４３と同様の機能を有する。
　記憶部８４は、ハードディスクやメモリなどの記憶装置からなり、演算処理部８５における各種処理に必要な処理情報やプログラム８４１を記憶する機能を有している。プログラム８４１は、演算処理部８５に読み込まれて実行されることにより各種処理部を実現するプログラムであり、通信Ｉ／Ｆ部８１などのデータ入出力機能を介して外部装置（図示せず）や記憶媒体（図示せず）から予め読み込まれて記憶部８４に保存される。記憶部８４で記憶される主な処理情報として、画像８４２、８４３、内部パラメータ８４４、外部パラメータ８４５がある。
　画像８４２は、プロジェクタ９によって計測対象点７１を含む構造物７に投影される画像である。画像８４３は、プロジェクタ９によって計測対象点７１を含む構造物７に投影された画像をカメラ６により撮影した画像である。内部パラメータ８４４は、カメラ６およびプロジェクタ９の焦点距離、画像中心、画素サイズ、など、カメラおよびプロジェクタの結像の特性を示すパラメータである。外部パラメータ８４５は、カメラ６とプロジェクタ９の位置、姿勢などを示すパラメータである。
　演算処理部８５は、マイクロプロセッサとその周辺回路を有し、記憶部８４からプログラム８４１を読み込んで実行することにより、上記ハードウェアとプログラム８４１とを協働させて各種処理部を実現する機能を有している。演算処理部８５で実現される主な処理部として、対応点検出部８５１と最適画像座標推定部８５２と三次元座標計算部８５３とがある。
　対応点検出部８５１は、プロジェクタ９によって構造物７に投影される画像８４２を記憶部４４から読み出し、その画像８４２中の計測対象点７１に対応する画像座標を対応点の座標として検出する機能を有する。また対応点検出部８５１は、カメラ６で撮影された画像８４３を記憶部８４から読み出し、画像中で最も輝度の高い箇所の座標を対応点の座標として検出する機能を有している。勿論、プロジェクタ９によって投影される画像を既知として、カメラ６による撮影画像の中からテンプレートマッチングやＳＩＦＴといった公知の対応点探索技術を用いて対応点を検出してもよい。
　最適画像座標推定部８５２は、以下の機能を有する。すなわち、最適画像座標推定部８５２は、画像８４２、８４３における計測対象点７１に対応する対応点の座標を対応点検出部８５１から入力し、内部パラメータ８４４と外部パラメータ８４５とを記憶部８４から入力する。さらに、最適画像座標推定部８５２は、画像８４２、８４３における計測対象点７１に対応する補正された対応点の座標を出力する。最適画像座標推定部８５２は、対応点の座標の補正量または該補正量の逆数を変数とする所定の多項式に基づき、補正された対応点の座標が内部パラメータ８４４と外部パラメータ８４５とで構成されるエピポーラ方程式を満たす補正ベクトルを計算し、補正された対応点の座標を出力する。
　最適画像座標推定部８５２は、固有多項式計算部８５２１と補正ベクトル計算部８５２２とを有する。固有多項式計算部８５２１は、対応点の座標の補正量または該補正量の逆数を固有値とする補正ベクトルが満たすべき拘束条件を表す所定の拘束条件行列を計算し、その拘束条件行列の固有方程式を解いて最小の補正量を算出する機能を有する。なお、所定の拘束条件行列は、対応点検出部８５１から入力した対応点の座標と記憶部８４から入力した内部パラメータ８４４と外部パラメータ８４５とに基づいて計算される。また補正ベクトル計算部８５２２は、上記最小の補正量に対応する上記拘束条件行列の固有ベクトルを計算し、上記固有ベクトルに所定のスケーリングを施して補正ベクトルを計算し、上記補正ベクトルから補正された対応点の座標を計算し出力する機能を有する。なお、拘束条件行列の固有ベクトルは、対応点検出部８５１から入力した対応点の座標と固有多項式計算部８５２１で計算された拘束条件行列と最小の補正量とに基づいて計算される。
　三次元座標計算部８５３は、補正された対応点の座標を最適画像座標推定部８５２から入力し、内部パラメータ８４４と外部パラメータ８４５とを記憶部８４から入力し、計測対象点７１の三次元座標を計算する機能を有する。また、三次元座標計算部８５３は、計算した計測対象点７１の三次元座標を画面表示部８３に出力し、また通信Ｉ／Ｆ部８１を通じて図示しない他の端末へ送信する機能を有する。
　図１３は、本実施形態に係る三角測量装置８の動作の一例を示すフローチャートである。以下、図１３を参照して、本実施形態に係る三角測量装置８の動作について説明する。
　三角測量装置４の対応点検出部８５１は、記憶部８４から読み出した画像８４２を、通信Ｉ／Ｆ部８１を通じてプロジェクタ９に送信し、プロジェクタ９によって画像８４２を構造物７に投影する（ステップＳ３１）。次に対応点検出部８５１は、構造物７に投影された画像８４２を含む画像をカメラ６により撮影し、撮影して得られた画像を通信Ｉ／Ｆ部８１でカメラ６から受信し、記憶部８４に画像８４３として記憶する（ステップＳ３２）。次に、対応点検出部８５１は、画像８４２、８４３中の計測対象点７１に対応する対応点の画像座標を検出し、最適画像座標推定部８５２に伝達する（ステップＳ３３）。以降、図１１に示したステップＳ２３～Ｓ２８と同様の処理が行われる（ステップＳ３４～Ｓ３９）。
　このように本実施形態によれば、アクティブステレオ法を使用して計測対象点の三次元位置を算出する三角測量装置において、大域的に最適な三次元位置を少ない演算量で計算することができる。その理由は、第１の実施形態の効果における説明と同様である。例えば、６次方程式ではなく３次方程式を解けばよく、また３次方程式の解から再投影誤差が大域的に最小となる解を選択できる。
　以上、実施形態および実施例を参照して本発明を説明したが、本発明は上記実施形態および実施例に限定されるものではない。本発明の構成や詳細には、本発明のスコープ内で当業者が理解し得る様々な変更をすることができる。
　本発明によれば、複数画像からの三次元形状復元といった用途に適用できる。この出願は、２０１４年９月１８日に出願された日本出願特願２０１４－１９０２４２を基礎とする優先権を主張し、その開示の全てをここに取り込む。
　上記の実施形態の一部又は全部は、以下の付記のようにも記載され得るが、以下には限られない。
［付記１］
　計測対象点の三次元位置を、ステレオ法を使用して算出する三角測量装置であって、
　前記計測対象点の画像を有する２枚の画像と前記画像を生成した光学機器の内部パラメータおよび外部パラメータとを記憶する記憶部と、前記記憶部に接続されたプロセッサとを有し、
　前記プロセッサは、
　前記２枚の画像における前記計測対象点に対応する対応点の座標と前記内部パラメータと前記外部パラメータとを入力し、補正された前記対応点の座標が前記内部パラメータと前記外部パラメータとで構成されるエピポーラ方程式を満たす補正ベクトルを、前記対応点の座標の補正量または該補正量の逆数を変数とする３次多項式に基づき計算し、補正された前記対応点の座標を算出し、
　前記補正された対応点の座標と前記内部パラメータと前記外部パラメータとを入力とし、前記計測対象点の三次元座標を計算する
ようにプログラミングされている三角測量装置。
［付記２］
　前記プロセッサは、
　前記対応点の座標と前記内部パラメータと前記外部パラメータとを入力し、前記対応点の座標の補正量または該補正量の逆数を固有値とする補正ベクトルが満たすべき拘束条件を表す所定の拘束条件行列を計算し、前記拘束条件行列の固有方程式を解いて最小の補正量を算出し、
　前記対応点の座標と前記拘束条件行列と前記最小の補正量とを入力とし、前記最小の補正量に対応する前記拘束条件行列の固有ベクトルを計算し、前記固有ベクトルに所定のスケーリングを施して補正ベクトルを計算し、前記補正ベクトルから補正された前記対応点の座標を計算する
ようにプログラミングされている付記１に記載の三角測量装置。
［付記３］
　前記対応点の補正量は、補正ベクトルの大きさを目的関数とし、定数倍の不定性を持つエピポーラ方程式を制約条件とする制約付き非線形最適化問題におけるラグランジュ未定乗数を表す
付記１または２に記載の三角測量装置。
［付記４］
　前記ステレオ法は、それぞれ異なる位置から前記計測対象点の２つの画像をカメラによって撮影するパッシブステレオ法である
付記１乃至３の何れかに記載の三角測量装置。
［付記５］
　前記ステレオ法は、ある位置から前記計測対象点の画像を投影するプロジェクタと別の位置から前記計測対象点の画像をカメラにより撮影するアクティブステレオ法である
付記１乃至３の何れかに記載の三角測量装置。
［付記６］
　前記所定の拘束条件行列を次式、
Figure JPOXMLDOC01-appb-M000001
により計算する付記２に記載の三角測量装置。
［付記７］
　前記計測対象点の画像を有する２枚の画像と前記画像を生成した光学機器の内部パラメータおよび外部パラメータとを記憶する記憶部と、前記記憶部に接続されたプロセッサとを有し、前記計測対象点の三次元位置を、ステレオ法を使用して算出する三角測量装置が実行する三角測量方法であって、
　前記プロセッサが、
　前記２枚の画像における前記計測対象点に対応する対応点の座標と前記内部パラメータと前記外部パラメータとを入力し、補正された前記対応点の座標が前記内部パラメータと前記外部パラメータとで構成されるエピポーラ方程式を満たす補正ベクトルを、前記対応点の座標の補正量または該補正量の逆数を変数とする３次多項式に基づき計算し、補正された前記対応点の座標を算出し、
　前記補正された対応点の座標と前記内部パラメータと前記外部パラメータとを入力とし、前記計測対象点の三次元座標を計算する
三角測量方法。
［付記８］
　前記プロセッサは、
　前記対応点の座標と前記内部パラメータと前記外部パラメータとを入力し、前記対応点の座標の補正量または該補正量の逆数を固有値とする補正ベクトルが満たすべき拘束条件を表す所定の拘束条件行列を計算し、前記拘束条件行列の固有方程式を解いて最小の補正量を算出し、
　前記対応点の座標と前記拘束条件行列と前記最小の補正量とを入力とし、前記最小の補正量に対応する前記拘束条件行列の固有ベクトルを計算し、前記固有ベクトルに所定のスケーリングを施して補正ベクトルを計算し、前記補正ベクトルから補正された前記対応点の座標を計算する
付記７に記載の三角測量方法。
［付記９］
　前記対応点の補正量は、補正ベクトルの大きさを目的関数とし、定数倍の不定性を持つエピポーラ方程式を制約条件とする制約付き非線形最適化問題におけるラグランジュ未定乗数を表す
付記７または８に記載の三角測量方法。
［付記１０］
　前記ステレオ法は、それぞれ異なる位置から前記計測対象点の２つの画像をカメラによって撮影するパッシブステレオ法である
付記７乃至９の何れかに記載の三角測量方法。
［付記１１］
　前記ステレオ法は、ある位置から前記計測対象点の画像を投影するプロジェクタと別の位置から前記計測対象点の画像をカメラにより撮影するアクティブステレオ法である
付記７乃至９の何れかに記載の三角測量方法。
［付記１２］
　前記所定の拘束条件行列を次式、
Figure JPOXMLDOC01-appb-M000002
により計算する付記８に記載の三角測量方法。
［付記１３］
　計測対象点の画像を有する２枚の画像と前記画像を生成した光学機器の内部パラメータおよび外部パラメータとを記憶する記憶部と、前記記憶部に接続されたプロセッサとを有し、前記計測対象点の三次元位置を、ステレオ法を使用して算出する三角測量装置の前記プロセッサに、
　前記２枚の画像における前記計測対象点に対応する対応点の座標と前記内部パラメータと前記外部パラメータとを入力し、補正された前記対応点の座標が前記内部パラメータと前記外部パラメータとで構成されるエピポーラ方程式を満たす補正ベクトルを、前記対応点の座標の補正量または該補正量の逆数を変数とする３次多項式に基づき計算し、補正された対応点の座標を算出し、
　前記補正された対応点の座標と前記内部パラメータと前記外部パラメータとを入力とし、前記計測対象点の三次元座標を計算し
を行わせるためのプログラム。
［付記１４］
　前記プロセッサは、
　前記対応点の座標と前記内部パラメータと前記外部パラメータとを入力し、前記対応点の座標の補正量または該補正量の逆数を固有値とする補正ベクトルが満たすべき拘束条件を表す所定の拘束条件行列を計算し、前記拘束条件行列の固有方程式を解いて最小の補正量を計算し、
　前記対応点の座標と前記拘束条件行列と前記最小の補正量とを入力とし、前記最小の補正量に対応する前記拘束条件行列の固有ベクトルを計算し、前記固有ベクトルに所定のスケーリングを施して補正ベクトルを計算し、前記補正ベクトルから補正された前記対応点の座標を計算する
付記１３に記載のプログラム。
［付記１５］
　前記対応点の補正量は、補正ベクトルの大きさを目的関数とし、定数倍の不定性を持つエピポーラ方程式を制約条件とする制約付き非線形最適化問題におけるラグランジュ未定乗数を表す
付記１３または１４に記載のプログラム。
［付記１６］
　前記ステレオ法は、それぞれ異なる位置から前記計測対象点の２つの画像をカメラによって撮影するパッシブステレオ法である
付記１３乃至１５の何れかに記載のプログラム。
［付記１７］
　前記ステレオ法は、ある位置から前記計測対象点の画像を投影するプロジェクタと別の位置から前記計測対象点の画像をカメラにより撮影するアクティブステレオ法である
付記１３乃至１５の何れかに記載のプログラム。
［付記１８］
　前記所定の拘束条件行列を次式、
Figure JPOXMLDOC01-appb-M000003
により計算する付記１４に記載のプログラム。
１…最適対応点座標推定部
１１…固有多項式計算部
１２…補正ベクトル計算部
２…三次元座標計算部
３１…プロセッサ
３２…プログラムメモリ
３３…記憶媒体

Claims (10)

  1. 　ステレオ法を使用して計測対象点の三次元位置を算出する三角測量装置であって、
    　前記計測対象点の画像を有する２枚の画像における前記計測対象点に対応する対応点の座標と前記２枚の画像を生成した光学機器の内部パラメータと外部パラメータとに基づき、前記対応点の座標の補正量または該補正量の逆数を変数とする固有多項式を用いて、補正された前記対応点の座標が前記内部パラメータと前記外部パラメータとで構成されるエピポーラ方程式を満たす補正ベクトルを計算し、計算した補正ベクトルに基づいて、補正された前記対応点の座標を算出する最適画像座標推定手段と、
    　前記補正された対応点の座標と前記内部パラメータと前記外部パラメータとに基づいて、前記計測対象点の三次元座標を計算する三次元座標計算手段と、
    を有する三角測量装置。
  2. 　前記最適画像座標推定手段は、
    　前記対応点の座標と前記内部パラメータと前記外部パラメータとを入力し、前記対応点の座標の補正量または該補正量の逆数を固有値とする補正ベクトルが満たすべき拘束条件を表す所定の拘束条件行列を計算し、前記拘束条件行列の固有方程式を解いて最小の補正量を算出する固有多項式計算手段と、
    　前記対応点の座標と前記拘束条件行列と前記最小の補正量とを入力とし、前記最小の補正量に対応する前記拘束条件行列の固有ベクトルを計算し、前記固有ベクトルに所定のスケーリングを施して補正ベクトルを計算し、前記補正ベクトルに基づいて補正された前記対応点の座標を計算する補正ベクトル計算手段と
    を有する請求項１に記載の三角測量装置。
  3. 　前記対応点の補正量は、補正ベクトルの大きさを目的関数とし、定数倍の不定性を持つエピポーラ方程式を制約条件とする制約付き非線形最適化問題におけるラグランジュ未定乗数を表す
    請求項１または２に記載の三角測量装置。
  4. 　前記ステレオ法は、それぞれ異なる位置から前記計測対象点の２つの画像をカメラによって撮影するパッシブステレオ法である
    請求項１乃至３の何れか１つに記載の三角測量装置。
  5. 　前記ステレオ法は、ある位置から前記計測対象点の画像を投影するプロジェクタと別の位置から前記計測対象点の画像をカメラにより撮影するアクティブステレオ法である
    請求項１乃至３の何れか１つに記載の三角測量装置。
  6. 　三角測量装置が、
    　２つの異なる位置から計測対象を撮影した画像における計測対象点に対応する対応点の座標と、前記２つの画像を撮影した光学機器の内部パラメータ及び外部パラメータとを入力し、
    　前記対応点の座標の補正量または該補正量の逆数を変数とする固有多項式に基づいて、補正された前記対応点の座標が前記内部パラメータと前記外部パラメータとで構成されるエピポーラ方程式を満たす補正ベクトルを計算し、
    　計算した前記補正ベクトルに基づき、補正された前記対応点の座標を算出し、
    　前記補正された対応点の座標と前記内部パラメータと前記外部パラメータとを入力し、前記計測対象点の三次元座標を計算する、
    三角測量方法。
  7. 　前記三角測量装置が、
    　前記対応点の座標と前記内部パラメータと前記外部パラメータとを入力し、前記対応点の座標の補正量または該補正量の逆数を固有値とする補正ベクトルが満たすべき拘束条件を表す所定の拘束条件行列を計算し、前記拘束条件行列の固有方程式を解いて最小の補正量を算出し、
    　前記対応点の座標と前記拘束条件行列と前記最小の補正量とに基づいて、前記最小の補正量に対応する前記拘束条件行列の固有ベクトルを計算し、前記固有ベクトルに所定のスケーリングを施して補正ベクトルを計算し、前記補正ベクトルに基づいて補正された前記対応点の座標を計算する
    請求項６に記載の三角測量方法。
  8. 　前記対応点の補正量は、補正ベクトルの大きさを目的関数とし、定数倍の不定性を持つエピポーラ方程式を制約条件とする制約付き非線形最適化問題におけるラグランジュ未定乗数を表す
    請求項６または７に記載の三角測量方法。
  9. 　コンピュータを、
    　前記計測対象点の画像を有する２枚の画像における前記計測対象点に対応する対応点の座標と前記２枚の画像を生成した光学機器の内部パラメータと外部パラメータとを入力し、前記対応点の座標の補正量または該補正量の逆数を変数とする固有多項式に基づいて、補正された前記対応点の座標が前記内部パラメータと前記外部パラメータとで構成されるエピポーラ方程式を満たす補正ベクトルを計算し、計算した補正ベクトルに基づいて補正された前記対応点の座標を算出する最適画像座標推定手段と、
    　前記補正された対応点の座標と前記内部パラメータと前記外部パラメータとに基づいて、前記計測対象点の三次元座標を計算する三次元座標計算手段と
    して機能させるためのプログラムを記録した記録媒体。
  10. 　前記最適画像座標推定手段は、
    　前記対応点の座標と前記内部パラメータと前記外部パラメータとに基づき、前記対応点の座標の補正量または該補正量の逆数を固有値とする補正ベクトルが満たすべき拘束条件を表す所定の拘束条件行列を計算し、前記拘束条件行列の固有方程式を解いて最小の補正量を算出する固有多項式計算手段と、
    　前記対応点の座標と前記拘束条件行列と前記最小の補正量とに基づいて、前記最小の補正量に対応する前記拘束条件行列の固有ベクトルを計算し、前記固有ベクトルに所定のスケーリングを施して補正ベクトルを計算し、前記補正ベクトルに基づいて補正された前記対応点の座標を計算する補正ベクトル計算手段と
    を有する請求項９に記載のプログラムを記録した記録媒体。
PCT/JP2015/004798 2014-09-18 2015-09-18 三角測量装置、三角測量方法およびそのプログラムを記録した記録媒体 WO2016042779A1 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
US15/510,737 US10126115B2 (en) 2014-09-18 2015-09-18 Triangulation device, triangulation method, and recording medium recording program therefor
JP2016548565A JP6547754B2 (ja) 2014-09-18 2015-09-18 三角測量装置、三角測量方法およびそのプログラム

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2014-190242 2014-09-18
JP2014190242 2014-09-18

Publications (1)

Publication Number Publication Date
WO2016042779A1 true WO2016042779A1 (ja) 2016-03-24

Family

ID=55532831

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2015/004798 WO2016042779A1 (ja) 2014-09-18 2015-09-18 三角測量装置、三角測量方法およびそのプログラムを記録した記録媒体

Country Status (3)

Country Link
US (1) US10126115B2 (ja)
JP (1) JP6547754B2 (ja)
WO (1) WO2016042779A1 (ja)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106871787A (zh) * 2017-04-13 2017-06-20 中国航空工业集团公司北京长城航空测控技术研究所 大空间线扫描成像三维测量方法
JP2018146363A (ja) * 2017-03-03 2018-09-20 三菱航空機株式会社 ３次元位置計測システム及び方法
CN110146071A (zh) * 2019-05-16 2019-08-20 成都精位科技有限公司 定位维度提升方法及装置
CN115690334A (zh) * 2023-01-03 2023-02-03 中国人民解放军63921部队 一种基于共面约束的三维重建整体优化方法
WO2023084601A1 (ja) * 2021-11-09 2023-05-19 日本電気株式会社 算出装置、算出方法、カメラパラメータ推定装置、及び非一時的なコンピュータ可読媒体

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10546385B2 (en) * 2016-02-25 2020-01-28 Technion Research & Development Foundation Limited System and method for image capture device pose estimation
US10225535B2 (en) * 2016-08-11 2019-03-05 Rabin Esrail Self-adjusting portable modular 360-degree projection and recording computer system
JP2020085596A (ja) * 2018-11-21 2020-06-04 三菱重工業株式会社 位置測定システム及び位置測定方法
WO2020208686A1 (ja) * 2019-04-08 2020-10-15 日本電気株式会社 カメラ校正装置、カメラ校正方法、及びプログラムが格納された非一時的なコンピュータ可読媒体
CN111504276B (zh) * 2020-04-30 2022-04-19 哈尔滨博觉科技有限公司 一种基于视觉投影尺度因子集的联合目标函数多推进器姿态角获取方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2014013147A (ja) * 2012-07-03 2014-01-23 Canon Inc ３次元計測装置及びロボット装置
JP2014096761A (ja) * 2012-11-12 2014-05-22 Canon Inc 画像処理装置、その制御方法、および制御プログラム

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090238426A1 (en) * 2008-03-19 2009-09-24 Uti Limited Partnership System and Methods for Identifying an Object within a Complex Environment
US9342886B2 (en) * 2011-04-29 2016-05-17 Qualcomm Incorporated Devices, methods, and apparatuses for homography evaluation involving a mobile device

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2014013147A (ja) * 2012-07-03 2014-01-23 Canon Inc ３次元計測装置及びロボット装置
JP2014096761A (ja) * 2012-11-12 2014-05-22 Canon Inc 画像処理装置、その制御方法、および制御プログラム

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
KEN'ICHI KANATANI ET AL.: "Triangulation from Two Views : Hartley-Sturm vs. Optimal Correction", IEICE TECHNICAL REPORT, vol. 107, no. 539, 3 March 2008 (2008-03-03), pages 335 - 342 *
KENTO YAMADA ET AL.: "Latest Algorithm for 3-D Reconstruction from Two Views", IEICE TECHNICAL REPORT, vol. 109, no. 182, 24 August 2009 (2009-08-24), pages 95 - 102 *
RICHARD I.HARTLEY: "Computer Vision and Image Understanding", TRIANGULATION, vol. 68, no. 2, December 1997 (1997-12-01), pages 146 - 157 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2018146363A (ja) * 2017-03-03 2018-09-20 三菱航空機株式会社 ３次元位置計測システム及び方法
CN106871787A (zh) * 2017-04-13 2017-06-20 中国航空工业集团公司北京长城航空测控技术研究所 大空间线扫描成像三维测量方法
CN106871787B (zh) * 2017-04-13 2019-02-22 中国航空工业集团公司北京长城航空测控技术研究所 大空间线扫描成像三维测量方法
CN110146071A (zh) * 2019-05-16 2019-08-20 成都精位科技有限公司 定位维度提升方法及装置
WO2023084601A1 (ja) * 2021-11-09 2023-05-19 日本電気株式会社 算出装置、算出方法、カメラパラメータ推定装置、及び非一時的なコンピュータ可読媒体
CN115690334A (zh) * 2023-01-03 2023-02-03 中国人民解放军63921部队 一种基于共面约束的三维重建整体优化方法

Also Published As

Publication number Publication date
US10126115B2 (en) 2018-11-13
JPWO2016042779A1 (ja) 2017-08-03
JP6547754B2 (ja) 2019-07-24
US20170254637A1 (en) 2017-09-07

Similar Documents

Publication Publication Date Title
WO2016042779A1 (ja) 三角測量装置、三角測量方法およびそのプログラムを記録した記録媒体
US9759548B2 (en) Image processing apparatus, projector and projector system including image processing apparatus, image processing method
JP6079333B2 (ja) 校正装置、方法及びプログラム
US9639942B2 (en) Information processing apparatus, information processing method, and storage medium
US9595106B2 (en) Calibration apparatus, calibration method, and program
US11488354B2 (en) Information processing apparatus and information processing method
KR20080029080A (ko) 단안 줌 카메라를 이용한 이동로봇의 자기위치 추정 시스템및 방법
JP6332951B2 (ja) 画像処理装置および画像処理方法、およびプログラム
JPWO2018168255A1 (ja) カメラパラメータ推定装置、カメラパラメータ推定方法、及びプログラム
JP7173285B2 (ja) カメラ校正装置、カメラ校正方法、及びプログラム
JP6641729B2 (ja) ラインセンサカメラのキャリブレーション装置及び方法
JP6694281B2 (ja) ステレオカメラおよび撮像システム
US11145072B2 (en) Methods, devices and computer program products for 3D mapping and pose estimation of 3D images
JP7044016B2 (ja) ラインセンサカメラのキャリブレーション装置及び方法
JP6573196B2 (ja) 距離情報補正装置、距離情報補正方法及び距離情報補正プログラム
JP6065670B2 (ja) ３次元計測システム、プログラム及び方法。
JP7121936B2 (ja) カメラ校正情報取得装置、画像処理装置、カメラ校正情報取得方法およびプログラム
Bingwei et al. Accurate geometric camera calibration technique using multi-views of a non-metric planar grid
WO2020144757A1 (ja) カメラパラメータ推定装置、カメラパラメータ推定方法、及びコンピュータ読み取り可能な記録媒体
JP5728399B2 (ja) 計測装置、方法及びプログラム
Pai et al. High-fidelity camera-based method for noncontact vibration testing of structures
WO2019058487A1 (ja) ３次元復元画像処理装置、３次元復元画像処理方法及び３次元復元画像処理プログラムを記憶したコンピュータ読み取り可能な記憶媒体
EP3136340A1 (en) Method for estimating a unified projection model for a camera and apparatus performing the same
KR101239671B1 (ko) 렌즈에 의한 왜곡 영상 보정방법 및 그 장치
JP2017163386A (ja) カメラパラメータ推定装置、カメラパラメータ推定方法、及びプログラム

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 15841276

Country of ref document: EP

Kind code of ref document: A1

WWE Wipo information: entry into national phase

Ref document number: 15510737

Country of ref document: US

ENP Entry into the national phase

Ref document number: 2016548565

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 15841276

Country of ref document: EP

Kind code of ref document: A1