Image synthesis method, image synthesizer, and recording medium on which image synthesis program is recorded

[0002] 1. Field of the Invention

[0003] The present invention relates to an image synthesis method, an image synthesizer, and a recording medium on which an image synthesis program has been recorded.

[0004] 2. Description of the Prior Art

[0005] A technique for calculating an optical flow from two images, and registering the two images on the basis of the obtained optical flow has been known. Description is made of a conventional method of calculating an optical flow.

[0006] (1) Lucas-Kanade Method

[0007] A large number of methods of calculating the apparent optical flow of a moving object in a moving image have been conventionally proposed. The Lucas-Kanade method which is a local gradient method out of the methods is one of the best methods. The reason for this is that the speed of processing is high, implementing is easy, and the result has confidence.

[0008] As for the details of the Lucas-Kanade method, see an article: B. Lucas and T. Kanade, “An Iterative Image Registration Technique with an Application to Stereo Vision,” In Seventh International Joint Conference on Artificial Intelligence (IJCAI-81), pp.674-979, 1981.

[0009] The outline of the Lucas-Kanade method will be described.

[0010] When a gray scale pattern I=(x, y, t) of image coordinates p=(x, y) at time t is moved to coordinates (x+δx, y+δy) with its gradation distribution kept constant after a very short time period (δt), the following optical flow constraint equation (1) holds: $\begin{array}{cc}\frac{\partial I}{\partial x}\frac{\delta x}{\delta t}+\frac{\partial I}{\partial y}\frac{\mathrm{\delta y}}{\delta t}+\frac{\partial I}{\partial t}=0& \left(1\right)\end{array}$

[0011] In order to calculate an optical flow {v=(δx/δt, δy/δt)=(u, v)} in a two-dimensional image, another constraint equation is required because the number of unknown parameters is two. Lucas and Kanade have assumed that the optical flow is constant in a local region of an object.

[0012] Suppose the optical flow is constant within a local region ω on an image, for example. In this case, a squared error E of a gray scale pattern which is to be minimized can be defined by the following equation (2) when the following substitutions are made: $\begin{array}{cc}{I}_0\left(P\right)=I\left(x,y,t\right),I_1\left(p+v\right)=I\left(x+u,y+v,t+\delta t\right)E=\sum _{\omega }{\left[I_1\left(p+v\right)-I_0\left(p\right)\right]}^2& \left(2\right)\end{array}$

[0013] When v is very small, the terms of second and higher degrees in Taylor's expansion can be ignored, so that a relationship expressed by the following equation (3) holds:

[0014] I₁(p+v)=I₁(p)+g(p)v  (3)

[0015] where g(p) is a linear differential of I₁(p).

[0016] The error E is minimized when the derivative of E with respect to v is zero, so that a relationship expressed by the following equation (4) holds: $\begin{array}{cc}\begin{array}{c}0=\frac{\partial }{\partial v}E\\ \approx \frac{\partial }{\partial v}\sum _{\omega }{\left[I_1\left(p\right)+g\left(p\right)v-I_0\left(p\right)\right]}^2\\ =\sum _{\omega }2g\left(p\right)\left[I_1\left(p\right)+g\left(p\right)v-I_0\left(p\right)\right]\end{array}& \left(4\right)\end{array}$

[0017] Therefore, the optical flow v is found by the following equation (5): $\begin{array}{cc}v\approx \frac{\sum _{\omega }g\left(p\right)\left[I_0\left(p\right)-I_1\left(p\right)\right]}{\sum _{\omega }{g\left(p\right)}^2}& \left(5\right)\end{array}$

[0018] Furthermore, the optical flow can be found with high precision by Newton-Raphson iteration, as expressed by the following equation (6): $\begin{array}{cc}{v}_{k+1}=v_k+\frac{\sum g_k\left[I_0-I_1^k\right]}{\sum {\left(g_k\right)}^2}{I}_1^k=I_1\left(p+v_k\right),g_k=g\left(p+v_k\right),I_0=I_0\left(p\right)& \left(6\right)\end{array}$

[0019] (2) Hierarchical Estimation Method

[0020] The largest problem of the gradient methods, including the Lucas-Kanade method, is that they cannot be applied to a large motion because a good initial value is required. Therefore, a method of producing images respectively having resolutions which differ at several levels like a pyramid hierarchical structure to solve the problem has been conventionally proposed.

[0021] Images having resolutions which differ at several levels are first previously produced from each of two consecutive images. An approximate optical flow is then calculated between the images having the lowest resolution. A more precise optical flow is calculated between the images having resolutions which are higher by one level. The processing is successively repeated until the optical flow is calculated between the images having the highest resolution.

[0022] FIG. 4, FIG. 3, FIG. 2, and FIG. 1 respectively illustrate an original image, an image having a lower resolution than that of the original image shown in FIG. 4, an image having a lower resolution than that of the image having a low resolution shown in FIG. 3, and an image having a lower resolution than that of the image having a low resolution shown in FIG. 2. In FIGS. 1 to 4, S indicates one patch.

[0023] An optical flow is gradually found from the image shown in FIG. 1 (an image in a hierarchy 1) the image shown in FIG. 2 (an a hierarchy 2 the image shown in FIG. 3 (an image in a hierarchy 3), and the image shown in FIG. 4 (an image in a hierarchy 4) in this order. In FIGS. 1 to 4, an arrow indicates an optical flow vector found for each patch.

[0024] However, the problem is that in a real image, there are few regions containing sufficient texture, so that a reliable optical flow is not obtained.

[0025] A technique for affixing a plurality of images to one another and synthesizing the images to obtain a seamless image having a wide field of view and having a high resolution (image mosaicking) has been conventionally actively studied. Classical applications include synthesis of aerial photographs or satellite photographs. In recent years, a method of synthesizing a plurality of digital images to obtain a seamless panoramic image, and constructing a virtual reality environment has been paid attention to.

[0026] The following two methods have been known as a technique for obtaining a panoramic image by synthesis.

[0027] The first method is a method of first translating a camera to previously pick up a plurality of images. The plurality of images obtained are simultaneously displayed on a monitor by a personal computer. A user designates corresponding points between the two images, to synthesize the two images.

[0028] In the first method, the motion of the camera is restricted to translation. In the first method, the user must designate the corresponding points.

[0029] The second method is a method of fixing a camera to a tripod and restricting the motion of the camera to only rotation on a horizontal surface, to pick up a plurality of images. The plurality of images obtained are projected on a cylindrical surface, to synthesize the images (see U.S. Pat. No. 5,396,583).

[0030] In the second method, the, motion of the camera must be restricted to only rotation on a horizontal surface. Further, the focal length or the angle of view of the camera must be measured.

[0031] An object of the present invention is to provide an image synthesis method in which highly precise positioning is performed even when a scene having a large depth is obtained by synthesis, an image synthesizer,and a recording medium on which an image synthesis program has been recorded.

[0032] Another object of the present invention is to provide an image synthesis method in which a seamless panoramic image can be obtained from a plurality of images, a camera for picking up the plurality of images is allowed to freely move, and the focal length need not be measured, an image synthesizer, and a recording medium on which an image synthesis program has been recorded.

[0033] A first image synthesis method according to the present invention is an image synthesis method for synthesizing a first image and a second image without a seam on the basis of coordinate values, on the first image, respectively corresponding to a plurality of feature points which are set in a portion, which is overlapped with the second image, of the first image and coordinate values, on the second image, at positions, on the second image, respectively corresponding to the feature points, characterized by comprising the first step of dividing both the images for each triangular patch having the three adjacent feature points as vertexes, and finding a transform equation for finding a position, on the first image, corresponding to a point on the second image for the triangular patch; the second step of finding the coordinate value, on the first image, corresponding to each of pixels in each of the triangular patches, on the second image using the transform equation calculated for the triangular patch; the third step of synthesizing the first image and the second image using t he coordinate value, on the first image, corresponding to each of the pixels on the second image, which is found at the second step.

[0034] An example of the transform equation found for each triangular patch is a planar projective transform matrix. The planar projective transform matrix is found on the basis of corresponding points at the three vertexes of the triangular patch and a fundamental matrix representing epipolar constraint conditions between the first image and the second image, for example. The planar projective transform matrix may be found on the basis of a correspondence between corresponding points at the three vertexes of the triangular patch and epipoles representing epipolar constraint conditions between the first image and the second image.

[0035] An example of the transform equation found for each triangular patch is an affine transform matrix. The affine transform matrix is found on the basis of the corresponding points at the three vertexes of the triangular patch.

[0036] An image synthesizer according to the present invention is an image synthesizer for synthesizing a first image and a second image without a seam on the basis of coordinate values, on the first image, respectively corresponding to a plurality of feature points which are set in a portion, which is overlapped with the second image, of the first image and coordinate values, on the second image, at positions, on the second image, respectively corresponding to the feature points, characterized by comprising first means for dividing both the images for each triangular patch having three adjacent feature points as vertexes, and finding a transform equation for finding a position, on the first image, corresponding to a point on the second image for the triangular patch; second means for finding the coordinate value, on the first image, corresponding to each of pixels in each of the triangular patches, on the second image using the transform equation calculated for the triangular patch; third means for synthesizing the first image and the second image using the coordinate value, on the first image, corresponding to each of the pixels on the second image which is found by the second means.

[0037] A computer readable recording medium on which an image synthesis program has been recorded according to the present invention is a computer readable recording medium on which an image synthesis program for synthesizing a first image and a second image without a seam on the basis of coordinate values, on the first image, respectively corresponding to a plurality of feature points which are set in a portion, which is overlapped with the second image, of the first image and coordinate values, on the second image, at positions, on the second image, respectively corresponding to the feature points has been recorded, characterized by comprising the first step of dividing both the images for each triangular patch having three adjacent feature points as vertexes, and finding a transform equation for finding a position, on the first image, corresponding to a point on the second image for the triangular patch; the second step of finding the coordinate value, on the first image, corresponding to each of pixels in each of the triangular patches, on the second image using the transform equation calculated for the triangular patch; and the third step of synthesizing the first image and the second image using the coordinate value, on the first image, corresponding to each of the pixels on the second image found at the second step.

[0038] A second image synthesis method according to the present invention is an image synthesis method for synthesizing a first image and a second image, characterized by comprising the first step of extracting from the first image a plurality of partial images which are effective for tracking by an optical flow between the first image and the second image as feature points; the second step of tracking a point, on the second image, corresponding to each of the feature points on the first image on the basis of the optical flow between both the images; the third step of finding a coordinate value, on the first image, corresponding to each of pixels on the second image on the basis of each of the feature points on the first image and the corresponding points on the second image which is found at the second step; and the fourth step of synthesizing the first image and the second image on the basis of the coordinate value, on the first image, corresponding to each of the pixels on the second image which is obtained at the third step, the third step comprising the step of dividing both the images for each triangular patch having three adjacent feature points as vertexes, and finding a transform equation for finding a position, on the first image, corresponding to the point on the second image for the triangular patch, and the step of finding the coordinate value, on the first image, corresponding to each of the pixels in the triangular patch on the second image using the transform equation calculated for the triangular patch.

[0039] An example of the first step is one comprising the steps of extracting an overlapped portion between the first image and the second image, and extracting from the portion, which is overlapped with the second image, of the first image a plurality of partial images which are effective for tracking by the optical flow between both the images as feature points.

[0040] At the step of extracting the overlapped portion between the first image and the second image, the overlapped portion between the first image and the second image is extracted on the basis of either one of the SSD method and the normalized cross-correlation method.

[0041] An example of the second step is one comprising the step of finding the optical flow from the first image and the second image on the basis of an optical flow estimation method using the hierarchically structured Lucas-Kanade method for interpolating the optical flow having low confidence out of the optical flows obtained in the respective stages of the optical flow estimation method using the optical flows in its surrounding regions, and the step of tracking the point, on the second image, corresponding to each of the feature points on the first image on the basis of the obtained optical flow.

[0042] The foregoing and other objects, features, aspects and advantages of the present invention will become more apparent from the following detailed description of the present invention when taken in conjunction with the accompanying drawings.

[0043] FIG. 1 is a schematic view for explaining a hierarchical estimation method, which illustrates an image in a hierarchy 1;

[0044] FIG. 2 is a schematic view for explaining a hierarchical estimation method, which illustrates an image in a hierarchy 2;

[0045] FIG. 3 is a schematic view for explaining a hierarchical estimation method, which illustrates an image in a hierarchy 3; FIG. 4 is a schematic view for explaining a hierarchical estimation method, which illustrates an image in a hierarchy 4;

[0046] FIG. 5 is a schematic view for explaining dilation processing performed in an optical flow estimation method employed in the present embodiment;

[0047] FIG. 6ais a schematic view showing an example of an optical flow which has been subjected to threshold processing for an image in a hierarchy;

[0048] FIG. 6bis a schematic view showing optical flow after filling;

[0049] FIG. 7 is a block diagram showing the configuration of a panoramic image synthesizer;

[0050] FIG. 8 is a flow chart showing the procedure for panoramic image synthesis processing;

[0051] FIG. 9 is an explanatory view for explaining a planar projective transform matrix;

[0052] FIG. 10 is a schematic view showing a triangular patch;

[0053] FIG. 11 is a schematic view showing epipolar constraint conditions;

[0054] FIG. 12 is a schematic view showing a triangular patch having N, M, and P as vertexes and epipoles e₁and e₂; and

[0055] FIG. 13 is an explanatory view for explaining an example of a pixel value blending method.

[0056] Description is now made of an embodiment in a case where the present invention is applied to a panoramic image synthesizer.

[0057] An optical flow estimation method employed in the panoramic image synthesizer will be described before describing the panoramic image synthesizer.

[0058] [1] Description of Optical Flow Estimation Method

[0059] An optical flow estimation method employed in the present invention presupposes hierarchical prediction for producing images having resolutions which differ at several levels like a pyramid hierarchical structure to gradually calculate an optical flow. A method of calculating an optical flow conforms to a gradient method such as the Lucas-Kanade method. That is, it presupposes an optical flow estimation method using a hierarchically structured gradient method. The Lucas-Kanade method is used as the gradient method.

[0060] The optical flow estimation method employed in the present embodiment is characterized in that an optical flow obtained in each of stages of the optical flow estimation method using the hierarchically structured Lucas-Kanade method is filled by dilation processing. This will be described in detail below.

[0061] One of advantages of the Lucas-Kanade method is that the result of trading has confidence. Tomasi and Kanade have showed that the trackability of a region can be calculated from image derivatives as follows (C. Tomasi and Kanade, “Shape and Motion from Image Streams: a Factorization Method-Part 3 Detection and Tracking of Point Features,” CMU-CS-91-132, Carnegie Mellon University, 1991).

[0062] From a 2×2 matrix of coefficients G in the following equation (7) having as elements the squares of differentials in the vertical and horizontal directions of an image ω in a region, its eigenvalues are calculated, thereby making it possible to determine the trackability of the region: $\begin{array}{cc}G=\sum _{p\in \omega }g\left(p\right){g\left(p\right)}^T& \left(7\right)\end{array}$

[0063] When both the eigenvalues of the matrix G are large, the region changes in orthogonal directions, and can be uniquely positioned. Consequently, the confidence γ of the result of tracking can be obtained by the following equation (8) from the smaller eigenvalue λ_minand a gray scale residual E between the regions after the tracking: $\begin{array}{cc}\gamma =\frac{{\lambda }_{\min }}{E}& \left(8\right)\end{array}$

[0064] The inventors and others of the present invention have developed a method of interpolating a region having low confidence using the result having high confidence in the same hierarchy in an optical flow. This uses the result in a hierarchy at the next coarser resolution level for only an initial value for tracking and does not utilize the result in a hierarchy at the current resolution level which is paid attention to. It is assumed that an optical flow in a region which is hardly textured has a value close to optical flows in its surrounding regions, to fill a flow field by morphology processing.

[0065] FIG. 5 shows how flow vector dilation processing is performed.

[0066] The left drawing represents a map of the confidence of a flow vector on a gray scale. It is assumed that the blacker the map is, the higher the confidence is.

[0067] Obtained flow is first subjected to threshold processing. A white portion is subjected to threshold processing because the result has low confidence.

[0068] The dilation processing of the result in the flow field is then performed as follows in imitation of hole filling processing using a morphology operation in a binary image. A flow vector u(i, j) in coordinates i,j of a region can be calculated, as expressed by the following equation (9), upon being weighted depending on confidence γ from flow vectors in its four vicinities: $\begin{array}{cc}u\left(i,j\right)=\sum _{p,q}^{}\frac{\gamma \left(i+p,j+q\right)\times u\left(i+p,j+q\right)}{{\gamma }_A}\left(p,q\right)=\left(0,1\right),\left(0,-1\right),\left(-1,0\right),\left(1,0\right){\gamma }_A=\sum _{p,q}^{}\gamma \left(i+p,j+q\right)& \left(9\right)\end{array}$

[0069] The processing is repeated until all regions having low confidence which have been subjected to threshold processing are filled. The filling processing is performed in each hierarchy. The flow vector u(i, j) in the coordinates i,j of the region may be calculated upon being weighted depending on confidence γ from flow vectors in its eight vicinities.

[0070] FIG. 6ashows optical flow which has been subjected to threshold processing for an image in a hierarchy, and FIG. 6bshows optical flow after filling. In FIG. 6a, arrows are optical flow vectors whose confidence is judged to be high by the threshold processing, and × marks indicate portions whose confidence is judged to be low.

[0071] [2] Description of Panoramic Image Synthesizer

[0072] Description is now made of a panoramic image synthesizer for obtaining a seamless panoramic image from a plurality of images.

[0073] FIG. 7 illustrates the configuration of a panoramic image synthesizer.

[0074] A display 21, a mouse 22 and a keyboard 23 are connected to a personal computer 10. The personal computer 10 comprises a CPU 11, a memory 12, a hard disk 13, and a drive (disk drive) 14 of a removable disk such as a CD-ROM.

[0075] The hard disk 13 stores a panoramic image synthesis program in addition to an OS (Operating System). The panoramic image synthesis program is installed in the hard disk 13 using a CD-ROM 20 storing the program. It is assumed that a plurality of images picked up by a digital camera are previously stored in the hard disk 13.

[0076] FIG. 8 shows the procedure for panoramic image synthesis processing performed by the CPU 11 in a case where panoramic image synthesis software is started.

[0077] For convenience of illustration, description is made of a case where two images are synthesized.

[0078] (i) Two images (a first image and a second image) which are designated by a user are read in the memory 12 (step 1).

[0079] (ii) Processing for extracting an overlapped portion of the two images is then performed (step 2). The processing for extracting the overlapped portion is performed on the basis of the SSD (Sum of Squared Difference) method and the normalized cross-correlation method, for example. When a change between the two images, for example, two successive images picked up by a video camera, is small, the processing for extracting the overlapped portion can be omitted.

[0080] (a) Description of SSD Method

[0081] In the SSD method, with respect to two original images, images I₁and I₂having lower resolutions than those of the original images are first respectively produced. An overlapped portion ω (size M×N) of the two images I₁and I₂having low resolutions is found using a squared error (E) per pixel, as expressed by the following equation (10). An amount of movement (d) between the images is changed in an allowable range, and the overlapped portion is extracted from the amount of movement (d) in a case where E is the smallest. $\begin{array}{cc}E\left(d\right)=\frac{\sum _{\omega }{\left(I_1\left(x\right)-I_2\left(x+d\right)\right)}^2}{M\times N}& \left(10\right)\end{array}$

[0082] (b) Description of Normalized Cross-correlation Method

[0083] In the normalized cross-correlation method, with respect to two original images, images I₁and I₂having lower resolutions than those of the original images are first respectively produced. An overlapped portion ω (size :M×N) of the two images I₁and I₂having low resolutions is found using a normalized cross-correlation coefficient C, as expressed by the following equation (11). An amount of movement (d) between the images is changed in an allowable range, and the overlapped portion is extracted from the amount of movement (d) in a case where C is the largest. $\begin{array}{cc}C\left(d\right)=\frac{\sum _{\omega }\left(I_1\left(x\right)-{\stackrel{\_}{I}}_1\right)\left(I_2\left(x+d\right)-{\stackrel{\_}{I}}_2\right)}{\sqrt{{\sigma }_1·{\sigma }_2}}& \left(11\right)\end{array}$

[0084] In the equation (11), {overscore (I)}₁and {overscore (I)}₂are averages of respective gray scale values of a first image I₁and a second image I₂in the overlapped portion of both the images in a case where the first image is fixed and the second image is moved by d. δ₁and δ₂are respective variances of the gray scale values of the first image I₁and the second image I₂in the overlapped portion of both the images in a case where the first image I₁is fixed and the second image I₂is moved by d.

[0085] (iii) Feature points are then extracted (step 3). That is, a plurality of partial images (rectangular regions) which are effective for tracking are extracted as feature points from the overlapped portion of the first image. In this case, the feature points are extracted such that they are not overlapped with each other. Specifically, the above-mentioned portion having a high eigenvalue λ_min(see FIG. 8) is extracted as the feature point.

[0086] iv) The feature points are then subjected to tracking processing (step 4). That is, a position, on the second image, of the extracted feature point on the first image is tracked.

[0087] Specifically, an optical flow vector for each patch of suitable size (13×13, for example) is first found by the optical flow estimation method described in the item [1]. The position, on the second image, corresponding to the feature point on the first image is found in sub-pixel precision by linear interpolation from flow vectors in the patches in four vicinities of the feature point on the first image.

[0088] (v) A planar projective transform matrix for each triangular patch is then calculated (step 5). Description is first made of a method of finding a planar projective transform matrix in a case where it is assumed that a scene is a single plane.

[0089] It is known in projective geometry that when a point M on a single plane in a three-dimensional space is observed from two different view points C and C′, as shown in FIG. 9, transformation between coordinates m and m′ on their respective image surfaces is linear, which is referred to as homography (O. Faugeras, “Three-Dimension Computer Vision: a Geometric Viewpoint”, MIT press, 1993.).

[0090] That is, in a case where the coordinates of an image are represented by homogeneous coordinates, a point m=(x, y, 1)^ton the first image has a corresponding point m′=(x′, y′, 1)^ton the second image, and the relationship therebetween is defined by the following equation (12) and more specifically, the following equation (13): $\begin{array}{cc}{m}^{\prime }=\mathrm{Hm}& \left(12\right)\\ (\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ 1\end{array})=(\begin{array}{ccc}{h}_1& h_2& h_3\\ h_4& h_5& h_6\\ h_7& h_8& 1\end{array})(\begin{array}{c}x\\ y\\ 1\end{array})& \left(13\right)\end{array}$

[0091] The transform matrix can be rewritten, as expressed by the following equation (14): $\begin{array}{cc}{x}^{\prime }=\frac{h_{1}x+h_{2}y+h_3}{h_{7}x+h_{8}y+1}{y}^{\prime }=\frac{h_{4}x+h_{5}y+h_6}{h_{7}x+h_{8}y+1}& \left(14\right)\end{array}$

[0092] Furthermore, the equation (14) can be rewritten, as expressed by the following equation (15).

[0093] h₁x+h₂y+h₃+∅h₄+∅h₅+∅h₆−x′xh₇x′yh₈=x′

[0094] xh₄+yh₅+h₆−y′xh₇−y′yh₈=y′  (15)

[0095] The number of unknown parameters in the planar projective transform matrix H is eight, and a set of corresponding points gives two equations, as expressed by the equation (15). If there are four or more sets of corresponding points, therefore, the planar projective transform matrix H can be found by a method of least squares.

[0096] Description was made of a method of finding the planar projective transform matrix H in a case where it is assumed that a scene is a single plane. In the case of a scene of the interior of a room, for example, however, an error is too large to approximate the scene by the single plane. The reason for this is that the depth of the scene is larger than the distance from a camera to an object. As a result, when a first image and a second image are synthesized, lines to be overlapped are doubled or blurred in an overlapped portion of both the images.

[0097] In the present embodiment, a method of acquiring planar projective transform matrices Hs where a scene to be synthesized has a large depth, will be described.

[0098] Specifically, a scene is divided into triangular patches on the basis of corresponding image feature points m₁to m₈in a first image 100 and a second image 200, to calculate a planar projective transform matrix for each of the patches, as shown in FIG. 10. Each of the triangular patches is a region having three adjacent feature points as vertexes.

[0099] In order to calculate the planar projective transform matrix for each of the patches, epipolar constraint conditions are utilized. Examples of a method of calculating the planar projective transform matrix for each patch include a method utilizing a fundamental matrix representing epipolar constraint conditions and a method utilizing a correspondence of epipoles representing epipolar constraint conditions.

[0100] (a) Method Utilizing Fundamental Matrix

[0101] Description is made of,a method utilizing a fundamental matrix representing epipolar constraint conditions between cameras. Since the degree of freedom of the planar projective transform matrix is eight, four sets of corresponding points are required. Since the number of vertexes of a triangular patch is three, however, the plane projective transform matrix is not generally solved.

[0102] Therefore, the inventors of the present invention have developed a method of calculating a planar projective transform matrix for each triangular patch using the fundamental matrix F representing epipolar constraint conditions between cameras and three sets of corresponding points.

[0103] A method of calculating a planar projective transform matrix for each triangular patch will be described in detail.

[0104] FIG. 11 illustrates the epipolar constraint conditions.

[0105] When a point M in a three-dimensional space is observed from different view points C₁and C₂, coordinates ml and m₂on respective imaging surfaces S₁and S₂of cameras must be on a plane composed of M, C₁and C₂. This is referred to as epipolar constraint conditions. Further, respective points of intersection e₁and e₂of a line connecting the view points C₁and C₂and the imaging surfaces S₁and S₂are referred to as epipoles. It is known that a fundamental matrix F representing the epipolar constraint conditions and the corresponding epipoles e₁and e₂can be calculated if there are seven or more sets of corresponding points between two images.

[0106] The epipolar constraint conditions between the two images are expressed by the following equation (16) using the fundamental matrix F and corresponding points m and m′:

[0107] m′^TFm=0  (16)

[0108] The fundamental matrix F is found from seven or more sets of, for example, eight sets of corresponding points. A method of finding the fundamental matrix F has been well known and hence, the description thereof is not repeated (see “Three-dimensional Vision” issued by Kyoritz Publishing Co., Ltd. on Apr. 20, 1998, for example). In recent years, a method of satisfactorily acquiring a fundamental matrix of a stereoscopic camera, which is not calibrated, from a pair of feature points, which correspond to each other, has been proposed. Therefore, it is preferable that the method is used (Z. Zhang, “Determining the Epipolar Geometry and its Uncertainty; A Review”, Int. Journal of Computer Vision, Vol. 27, No. 2, pp. 161-195, 1988).

[0109] Since m′=Hm, as expressed by the foregoing equation (12), the equation (12) is substituted in the equation (16), to obtain the following equation (17):

[0110] (Hm)^TFm=m^TH^TFm=0  (17)

[0111] A matrix H^TF in the foregoing equation (17) must be a skew-symmetric matrix, as expressed by the following equation (18), because it is an outer product of a vector m: $\begin{array}{cc}{H}^{T}F=(\begin{array}{ccc}0& -a_1& a_2\\ a_1& 0& -a_3\\ -a_2& a_3& 0\end{array})& \left(18\right)\end{array}$

[0112] Six equations are obtained from the foregoing equation (15) on the basis of three sets of corresponding points. Further, six equations are obtained from the foregoing equation (18). That is, 12 equations are obtained. The number of unknowns in a planar projective transform matrix H is eight, thereby making it possible to calculate the planar projective transform matrix, H for each triangular patch by a method of least squares using the 12 equations.

[0113] Specifically, when the three sets of corresponding points are taken as {(m₁, m₁′), (m₂, m₂′), (m₃, m₃′)}, six equations are obtained, as expressed by the following equation (19), on the basis of the foregoing equation (15): $\begin{array}{cc}(\begin{array}{cccccccc}{x}_1& y_1& 1& 0& 0& 0& -x_1^{\prime }{x}_1& -x_1^{\prime }{y}_1\\ 0& 0& 0& x_1& y_1& 1& -y_1^{\prime }{x}_1& -y_1^{\prime }{y}_1\\ x_2& y_2& 1& 0& 0& 0& -x_2^{\prime }{x}_2& -x_2^{\prime }{y}_2\\ 0& 0& 0& x_2& y_2& 1& -y_2^{\prime }{x}_2& -y_2^{\prime }{y}_2\\ x_3& y_3& 1& 0& 0& 0& -x_3^{\prime }{x}_3& -x_3^{\prime }{y}_3\\ 0& 0& 0& x_3& y_3& 1& -y_3^{\prime }{x}_3& -y_3^{\prime }{y}_3\end{array})(\begin{array}{c}{h}_1\\ h_2\\ h_3\\ h_4\\ h_5\\ h_6\\ h_7\\ h_8\end{array})=(\begin{array}{c}{x}_1^{\prime }\\ y_1^{\prime }\\ x_2^{\prime }\\ y_2^{\prime }\\ x_3^{\prime }\\ y_3^{\prime }\end{array})& \left(19\right)\end{array}$

[0114] The foregoing equation is expressed by the following equation (20): $\begin{array}{cc}\begin{array}{c}{H}^{T}F=(\begin{array}{ccc}{h}_1& h_4& h_7\\ h_2& h_5& h_8\\ h_3& h_6& 1\end{array})(\begin{array}{ccc}{f}_1& f_2& f_3\\ f_4& f_5& f_6\\ f_7& f_8& f_9\end{array})=(\begin{array}{ccc}0& -a_1& a_2\\ a_1& 0& -a_3\\ -a_2& a_3& 0\end{array})\\ =(\begin{array}{ccc}{h}_1{f}_1+h_4{f}_4+h_7{f}_7& h_1{f}_2+h_4{f}_5+h_7{f}_8& h_1{f}_3+h_4{f}_6+h_7{f}_9\\ h_2{f}_1+h_5{f}_4+h_8{f}_7& h_2{f}_2+h_5{f}_5+h_8{f}_8& h_2{f}_3+h_5{f}_6+h_8{f}_9\\ h_3{f}_1+h_6{f}_4+f_7& h_3{f}_2+h_6{f}_5+f_8& h_3{f}_3+h_6{f}_6+f_9\end{array})\end{array}& \left(20\right)\end{array}$

[0115] Three equations are obtained, as expressed by the following equation (21) from conditions under which a diagonal matrix is zero:

[0116] h₁f₁+h₄f₄+h₇f₇=0

[0117] h₂f₂+h₅f₅+h₈f₈=0

[0118] h₃f₃+h₆f₆+f₉f₉=0  (21)

[0119] One equation is obtained, as expressed by the following equation (22), from conditions under which −a1+a1=0 in the skew-symmetric matrix in the foregoing equation (20):

[0120] f₂h₁+f₁h₂+f₅h₄+f₄h₅+f₈h₇+f₇h₈=0  (22)

[0121] One equation is obtained, as expressed by the following equation (23), from conditions under which a2−a2=0 in the skew-symmetric matrix in the foregoing equation (20):

[0122] f₃h₁+f₁h₃+f₆h₄+f₄h₆+f₉h₇=−f₇  (23)

[0123] One equation is obtained, as expressed by the following equation (24), from:conditions under which −a3+a3=0 in the skew-symmetric matrix in the foregoing equation (20):

[0124] f₃h₂+f₂h₃+f₆h₅+f₅h₆+f₉h₈=−f₈  (24)

[0125] Six equations are obtain,ed, as expressed by the following equation (25), from the foregoing equations (21) to (24): $\begin{array}{cc}(\begin{array}{cccccccc}{f}_1& 0& 0& f_4& 0& 0& f_7& 0\\ 0& f_2& 0& 0& f_5& 0& 0& f_8\\ 0& 0& f_3& 0& 0& f_6& 0& 0\\ f_2& f_1& 0& f_5& f_4& 0& f_8& f_7\\ f_3& 0& f_1& f_6& 0& f_4& f_9& 0\\ 0& f_3& f_2& 0& f_6& f_5& 0& f_9\end{array})(\begin{array}{c}{h}_1\\ h_2\\ h_3\\ h_4\\ h_5\\ h_6\\ h_7\\ h_8\end{array})=(\begin{array}{c}0\\ 0\\ -f_9\\ 0\\ -f_7\\ -f_8\end{array})& \left(25\right)\end{array}$

[0126] When the foregoing equations (19) and (25) are made simultaneous, 12 equations are obtained, as expressed by the following equation (26): $\begin{array}{cc}(\begin{array}{cccccccc}{x}_1& y_1& 1& 0& 0& 0& -x_1^{\prime }{x}_1& -x_1^{\prime }{y}_1\\ 0& 0& 0& x_1& y_1& 1& -y_1^{\prime }{x}_1& -y_1^{\prime }{y}_1\\ x_2& y_2& 1& 0& 0& 0& -x_2^{\prime }{x}_2& -x_2^{\prime }{y}_2\\ 0& 0& 0& x_2& y_2& 1& -y_2^{\prime }{x}_2& -y_2^{\prime }{y}_2\\ x_3& y_3& 1& 0& 0& 0& -x_3^{\prime }{x}_3& -x_3^{\prime }{y}_3\\ 0& 0& 0& x_3& y_3& 1& -y_3^{\prime }{x}_3& -y_3^{\prime }{y}_3\\ f_1& 0& 0& f_4& 0& 0& f_7& 0\\ 0& f_2& 0& 0& f_5& 0& 0& f_8\\ 0& 0& f_3& 0& 0& f_6& 0& 0\\ f_2& f_1& 0& f_5& f_4& 0& f_8& f_7\\ f_3& 0& f_1& f_6& 0& f_4& f_9& 0\\ 0& f_3& f_2& 0& f_6& f_5& 0& f_9\end{array})(\begin{array}{c}{h}_1\\ h_2\\ h_3\\ h_4\\ h_5\\ h_6\\ h_7\\ h_8\end{array})=(\begin{array}{c}{x}_1^{\prime }\\ y_1^{\prime }\\ x_2^{\prime }\\ y_2^{\prime }\\ x_3^{\prime }\\ y_3^{\prime }\\ 0\\ 0\\ -f_9\\ 0\\ -f_7\\ -f_8\end{array})& \left(26\right)\end{array}$

[0127] Since the number of unknowns in the planar projective transform matrix H is eight, the planar projective transform matrix H corresponding to the triangular patch is calculated using the method of least squares.

[0128] (b) Method Utilizing Correspondence Between Epipoles

[0129] Description is made of a method utilizing a correspondence between epipoles. The correspondence between epipoles is found from a fundamental matrix F. This method is a method of finding a planar projective transform matrix for each triangular patch from:sets of corresponding points at three vertexes of a triangular patch between two images and corresponding points of the epipoles between the images.

[0130] FIG. 12 illustrates a triangular patch having N, M, and P as vertexes, and epipoles e₁and e₂.

[0131] The epipoles e₁and e₂are always the same irrespective of the position and the posture of the triangular patch. A total of four sets of corresponding points are found from the sets of corresponding points at the three vertexes of the triangular patch between the two images (three sets of corresponding points) and corresponding points of the eipoles between the two images (one set of corresponding points).

[0132] Eight equations are obtained, expressed by the following equation (27), on the basis of the four sets of corresponding points (x_i, yi) and (x_i′, y_i′) (i=1, 2, 3, e) and the foregoing equation (15); $\begin{array}{cc}(\begin{array}{cccccccc}{x}_1& y_1& 1& 0& 0& 0& -x_1^{\prime }{x}_1& -x_1^{\prime }{y}_1\\ 0& 0& 0& x_1& y_1& 1& -y_1^{\prime }{x}_1& -y_1^{\prime }{y}_1\\ x_2& y_2& 1& 0& 0& 0& -x_2^{\prime }{x}_2& -x_2^{\prime }{y}_2\\ 0& 0& 0& x_2& y_2& 1& -y_2^{\prime }{x}_2& -y_2^{\prime }{y}_2\\ x_3& y_3& 1& 0& 0& 0& -x_3^{\prime }{x}_3& -x_3^{\prime }{y}_3\\ 0& 0& 0& x_3& y_3& 1& -y_3^{\prime }{x}_3& -y_3^{\prime }{y}_3\\ x_e& y_e& 1& 0& 0& 0& -x_e^{\prime }{x}_e& -x_e^{\prime }{y}_e\\ 0& 0& 0& x_e& y_e& 1& -y_e^{\prime }{x}_e^{\prime }& -y_e^{\prime }{y}_e\end{array})(\begin{array}{c}{h}_1\\ h_2\\ h_3\\ h_4\\ h_5\\ h_6\\ h_7\\ h_8\end{array})=(\begin{array}{c}{x}_1^{\prime }\\ y_1^{\prime }\\ x_2^{\prime }\\ y_2^{\prime }\\ x_3^{\prime }\\ y_3^{\prime }\\ x_e^{\prime }\\ y_e^{\prime }\end{array})& \left(27\right)\end{array}$

[0133] Since the number of elements (h₁to h₈) in the plane projective transform matrix is eight, the elements h₁to h₈are found from the eight equations.

[0134] (vi) An image is produced by overlapping both the two images with each other such that the positions of the two images coincide with each other on the basis of the found planar r projective transform matrix for each triangular patch (step 6).

[0135] Specifically, a coordinate value, on the first image, corresponding to each of pixels in each of the triangular patches on the second image is found using a planar projective transform matrix H calculated for the triangular patch. On the other hand, a coordinate value, on the first image, corresponding to each of pixels which do not belong to any of the triangular patches on the second image is found using a planar projective transform matrix H calculated for the triangular patch closest to the pixel. The coordinate value, on the first image, corresponding to each of the pixels which do not belong to any of the triangular patches on the second image may be found using a planar projective transform matrix calculated using all corresponding points in an overlapped portion of both the images without dividing a scene into triangles. The planar projective transform matrix calculated using all the corresponding points in the overlapped portion of both the images is a planar projective transform matrix calculated assuming that the scene is a plane, and corresponds to an average of the plane projective transformation matrices for the triangular patches.

[0136] The first image and the second image are synthesized using the coordinate value, on the first image, corresponding to each of the pixels on the second image.

[0137] At the foregoing step 5, a transform matrix for transforming a point on the first image into a position on the second image is found. Therefore, an image obtained by synthesizing the first image and the second image is obtained using the second image as a reference image.

[0138] At the foregoing step 5, a transform matrix for transforming a point, on the second image, corresponding to each of feature points on the first image into a position on the first image may be found on the basis of the point on the second image. When the transform matrix for transforming the point on the second image into the position on the first image is found, an image obtained by synthesizing the first image and the second image is obtained using the first image as a reference image.

[0139] (vii) Thereafter, pixel values are blended with respect to an overlapped portion of the two images (step 7). That is, the first image and the second image are not picked up under the same conditions, so that one of the images may, in some cases, be darker than the other image. Therefore, the pixel values of both the images are blended in the overlapped portion of the images.

[0140] An example of a method of blending pixel values will be described. FIG. 13 illustrates an image produced at the step 6. A hatched portion shown in FIG. 13 is a portion where a first image 31 and a second image 32 are overlapped with each other (an overlapped portion). In FIG. 13, G1 and G2 respectively indicate the position of the center of gravity of the first image and the position of the center of gravity of the second image.

[0141] Letting I₁be a pixel value, on the first image 31, of an arbitrary point P in the overlapped portion, I₂be a pixel value, on the second image 32, of the point P, d₁be the distance between the point P and the position of the center of gravity G₁of the first image, and d₂be the distance between the point P and the position of the center of gravity G₂of the second image, a pixel value Blen, after the blending, of the point P is expressed by the following equation (28):

[0142] $\begin{array}{cc}\mathrm{Blen}=\frac{d_2·I_1+d_1·I_2}{d_1+d_2}& \left(28\right)\end{array}$

[0143] Although in the above-mentioned embodiment, the planar projective transform matrix H is found for each of the triangular patches, and the coordinate value, on the first image, corresponding to each of the pixels in the triangular patch on the second image is found using the planar projective transform matrix H calculated for the triangular patch, an affine transform matrix A may be found for each of the triangular patches, and the coordinate value, on the first image, corresponding to each of the pixels in the triangular patch on the second image may be found using the affine transform matrix A calculated for the triangular patch.

[0144] The affine transform matrix A is expressed by the following equation (29). Since the number of elements h₁to h₆in the affine transform matrix A is six, the elements h₁to h₆can be calculated from only corresponding points at the three vertexes of the triangular patch between the two images: $\begin{array}{cc}A=(\begin{array}{ccc}{h}_1& h_2& h_3\\ h_4& h_5& h_6\\ 0& 0& 1\end{array})& \left(29\right)\end{array}$

[0145] The relationship among the elements h₁to h₆in the transform matrix A and corresponding points (x_i, y_i) and (x_i′, y_i′) (i=1 , 2 , 3) at the three vertexes of the triangular patch between the two images is expressed by the following equation (30): $\begin{array}{cc}(\begin{array}{cccccc}{x}_1& y_1& 1& 0& 0& 0\\ 0& 0& 0& x_1& y_1& 1\\ x_2& y_2& 1& 0& 0& 0\\ 0& 0& 0& x_2& y_2& 1\\ x_3& y_3& 1& 0& 0& 0\\ 0& 0& 0& x_3& y_3& 1\end{array})(\begin{array}{c}{h}_1\\ h_2\\ h_3\\ h_4\\ h_5\\ h_6\end{array})=(\begin{array}{c}{x}_1^{\prime }\\ y_1^{\prime }\\ x_2^{\prime }\\ y_2^{\prime }\\ x_3^{\prime }\\ y_3^{\prime }\end{array})& \left(30\right)\end{array}$

[0146] Although the present invention has been described and illustrated in detail, it is clearly understood that the same is by way of illustration and example only and is not to be taken by way of limitation, the spirit and scope of the present invention being limited only by the terms of the appended claims.