WO2019058487A1 - ３次元復元画像処理装置、３次元復元画像処理方法及び３次元復元画像処理プログラムを記憶したコンピュータ読み取り可能な記憶媒体 - Google Patents

３次元復元画像処理装置、３次元復元画像処理方法及び３次元復元画像処理プログラムを記憶したコンピュータ読み取り可能な記憶媒体 Download PDF

Info

Publication number
WO2019058487A1
WO2019058487A1 PCT/JP2017/034141 JP2017034141W WO2019058487A1 WO 2019058487 A1 WO2019058487 A1 WO 2019058487A1 JP 2017034141 W JP2017034141 W JP 2017034141W WO 2019058487 A1 WO2019058487 A1 WO 2019058487A1
Authority
WO
WIPO (PCT)
Prior art keywords
images
image
dimensional
image processing
calculated
Prior art date
Application number
PCT/JP2017/034141
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 PCT/JP2017/034141 priority Critical patent/WO2019058487A1/ja
Publication of WO2019058487A1 publication Critical patent/WO2019058487A1/ja
Priority to US16/746,006 priority patent/US11176702B2/en

Links

Images

Classifications

    • 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
    • G06T7/74—Determining position or orientation of objects or cameras using feature-based methods involving reference images or patches
    • 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/579—Depth or shape recovery from multiple images from motion
    • 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
    • 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/10016—Video; Image sequence
    • 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/10068—Endoscopic image
    • 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/30—Subject of image; Context of image processing
    • G06T2207/30244—Camera pose

Definitions

  • the present invention relates to a three-dimensional restored image processing apparatus, a three-dimensional restored image processing method, and a three-dimensional restored image processing program for restoring a three-dimensional shape of an object from a plurality of images obtained by imaging the object from different viewpoints.
  • a computer readable storage medium storing
  • SfM Structure from Motion
  • an imaging unit that picks up an image with three-dimensional coordinates of a feature point that is consistent with the coordinates of the observed feature point from the coordinates of a common feature point of a subject in a plurality of images by a technique called bundle adjustment
  • the posture of the is calculated.
  • the difference between the coordinates at which a feature point appears on the image and the coordinates on the image of the feature point actually observed when the feature point at the estimated three-dimensional coordinates is imaged in the posture of an imaging unit Estimating the three-dimensional coordinates of the feature point and the attitude of the imaging unit are performed so that the cost function that increases with the size decreases.
  • the cost function is a non-linear function related to the three-dimensional coordinates of the feature point and the attitude of the imaging unit, and it is difficult to analytically find a solution that minimizes the cost function. Therefore, the gradient method is to search for a solution that minimizes the cost function.
  • the posture of the imaging unit in the reference coordinate system is calculated as an initial value by combining a plurality of relative postures calculated between two images. Methods are widely used.
  • the relative orientation between the two images is calculated from the coordinates of the feature points associated between the two images.
  • an 8-point algorithm, a 5-point algorithm, etc. are known.
  • the eight-point algorithm is an algorithm for obtaining a basic matrix (F matrix) using coordinates of eight sets of feature points common to two images, and calculating a relative attitude from the F matrix.
  • the five-point algorithm obtains a basic matrix (E matrix) using coordinates of five sets of feature points common to two images and camera internal parameters representing the characteristics of the imaging unit, and calculates a relative attitude from the E matrix Algorithm.
  • the relative attitude may not be uniquely calculated only from the coordinates of the feature point. That is, when the object is a plane, in the 5-point algorithm, one of the two relative postures obtained by solving the equation may be a relative posture which does not match the actual camera operation. As described above, when bundle adjustment is performed using the initial value obtained from the relative posture calculated from the two images in which the relative posture can not be calculated uniquely, the processing time until the result by the gradient method converges will be long. And may converge to an incorrect value.
  • An object of the present invention is to provide an image processing apparatus, a three-dimensional restored image processing method, and a computer readable storage medium storing a three-dimensional restored image processing program.
  • a three-dimensional restored image processing apparatus includes an acquiring unit for acquiring a plurality of images captured from different viewpoints, and an image pair of two images from the acquired plurality of images. And an attitude calculation unit for calculating a relative attitude between the two images constituting each of the selected image pairs, and the relative attitude between the two images constituting the selected image pair being Based on the determination result of the uniqueness determination unit that determines whether or not it can calculate uniquely, and the uniqueness determination unit, the relative attitude used to calculate the initial value of the bundle adjustment is determined, and the determined relative attitude is determined. And a bundle adjustment unit that performs bundle adjustment from the initial value of the bundle adjustment calculated based on the calculated bundle adjustment value.
  • a three-dimensional restored image processing method includes acquiring a plurality of images captured from different viewpoints, and selecting two images from the acquired plurality of images as an image pair. Calculating the relative attitude between the two images constituting each of the selected image pair, and uniquely calculating the relative attitude between the two images constituting the selected image pair And the relative attitude used to calculate the initial value of the bundle adjustment based on the result of the determination, and the initial value of the bundle adjustment calculated based on the determined relative attitude And performing more bundle adjustments.
  • a computer-readable storage medium storing the three-dimensional restored image processing program according to the third aspect of the present invention comprises: acquiring a plurality of images captured from different viewpoints; and 2 from the acquired plurality of images Selecting images one by one as an image pair, calculating a relative attitude between the two images constituting each of the selected image pairs, and between the two images constituting the selected image pair Determining whether or not the relative posture can be uniquely calculated, and determining the relative posture used to calculate the initial value of the bundle adjustment based on the result of the determination, and based on the determined relative posture
  • a three-dimensional restored image processing program is stored that causes a computer to execute bundle adjustment from the calculated initial value of bundle adjustment.
  • FIG. 1 is a block diagram showing the configuration of an imaging apparatus according to each embodiment of the present invention.
  • FIG. 2 is a block diagram showing the configuration of the three-dimensional restored image processing apparatus in the first embodiment.
  • FIG. 3 is a flowchart showing the operation of the three-dimensional restored image processing apparatus in the first embodiment.
  • FIG. 4 is a diagram showing an example of selection of an image pair.
  • FIG. 5 is a diagram showing the concept of a transformation matrix.
  • FIG. 6A is a diagram showing the relationship between the imaging unit and two images when the relative attitude can be calculated uniquely.
  • FIG. 6B is a diagram showing the relationship between the imaging unit and two images when the relative attitude can not be calculated uniquely.
  • FIG. 6A is a diagram showing the relationship between the imaging unit and two images when the relative attitude can not be calculated uniquely.
  • FIG. 7 is a flowchart showing the operation of the three-dimensional restored image processing apparatus of the first modification of the first embodiment.
  • FIG. 8 is a block diagram showing the configuration of a three-dimensional restored image processing apparatus according to the second embodiment.
  • FIG. 9 is a flowchart showing the operation of the three-dimensional restored image processing apparatus according to the second embodiment.
  • FIG. 10 is a flowchart showing the operation of the three-dimensional restored image processing apparatus of the modified example 1 of the second embodiment.
  • FIG. 1 is a block diagram showing the configuration of an imaging apparatus according to each embodiment of the present invention.
  • the imaging device 1 illustrated in FIG. 1 includes a CPU 2, a RAM 3, a ROM 4, an imaging unit 5, a display unit 6, and a three-dimensional restored image processing device 7.
  • the imaging apparatus 1 is, for example, an endoscope apparatus having an imaging function.
  • the imaging device 1 may not be an endoscope device, as long as it is a device having an imaging function, such as a digital camera, a webcam, a smartphone, and an on-vehicle camera.
  • the CPU 2 is a control circuit that controls each part of the imaging device 1.
  • the CPU 2 performs control in accordance with the program stored in the ROM 4.
  • the CPU 2 controls an operation of imaging by the imaging unit 5 or controls an operation of display on the display unit 6.
  • the CPU 2 also inputs camera internal parameters indicating the characteristics of the imaging unit 5 to the three-dimensional restored image processing device 7.
  • the camera internal parameters include the focal length of the optical system of the imaging unit 5 and the coordinates of the image center which is the coordinates at which the optical axis of the optical system intersects with the imaging surface of the imaging device.
  • FIG. 1 shows the CPU, the control circuit of the imaging device 1 may be configured by an ASIC, an FPGA, or the like.
  • the RAM 3 temporarily stores various data such as data calculated by the CPU 2, an image obtained by the imaging unit 5, and a processing result of the three-dimensional restored image processing apparatus 7.
  • the RAM 3 is configured of volatile memory such as DRAM and SDRAM.
  • the ROM 4 stores various data such as programs executed by the CPU 2 and camera internal parameters.
  • the ROM 4 is configured by a non-volatile memory such as a flash ROM.
  • the ROM 4 may be replaced with another storage medium such as a hard disk drive (HDD) and a solid state drive (SSD).
  • HDD hard disk drive
  • SSD solid state drive
  • the imaging unit 5 includes an optical system and an imaging device, and captures an image of a subject (not shown) to obtain an image of the subject.
  • the optical system forms an image of a subject (not shown) on the imaging surface of the imaging device.
  • the optical system may include a focus lens and a zoom lens.
  • the imaging device is, for example, a CMOS image sensor, and converts an image of a subject formed on an imaging surface into an image signal (hereinafter, simply referred to as an image) which is an electrical signal.
  • the image generated by the imaging device is stored, for example, in the RAM 3.
  • the display unit 6 displays various images such as an image obtained by the imaging unit 5.
  • the display unit 6 is a liquid crystal display, an organic EL display, or the like.
  • the three-dimensional restored image processing device 7 performs processing to restore the three-dimensional shape of the subject from the image obtained by the imaging unit 5, for example.
  • the three-dimensional restored image processing device 7 is configured by a processor such as a DSP.
  • the three-dimensional restored image processing apparatus 7 will be further described with reference to FIG.
  • FIG. 2 is a block diagram showing the configuration of the three-dimensional restored image processing apparatus 7 in the first embodiment.
  • the three-dimensional restored image processing device 7 includes an acquisition unit 101, a posture calculation unit 102, a uniqueness determination unit 104, a bundle adjustment unit 105, and a three-dimensional coordinate calculation unit 106.
  • the functions of these blocks are realized by, for example, software processing.
  • the acquisition unit 101, the posture calculation unit 102, the uniqueness determination unit 104, the bundle adjustment unit 105, and the three-dimensional coordinate calculation unit 106 may be configured by hardware.
  • the acquisition unit 101 acquires images of a plurality of viewpoints. Then, the acquisition unit 101 inputs the acquired image to the posture calculation unit 102.
  • the images of a plurality of viewpoints acquired by the acquisition unit 101 may be, for example, moving images obtained by continuously imaging the same subject while changing the posture of the imaging unit 5 or the same subject It may be a plurality of individual still images obtained by imaging a plurality of times with different orientations of the imaging unit 5. Furthermore, the acquisition unit 101 is not limited to one that acquires an image directly from the imaging unit 5.
  • the acquisition unit 101 stores the captured image in the storage device It may be obtained from In such a configuration, the three-dimensional restored image processing device 7 may not necessarily be mounted on the imaging device 1, and may be, for example, a personal computer or the like that is communicatively connected to the imaging device 1.
  • the posture calculation unit 102 selects two of the plurality of images input from the acquisition unit 101 as an image pair, and the posture of the imaging unit 5 when capturing the first image forming the selected image pair
  • the calculation of the relative attitude representing the change in the attitude of the imaging unit 5 when the second image is captured is repeated.
  • the relative attitude is calculated by a five-point algorithm.
  • the five-point algorithm finds the basic matrix (E matrix) using the coordinates of five sets of feature points (corresponding feature points) common to two images and camera internal parameters, and calculates the relative attitude from this E matrix Algorithm.
  • the camera internal parameters include the focal length of the optical system of the imaging unit 5 and the coordinates of the image center which is the coordinates at which the optical axis of the optical system intersects with the imaging surface of the imaging device.
  • the uniqueness determination unit 104 determines the relative posture of the image pair selected by the posture calculation unit 102 from the distribution of the feature points included in the two images constituting the image pair selected by the posture calculation unit 102 and the relative posture. It is determined whether or not it can be uniquely identified. Although details will be described later, when the subject is a plane, the relative posture may not be identified uniquely depending on the selection of the image pair. The uniqueness determination unit 104 determines, for each image pair, whether it is possible to uniquely identify such a relative posture. Then, the uniqueness determination unit 104 sends information indicating whether or not the relative posture can be uniquely identified to the bundle adjustment unit 105 together with the relative posture calculated by the posture calculation unit 102.
  • the bundle adjustment unit 105 performs bundle adjustment for estimating the three-dimensional coordinates of the feature points in the image and the posture of the imaging unit 5 using the relative posture calculated for each image pair by the posture calculation unit 102. At this time, the bundle adjustment unit 105 is used to calculate the initial value of the bundle adjustment for the relative orientation calculated from the image pair determined to be uniquely identifiable by the uniqueness determination unit 104. The relative posture calculated from the image pair determined not to uniquely identify the posture is not used to calculate the initial value of the bundle adjustment. Details of bundle adjustment will be described later.
  • the three-dimensional coordinate calculation unit 106 calculates three-dimensional coordinates (three-dimensional shape of the subject) of all points in the image from the three-dimensional coordinates of the feature points calculated by the bundle adjustment unit 105 and the camera posture of the imaging unit 5 .
  • FIG. 3 is a flowchart showing the operation of the three-dimensional restored image processing apparatus 7 in the first embodiment.
  • step S101 the acquisition unit 101 acquires, for example, a plurality of images captured by the imaging unit 5 and stored in the RAM 3, and inputs the acquired plurality of images to the posture calculation unit 102.
  • the plurality of images acquired by the acquisition unit 101 may be, for example, moving images obtained by continuously imaging the same subject while making the orientation of the imaging unit 5 different. It may be a plurality of still images not continuous in time obtained by imaging the same subject a plurality of times by changing the posture of the imaging unit 5 differently.
  • step S102 the posture calculation unit 102 selects two of the plurality of input images as an image pair.
  • the combination of the two images selected as the image pair may be different from the combinations in the other image pairs. However, combinations that are considered to have no correlation between the two images (combinations in which there is a high possibility that there is no corresponding feature point between the two images) may be excluded from the selection.
  • FIG. 4 An example of selecting an image pair will be described with reference to FIG.
  • four images of a key frame 1, a key frame 2, a key frame 3, and a key frame 4 are input to the attitude calculation unit 102 as a plurality of images.
  • Key frames 1, 2, 3, and 4 are images obtained by imaging in this order.
  • the combination of image pairs in the case of four sheets is a combination of key frames 1 and 2, a combination of key frames 1 and 3, a combination of key frames 1 and 4, a combination of key frames 2 and 3, a combination of key frames 2 and 4 , And 6 combinations of key frames 3 and 4.
  • the key frame 1 and the key frame 4 may be excluded from the selection of the image pair because the time lag is large.
  • step S103 the posture calculation unit 102 calculates the relative posture between the two images from the coordinates of the corresponding feature point between the two images.
  • the process of step S103 will be described.
  • the relative attitude a conversion matrix that represents a change in the coordinate system of the imaging unit 5 that has captured two images is calculated.
  • FIG. 5 is a diagram showing the concept of a transformation matrix. As shown in FIG.
  • the transformation matrix for an image pair including X is a matrix that transforms the coordinate system of the image imaged at coordinate O into the coordinate system of the image imaged at coordinate O ′.
  • the transformation matrix represents a 3 ⁇ 3 rotation matrix R that represents rotation of the coordinate system of the imaging unit 5 and a translation vector t that is a three-dimensional vector that represents translational movement of the coordinate system of the imaging unit 5.
  • a method of calculating a transformation matrix there are known an 8-point algorithm, a 5-point algorithm, and the like which solve equations by linear solution from coordinates of corresponding feature points (corresponding feature points) between two images.
  • the 8-point algorithm can not solve the equations because the number of independent equations is smaller than the number of unknowns of the equations. Therefore, in the present embodiment, a five-point algorithm is used which is a method capable of calculating the relative posture between two images even when the subject is a plane.
  • a conversion matrix representing a relative posture between two images is calculated by solving a non-linear equation created from the coordinates of five sets of corresponding feature points.
  • a 3 ⁇ 3 matrix called E matrix determined by the relative posture between two images is calculated, and further, the E matrix is subjected to singular value decomposition to calculate a transformation matrix.
  • the E matrix calculated from the combination of key frames 1 and 2 is E 12
  • the E matrix calculated from the combination of key frames 1 and 3 is calculated from the combination of E 13 and key frames 2 and 3.
  • the E matrix calculated from the combination of E 23 and key frames 2 and 4 is shown as E 24 and the E matrix calculated from the combination of key frames 3 and 4 as E 34 .
  • the coordinates of the corresponding feature points may be extracted by tracking, for each frame, easily trackable feature points extracted in a certain frame in the moving image. it can.
  • a method of extracting feature points in an image methods such as FAST corner detection and Harris corner detection are known.
  • a method of tracking feature points in a moving image a block matching method, Lukas-Kanade method, etc. are known.
  • the plurality of images acquired by the acquisition unit 101 are not temporally continuous images, the coordinates of the corresponding feature points are acquired for each image, and feature points that are easy to be associated with are obtained between the acquired feature points. The similarity can be extracted by comparing for each image.
  • Scale-Invariant Feature Transform SIFT
  • SURF Speeded Up Robust Features
  • RANSAC adopts a relative posture in which the number of feature points whose errors are equal to or less than a threshold is the smallest.
  • LMedS a relative attitude with the smallest median value of errors is adopted.
  • these errors for example, a Sampson error is used.
  • step S107 after the relative attitude is calculated, the uniqueness determination unit 104 determines the uniqueness to determine whether or not the relative attitude between the two images selected as the image pair can be calculated uniquely. I do.
  • the uniqueness determination in step S107 will be described.
  • a transformation matrix H representing this projective transformation is a rotation matrix R from a coordinate system of an image imaged at a plane normal vector n, a distance to the plane d, and a coordinate O to a coordinate system of an image imaged at a coordinate O '.
  • the translation vector t is a rotation matrix R from a coordinate system of an image imaged at a plane normal vector n, a distance to the plane d, and a coordinate O to a coordinate system of an image imaged at a coordinate O '.
  • the uniqueness determination unit 104 determines, out of two images, the relative orientation between the two images calculated by the orientation calculation unit 102 and the three-dimensional coordinates of the feature point calculated using the relative orientation.
  • the number n1 of feature points whose three-dimensional distance is close to the optical center when the first image is captured and the number n2 of feature points whose three-dimensional distance is near the optical center when the second image is captured Ask for Then, uniqueness determination unit 104 determines that the posture change can be calculated uniquely when min (n1, n2) / max (n1, n2) is equal to or greater than the threshold value, that is, the difference between n1 and n2 is small.
  • min (n1, n2) is the smaller value of n1 and n2
  • max (n1, n2) is the larger value of n1 and n2.
  • DLT Direct Linear Transformation
  • the midpoint method can be used.
  • the three-dimensional coordinates of the feature point based on the coordinate system of the first image is a three-dimensional vector X
  • the three-dimensional coordinates of the feature point based on the coordinate system of the second image is a three-dimensional vector
  • Two-dimensional vectors x ' 1 and x' 2 for coordinates of feature points in X ', the first image and the second image, f for the focal length of the optical system of the imaging unit 5 for imaging the object
  • the coordinate O 1 of the optical center in the first image the second image
  • the number n1 and n2 of feature points can be calculated by calculating the distance between the three-dimensional coordinates X and the optical center coordinates O 1 and O 2 in the coordinate system of the first image of all feature points. It is calculated. Uniqueness determination is performed as described above from the numbers n1 and n2 of these feature points.
  • step S108 it is determined whether or not the bundle adjustment unit 105 can uniquely calculate the relative orientation between the two images selected as the image pair by the uniqueness determination unit 104. Determine When it is determined in step S108 that the relative attitude between the two images selected as the image pair can be uniquely calculated, the process proceeds to step S109. If it is determined in step S108 that the relative orientation between the two images selected as the image pair can not be calculated uniquely, the process proceeds to step S110. That is, in the present embodiment, when it is determined that the relative posture between the two images selected as the image pair can not be calculated uniquely, the calculation of the initial value described in step S109 is not performed.
  • step S109 the bundle adjustment unit 105 calculates an image posture in the reference coordinate system as an initial value of bundle adjustment using the relative posture between the two images selected as the image pair.
  • an example of the calculation method of the initial value will be described.
  • the image attitude in the reference coordinate system as an initial value includes an initial value R j 0 of a rotation matrix representing the rotation of the image j with respect to the reference coordinate system, and a translation vector t j 0 representing translation of the image j with respect to the reference coordinate system. .
  • the Jbase-th coordinate system is first selected as a reference coordinate image.
  • the initial value R Jbase 0 of the rotation matrix of the reference coordinate image and the initial value t Jbase 0 of the translation vector are respectively set as follows.
  • R j 0 is calculated as a rotation matrix of image j with respect to R jbase 0
  • t j 0 as a translation vector of image j with respect to t jbase 0 .
  • the image orientation in the reference coordinate system of the image J whose orientation in the reference coordinate system is known is R J 0 , t J 0
  • the relative orientation of the image j with respect to the image J is R J, j , t J , j (that is, between the three-dimensional coordinates P ′ ij of the feature point i in the image j and the three-dimensional coordinates P ′ iJ of the feature point i in the image J
  • P ′ ij R J, j T
  • the initial value P i 0 of the three-dimensional coordinates of the feature point with respect to the reference coordinate system is, for example, the same DLT as when the three-dimensional coordinates of the feature point are calculated from the relative posture between two images in the uniqueness determination described above. Calculated using the method.
  • M be the number of images for which the relative coordinates with respect to the reference coordinate system have been calculated be M
  • the coordinates in the image j be x ' and ij
  • the three-dimensional coordinates of the feature point i in the coordinate system of the image j and X ij R j 0T (X i -t j 0).
  • the following two-dimensional vector x ij is defined from the focal length f of the optical system and the coordinates a of the image center.
  • step S110 after the calculation of the initial value, the bundle adjustment unit 105 determines whether or not the image orientation in the reference coordinate system as the initial value of the bundle adjustment has been calculated for all the images. If it is determined in step S110 that the initial values for all the images have not been calculated, the process returns to step S102. In this case, the calculation of the relative posture, the determination of the uniqueness, and the calculation of the initial value described above are performed for another image pair. When it is determined in step S110 that the initial values for all the images have been calculated, the process proceeds to step S111.
  • step S111 the bundle adjustment unit 105 performs bundle adjustment to calculate the three-dimensional coordinates of all feature points and the posture of the imaging unit 5 at the time of capturing an image.
  • the attitude of the imaging unit 5 includes the three-dimensional coordinates of the optical center, the direction of the optical axis, and the rotation angle of the imaging element around the optical axis direction. The following describes bundle adjustment.
  • the bundle adjustment is a process of calculating three-dimensional coordinates of feature points and an attitude of an imaging unit which are in good agreement with the coordinates of the observed feature points on the image.
  • three-dimensional coordinates of a feature point that minimizes the cost C represented by (Expression 8) below and the attitude of the imaging unit are calculated. Therefore, assuming that the three-dimensional coordinates of the feature point that minimizes the cost C is P I m and the orientation of the imaging unit is R J m and t J m , P I m , R J m and t J m are It is calculated from equation 9).
  • P i in (Expression 9) is a three-dimensional vector representing three-dimensional coordinates in the reference coordinate system of the i-th feature point (feature point i).
  • t j is a translation vector representing coordinate conversion between the coordinate system of the imaging unit 5 and the reference coordinate system when the j-th image (image j) is captured.
  • R j is a rotation matrix representing coordinate conversion between the coordinate system of the imaging unit 5 and the reference coordinate system when the j-th image (image j) is captured, and 3 ⁇ the determinant is 1 It is a 3-orthogonal matrix.
  • q ij is a two-dimensional vector obtained by normalizing the coordinates of the feature point i in the image j.
  • q ij is a normalized coordinate in the image j of the feature point i which can be expressed by the following (Expression 10).
  • p ij (P i , R j , t j ) should be reflected in the image j captured in the coordinate system in which the feature point i at the three-dimensional coordinate P i is indicated by t j and R j Coordinates of p ij (P i , R j , t j ) can be expressed by the following (Expression 11).
  • q ij -p ij (P i , R j , t j ) is the feature point i on the image j
  • the normalized coordinates q ij detected from the image and the orientations of the image j are t j and R j
  • This q ij -p ij (P i , R j , t j ) is called a reprojection error.
  • the function ⁇ ( ⁇ ) of (Equation 8) is a function with the error ⁇ as an argument.
  • the function ⁇ ( ⁇ ) it is possible to use, for example, a function of a square error shown in the following (Equation 12) or a Huber loss function shown in the following (Equation 13).
  • the function ⁇ ( ⁇ ) for calculating the cost C is integrated for all feature points i for which coordinates on the image have been calculated, in all the images j for which the postures R j and t j have been calculated.
  • the cost C is a nonlinear function of the P i, R j, t j, it is difficult to obtain P i that minimizes the cost C, R j, a t j analytically. Therefore, in the present embodiment, P i , R j and t j which minimize the cost C are obtained by the gradient method from the above-described initial values P i 0 , R j 0 and t j 0 .
  • the gradient method, P i in the n th step, R j is t j P i n, R j n, from t j n P i, R j, when the minutely varying the t j cost C update amount [Delta] P i n obtained from a change in, ⁇ R j n, ⁇ t n based on the j (n + 1)
  • P i in th step, R j is t j P i n + 1
  • the initial value of bundle adjustment using the relative orientation calculated from the image pair There is no calculation of The accuracy of bundle adjustment can be improved by not performing the calculation of the initial value using such an ambiguous relative attitude.
  • the three-dimensional restored image processing apparatus After bundle adjustment, the three-dimensional restored image processing apparatus outputs a bundle adjustment result. Thereafter, the imaging device 1 records the bundle adjustment result in the RAM 3 or outputs the result to the display unit 6.
  • the process of FIG. 3 ends.
  • calculation of three-dimensional coordinates (three-dimensional shape of the subject) of all points in the image by the three-dimensional coordinate calculation unit 106 may be performed.
  • P I m , R J m and t J m obtained as a result of bundle adjustment are used to calculate three-dimensional coordinates of all points in this image.
  • a three-dimensional image of the restored object is displayed on the display unit 6 as necessary.
  • the calculation of the initial value described in step S109 is a row. I can not do it. As a result, even when the subject is flat, bundle adjustment can be performed in a short time and with high accuracy.
  • FIG. 7 is a flowchart showing the operation of the three-dimensional restored image processing apparatus 7 of the first modification.
  • an initial value is calculated each time an image pair is selected.
  • the bundle adjustment is performed after initial values for all the images that can uniquely calculate the relative posture are calculated.
  • the process proceeds to step S111 every time the initial value of a certain image pair is calculated in step S109, and the initial values of all the images whose initial values have been calculated Bundle adjustment is performed using. Thereafter, the process proceeds to step S112, and the process ends when bundle adjustment using all the images determined to be able to uniquely calculate the relative attitude is performed.
  • Bundle adjustment may be performed sequentially as in the first modification.
  • an image pair is selected among all the images other than ones among all the images, and the relative orientation is unique among the N (N-1) / 2 image pairs created for N images.
  • the postures R j 0 and t j 0 of the images relative to the reference coordinate system as initial values and the three-dimensional coordinates P i 0 of the feature points may be calculated from the relative postures of the image pairs determined to be able to be calculated.
  • the posture calculation unit 102 sets a solution (solution 1) in which the subject (feature point) is considered to be in front of the imaging unit 5 among the four solutions, as the final relative posture.
  • the posture calculation unit 102 may calculate two types of relative postures. When the subject is a plane and the relative posture can not be determined uniquely, that is, when all feature points for both solution 1 and solution 2 are in front of the imaging unit 5, these two solutions are used.
  • the uniqueness determination it is determined that the relative attitude can not be calculated uniquely regardless of which relative attitude is used. Therefore, even when the attitude calculation unit 102 calculates the relative attitude based on two solutions of the solution 1 and the solution 2, the uniqueness determination unit 104 uses only one of the relative attitudes calculated by the attitude calculation unit 102. Uniqueness determination may be performed.
  • FIG. 8 is a block diagram showing the configuration of the three-dimensional restored image processing apparatus 7 in the second embodiment.
  • the same components as in FIG. 2 are denoted by the same reference numerals as in FIG. 2.
  • the description of the components given the same reference numerals is omitted.
  • the three-dimensional restored image processing device 7 of the second embodiment further includes a plane determination processing unit 103.
  • the plane determination processing unit 103 is provided between the posture calculation unit 102 and the uniqueness determination unit 104, and determines whether the subject is a plane.
  • uniqueness determination is performed based only on the distribution of corresponding feature points between two images regardless of the shape of the subject.
  • the relative posture can be uniquely calculated regardless of the distribution of corresponding feature points between the two images. Therefore, in the second embodiment, it is determined whether or not the subject is a plane prior to the uniqueness determination, and the uniqueness determination is performed only when the subject is a plane.
  • FIG. 9 is a flowchart showing the operation of the three-dimensional restored image processing apparatus 7 in the second embodiment.
  • the same processing as that of FIG. 3 is given the same reference numeral as that of FIG. The description of the processes given the same reference numerals will be omitted. That is, in FIG. 9, after the relative orientation between the two images is calculated in step S103, the process proceeds to step S105.
  • the plane determination processing unit 103 determines whether the subject is a plane. For example, the plane determination processing unit 103 determines whether the subject is a plane based on whether the correspondence of coordinates of feature points between two images can be represented by a relation of projective transformation.
  • the projective transformation matrix H between the images is obtained from the coordinates of the corresponding feature points in the first image and the second image, and the coordinates of the corresponding feature points in the first image are transformed with the projective transformation matrix H It is possible to determine whether the object is a plane by comparing the error with the threshold value, using the difference between the coordinates of the second point and the coordinates of the second corresponding feature point as the error of projection conversion.
  • a threshold value for determining whether the subject is a plane or not from the error of projection conversion it is possible to use a GRIC (Geometric Robust Information Criterion) value.
  • GRIC Garnier Robust Information Criterion
  • step S106 the uniqueness determination unit 104 determines whether the subject is a plane based on the information from the plane determination processing unit 103. If it is determined in step S106 that the subject is flat, the process proceeds to step S107. Thereafter, the same processing as in the first embodiment is performed. On the other hand, when it is determined in step S106 that the subject is not flat, the process proceeds to step S109. That is, the uniqueness determination is omitted.
  • unnecessary determination of uniqueness can be omitted by determining whether or not the subject is a plane prior to determination of uniqueness, leading to efficient processing. .
  • FIG. 10 is a flowchart showing the operation of the three-dimensional restored image processing apparatus 7 of the first modification.
  • the first modification is a process that can be applied when a plurality of images acquired by the acquisition unit 101 are images obtained by imaging the same subject so as to have a common visual field.
  • the result of plane determination for one image pair can also be used for plane determination for another image pair. That is, plane determination may be performed only once.
  • step S104 the plane determination processing unit 103 determines whether the plane determination of the subject has been performed. Then, in step S104, the process proceeds to step S105 only when it is determined that the subject's plane determination has not been performed. When it is determined in step S104 that the plane determination of the subject has been performed, the process of step S105 is skipped.
  • the plane determination is performed only once. This leads to more efficient processing.
  • three-dimensional coordinates of the subject can be calculated by the DLT method. That is, three-dimensional coordinates for all corresponding feature points whose coordinates are commonly calculated between two images are calculated, and a first principal component vector and a second principal component vector obtained by PCA (principal component decomposition) of three-dimensional coordinates are stretched. By obtaining the surface, it is possible to determine the main plane on which the subject is stretched. When the ratio of the number of feature points whose main plane three-dimensional distance is equal to or less than the threshold value is sufficiently large among all corresponding feature points, it is assumed that all corresponding feature points are on the plane. At this time, it is determined that the subject is a plane.
  • the threshold value may be determined on the basis of the optical center of the optical system that has captured the image and the distance to the subject, or the threshold value for the first principal component or the second principal component calculated by PCA. It may be determined on the basis of the unique value.

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Image Analysis (AREA)
  • Image Processing (AREA)

Abstract

３次元画像処理装置（７）は、異なる視点で撮像された複数枚の画像を取得する取得部（１０１）と、取得された複数枚の画像から２枚の画像を画像ペアとして選択し、選択した画像ペアを構成する２枚の画像の間の相対姿勢を算出する姿勢算出部（１０２）と、選択された画像ペアを構成する２枚の画像の間の相対姿勢が一意に算出できるか否かを判断する一意性判断部（１０４）と、一意性判断部（１０４）の判断結果に基づいてバンドル調整の初期値の算出に用いる相対姿勢を決定し、決定された相対姿勢に基づいて算出されたバンドル調整の初期値よりバンドル調整を行うバンドル調整部（１０５）とを備える。

Description

３次元復元画像処理装置、３次元復元画像処理方法及び３次元復元画像処理プログラムを記憶したコンピュータ読み取り可能な記憶媒体
　本発明は、被写体を異なる複数の視点から撮像して得られた複数枚の画像から被写体の３次元形状を復元する３次元復元画像処理装置、３次元復元画像処理方法及び３次元復元画像処理プログラムを記憶したコンピュータ読み取り可能な記憶媒体に関する。
　被写体を異なる複数の視点から撮像して得られた複数枚の画像から、画像を撮像した撮像部の姿勢と被写体の３次元形状とを復元するSfM(Structure from Motion)という技術が知られている。SfMでは、複数枚の画像における被写体の共通の特徴点の座標から、バンドル調整という技術により、観測された特徴点の座標と整合性が取れた特徴点の３次元座標と画像を撮像した撮像部の姿勢とが算出される。SfMに関する技術として、特開２０１５－１１１８２２号公報、D. Nister, “An Efficient Solution to the Five-Point Relative Pose Problem,” TPAMI (2004)、J.L. Schonberger, and J.-M. Frahm, "Structure-from-Motion Revisited," CVPR （2016）及びD. Martinec and T. Pajdla, “Robust Rotation and Translation Estimation in Multiview Reconstruction,” CVPR (2007)において開示されている技術が知られている。
　バンドル調整では、推定された３次元座標にある特徴点をある撮像部の姿勢で撮像したときに特徴点が画像上に現れる座標と実際に観測された特徴点の画像上の座標との差の大きさにつれて大きくなるコスト関数が小さくなるよう、特徴点の３次元座標と撮像部の姿勢を推定することが行われる。コスト関数は、特徴点の３次元座標と撮像部の姿勢とに関する非線形関数であり、コスト関数が最小化する解を解析的に求めることは困難である。そこで、勾配法によりコスト関数が最小化する解を探索することになる。しかしながら、勾配法を行うためには、ある基準座標系での特徴点の３次元座標と撮像部の姿勢についての初期値が必要である。コスト関数は複雑な形状をしており、高精度な３次元復元を行うためには、正しい解に十分に近い初期値から勾配法の処理を開始する必要がある。
　バンドル調整の初期値を算出する手法としては、２枚の画像間で算出した相対姿勢を複数組み合わせることで、全画像の共通の座標系である基準座標系における撮像部の姿勢を初期値として算出する手法が広く行われている。２枚の画像間の相対姿勢は、２枚の画像間で対応付けた特徴点の座標から算出される。相対姿勢を算出するためのアルゴリズムとして、８点アルゴリズム及び５点アルゴリズム等が知られている。８点アルゴリズムは、２枚の画像間で共通する８組の特徴点の座標を用いて基礎行列（Ｆ行列）を求め、Ｆ行列から相対姿勢を算出するアルゴリズムである。５点アルゴリズムは、２枚の画像間で共通する５組の特徴点の座標と撮像部の特性を表すカメラ内部パラメータとを用いて基本行列（Ｅ行列）を求め、Ｅ行列から相対姿勢を算出するアルゴリズムである。
　ここで、被写体が平面であって画面上の全特徴点が同一の平面に乗っている場合、８点アルゴリズムでは相対姿勢を算出することができないため５点アルゴリズムを使用する必要がある。ただし、５点アルゴリズムを使用した場合でも、被写体が平面である場合は、原理的に、特徴点の座標だけからでは相対姿勢を一意に算出できないことがある。つまり、被写体が平面であるときには、５点アルゴリズムでは方程式を解いて得られる２通りの相対姿勢のうちの片方の相対姿勢が実際のカメラ動作と一致しない相対姿勢であることが起こりえる。このように、相対姿勢を一意に算出できない２枚の画像によって算出した相対姿勢から求めた初期値を用いてバンドル調整が行われると、勾配法による結果が収束するまでの処理時間が長くなったり、誤った値に収束したりする可能性がある。
　本発明は、前記の事情に鑑みてなされたものであり、被写体が平面であるときでも、短時間で、かつ、精度よくバンドル調整を行うことができる初期値を算出することができる３次元復元画像処理装置、３次元復元画像処理方法及び３次元復元画像処理プログラムを記憶したコンピュータ読み取り可能な記憶媒体を提供することを目的とする。
　本発明の第１の態様の３次元復元画像処理装置は、異なる視点で撮像された複数枚の画像を取得する取得部と、取得された前記複数枚の画像から２枚ずつの画像を画像ペアとして選択し、選択した前記画像ペアそれぞれを構成する２枚の画像の間の相対姿勢を算出する姿勢算出部と、選択された前記画像ペアを構成する２枚の画像の間の前記相対姿勢が一意に算出できるか否かを判断する一意性判断部と、前記一意性判断部の判断結果に基づいてバンドル調整の初期値の算出に用いる前記相対姿勢を決定し、決定された前記相対姿勢に基づいて算出された前記バンドル調整の初期値よりバンドル調整を行うバンドル調整部とを具備する。
　本発明の第２の態様の３次元復元画像処理方法は、異なる視点で撮像された複数枚の画像を取得することと、取得した前記複数枚の画像から２枚ずつの画像を画像ペアとして選択し、選択した前記画像ペアそれぞれを構成する２枚の画像の間の相対姿勢を算出することと、選択された前記画像ペアを構成する２枚の画像の間の前記相対姿勢が一意に算出できるか否かを判断することと、前記判断の結果に基づいてバンドル調整の初期値の算出に用いる前記相対姿勢を決定し、決定された前記相対姿勢に基づいて算出された前記バンドル調整の初期値よりバンドル調整を行うこととを具備する。
　本発明の第３の態様の３次元復元画像処理プログラムを記憶したコンピュータ読み取り可能な記憶媒体は、異なる視点で撮像された複数枚の画像を取得することと、取得した前記複数枚の画像から２枚ずつの画像を画像ペアとして選択し、選択した前記画像ペアそれぞれを構成する２枚の画像の間の相対姿勢を算出することと、選択された前記画像ペアを構成する２枚の画像の間の前記相対姿勢が一意に算出できるか否かを判断することと、前記判断の結果に基づいてバンドル調整の初期値の算出に用いる前記相対姿勢を決定し、決定された前記相対姿勢に基づいて算出された前記バンドル調整の初期値よりバンドル調整を行うこととをコンピュータに実行させる３次元復元画像処理プログラムを記憶している。
図１は、本発明の各実施形態に係る撮像装置の構成を示すブロック図である。 図２は、第１の実施形態における３次元復元画像処理装置の構成を示すブロック図である。 図３は、第１の実施形態における３次元復元画像処理装置の動作を示すフローチャートである。 図４は、画像ペアの選択例を示す図である。 図５は、変換行列の概念を示す図である。 図６Ａは、相対姿勢を一意に算出できるときの撮像部と２枚の画像との関係を示す図である。 図６Ｂは、相対姿勢を一意に算出できないときの撮像部と２枚の画像との関係を示す図である。 図７は、第１の実施形態の変形例１の３次元復元画像処理装置の動作を示すフローチャートである。 図８は、第２の実施形態における３次元復元画像処理装置の構成を示すブロック図である。 図９は、第２の実施形態における３次元復元画像処理装置の動作を示すフローチャートである。 図１０は、第２の実施形態の変形例１の３次元復元画像処理装置の動作を示すフローチャートである。
　以下、図面を参照して本発明の実施形態を説明する。　
　［第１の実施形態］
　まず、本発明の第１の実施形態を説明する。図１は、本発明の各実施形態に係る撮像装置の構成を示すブロック図である。図１に示す撮像装置１は、ＣＰＵ２と、ＲＡＭ３と、ＲＯＭ４と、撮像部５と、表示部６と、３次元復元画像処理装置７とを有している。この撮像装置１は、例えば撮像機能を有する内視鏡装置である。撮像装置１は、例えば、デジタルカメラ、ウェブカム、スマートフォン、車載カメラ等、撮像機能を有する装置であれば、内視鏡装置でなくてもよい。
　ＣＰＵ２は、撮像装置１の各部の制御をする制御回路である。ＣＰＵ２は、ＲＯＭ４に記録されているプログラムに従って制御を行う。例えば、ＣＰＵ２は、撮像部５による撮像の動作を制御したり、表示部６における表示の動作を制御したりする。また、ＣＰＵ２は、３次元復元画像処理装置７に撮像部５の特性を示すカメラ内部パラメータを入力する。カメラ内部パラメータは、撮像部５の光学系の焦点距離と、光学系の光軸が撮像素子の撮像面と交わる座標である画像中心の座標とを含む。ここで、図１では、ＣＰＵであるが、撮像装置１の制御回路は、ＡＳＩＣ、ＦＰＧＡ等で構成されてもよい。
　ＲＡＭ３は、ＣＰＵ２で計算されたデータ、撮像部５で得られた画像、３次元復元画像処理装置７における処理結果といった各種のデータを一時的に記憶する。ＲＡＭ３は、例えばＤＲＡＭ、ＳＤＲＡＭといった揮発性メモリによって構成される。
　ＲＯＭ４は、ＣＰＵ２によって実行されるプログラム、カメラ内部パラメータといった各種のデータを記憶している。ＲＯＭ４は、フラッシュＲＯＭ等の不揮発性メモリによって構成される。ＲＯＭ４は、ハードディスクドライブ（ＨＤＤ）及びソリッドステートドライブ（ＳＳＤ）といった他の記憶媒体に置き替えられてもよい。
　撮像部５は、光学系と撮像素子とを含み、図示しない被写体を撮像して被写体に係る画像を取得する。光学系は、図示しない被写体の像を撮像素子の撮像面に結像させる。光学系は、フォーカスレンズ、ズームレンズを含んでいてもよい。撮像素子は、例えばＣＭＯＳイメージセンサであり、撮像面に結像した被写体の像を電気信号である画像信号（以下、単に画像という）に変換する。撮像素子で生成された画像は、例えばＲＡＭ３に記憶される。
　表示部６は、撮像部５で得られた画像等の各種の画像を表示する。表示部６は、液晶ディスプレイ及び有機ＥＬディスプレイ等である。
　３次元復元画像処理装置７は、例えば撮像部５で得られた画像から被写体の３次元形状を復元する処理を行う。３次元復元画像処理装置７は、例えばＤＳＰといったプロセッサによって構成される。以下、図２を参照して３次元復元画像処理装置７についてさらに説明する。図２は、第１の実施形態における３次元復元画像処理装置７の構成を示すブロック図である。図２に示すように、３次元復元画像処理装置７は、取得部１０１と、姿勢算出部１０２と、一意性判断部１０４と、バンドル調整部１０５と、３次元座標算出部１０６とを有する。これらのブロックの機能は、例えばソフトウェア処理によって実現される。勿論、取得部１０１と、姿勢算出部１０２と、一意性判断部１０４と、バンドル調整部１０５と、３次元座標算出部１０６は、ハードウェアによって構成されてもよい。
　取得部１０１は、複数の視点の画像を取得する。そして、取得部１０１は、取得した画像を姿勢算出部１０２に入力する。取得部１０１によって取得される複数の視点の画像は、例えば同一の被写体を撮像部５の姿勢を異ならせながら連続的に撮像することによって得られた動画像であってもよいし、同一の被写体を撮像部５の姿勢を異ならせて複数回の撮像をすることによって得られた個別の複数の静止画であってもよい。また、取得部１０１は、撮像部５から直接的に画像を取得するものに限らない。例えば、撮像装置１の内部の記憶装置（例えばＲＯＭ４）又は撮像装置１の外部の記憶装置に撮影済みの画像が記憶されている場合には、取得部１０１は、その撮影済みの画像を記憶装置から取得してもよい。このような構成の場合、３次元復元画像処理装置７は、必ずしも撮像装置１に搭載されていなくてもよく、例えば、撮像装置１と通信接続して使用するパソコン等でもよい。
　姿勢算出部１０２は、取得部１０１から入力された複数の画像のうちの２枚を画像ペアとして選択し、選択した画像ペアを構成する１枚目の画像を撮像したときの撮像部５の姿勢と２枚目の画像を撮像したときの撮像部５の姿勢との変化を表す相対姿勢を算出することを繰り返す。詳細は後で説明するが、本実施形態では、相対姿勢は、５点アルゴリズムによって算出される。５点アルゴリズムは、２枚の画像間で共通の５組の特徴点（対応特徴点）の座標とカメラ内部パラメータとを用いて基本行列（Ｅ行列）を求め、このＥ行列から相対姿勢を算出するアルゴリズムである。前述したように、カメラ内部パラメータは、撮像部５の光学系の焦点距離と、光学系の光軸が撮像素子の撮像面と交わる座標である画像中心の座標とを含む。
　一意性判断部１０４は、姿勢算出部１０２で選択された画像ペアを構成する２枚の画像に含まれる特徴点の分布と相対姿勢とから、姿勢算出部１０２で選択された画像ペアで相対姿勢を一意に特定できるか否かを判断する。詳細は、後で説明するが、被写体が平面である場合、画像ペアの選択の仕方によっては相対姿勢を一意に特定できないことがある。一意性判断部１０４は、このような相対姿勢を一意に特定できるか否かの判断を画像ペア毎に行う。そして、一意性判断部１０４は、相対姿勢を一意に特定できるか否かを示す情報を姿勢算出部１０２で算出された相対姿勢とともにバンドル調整部１０５に送る。
　バンドル調整部１０５は、姿勢算出部１０２で画像ペア毎に算出された相対姿勢を用いて画像内の特徴点の３次元座標と撮像部５の姿勢とを推定するためのバンドル調整を行う。このとき、バンドル調整部１０５は、一意性判断部１０４で相対姿勢を一意に特定できると判断された画像ペアから算出された相対姿勢についてはバンドル調整の初期値を算出するために使用し、相対姿勢を一意に特定できないと判断された画像ペアから算出された相対姿勢についてはバンドル調整の初期値を算出するために使用しない。バンドル調整の詳細は後で説明する。
　３次元座標算出部１０６は、バンドル調整部１０５で算出された特徴点の３次元座標と撮像部５のカメラ姿勢とから画像内の全点の３次元座標（被写体の３次元形状）を算出する。
　次に、３次元復元画像処理装置７の動作を説明する。図３は、第１の実施形態における３次元復元画像処理装置７の動作を示すフローチャートである。
　ステップＳ１０１において、取得部１０１は、例えば撮像部５によって撮像されてＲＡＭ３に記憶された複数の画像を取得し、取得した複数の画像を姿勢算出部１０２に入力する。前述したように、取得部１０１によって取得される複数の画像は、例えば同一の被写体を撮像部５の姿勢を異ならせながら連続的に撮像することによって得られた動画像であってもよいし、同一の被写体を撮像部５の姿勢を異ならせて複数回撮像することによって得られた時間的に連続していない複数の静止画であってもよい。
　ステップＳ１０２において、姿勢算出部１０２は、入力された複数の画像のうちの２枚を画像ペアとして選択する。画像ペアとして選択する２枚の組み合わせは、他の画像ペアにおける組み合わせと異なっていればよい。ただし、２枚の画像間の相関がないと考えられる組み合わせ（２枚の画像間で対応する特徴点がない可能性が高い組み合わせ）は選択から除外するようにしてもよい。
　図４を参照して画像ペアの選択例を説明する。図４では、複数の画像としてキーフレーム１、キーフレーム２、キーフレーム３、キーフレーム４の４枚の画像が姿勢算出部１０２に入力されている。キーフレーム１、２、３、４はこの順の撮像によって得られた画像である。４枚の場合の画像ペアの組み合わせは、キーフレーム１、２の組み合わせ、キーフレーム１、３の組み合わせ、キーフレーム１、４の組み合わせ、キーフレーム２、３の組み合わせ、キーフレーム２、４の組み合わせ、キーフレーム３、４の組み合わせの６通りである。しかしながら、キーフレーム１とキーフレーム４とは時間的なずれが大きいために画像ペアの選択からは除外されてよい。
　ステップＳ１０３において、姿勢算出部１０２は、２枚の画像間で対応する特徴点の座標から２枚の画像間の相対姿勢を算出する。以下、ステップＳ１０３の処理を説明する。本実施形態では、相対姿勢として、２枚の画像を撮像した撮像部５の座標系の変化を表す変換行列を算出する。図５は、変換行列の概念を示す図である。図５に示すように、被写体の任意の点Pを座標Oから撮像部５によって撮像したときに得られる画像と被写体の点Pを座標O´から撮像部５によって撮像したときに得られる画像とを含む画像ペアについての変換行列は、座標Oで撮像した画像の座標系を座標O´で撮像した画像の座標系に変換する行列である。この変換行列は、撮像部５の座標系の回転を表す３×３回転行列Rと撮像部５の座標系の並進移動を表す３次元ベクトルである並進ベクトルtとを表す。
　変換行列を算出する手法としては、２枚の画像間で対応する特徴点（対応特徴点）の座標から線型解法により方程式を解く８点アルゴリズム及び５点アルゴリズム等が知られている。しかしながら、被写体が平面であるときには、８点アルゴリズムでは方程式の未知数の個数より独立な式の個数が小さくなるために方程式を解くことができない。そこで、本実施形態では、被写体が平面であるときにも２枚の画像間の相対的な姿勢を算出することができる手法である５点アルゴリズムが用いられる。５点アルゴリズムは、カメラ内部パラメータである、撮像部５の光学系の焦点距離fと、光学系の光軸が撮像素子の撮像面と交わる座標である画像中心の座標a=(a,b)Tとが既知であるときに用いることができる手法である。この５点アルゴリズムでは、５組の対応特徴点の座標から作成された非線形方程式を解くことにより、２枚の画像間の相対的な姿勢を表す変換行列が算出される。５点アルゴリズムでは、２枚の画像間の相対的な姿勢で決まるＥ行列とよばれる３×３行列を算出し、さらにＥ行列を特異値分解して変換行列を算出する。なお、図４において、キーフレーム１、２の組み合わせから算出されるＥ行列がE12、キーフレーム１、３の組み合わせから算出されるＥ行列がE13、キーフレーム２、３の組み合わせから算出されるＥ行列がE23、キーフレーム２、４の組み合わせから算出されるＥ行列がE24、キーフレーム３、４の組み合わせから算出されるＥ行列がE34として示されている。
　対応特徴点の座標は、取得部１０１によって取得された複数枚の画像が動画像であるときには、動画中のあるフレームで抽出した追跡しやすい特徴点をフレーム毎に追跡することで抽出することができる。画像の中の特徴点を抽出する手法としては、FASTコーナ検出及びHarrisコーナ検出といった方法が知られている。また、動画において特徴点を追跡する手法としては、ブロックマッチング法及びLukas-Kanade法等が知られている。また、対応特徴点の座標は、取得部１０１によって取得された複数枚の画像が時間的に連続した画像でないときには、対応を付けやすい特徴点を画像毎に取得し、取得した特徴点の間の類似度を画像毎に比較することで抽出することができる。対応を付けやすい特徴点としては、例えば、SIFT（Scale-Invariant Feature Transform）やSURF（Speeded Up Robust Features）がよく知られている。なお、対応特徴点は、自動的に抽出するだけに限らず、複数枚の画像の対応する点を作業者が１点ずつ指定することで抽出されてもよい。
　ここで、５点アルゴリズムのような対応特徴点の座標から作成された方程式を解くことで画像間の相対姿勢を算出する手法では、対応特徴点の座標の小さな誤差が、算出された相対姿勢の大きな誤差につながることがある。そのため、本実施形態では、相対姿勢の算出の際に、RANSAC（RANdom SAmple Consensus）及びLMedS（Least Median of Squares：最小メジアン法）といったロバスト推定、すなわち方程式を作成するのに必要な最小点数だけの対応特徴点を対応特徴点から抜き出して方程式を解く処理を繰り返し、全対応特徴点の座標を最もよく説明する相対姿勢を採用する処理が行われる。例えばRANSACでは、誤差がしきい値以下である特徴点の数が最も小さい相対姿勢を採用する。また、LMedSでは、誤差の２乗の中央値が最も小さい相対姿勢を採用する。これらの誤差としては、例えば、サンプソン誤差が使用される。
　ここで、図３の説明に戻る。相対姿勢が算出された後のステップＳ１０７において、一意性判断部１０４は、画像ペアとして選択された２枚の画像の間の相対姿勢を一意に算出できるか否かを判断するための一意性判断を行う。以下、ステップＳ１０７の一意性判断について説明する。
　D. Nister, “An Efficient Solution to the Five-Point Relative Pose Problem,” TPAMI (2004)には、被写体が平面であるときには、「画像上の全ての点が片方の画像を撮像した視点位置に近ければ」、相対姿勢を一意に算出できないことが開示されている。すなわち、画像ペアを構成する２枚の画像のうちの１枚目の画像を撮像した撮像部５の光学中心に近い特徴点と２枚目の画像を撮像した撮像部５の光学系の光学中心に近い特徴点の両方が２枚の画像の両方に映っており、かつ、それらの特徴点が対応特徴点であるときには、被写体が平面であっても相対姿勢を一意に算出することができる。この状態は、図６Ａに示すように、撮像部５の視点が座標Oから座標O´に移動したときに、OO´の垂直二等分面ＰＬの両側に対応特徴点がある状態と表現することもできる。一方で、全ての対応特徴点が１枚目の画像を撮像した撮像部５の光学系の光学中心と２枚目の画像を撮像した撮像部５の光学系の光学中心との何れかにだけ近いときには、相対姿勢を一意に算出できない、この状態は、図６Ｂに示すように、撮像部５の視点が座標Oから座標O´に移動したときに、OO´の垂直二等分面ＰＬの片側にだけ対応特徴点がある状態と表現することもできる。
　図５にも示すように、被写体が平面であるとき、２枚の画像の間の対応特徴点の間には射影変換の関係がある。この射影変換を表す変換行列Hは、平面の法線ベクトルnと、平面までの距離dと、座標Oで撮像した画像の座標系から座標O´で撮像した画像の座標系への回転行列Rと並進ベクトルtとから算出することができる。
　ここで、HからR、t、d、nを求めると、数式上、正しい解｛R,t,d,n｝が４通り算出される。４通りの解によって表される姿勢変化のうち、実際の姿勢変化は、被写体（特徴点）が撮像部５の前方にある姿勢変化である。このため、相対姿勢を一意に決めることができるのは、４通りの解が以下の関係になるときである。なお、（解２）、（解３）の状態において、撮像部５の前方にある特徴点と後方にある特徴点とは図６Ａで示したOO´の垂直二等分面ＰＬを境にして別れる。　
　　　（解１）：全特徴点が撮像部５の前方にある
　　　（解２）、（解３）:一部の特徴点が撮像部５の前方、一部の特徴点が撮像部５の後方にある
　　　（解４）:全特徴点が撮像部５の後方にある
　一方、４通りの解が以下の関係になるときには、相対姿勢を一意に決めることができない。これは、図６Ｂで示す状態しかないときである。この場合、（解１）と（解２）の何れが実際の姿勢変化を表しているかを数式だけでは特定することができない。　
　　　（解１）、（解２）：全特徴点が撮像部５の前方にある
　　　（解３）、（解４）：全特徴点が撮像部５の後方にある
　このような原理から、画像ペアを構成する２枚の画像のうちの１枚目の画像を撮像した撮像部５の光学中心に近い特徴点と２枚目の画像を撮像した撮像部５の光学系の光学中心に近い特徴点の両方が２枚の画像の両方に映っており、かつ、それらの特徴点が対応特徴点であるときには、被写体が平面であっても相対姿勢を一意に算出することができ、そうでないときには、相対姿勢を一意に算出することができないことが分かる。
　以下、一意性判断の具体的な一例を説明する。まず、一意性判断部１０４は、姿勢算出部１０２で算出された２枚の画像間の相対姿勢とその相対姿勢を用いて算出した特徴点の３次元座標とから、２枚の画像のうちの１枚目の画像が撮像されたときの光学中心に３次元距離が近い特徴点の数n1と２枚目の画像が撮像されたときの光学中心に３次元距離が近い特徴点の数n2とを求める。そして、一意性判断部１０４は、min(n1,n2)/max(n1,n2)がしきい値以上である、すなわちn1とn2との差が小さいとき、姿勢変化を一意に算出できると判断する。ここで、min(n1,n2)は、n1とn2のうちの小さい方の値であり、max(n1,n2)はn1とn2のうち大きい方の値である。
　なお、一意性判断に際して特徴点の３次元座標を求めるための手法として、２枚の画像間の相対姿勢と対応特徴点の画像上での座標とが既知であるときにはDLT（Direct Linear Transformation）法又は中点法を用いることができる。DLT法では、１枚目の画像の座標系を基準とした特徴点の３次元座標を３次元ベクトルX、２枚目の画像の座標系を基準とした特徴点の３次元座標を３次元ベクトルX’、１枚目の画像及び２枚目の画像での特徴点の座標をそれぞれ２次元ベクトルx´1及びx´2、被写体を撮像した撮像部５の光学系の焦点距離をf、画像中心の座標を２次元ベクトルaとし、行列式が1である３×３直交行列Rと、３次元ベクトルtとを用いてX’=RX＋tと表わすことができるとき、以下の条件式（式１）及び（式２）を満足するようにして特徴点の３次元座標を表す３次元ベクトルXが算出される。
Figure JPOXMLDOC01-appb-M000001
　以下、より具体的に説明する。３次元ベクトルXをX＝（X Y Z）Tとおくと、（式１）及び（式２）のzj×Xj＝0 (j=1,2)の左辺は３次元ベクトルになる。そして、zj×Xjによって表される３次元ベクトルの各要素は、X、Y、Zに関する１次の項と定数項とを持つ。このため、zj×Xjを以下の（式３）のように書き換えることができる。
Figure JPOXMLDOC01-appb-M000002
 （式１）及び（式２）より、（式３）によって書き換えられたベクトルは零ベクトル0に等しい。そこで、以下の（式４）で表される行列Aを定義する。
Figure JPOXMLDOC01-appb-M000003
 このとき、以下の（式５）の関係が成り立つ。
Figure JPOXMLDOC01-appb-M000004
 （式５）で示すように、（X Y Z 1）Tは４×４行列Aの零空間となるはずである。実際には、x´1及びx´2が誤差を含むときには、必ずしも行列Aに零空間は存在しない。この場合、行列Aの最も小さい特異値に対する右特異ベクトルが張る空間を行列Aの零空間と考え、最も小さい特異値に対する右特異ベクトルの第４成分が1となるように正規化することで（X Y Z 1）Tを得ることができる。
　１枚目の画像の座標系でXにある点を２枚目の画像の座標系でRX＋tと表すことができるときでは、１枚目の画像における光学中心の座標O1、２枚目の画像における光学中心の座標O2は１枚目の画像の座標系でそれぞれO1=0、O2=－RTtと表される。この関係から、全特徴点の１枚目の画像の座標系での３次元座標Xと光学中心の座標O1、座標O2との距離を算出することで、特徴点の数n1及びn2が算出される。これらの特徴点の数n1及びn2から前述したようにして一意性判断が行われる。
　ここで、図３の説明に戻る。一意性判断の後のステップＳ１０８において、バンドル調整部１０５は、一意性判断部１０４により、画像ペアとして選択された２枚の画像の間の相対姿勢を一意に算出できると判断されているか否かを判定する。ステップＳ１０８において、画像ペアとして選択された２枚の画像の間の相対姿勢を一意に算出できると判断されているときに、処理はステップＳ１０９に移行する。ステップＳ１０８において、画像ペアとして選択された２枚の画像の間の相対姿勢を一意に算出できないと判断されているときに、処理はステップＳ１１０に移行する。すなわち、本実施形態では、画像ペアとして選択された２枚の画像の間の相対姿勢を一意に算出できないと判断されているときには、ステップＳ１０９で説明する初期値の算出は行われない。
　ステップＳ１０９において、バンドル調整部１０５は、画像ペアとして選択された２枚の画像の間の相対姿勢を用いてバンドル調整の初期値としての基準座標系での画像姿勢を算出する。以下、初期値の算出手法の例を説明する。
　初期値としての基準座標系での画像姿勢は、基準座標系に対する画像jの回転を表す回転行列の初期値Rj 0、基準座標系に対する画像jの並進移動を表す並進ベクトルtj 0を含む。これらを算出するために、例えば、まず、Jbase番目の座標系が基準座標画像として選択される。基準座標画像の回転行列の初期値RJbase 0、並進ベクトルの初期値tJbase 0をそれぞれ以下のように設定する。Rj 0はRJbase 0に対する画像jの回転行列、tj 0はtjbase 0に対する画像jの並進ベクトルとして算出される。
Figure JPOXMLDOC01-appb-M000005
 なお、基準座標系での姿勢が既知である画像Jの、基準座標系での画像姿勢がRJ 0 、tJ 0であり、画像Jに対する画像jの相対姿勢がRJ,j、tJ,jであるとき(すなわち、特徴点iの画像jでの３次元座標P´ijと特徴点iの画像Jでの３次元座標P´iJとの間にP´ij＝RJ,j T(P´iJ－tJ,j)の関係があるとき)、基準座標系での画像jの姿勢Rj 0、tj 0は、それぞれ、Rj 0＝RJ,jRJ 0T、tj 0＝Rj 0tJ,j＋tj 0により算出される。
　基準座標系に対する特徴点の３次元座標の初期値Pi 0は、例えば、前述した一意性判断において２枚の画像の間の相対姿勢から特徴点の３次元座標を算出したときと同様のDLT法を用いて算出される。基準座標系に対する相対的な座標を算出した画像の枚数をM枚とし、i番目の特徴点の３次元座標をXi＝（Xi Yi Zi）T、画像jでの座標をx´ijとし、画像jの座標系での特徴点iの３次元座標をXij＝Rj 0T(Xi－tj 0)とする。また、光学系の焦点距離fと画像中心の座標aとから以下の２次元ベクトルxijを定義する。
Figure JPOXMLDOC01-appb-M000006
 ２次元ベクトルxijの第１要素を第１要素とし、２次元ベクトルxijの第２要素を第２要素とし、第３要素を1とする３次元ベクトルをzijとすると、以下の（式６）の関係を満たすaij,kl (k=1,2,3、l=1,2,3,4)が存在する。
Figure JPOXMLDOC01-appb-M000007
 （式６）のうち、k=1,2の aij,klを用いて以下の（式７）で表される行列Aiを定義する。そして、Ai (Xi Yi Zi 1)T = 0を満足するXi、Yi、Ziを求めることで算出されたXiをバンドル調整での特徴点iの３次元座標初期値Pi 0とする。
Figure JPOXMLDOC01-appb-M000008
 実際には、x´ij、Rj 0、tj 0 (j=1～M)は誤差を含む。このため、必ずしも行列Aiに零空間は存在しない。この場合、最も小さい特異値に対する右特異ベクトルを零空間と考え、第４成分が1となるよう正規化することで（Xi Yi Zi）Tを得ることができる。
　ここで、図３の説明に戻る。初期値の算出の後のステップＳ１１０において、バンドル調整部１０５は、全画像について、バンドル調整の初期値としての基準座標系での画像姿勢を算出したか否かを判定する。ステップＳ１１０において、全画像についての初期値を算出していないと判定されたときに、処理はステップＳ１０２に戻る。この場合、別の画像ペアに対して前述した、相対姿勢の算出、一意性判断、初期値の算出が行われる。そして、ステップＳ１１０において、全画像についての初期値を算出したと判定されたときに、処理はステップＳ１１１に移行する。
　ステップＳ１１１において、バンドル調整部１０５は、バンドル調整を行って全特徴点の３次元座標と画像を撮像した時点の撮像部５の姿勢とを算出する。ここでの撮像部５の姿勢は、光学中心の３次元座標と光軸の方向及び光軸方向を軸とした撮像素子の回転角度を含む。以下、バンドル調整について説明する。
　バンドル調整は、観測された特徴点の画像上の座標とよく一致する、特徴点の３次元座標と撮像部の姿勢とを算出する処理である。バンドル調整では、例えば、以下の（式８）で表されるコストCを最小化する特徴点の３次元座標と撮像部の姿勢とが算出される。
Figure JPOXMLDOC01-appb-M000009
 したがって、コストCを最小化する特徴点の３次元座標をPI mとし、撮像部の姿勢をRJ m、tJ mとすると、PI m、RJ m、tJ mは以下の（式９）から算出される。
Figure JPOXMLDOC01-appb-M000010
 ここで、（式９）のPiは、i番目の特徴点（特徴点i）の基準座標系における３次元座標を表す３次元ベクトルである。tjは、j番目の画像（画像j）が撮像されたときの撮像部５の座標系と基準座標系との間の座標変換を表す並進ベクトルである。Rjは、j番目の画像（画像j）が撮像されたときの撮像部５の座標系と基準座標系との間の座標変換を表す回転行列であって、行列式が1である３×３直交行列である。このとき、画像jを撮像した時点での撮像部５の座標系での特徴点iの３次元座標P´ijは、P´ij＝Rj T(Pi－tj)と表される。qijは特徴点iの画像jでの座標を正規化した2次元ベクトルである。qijは、以下の（式１０）で表すことができる特徴点iの画像jでの正規化座標である。
Figure JPOXMLDOC01-appb-M000011
 ここで、（式１０）におけるQij,x、Qij,yは、特徴点iの画像jでの座標Qij=(Qij,x Qij,y)Tのｘ座標及びｙ座標である。また、pij(Pi,Rj,tj)は、３次元座標Piにある特徴点iがtj及びRjで示される座標系で撮像が行われた画像jにおいて写っているはずの座標である。pij(Pi,Rj,tj)は、以下の（式１１）で表すことができる。
Figure JPOXMLDOC01-appb-M000012
 なお、qij－pij(Pi,Rj,tj)は、特徴点iが画像j上で、画像から検出された正規化座標qijと画像jの姿勢がtj及びRjであった場合の特徴点iの画像上の座標pij(Pi,Rj,tj)との差である。このqij－pij(Pi,Rj,tj)は、再投影誤差とよばれている。また、（式８）の関数ρ（δ）は誤差δを引数とする関数である。関数ρ（δ）としては、例えば、以下の（式１２）で示される２乗誤差の関数や以下の（式１３）で示されるHuberの損失（ロス）関数等を用いることができる。コストCの算出のための関数ρ（δ）は、姿勢Rj,tjが算出されたすべての画像j中で、画像上の座標が算出されたすべての特徴点iについて積算する。
Figure JPOXMLDOC01-appb-M000013
　通常、コストCは、Pi,Rj,tjについての非線形関数であり、コストCを最小化するPi,Rj,tjを解析的に求めることが困難である。したがって、本実施形態では、コストCを最小化するPi,Rj,tjを、前述した初期値Pi 0、Rj 0、tj 0からの勾配法によって求める。ここで、勾配法は、nステップ目でのPi,Rj,tjであるPi n、Rj n、tj nからPi,Rj,tjを微小に変化させたときのコストCの変化から求められる更新量ΔPi n、ΔRj n、Δtn jに基づいて(n+1)ステップ目でのPi,Rj,tjであるPi n+1、Rj n+1、tj n+1をそれぞれPi n+1＝Pi n＋ΔPi n、Rj n+1＝Rj n+1＋ΔRj n、tj n+1＝tj n+1＋Δtj nより算出していき、Pi,Rj,tjが収束した値をPI m、RJ m、tJ mとする手法である。勾配法には、最急降下法等の様々な方法が存在する。一般に大規模な非線形関数を最小化する際の勾配法では、より収束が速い方法としてLevenberg-Marquardt法が使用される。
　なお、コストCはPi,Rj,tjについて複雑な形状をしているので、十分に真の解に近い初期値から最適化を開始しない場合には、収束までの処理時間が長くなったり、誤差が大きい誤った値に収束したりする。したがって、高精度に３次元復元を行うためには、よい初期値から処理を開始する必要がある。本実施形態では、画像ペアとして選択された２枚の画像の間の相対姿勢を一意に算出できないと判断されているときには、その画像ペアから算出された相対姿勢を用いてのバンドル調整の初期値の算出は行われない。このような曖昧な相対姿勢を用いての初期値の算出を行わないようにすることでバンドル調整の精度を向上させることができる。
　バンドル調整後、３次元復元画像処理装置は、バンドル調整結果を出力する。この後、撮像装置１は、バンドル調整結果をＲＡＭ３に記録したり、表示部６に出力したりする。
　以上のようなバンドル調整の後、図３の処理は終了する。バンドル調整の後で、３次元座標算出部１０６による画像内の全点の３次元座標（被写体の３次元形状）の算出が行ってもよい。この画像内の全点の３次元座標の算出にバンドル調整の結果として得られたPI m、RJ m、tJ mが用いられる。そして、画像内の全点の３次元座標の算出の後、必要に応じて表示部６に復元された被写体の３次元画像が表示される。
　以上説明したように本実施形態によれば、画像ペアとして選択された２枚の画像の間の相対姿勢を一意に算出できないと判断されているときには、ステップＳ１０９で説明する初期値の算出が行われない。これにより、被写体が平面であるときでも、短時間で、かつ、精度よくバンドル調整を行うことができる。
　［変形例１］
　以下、第１の実施形態の変形例を説明する。図７は、変形例１の３次元復元画像処理装置７の動作を示すフローチャートである。図３の処理では、画像ペアが選択される毎に初期値が算出される。そして、相対姿勢を一意に算出できる全画像についての初期値が算出された後でバンドル調整が行われる。これに対し、図７に示す変形例ではステップＳ１０９においてある画像ペアについての初期値が算出される毎に処理がステップＳ１１１に移行し、それまでに初期値が算出済みの全画像の初期値を用いてバンドル調整が行われる。その後、処理がステップＳ１１２に移行し、相対姿勢を一意に算出できると判断された全画像を用いてのバンドル調整が行われたときに処理が終了する。このような変形例１のように逐次的にバンドル調整が行われてもよい。
　［変形例２］
　また、全画像の間で自分以外の全画像との間で画像ペアを選択し、N枚の画像に対して作成したN(N-1)/2組の画像ペアのうち、相対姿勢を一意に算出できると判断された画像ペアの相対姿勢から、初期値としての基準座標系に対する各画像の姿勢Rj 0、tj 0及び特徴点の３次元座標Pi 0が算出されてもよい。
　［変形例３］
　前述した実施形態では、一意性判断において、２枚の画像の間の相対姿勢を一意に算出できないと判断されたときには必ずその相対姿勢は初期値の算出には用いられない。しかしながら、必ずしも全ての相対姿勢が用いられないようにしなくてもよい。
　［変形例４］
　前述したように、２枚の画像の間の相対姿勢には４通りの解が存在する。姿勢算出部１０２は、この４通りの解のうち、被写体（特徴点）が撮像部５の前方にあると考えられる解（解１）を最終的な相対姿勢とする。しかしながら、姿勢算出部１０２の構成によっては、姿勢算出部１０２は、２通りの相対姿勢を算出することもある。被写体が平面であり、相対姿勢を一意に判断することができない場合、すなわち、解１と解２の両方について全特徴点が撮像部５の前方にある場合、この２通りの解を用いて行った一意性判断を行うと、いずれの相対姿勢を用いた場合でも、一意に相対姿勢を算出することができないと判断される。そのため、姿勢算出部１０２は、解１と解２の２通りの解による相対姿勢を算出する場合でも、一意性判断部１０４は、姿勢算出部１０２によって算出された何れかの相対姿勢だけを用いて一意性判断を行ってもよい。
　［第２の実施形態］
　次に、本発明の第２の実施形態を説明する。図８は、第２の実施形態における３次元復元画像処理装置７の構成を示すブロック図である。ここで、図８において図２と同一の構成については図２と同一の参照符号を付している。この同一の参照符号を付した構成については説明を省略する。
　第２の実施形態の３次元復元画像処理装置７は、さらに、平面判定処理部１０３を有している。平面判定処理部１０３は、姿勢算出部１０２と一意性判断部１０４との間に設けられ、被写体が平面であるか否かを判定する。第１の実施形態では、被写体の形状によらずに、２枚の画像の間の対応特徴点の分布の仕方のみから一意性判断が行われる。しかしながら、被写体が平面でない場合には、２枚の画像の間の対応特徴点の分布の仕方によらずに相対姿勢は一意に算出できる。そこで、第２の実施形態は、一意性判断に先立って被写体が平面であるか否かを判定し、被写体が平面であるときだけ一意性判断を行うものである。
　図９は、第２の実施形態における３次元復元画像処理装置７の動作を示すフローチャートである。ここで、図９において図３と同一の処理については図３と同一の符号を付している。この同一の符号を付した処理については説明を省略する。すなわち、図９では、ステップＳ１０３において２枚の画像の間の相対姿勢が算出された後で処理がステップＳ１０５に移行する。ステップＳ１０５において、平面判定処理部１０３は、被写体が平面であるか否かを判定する。例えば、平面判定処理部１０３は、２枚の画像の間の特徴点の座標の対応を射影変換の関係で表すことができるかどうかから被写体が平面であるか否かを判定する。被写体が平面である場合、１枚目の画像と２枚目の画像での対応特徴点の座標をそれぞれ２次元ベクトルx´1＝（x´1 y´1）T、x´2＝（x´2 y´2）Tとするとき以下の（式１４）で示す条件を満たす３×３行列Hが存在する。
Figure JPOXMLDOC01-appb-M000014
 そこで、１枚目の画像と２枚目の画像での対応特徴点の座標から画像間の射影変換行列Hを求め、１枚目の画像の対応特徴点の座標を射影変換行列Hにて変換した座標と、２枚目の対応特徴点の座標との差を射影変換の誤差とし、この誤差をしきい値と比較することで被写体が平面であるかを判定することができる。射影変換の誤差から被写体が平面であるか判定するしきい値としては、GRIC(Geometric Robust Information Criterion)値を用いることができる。このような処理の後、平面判定処理部１０３は、被写体が平面であるか否かを示す情報を一意性判断部１０４に出力する。その後、処理はステップＳ１０６に移行する。
　ステップＳ１０６において、一意性判断部１０４は、平面判定処理部１０３からの情報に基づいて被写体が平面であるか否かを判定する。ステップＳ１０６において被写体が平面であると判定されたときには、処理はステップＳ１０７に移行する。これ以後、第１の実施形態と同様の処理が行われる。一方、ステップＳ１０６において被写体が平面でないと判定されたときには、処理はステップＳ１０９に移行する。すなわち、一意性判断が省略される。
　以上説明したように本実施形態によれば、一意性判断に先立って被写体が平面であるか否かを判定することで、不要な一意性判断を省略することができ、処理の効率化につながる。
　なお、第２の実施形態において、第１の実施形態で説明した各種の変形例を適用してもよい。
　［変形例１］
　以下、第２の実施形態の変形例を説明する。図１０は、変形例１の３次元復元画像処理装置７の動作を示すフローチャートである。変形例１は、取得部１０１によって取得された複数枚の画像が、共通の視野を持つように同一の被写体が撮像されて得られた画像等のときに適用できる処理である。この場合、１つの画像ペアについての平面判定の結果を他の画像ペアについての平面判定にも利用できる。すなわち、平面判定は１回だけ行われればよい。
　したがって、変形例１では、ステップＳ１０３において２枚の画像の間の相対姿勢が算出された後、処理がステップＳ１０４に移行する。ステップＳ１０４において、平面判定処理部１０３は、被写体の平面判定を実行済みであるか否かを判定する。そして、ステップＳ１０４において、被写体の平面判定を実行済みでないと判定されたときだけ処理がステップＳ１０５に移行する。ステップＳ１０４において、被写体の平面判定を実行済みであると判定されたときにはステップＳ１０５の処理がスキップされる。
　このように変形例１では、１つの画像ペアについての平面判定の結果を他の画像ペアについての平面判定にも利用できるときには、平面判定が１回だけ行われる。これにより、処理の効率化につながる。
　［変形例２］
　前述した平面判定では、２枚の画像の間の特徴点の座標の対応を射影変換の関係で表すことができるかどうかから被写体が平面であるか否かが判定される。しかしながら、平面判定の手法は、これに限るものではない。例えば、姿勢算出部１０２で算出された２枚の画像の間の相対姿勢から算出される被写体の３次元座標から被写体が平面であるか否かが判定されてもよい。
　２枚の画像の間の相対姿勢が既知であるとき、例えばDLT法にて被写体の３次元座標を算出することができる。つまり、２枚の画像に共通で座標を算出した全ての対応特徴点に対する３次元座標を算出し、３次元座標をPCA(主成分分解)した第１主成分ベクトルと第２主成分ベクトルが張る面を得ることで、被写体が張る主要な平面を求めることができる。全ての対応特徴点のうちで主要平面の３次元距離がしきい値以下の特徴点数の割合が十分に大きいときには、全ての対応特徴点が平面上にあるとする。このとき、被写体は平面であると判定される。ここで、しきい値は、画像を撮像した光学系の光学中心と被写体までの距離とを基準として決められてもよいし、PCAで算出された第１主成分に対する固有値や第２主成分に対する固有値を基準として決められてもよい。

Claims (9)

  1. 　異なる視点で撮像された複数枚の画像を取得する取得部と、
    　取得された前記複数枚の画像から２枚ずつの画像を画像ペアとして選択し、選択した前記画像ペアそれぞれを構成する２枚の画像の間の相対姿勢を算出する姿勢算出部と、
    　選択された前記画像ペアを構成する２枚の画像の間の前記相対姿勢が一意に算出できるか否かを判断する一意性判断部と、
    　前記一意性判断部の判断結果に基づいてバンドル調整の初期値の算出に用いる前記相対姿勢を決定し、決定された前記相対姿勢に基づいて算出された前記バンドル調整の初期値よりバンドル調整を行うバンドル調整部と、
    　を具備する３次元復元画像処理装置。
  2. 　前記バンドル調整部は、２枚の画像の間の前記相対姿勢が一意に算出できると判断された前記相対姿勢のみを前記バンドル調整の初期値の算出に用いる前記相対姿勢として決定する請求項１に記載の３次元復元画像処理装置。
  3. 　前記一意性判断部は、前記姿勢算出部によって算出された１つ又は２つの前記相対姿勢のうちの何れかを用いて前記相対姿勢が一意に算出できるか否かを判断する請求項１又は２に記載の３次元復元画像処理装置。
  4. 　前記画像ペアを構成する２枚の画像から、被写体が平面であるか否かを判定する平面判定処理部をさらに具備し、
    　前記一意性判断部は、前記被写体が平面であると判定されたときに前記相対姿勢が一意に算出できる否かを判断し、前記被写体が平面でないと判定されたときに前記相対姿勢が一意に算出できると判断する請求項１乃至３の何れか１項に記載の３次元復元画像処理装置。
  5. 　前記画像ペアを構成する２枚の画像から、被写体が平面であるか否かを判定する平面判定処理部をさらに具備し、
    　前記一意性判断部は、１つの前記画像ペアについて前記平面判定処理部よって判定された前記被写体が平面であるか否かの判定結果を他の全ての画像ペアについての前記被写体が平面であるか否かの判定結果としても利用する請求項１乃至３の何れか１項に記載の３次元復元画像処理装置。
  6. 　前記平面判定処理部は、前記画像ペアを構成する２枚の画像のそれぞれの特徴点の座標を算出し、前記画像ペアを構成する２枚の画像における前記特徴点の座標の関係が射影変換の関係で表すことできるときに前記被写体が平面であると判定する請求項４又は５に記載の３次元復元画像処理装置。
  7. 　前記平面判定処理部は、前記画像ペアについて算出された前記相対姿勢から前記画像ペアを構成する２枚の画像における特徴点の３次元座標を算出し、算出した前記３次元座標を有する前記特徴点が同一平面上に存在するときに前記被写体が平面であると判定する請求項４又は５に記載の３次元復元画像処理装置。
  8. 　異なる視点で撮像された複数枚の画像を取得することと、
    　取得した前記複数枚の画像から２枚ずつの画像を画像ペアとして選択し、選択した前記画像ペアそれぞれを構成する２枚の画像の間の相対姿勢を算出することと、
    　選択された前記画像ペアを構成する２枚の画像の間の前記相対姿勢が一意に算出できるか否かを判断することと、
    　前記判断の結果に基づいてバンドル調整の初期値の算出に用いる前記相対姿勢を決定し、決定された前記相対姿勢に基づいて算出された前記バンドル調整の初期値よりバンドル調整を行うことと、
    　を具備する３次元復元画像処理方法。
  9. 　異なる視点で撮像された複数枚の画像を取得することと、
    　取得した前記複数枚の画像から２枚ずつの画像を画像ペアとして選択し、選択した前記画像ペアそれぞれを構成する２枚の画像の間の相対姿勢を算出することと、
    　選択された前記画像ペアを構成する２枚の画像の間の前記相対姿勢が一意に算出できるか否かを判断することと、
    　前記判断の結果に基づいてバンドル調整の初期値の算出に用いる前記相対姿勢を決定し、決定された前記相対姿勢に基づいて算出された前記バンドル調整の初期値よりバンドル調整を行うことと、
    　をコンピュータに実行させるための３次元復元画像処理プログラムを記憶したコンピュータ読み取り可能な記憶媒体。
PCT/JP2017/034141 2017-09-21 2017-09-21 ３次元復元画像処理装置、３次元復元画像処理方法及び３次元復元画像処理プログラムを記憶したコンピュータ読み取り可能な記憶媒体 WO2019058487A1 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
PCT/JP2017/034141 WO2019058487A1 (ja) 2017-09-21 2017-09-21 ３次元復元画像処理装置、３次元復元画像処理方法及び３次元復元画像処理プログラムを記憶したコンピュータ読み取り可能な記憶媒体
US16/746,006 US11176702B2 (en) 2017-09-21 2020-01-17 3D image reconstruction processing apparatus, 3D image reconstruction processing method and computer-readable storage medium storing 3D image reconstruction processing program

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2017/034141 WO2019058487A1 (ja) 2017-09-21 2017-09-21 ３次元復元画像処理装置、３次元復元画像処理方法及び３次元復元画像処理プログラムを記憶したコンピュータ読み取り可能な記憶媒体

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US16/746,006 Continuation US11176702B2 (en) 2017-09-21 2020-01-17 3D image reconstruction processing apparatus, 3D image reconstruction processing method and computer-readable storage medium storing 3D image reconstruction processing program

Publications (1)

Publication Number Publication Date
WO2019058487A1 true WO2019058487A1 (ja) 2019-03-28

Family

ID=65810325

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2017/034141 WO2019058487A1 (ja) 2017-09-21 2017-09-21 ３次元復元画像処理装置、３次元復元画像処理方法及び３次元復元画像処理プログラムを記憶したコンピュータ読み取り可能な記憶媒体

Country Status (2)

Country Link
US (1) US11176702B2 (ja)
WO (1) WO2019058487A1 (ja)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111489438A (zh) * 2020-04-09 2020-08-04 深圳趣途科技有限责任公司 重建三维模型方法、重建三维模型***及计算机存储介质
JP2021192226A (ja) * 2020-06-30 2021-12-16 ベイジン バイドゥ ネットコム サイエンス テクノロジー カンパニー リミテッド 画像を処理するための方法と装置、電子機器、コンピュータ可読記憶媒体及びコンピュータープログラム

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004334819A (ja) * 2003-03-13 2004-11-25 Toshiba Corp ステレオキャリブレーション装置とそれを用いたステレオ画像監視装置
JP2009014628A (ja) * 2007-07-06 2009-01-22 Topcon Corp 位置計測装置及び位置計測方法

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2007137388A1 (en) * 2006-05-26 2007-12-06 Corporation Spg Data3D Photogrammetric system and techniques for 3d acquisition
CN104660900B (zh) 2013-10-30 2018-03-02 株式会社摩如富 图像处理装置及图像处理方法
US9460513B1 (en) * 2015-06-17 2016-10-04 Mitsubishi Electric Research Laboratories, Inc. Method for reconstructing a 3D scene as a 3D model using images acquired by 3D sensors and omnidirectional cameras
CN108537876B (zh) * 2018-03-05 2020-10-16 清华-伯克利深圳学院筹备办公室 三维重建方法、装置、设备及存储介质

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004334819A (ja) * 2003-03-13 2004-11-25 Toshiba Corp ステレオキャリブレーション装置とそれを用いたステレオ画像監視装置
JP2009014628A (ja) * 2007-07-06 2009-01-22 Topcon Corp 位置計測装置及び位置計測方法

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111489438A (zh) * 2020-04-09 2020-08-04 深圳趣途科技有限责任公司 重建三维模型方法、重建三维模型***及计算机存储介质
JP2021192226A (ja) * 2020-06-30 2021-12-16 ベイジン バイドゥ ネットコム サイエンス テクノロジー カンパニー リミテッド 画像を処理するための方法と装置、電子機器、コンピュータ可読記憶媒体及びコンピュータープログラム
JP7246427B2 (ja) 2020-06-30 2023-03-27 ベイジン バイドゥ ネットコム サイエンス テクノロジー カンパニー リミテッド 画像を処理するための方法と装置、電子機器、コンピュータ可読記憶媒体及びコンピュータープログラム
US11704896B2 (en) 2020-06-30 2023-07-18 Beijing Baidu Netcom Science And Technology Co., Ltd. Method, apparatus, device and storage medium for image processing

Also Published As

Publication number Publication date
US11176702B2 (en) 2021-11-16
US20200202563A1 (en) 2020-06-25

Similar Documents

Publication Publication Date Title
JP6330987B2 (ja) 画像処理装置、画像処理方法、及び記憶媒体
Moulon et al. Adaptive structure from motion with a contrario model estimation
CN101563709B (zh) 校准照相机***
KR101791590B1 (ko) 물체 자세 인식장치 및 이를 이용한 물체 자세 인식방법
JP6324025B2 (ja) 情報処理装置、情報処理方法
US10853960B2 (en) Stereo matching method and apparatus
US10460471B2 (en) Camera pose estimating method and system
JP2008506953A5 (ja)
JP6192507B2 (ja) 画像処理装置、その制御方法、および制御プログラム、並びに撮像装置
JP6936974B2 (ja) 位置姿勢推定装置、位置姿勢推定方法及びプログラム
US10346949B1 (en) Image registration
US11176702B2 (en) 3D image reconstruction processing apparatus, 3D image reconstruction processing method and computer-readable storage medium storing 3D image reconstruction processing program
JP6359985B2 (ja) デプス推定モデル生成装置及びデプス推定装置
JP5937977B2 (ja) 変換行列推定装置、変換行列推定方法、及びプログラム
CN109859313B (zh) 3d点云数据获取方法、装置、3d数据生成方法及***
CN109902695B (zh) 一种面向像对直线特征匹配的线特征矫正与提纯方法
CN111089579B (zh) 异构双目slam方法、装置及电子设备
JP6080424B2 (ja) 対応点探索装置、そのプログラムおよびカメラパラメータ推定装置
JP2011022084A (ja) 三次元姿勢測定装置および三次元姿勢測定方法
JP6843552B2 (ja) 画像処理装置、画像処理方法およびプログラム。
Georgiev et al. A fast and accurate re-calibration technique for misaligned stereo cameras
JP6584139B2 (ja) 情報処理装置、情報処理方法及びプログラム
KR20150119770A (ko) 카메라를 사용한 3차원 좌표 측정 장치 및 방법
JP2007034964A (ja) カメラ視点運動並びに３次元情報の復元及びレンズ歪パラメータの推定方法、装置、カメラ視点運動並びに３次元情報の復元及びレンズ歪パラメータの推定プログラム
JP2000353244A (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: 17926119

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 17926119

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: JP