JP2017059048A - Arithmetic unit and operation method - Google Patents

Abstract

PROBLEM TO BE SOLVED: To make it possible to quickly and highly accurately estimate the least square solution of a projective transformation matrix.SOLUTION: Provided is an arithmetic unit which derives an M-th dimensional projective transformation matrix by solving simultaneous equations set based on coordinate point pairs before and after the projective transformation, the arithmetic unit comprising: determination means for determining 2N simultaneous equations set based on N coordinate point pairs, the simultaneous equations being represented as Ax=b, by a coefficient matrix A, a solution vector b, and a variable vector x constituted of the elements of the transformation matrix; breakdown means for breaking down the coefficient matrix A into a partial orthogonal matrix Qwith inter-vector orthogonalization partially performed in the coefficient matrix A and an upper triangular matrix R; and derivation means for determining the variable vector x by processing a deformed simultaneous equations obtained by multiplying both sides of the simultaneous equations by a transposition matrix Qof the partial orthogonal matrix Qto derive the elements of the projective transformation matrix.SELECTED DRAWING: Figure 4

Description

  The present invention relates to a technique for estimating a projective transformation matrix between two images.

  The affine transformation matrix having the degree of freedom “6” can enlarge, reduce, rotate, and translate a two-dimensional image. A two-dimensional projective transformation matrix is obtained by adding two degrees of freedom to express the depth (perspective) in the horizontal direction and the vertical direction. The affine transformation matrix in which the transformation target is expanded to a three-dimensional space has a degree of freedom of “12” and the projective transformation matrix has a degree of freedom of “15”, and transformation processing using these is used for computer graphics (CG) generation. .

  On the other hand, in the field of computer vision, a transformation matrix is also estimated from a plurality of images. For example, for the purpose of camera shake correction of a moving image shooting camera, a transformation matrix is estimated from a plurality of images shot by a camera or the like and correction is performed. A function for obtaining a two-dimensional projective transformation matrix from coordinate data of corresponding points in two images is used in many applications, and this transformation matrix is called a homography matrix. The homography matrix can be used in a wide range of applications such as panoramic image generation, projector projection position correction, camera calibration using a plane pattern, object shape estimation / restoration, and obstacle detection on a road surface.

  However, the processing for obtaining a homography matrix by solving simultaneous equations based on projective transformation has a large amount of processing. Therefore, Patent Document 1 discloses a technique for selectively using two types, a local homography matrix obtained by an approximation operation effective when the correlation between images is high, and a global homography matrix obtained by solving simultaneous equations. Specifically, when the correlation between images is high, a local homography matrix capable of high-speed calculation is calculated, and only when the correlation is low, simultaneous equations are solved to obtain a global homography matrix.

JP2012-86285A

However, if there are many corresponding points when creating simultaneous equations, numerical calculation errors tend to accumulate in the product-sum operation when calculating the matrix product (A ^T A) of the coefficient matrix A. On the other hand, as a technique for avoiding error accumulation, there is a method using QR decomposition of a matrix, but it is necessary to orthogonalize the column vectors of the coefficient matrix A to each other, which is very heavy processing.

  The present invention has been made in view of such a problem, and an object of the present invention is to provide a technique capable of estimating a least-squares solution of a projective transformation matrix with high speed and accuracy.

In order to solve the above-described problems, the arithmetic device according to the present invention has the following configuration. That is, in an arithmetic device that derives a transformation matrix of M-dimensional projective transformation by solving simultaneous equations formed based on coordinate point pairs before and after transformation by the projective transformation, based on N coordinate point pairs. 2N simultaneous equations (where 2N> (M + 1) ² −1), and Ax = b by a variable vector x composed of a coefficient matrix A, a solution vector b, and each element of the transformation matrix A determination means for determining simultaneous equations expressed as follows: the coefficient matrix A is decomposed into a partial orthogonal matrix Q _p and an upper triangular matrix R _p partially orthogonalized between column vectors in the coefficient matrix A And determining the variable vector x by processing the modified simultaneous equations obtained by multiplying both sides of the simultaneous equations by the transposed matrix Q _p ^T of the partial orthogonal matrix Q _p . Deriving means for deriving each element of the transformation matrix, and at least one set of column vectors in the partial orthogonal matrix Q _p are not orthogonal to each other.

Alternatively, in an arithmetic unit that derives a transformation matrix of M-dimensional projective transformation by solving simultaneous equations that are formulated based on coordinate point pairs before and after transformation by the projective transformation, based on N coordinate point pairs 2N (2N> (M + 1) ² −1) simultaneous equations that are formulated, and N rows ((M + 1) ² −1) columns of coefficient matrix A _h and N rows and M columns of solution matrix B , A determining means for determining a simultaneous equation expressed as A _h X = B by a variable matrix X of ((M + 1) ² −1) rows and M columns composed of each element of the transformation matrix, and the coefficient matrix A _h and a decomposing means for decomposing into a orthogonalizing partially portion orthogonal matrix was performed Q _p and an upper triangular matrix R _p between column vectors in the coefficient matrix a _h, transpose matrix Q _p of the partial orthogonal matrix Q _p deformation communicating obtained by multiplying the ^T on both sides of the simultaneous equations The solution matrix B determined by processing an equation, orthogonal the a derivation means for deriving the respective elements of the transformation matrix of the projective transformation, has at least one set of row vectors in the partial orthogonal matrix Q _p are each Not done.

  ADVANTAGE OF THE INVENTION According to this invention, the technique which makes it possible to estimate the least squares solution of a projective transformation matrix quickly and accurately can be provided.

It is a figure which shows the coefficient matrix and transposed matrix of simultaneous equations for calculating | requiring a homography matrix exemplarily. It is a figure which shows the coefficient matrix of a normal equation exemplarily. It is a figure explaining transition of a coefficient matrix in a 1st embodiment. It is a flowchart of the process in 1st Embodiment. It is a figure explaining transition of a coefficient matrix in a 2nd embodiment. It is a flowchart of the process in 2nd Embodiment. It is a figure explaining transition of a coefficient matrix in a 3rd embodiment. It is a flowchart of the process in 3rd Embodiment. It is an internal block diagram of the arithmetic unit comprised by PC.

  Hereinafter, preferred embodiments of the present invention will be described in detail with reference to the drawings. The following embodiments are merely examples, and are not intended to limit the scope of the present invention.

(First embodiment)
As a first embodiment of an arithmetic device according to the present invention, an arithmetic device that estimates a least square solution such as a homography matrix or a three-dimensional projective transformation matrix at high speed and with high accuracy will be described. Specifically, focusing on the structure of the coefficient matrix of the simultaneous equations, the least square estimation is performed by partially orthogonalizing the column vectors. This enables high-precision and high-speed estimation.

<Prerequisite technology>
First, the principle of camera shake correction of a moving image shooting camera realized using a homography matrix and how to obtain the homography matrix will be described.

  A blur between consecutive frames obtained by moving image shooting is detected as a geometric transformation matrix, and a matrix obtained by smoothing the geometric transformation matrix over a dozen frames is derived. By inversely correcting the image using the matrix derived in this way, blurring between frames and viewpoint movement can be smoothed, and a camera shake correction function during moving image shooting can be realized. In particular, when the distance from the camera to the subject is sufficiently far away and it can be considered that each point of interest of the subject exists on the same plane, the geometric transformation matrix that represents the blur between frames should be represented by a homography matrix. Can do. At that time, the coordinates of the pixel between the temporally preceding frame and the succeeding frame are associated by Expression (1).

  The matrix of 3 rows and 3 columns on the right side in Equation (1) is a homography matrix, and the pixel (coordinate point before conversion) of the preceding frame is the pixel (x, y) coordinate ( Converted to a coordinate point after conversion. Hereinafter, these two coordinate points before and after conversion are referred to as a coordinate point pair. Further, the expression (1) is expressed by a conversion expression using homogeneous coordinates in order to express the nonlinear projective transformation like a linear operation.

  When the 9th element in the lower right is normalized (normalized) to be “1”, the other 8 elements can take free values, and a matrix with 8 degrees of freedom is obtained. On the other hand, when the coordinates between the two images represented by the relationship of Expression (1) are developed, two equations shown in Expression (2) are obtained.

  Eight equations, that is, eight simultaneous equations can be determined and formulated from four coordinate point pairs, and one homography matrix can be obtained by solving these simultaneous equations. However, the four coordinate point pairs need to exist on the same plane in the subject space, and if a coordinate point pair that is not on the same plane is included, it becomes a meaningless homography matrix. Therefore, it is necessary to evaluate how much the solution is appropriate by converting the coordinate data of other coordinate point pairs.

On the other hand, when obtaining a homography matrix from more than 5 sets of coordinate points, it becomes an over-conditioning simultaneous equation with more than 8 equations, and there are no solutions that satisfy all equations, and the equations cannot be solved. May end up. In M-dimensional projective transformation, when N or more (2N> (M + 1) ² −1) coordinate point pairs are used, 2N simultaneous equations are formed to form over-conditional simultaneous equations.

In this case, the matrix ^AT obtained by transposing the coefficient matrix A of the equation is multiplied on both sides, that is, the coefficient matrix A and the solution vector b. As a result, the coefficient matrix and the solution vector are projected onto an eight-dimensional space and converted into an 8-ary simultaneous equation called a “normal equation” to calculate a solution. When the variable vector is represented by x, the normal equation is expressed as (A ^T A) x = A ^T b. Here, (A ^T A) is a coefficient matrix of the normal equation, and the size of the matrix is 8 × 8. A ^T b on the right side is an 8-dimensional vector.

  A solution obtained by solving a normal equation is a least square method solution (solution based on the least square method) for five or more pairs of coordinate points, and is a solution reflecting all the coordinate point pairs. The method of obtaining a least square solution using this normal equation is a general method often used in statistical processing.

<Device configuration>
FIG. 9 is a block diagram illustrating a hardware configuration of the arithmetic device. Here, since an example in which the arithmetic unit 101 is configured by a general personal computer (PC) is shown, the display 123, the communication I / F 124, the keyboard 125, and the like are included, but these are not essential configurations. .

  In addition to executing various calculations, the CPU 120 comprehensively controls the calculation device 101. The CPU 120 implements each process shown in FIG. 4 by executing a calculation program stored in, for example, the ROM 122 or the hard disk drive (HDD) 126.

  The HDD 126 stores, for example, a calculation program or various control programs used in the calculation device 101. Also, various information related to the arithmetic program or various control programs is stored. Here, various information includes information related to simultaneous equations. A RAM 121 is also used to temporarily store various information.

  The keyboard 125 is a functional unit that accepts input of various information used for calculation from the user. The display 123 is a functional unit that provides the user with various types of information such as the progress of calculation processing and calculation results. The communication interface (I / F) 124 is an interface for connecting to an external device (not shown), and is, for example, an interface for wired communication or wireless communication. Via the I / F 124, for example, an image of a plurality of frames to be processed can be input from an external imaging device or the like, or a pair of feature points can be input. Further, the obtained projective transformation matrix can be output to an external image correction device or the like via the I / F 124.

<Overview of homing matrix derivation>
First, a specific example of an overconditional simultaneous equation for N pairs (N is 5 or more) of coordinate point pairs and its normal equation will be shown. As described above, 2N 8-ary simultaneous equations can be determined from N coordinate point pairs. This simultaneous equation is expressed as Ax = b using a coefficient matrix A, a variable vector x, and a solution vector b.

FIG. 1 is a diagram exemplarily showing a coefficient matrix / solution vector of simultaneous equations for obtaining a homography matrix and a transposed matrix of the coefficient matrix. FIG. 1 (A) represents 2N simultaneous equations in matrix form. FIG. 1B shows the matrix of FIG. 1A rearranged based on the solution vector. Hereinafter, the matrix and the column vector shown in FIG. 1B are referred to as a coefficient matrix A and a solution vector b, respectively. FIG. 1C shows ^AT , which is a transposed matrix of the coefficient matrix A.

  FIG. 2 is a diagram exemplarily showing a coefficient matrix of a normal equation. The matrix shown in FIG. 2 (A) and the column vector shown in FIG. 2 (B) are obtained by multiplying each of the coefficient matrix A and the solution vector b by the transposed matrix of FIG. is there.

That is, each corresponds to a coefficient matrix (A ^T A) and a solution vector (A ^T b) existing on both sides of the normal equation. The symbol Σ represents the sum of subscripts i = 1, 2,. This coefficient matrix is a symmetric matrix. That is, in FIG. 2A, the lower left 2 rows and 6 columns are described using only the Σ symbol, but are the same as the elements obtained by transposing the upper right 2 columns and 6 rows.

  As can be seen from FIG. 2A, zero elements are regularly arranged in a normal equation for obtaining a homography matrix, and a redundant portion exists in the coefficient matrix.

In the QR decomposition method, the coefficient matrix A is decomposed into an orthogonal matrix Q and an upper triangular matrix R. On the other hand, in the first embodiment, the coefficient matrix A is decomposed into partially orthogonal matrices Q _p partially orthogonalized and an upper triangular matrix R _p corresponding to Q _p . Then, by multiplying Q _{p by} the transposed matrix of Q _p , the zero element of the obtained 8 × 8 matrix is further increased compared to the matrix of FIG. Thus, the equation is more easily solved than the normal equation shown in FIG.

  In particular, in the first embodiment, in the four column vectors excluding the first and fourth columns in the coefficient matrix shown in FIG. 1 (B), the average of each of the upper N and lower N is calculated. Calculate and subtract these averages from each element. There is no need to do anything where N "0" s are lined up. In the M-dimensional projective transformation, the same processing is performed on 2M column vectors including the values of the coordinate points before conversion included in each of the N coordinate point pairs.

Each column vector obtained by subtracting the average value is orthogonal to the first column and the fourth column. _Let Q _p be the matrix processed in this way. That is, column vectors in a general QR decomposition are orthogonal to each other. In contrast, the column vectors of Q _p in the first embodiment is only orthogonal to the particular column, there is a set of column vectors which are not orthogonal to each other.

FIG. 3 is a diagram for explaining the transition of the coefficient matrix in the first embodiment. Specifically, this is equivalent to a process of performing a similar process on both sides of a normal equation that is a simultaneous equation and generating a modified simultaneous equation. 3 (A) is a coefficient matrix A, a representation by the product of the above Q _p and an upper triangular matrix R _p. That is, the upper triangular matrix shown on the right side shown in FIG. 3 (A) is R _p. In FIG. 3, the result of subtracting the above average from the u coordinate and the v coordinate is expressed as du and dv.

Then, multiplying the transposed matrix _Q ^{p T} of _{Q p} into a column vector b of the left side of _{Q p} and the right side. At this time, the calculation result of Q _p ^T Q _p is an 8 × 8 symmetric matrix as shown in FIG. 3B. Since it is a symmetric matrix, the description of the lower left triangular matrix part is omitted. Then, the following processing is performed on the 8 × 8 coefficient matrix. Note that the following process corresponds to a process of multiplying an inverse matrix of Q _p ^T Q _p as a result.

  First, the reciprocal of N, that is, 1 / N is calculated and multiplied by the first and fourth rows. Next, in order to diagonalize the 2 × 2 region in the vicinity of the diagonal elements in the 2nd, 3rd and 5th and 6th rows, an inverse matrix of the 2 × 2 matrix is obtained, and using this inverse matrix, Each column in the third row and the fifth and sixth rows is converted as a secondary column vector. Note that the 2 × 2 region in the vicinity of the diagonal element is known to become a unit matrix when converted, and thus conversion is not necessary.

  That is, column vectors that must be actually converted by the inverse matrix described above exist only in the seventh column, the eighth column, and the solution vector. By the conversion by the 2 × 2 inverse matrix, the 8 × 8 coefficient matrix becomes a matrix obtained by unitizing the upper left 6 × 6 region.

  Next, the elements in the first to sixth columns of the seventh and eighth rows are erased to zero. The sum of products is calculated between the four significant elements in the seventh row and the four significant elements in the seventh column, and is subtracted from the element in the seventh row and the seventh column where the row and the column intersect. This sum-of-product / subtraction process is similarly performed between the seventh row and the eighth column, the eighth row and the seventh column, and the eighth row and the eighth column.

  With the above processing, the elements in the first to sixth columns in the seventh and eighth rows are zero-erased. Next, the remaining 2 × 2 elements in the seventh and eighth columns are diagonalized. Here, it is only necessary to transform the elements of the corresponding solution vector with the 2 × 2 inverse matrix described above. Thereby, each element of the variable vector x can be determined.

<Operation of arithmetic unit>
FIG. 4 is a flowchart of processing in the first embodiment. The following process is started, for example, by inputting a coordinate point pair group between a temporally preceding frame and a subsequent frame.

  In step S401, the CPU 120 selects N sets (N is 5 or more) of coordinate data from the coordinate point pair group, and generates a coefficient matrix A and a solution vector b in the format shown in FIG.

In step S402, the CPU 120 decomposes the 2N rows × 8 columns coefficient matrix A into a matrix Q _p and a matrix R _{p in the} format shown in FIG. As described above, Q _p is obtained by calculating and subtracting the average of the upper N pieces and the lower N pieces in the column vector excluding the first and fourth columns of the coefficient matrix A.

In step S403, CPU 120 is a transposed matrix _Q ^{p T} of _{Q p,} multiplies the _{Q p} and b.

  In step S404, the CPU 120 calculates an inverse of N and an inverse matrix of 2 × 2 matrix. In step S405, the CPU 120 multiplies the matrix shown in FIG. 3B by the reciprocal number of N in the 1st and 4th rows and the 2 × 2 inverse matrix in the 2nd, 3rd and 5th and 6th rows. In step S406, the CPU 120 zero-erases the first and sixth columns of the seventh and eighth rows with respect to the matrix shown in FIG.

In step S407, the CPU 120 solves the binary simultaneous equations from the remaining coefficients existing in the seventh and eighth rows, and calculates two parameters. In step S408, CPU 120 performs a backward substitution processing for deleting an element that remains in 7,8 rows of _{Q p.} In step S409, CPU 120 performs a backward substitution process of erasing the non-diagonal elements of the upper triangular matrix _{R p.}

Through this series of processes, one homography matrix estimated by least squares can be obtained. In the above-described processing, only significant coefficients that are not zero are processed efficiently, so that the processing steps are divided into fine steps. However, it is also possible to process Q _p ^T Q _p shown in FIG. 3B by a well-known method such as Gaussian elimination or Gauss-Jordan method. In that case, the above-described steps S404 to S408 are replaced with these methods. In the above description, the zero element is also processed without distinction. However, the zero element may be detected and the process may be skipped.

  Note that the coordinate point pair selected in S401 may include outliers that deviate from the premise. The outlier means a coordinate point pair having a difference (deviation) larger than a predetermined threshold between coordinates obtained as a result of transforming coordinates in a preceding frame by a homography matrix and coordinates in a corresponding subsequent frame. .

  Therefore, some or all of the coordinate data may be exchanged to obtain a number of least square estimation matrices. That is, the processing shown in FIG. 4 is performed for N coordinate point pairs of various combinations, the degree of deviation by the conversion is evaluated, the best parameter for evaluation is selected, and the final homography matrix is derived. Also good.

  As described above, according to the first embodiment, the least square solution such as a homography matrix or a three-dimensional projective transformation matrix can be obtained at high speed by orthogonalizing only a part of column vectors that can be orthogonalized easily. It becomes possible to estimate with high accuracy. That is, orthogonal processing is not performed between at least one set of column vectors. This makes it possible to evaluate a large number of parameters calculated within a unit time and to select a highly reliable conversion matrix. As a result, for example, when processing moving image data, many parameters can be evaluated within one frame time, so that a homography matrix with higher accuracy can be derived. In addition, it is possible to improve performances such as image composition, alignment, camera image stabilization, and three-dimensional shape estimation / restoration of an object.

(Second Embodiment)
In the second embodiment, as in the first embodiment, the coefficient matrix A is decomposed into a partially orthogonal matrix Q _p partially orthogonalized and an upper triangular matrix R _p corresponding to the Q _p . However, the calculation method of Q _p is different from that of the first embodiment.

Specifically, in the coefficient matrix A, after calculating and subtracting partial averages for the second and third columns to make the first column orthogonal, the second and third columns are further orthogonalized. , A partial orthogonal matrix Q _p . In this case, when calculating Q _p ^T Q _p for the 4th, 5th and 6th columns, it is sufficient to use the results of the 1st, 2nd and 3rd columns, so the processing can be omitted.

FIG. 5 is a diagram for explaining the transition of the coefficient matrix in the second embodiment. FIG. 5 (A), the coefficient matrix A, a representation by the product of the above Q _p and an upper triangular matrix R _p. As in FIG. 3, the result of subtracting the above-mentioned partial average from the u coordinate and the v coordinate is denoted as du and dv.

FIG. 5B shows a matrix obtained by Q _p ^T Q _p . Although the non-zero elements in the seventh and eighth columns are increasing, the upper left 6 × 6 submatrix is already a diagonal matrix immediately after the calculation of Q _p ^T Q _p . Therefore, in the matrix shown in FIG. 5B, this 6 × 6 submatrix can be unitized by dividing each of the first to sixth rows by the value of the diagonal element.

  Next, the elements in the first to sixth columns of the seventh and eighth rows of the matrix shown in FIG. First, the sum of products is calculated between the six elements in the seventh row and the first to sixth columns and the six elements in the seventh column and the first to sixth rows, and is subtracted from the element in the seventh row and the seventh column where the row and the column intersect. . This sum-of-product / subtraction process is similarly performed between the seventh row and the eighth column, the eighth row and the seventh column, and the eighth row and the eighth column. With the above processing, the elements in the first to sixth columns in the seventh and eighth rows are zero-erased. Since the subsequent processing is the same as that of the first embodiment, description thereof is omitted.

<Operation of arithmetic unit>
FIG. 6 is a flowchart of processing in the second embodiment. Similar to the first embodiment, for example, it is started by inputting a pair of coordinate points between a temporally preceding frame and a subsequent frame. In addition, since it is the same as that of 1st Embodiment except three steps of step S602, S604, and S605, description is abbreviate | omitted.

In step S602, the CPU 120 decomposes the coefficient matrix A of 2N rows × 8 columns into a matrix Q _p and a matrix R _{p in the} format shown in FIG. As described above, Q _p is obtained by calculating and subtracting partial averages for the second and third columns to make the first column orthogonal, and then making the second and third columns orthogonal. can get. As described above, it is sufficient for the fourth, fifth, and sixth columns to use the results of the first, second, and third columns.

In step S604, the CPU 120 calculates the reciprocal of the diagonal elements of the upper left 6 × 6 matrix which is a diagonal matrix in the matrix (result of Q _p ^T Q _p ) shown in FIG. Although there are a total of six diagonal elements, the diagonal elements of 4-6 rows have the same values as those of 1-3 rows, so in practice only three reciprocals need be calculated.

  In step S605, the CPU 120 multiplies the corresponding lines of 1 to 6 by the three reciprocals calculated in step S604. Actually, it is not necessary to multiply all the elements of the row, and it is only necessary to multiply the three elements of the seventh and eighth columns and the solution vector.

  In the second embodiment, the column vectors of the seventh and eighth columns are not orthogonalized. However, as in the first embodiment, these two columns are assigned to the first and fourth columns. It is also possible to make it orthogonal. In this case, the processing after orthogonalization is almost the same as in the first embodiment, except that the number of elements (= the number of products) when calculating the product sum is reduced from six to four.

  As described above, according to the second embodiment, the least square solution such as a homography matrix or a three-dimensional projective transformation matrix can be obtained at high speed by orthogonalizing only a part of column vectors that can be orthogonalized easily. It becomes possible to estimate with high accuracy.

(Third embodiment)
In the third embodiment, a mode will be described in which a least square solution such as a homography matrix or a three-dimensional projective transformation matrix is estimated at high speed and with a smaller matrix size. In particular, a method for holding and processing coefficient data when the matrix size is reduced to half will be described.

  As shown in FIG. 1B, the coefficient matrix A that is an object for orthogonalizing column vectors has a size of 2N rows × 8 columns. That is, when the number of N is 500, for example, the storage area for storing the coefficient matrix is 8K words, which is considerably large. For this reason, there is a possibility that appropriate processing cannot be performed when the memory size cannot be made very large. In particular, in a processor for embedded use and the like, it is assumed that the memory area is divided and used for other processing, so that the smaller the matrix size, the better.

FIG. 7 is a diagram for explaining the transition of the coefficient matrix in the third embodiment. FIG. 7A shows a coefficient matrix A _h and a solution matrix B. The coefficient matrix A _h is a matrix in which the number of rows is N rows × 8 columns, which is a half of FIG. 1B, and is generated by storing coordinate data in the matrix. In the case of M dimensions, the coefficient matrix A _h has N rows ((M + 1) ² −1) columns. Further, the N + 1 and subsequent rows of the solution vector b shown in FIG. 1B are moved to the second column to obtain a solution matrix B of N rows × 2 columns. In the case of M dimensions, the solution matrix B has N rows and M columns. Further, the variable matrix X has ((M + 1) ² −1) rows and M columns. At this time, the simultaneous equations are expressed as A _h X = B.

Here, each column excluding the first column and the fourth column of the coefficient matrix A _h is orthogonalized with the first column. That is, the partial orthogonal matrix Q _p by subtracting and removing from the data to calculate the average of the N data in each column. Then, a corresponding upper triangular matrix R _p of 8 rows × 8 columns is generated.

FIG. 7B shows the coefficient matrix A _h as a product of the above-described Q _p and the upper triangular matrix R _p . As in FIG. 3, the result of subtracting the above average from the u coordinate and v coordinate is denoted as du and dv.

FIG. 7C shows an 8 × 8 matrix obtained by Q _p ^T Q _p . Here, 17 elements with a Σ symbol in FIG. 7C are calculated. The element value “N” in the upper left corner is a definite value and is not included in the calculation target.

  Among the 17 elements, 7 elements including y in the 5th and 6th columns are moved to the 7th and 8th columns. Alternatively, a column vector having the same configuration as the seventh and eighth columns of the matrix shown in FIG. 3B is formed by adding to the moved previous element. However, the bottom element in the seventh column is ignored. Thereafter, in order to symmetrize the matrix, the upper triangular area of the matrix is folded back with a diagonal component and copied to the lower triangular area.

FIG. 7D is an 8-row × 2-column matrix obtained by multiplying the solution matrix B by Q _p ^T. In FIG. 7D, 10 elements with a Σ symbol are calculated, and movement and addition are performed in the same manner as described above to obtain the same result as the vector shown in FIG.

  The subsequent processing is the same as the processing after S404 of the processing flow in the first embodiment. This is for processing the same matrix in FIG.

<Operation of arithmetic unit>
FIG. 8 is a flowchart of processing in the third embodiment. As described above, S404 and subsequent steps are the same as those in the first embodiment, and steps S801 to S805 before that will be described.

In step S801, the CPU 120 selects N sets (N is 5 or more) of coordinate data from the coordinate point pair group, and generates a coefficient matrix A _h of N rows × 8 columns in the format shown in FIG. To do.

In step S802, the CPU 120 decomposes the coefficient matrix A _h into a matrix Q _p and a matrix R _{p in the} format shown in FIG. 7B. As described above, Q _p is obtained by subtracting the average value for each column from each column except the first column and the fourth column of the coefficient matrix A _h .

In step S803, CPU 120 is a transposed matrix _Q ^{p T} of _{Q p,} multiplies the _{Q p.} As described above, the formed matrix is an 8 × 8 matrix, but at this stage, only 17 elements are calculated.

  In step S804, the CPU 120 transforms the 8 × 8 matrix being formed so as to be similar to the matrix shown in FIG. For example, as described above, the seven elements including y in the fifth and sixth columns are moved to the seventh and eighth columns, or added to the moved previous elements. In step S805, the CPU 120 folds the upper triangular area of the matrix with a diagonal component and copies it to the lower triangular area to make the matrix symmetric.

  As a result of the above processing, a symmetric matrix of 8 rows × 8 columns shown in FIG. 3B is obtained. Therefore, if the processing of S404 to S409 is performed on this processing result, the same result as in the first embodiment can be obtained.

  As described above, according to the third embodiment, it is possible to reduce the storage area by reducing the size of the matrix storing the coefficient matrix. Note that the matrix size can be further reduced by doing the following.

In the coefficient matrix _Ah , all the data stored in the fourth column is “0”, so it can be considered that there are seven columns except for the fourth column. Further, since all the data stored in the first column is “1”, it can be considered that there are six columns excluding the first column.

  In that case, the matrix size obtained by multiplication with the transposed matrix is 6 rows × 6 columns. In order to obtain the 8-row × 8-column matrix shown in FIG. 3B, the elements need to be rearranged significantly, but the storage area for storing the coefficient matrix can be further reduced by 25%.

(Modification)
In the first to third embodiments described above, the description has been made on the homography matrix that is a two-dimensional projective transformation matrix. However, the present invention can also be applied to least square estimation of an M-dimensional, for example, three-dimensional projective transformation matrix.

Since the degree of freedom of the homography matrix is 8 (= 3 × 3-1), the size of the coefficient matrix A for N coordinate point pairs is 2N rows × 8 columns. The calculation result of Q _p ^T Q _p is a matrix of 8 rows × 8 columns.

On the other hand, since the degree of freedom of the three-dimensional projective transformation matrix is 15 (= 4 × 4-1), the size of the coefficient matrix is 2N rows × 15 columns. That is, the calculation result of the corresponding Q _p ^T Q _p is a matrix of 15 rows × 15 columns. Therefore, when all the column vectors of the 15-column coefficient matrix are completely orthogonalized using QR decomposition, the calculation amount is nearly four times that in the case of 8 columns.

On the other hand, similarly to the above-described embodiment, when partial orthogonalization centering on removal of the average value is performed, processing can be performed in a short time. That is, the result of Q _p ^T Q _p is a 15 × 15 matrix, but non-zero off-diagonal components can be effectively zeroed out.

(Other examples)
The present invention supplies a program that realizes one or more functions of the above-described embodiments to a system or apparatus via a network or a storage medium, and one or more processors in a computer of the system or apparatus read and execute the program This process can be realized. It can also be realized by a circuit (for example, ASIC) that realizes one or more functions.

  101 arithmetic unit; 120 CPU; 121 RAM; 122 ROM

Claims (10)

An arithmetic device for deriving a transformation matrix of M-dimensional projective transformation by solving simultaneous equations formed based on a pair of coordinate points before and after transformation by the projective transformation,
2N (2N> (M + 1) ² −1) simultaneous equations formulated based on N coordinate point pairs, each of which includes a coefficient matrix A, a solution vector b, and elements of the transformation matrix Determining means for determining simultaneous equations expressed as Ax = b by the variable vector x
Decomposition means for decomposing the coefficient matrix A into a partial orthogonal matrix Q _p and an upper triangular matrix R _p partially orthogonalized between column vectors in the coefficient matrix A;
The variable vector x is determined by processing a modified simultaneous equation obtained by multiplying both sides of the simultaneous equations by the transposed matrix Q _p ^T of the partial orthogonal matrix Q _p , and each element of the transformation matrix of the projective transformation is defined as Deriving means for deriving;
Have
The arithmetic unit, wherein at least one set of column vectors in the partial orthogonal matrix Q _p is not orthogonal to each other. 2. The decomposition means performs orthogonalization by subtracting an average value of significant elements of each column vector from the significant elements in a predetermined column vector in the coefficient matrix A. The computing device described. The computing device according to claim 2, wherein the predetermined column vector includes 2M column vectors including values of coordinate points before conversion included in the N coordinate point pairs. The derivation means performs a predetermined calculation on both sides of the simultaneous equations so that a non-diagonal component of Q _p ^T Q _p existing on the left side of the modified simultaneous equations becomes zero. 4. The arithmetic unit according to any one of items 1 to 3. The derivation means determines the variable vector x by multiplying both sides of the modified simultaneous equations by an inverse matrix of Q _p ^T Q _p and further multiplying both sides by an inverse matrix of R _p. The arithmetic device according to claim 4. The arithmetic unit according to claim 5, wherein the multiplication is realized by performing a backward substitution process. An arithmetic device for deriving a transformation matrix of M-dimensional projective transformation by solving simultaneous equations formed based on a pair of coordinate points before and after transformation by the projective transformation,
2N (2N> (M + 1) ² −1) simultaneous equations formulated based on N coordinate point pairs, and N rows ((M + 1) ² −1) columns of coefficient matrix A _h , N rows and M columns of solution matrix B, and ((M + 1) ² −1) rows and M columns of variable matrix X configured by each element of the transformation matrix determine simultaneous equations expressed as A _h X = B. A determination means;
The coefficient matrix A _h, and the coefficient matrix A partial orthogonal matrix was orthogonalization between column vectors partially in _h Q _p and an upper triangular matrix R _p and the decomposing decomposition means,
The solution matrix B is determined by processing a modified simultaneous equation obtained by multiplying both sides of the simultaneous equations by the transposed matrix Q _p ^T of the partial orthogonal matrix Q _p , and each element of the transformation matrix of the projective transformation is _defined as Deriving means for deriving;
Have
The arithmetic unit, wherein at least one set of column vectors in the partial orthogonal matrix Q _p is not orthogonal to each other. An arithmetic method in an arithmetic device for deriving a transformation matrix of M-dimensional projective transformation by solving simultaneous equations formed based on coordinate point pairs before and after transformation by the projective transformation,
2N (2N> (M + 1) ² −1) simultaneous equations formulated based on N coordinate point pairs, each of which includes a coefficient matrix A, a solution vector b, and elements of the transformation matrix A determination step of determining simultaneous equations expressed as Ax = b by a variable vector x
A decomposition step of decomposing the coefficient matrix A into a partial orthogonal matrix Q _p and an upper triangular matrix R _p partially orthogonalized between column vectors in the coefficient matrix A;
The variable vector x is determined by processing a modified simultaneous equation obtained by multiplying both sides of the simultaneous equations by the transposed matrix Q _p ^T of the partial orthogonal matrix Q _p , and each element of the transformation matrix of the projective transformation is defined as A derivation process to derive,
Including
The calculation method according to claim 1, wherein at least one set of column vectors in the partial orthogonal matrix Q _p is not orthogonal to each other. An arithmetic method in an arithmetic device for deriving a transformation matrix of M-dimensional projective transformation by solving simultaneous equations formed based on coordinate point pairs before and after transformation by the projective transformation,
2N (2N> (M + 1) ² −1) simultaneous equations formulated based on N coordinate point pairs, and N rows ((M + 1) ² −1) columns of coefficient matrix A _h , N rows and M columns of solution matrix B, and ((M + 1) ² −1) rows and M columns of variable matrix X configured by each element of the transformation matrix determine simultaneous equations expressed as A _h X = B. A decision process;
The coefficient matrix A _h, and the coefficient matrix A partial orthogonal matrix was orthogonalization between column vectors partially in _h Q _p and an upper triangular matrix R _p and the decomposing decomposition step,
The solution matrix B is determined by processing a modified simultaneous equation obtained by multiplying both sides of the simultaneous equations by the transposed matrix Q _p ^T of the partial orthogonal matrix Q _p , and each element of the transformation matrix of the projective transformation is _defined as A derivation process to derive,
Have
The arithmetic unit, wherein at least one set of column vectors in the partial orthogonal matrix Q _p is not orthogonal to each other.   The program for functioning a computer as each means of the arithmetic unit of any one of Claims 1 thru | or 7.

Images