JP2007257518A - Image matching device, image matching method and image matching program - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To provide an image matching device, an image matching method and an image matching program, capable of performing higher-accuracy image matching. <P>SOLUTION: This image matching device has: a movement clustering processing part 109 allocating a label to each area obtained by dividing a reference image; a movement parameter estimation part 110 estimating a first movement parameter of a movement model; a movement parameter correction part 113 correcting the first movement parameter such that potential energy between a point on the reference image and a point on a target image decided by the first movement parameter becomes minimum; a movement parameter allocation part 111 determining a second movement parameter from the first movement parameter, and allocating it to each second lattice point inside the area corresponding to the label; a correction item calculation part 112 calculating a correction item from the second movement parameter and a solution processing result of a first motion equation; and a second solution processing part 108 performing solution processing of a second motion equation wherein the correction item is added to the first motion equation. <P>COPYRIGHT: (C)2008,JPO&INPIT

Description

  The present invention relates to an image matching apparatus, an image matching method, and an image matching program that detect and associate corresponding points from two images of a target image and a reference image.

  Image matching technology that detects the corresponding point of each pixel of the other reference image from one target image and obtains the corresponding relationship includes, for example, motion detection in moving images, stereo matching, image morphing, image recognition, and moving image coding Are used in the field of image processing technology.

  Such image matching techniques can be classified into four main types: an optical flow method, a block-based method, a gradient method, and a Bayesian method (see, for example, Non-Patent Document 1).

  The optical flow method derives an optical flow equation that “a change in luminance is constant” and obtains a flow using the optical flow equation as a constraint. The block-based method is a method in which an image is divided into predetermined blocks and a motion is obtained by template matching for each block. The gradient method is a method for performing matching in a direction in which the luminance gradient of an image decreases. The Bayesian method is a technique for obtaining matching that is probabilistically plausible.

  In addition, as a conventional image matching technique, a plurality of multi-resolution image pyramids are generated using a plurality of multi-resolution filters, and the generated image pyramids are sequentially matched from the upper layer to the lower layer, thereby performing a large process. There is an image matching technique with high robustness capable of associating images from movement to small movement (see, for example, Patent Document 1).

A. MuratTekalp, “Digital Video Processing”, Prentice Hall, 1995 Japanese Patent No. 2927350

  However, such a conventional image matching technique has the following problems. The optical flow method has a problem that it is sensitive to noise and it is inherently difficult to cope with fast movement. In the block-based method, image matching is performed for each block in the image. However, there is a problem that it is difficult to essentially cope with the movement of the object in the image such as deformation or rotation. In the gradient method, since matching is performed in a direction in which the luminance gradient of the image decreases, there is a problem that it is difficult to stably search for the movement of the object. The Bayesian method has a problem that it is difficult to obtain a global optimum.

  On the other hand, in the technique disclosed in Patent Document 1, in principle, a plurality of multi-resolution filters are used to perform matching processing from the highest layer to the lowest layer of the multi-resolution image pyramid. It is a typical optimization process and is not necessarily optimized globally.

  The present invention has been made in view of the above, and an object thereof is to provide an image matching apparatus, an image matching method, and an image matching program capable of performing image matching with higher accuracy.

  In order to solve the above-described problems and achieve the object, the present invention is an image matching device for obtaining a correspondence relationship between a target image and a reference image, and a plurality of first lattice points on the target image. And a potential energy based on the correlation of the image between each of the second grid points corresponding to the first grid point on the reference image on a one-to-one basis, The image energy for calculating the image energy force received by the second lattice point by the gradient of the potential energy based on the position of the lattice point and the position of the first lattice point corresponding to the second lattice point Elastic energy force calculation for calculating a first elastic energy force received from an elastic energy between the force calculation unit and each second lattice point and another second lattice point adjacent to the second lattice point And A frictional force calculation unit for calculating a frictional force acting on each second lattice point, and a first relating to each second lattice point by the image energy force, the first elastic energy force, and the frictional force. A first solution processing unit that solves the equation of motion by numerical analysis, and a result of the solution processing by the first solution processing unit, the reference image is divided into a plurality of regions, and each region is divided into regions. A clustering processing unit for assigning a label for identifying the first motion, and a first motion which is a parameter for defining a motion model indicating a mapping relation between the second lattice point and the first lattice point in the region to which the label is assigned Between a point on the target image and a point on the reference image determined by the estimated first motion parameter and a motion parameter estimation unit that estimates a parameter by numerical analysis A motion parameter correcting unit that corrects the first motion parameter so that the potential energy is minimized, and the second lattice point in the region corresponding to the label from the corrected first motion parameter. A second motion parameter that is the first motion parameter that is optimal for the second motion parameter is determined, and the determined second motion parameter is assigned to each second grid point in the region corresponding to the label. Correction term calculation for calculating a correction term for correcting the first equation of motion based on a motion parameter assignment unit, the assigned second motion parameter, and a result of the solution processing by the first solution processing unit And second solution processing for solving the second motion equation relating to each of the second lattice points obtained by adding the correction term to the first motion equation by numerical analysis And a mapping processing unit for obtaining a correspondence relationship between the target image and the reference image from a result of the solution processing by the second solution processing unit.

  The present invention also provides an image matching method and an image matching program that can be executed by the image matching apparatus.

  According to the present invention, the reference image is divided into a plurality of regions, and the first motion parameter is the optimal motion parameter for classifying the motion model that defines the motion of the second grid point for each region. 2 motion parameters are obtained and solution processing is performed by numerical analysis of the equation of motion indicating the motion model to obtain the correspondence between the target image and the reference image. Uniformity can be expressed robustly, more accurate image matching can be performed, and as a result, an interpolated frame with higher image quality can be generated.

  Exemplary embodiments of an image matching apparatus, an image matching method, and an image matching program according to the present invention are explained in detail below with reference to the accompanying drawings.

The image matching apparatus according to the following embodiment is applied to a motion compensation apparatus that generates a compensated image by performing motion compensation of a moving image, stereo matching, image morphing,
The present invention can also be applied to image processing techniques such as image recognition and moving image encoding.

(Embodiment 1)
FIG. 1 is a block diagram illustrating a configuration of the image matching apparatus according to the first embodiment. As shown in FIG. 1, the image matching apparatus 100 according to the present embodiment includes a first solution processing unit 101, an image energy force calculation unit 103, an elastic energy force calculation unit 104, a friction force calculation unit 105, Second solution processing unit 108, motion clustering processing unit 109, motion parameter estimation unit 110, motion parameter correction unit 113, motion parameter allocation unit 111, correction term calculation unit 112, mapping processing unit 107, frame A memory 106 is mainly provided.

  In the image matching device 100 according to the first embodiment, each lattice point obtained by dividing each of the target image and the reference image into a predetermined lattice shape is treated as a dynamic mass point, and the first lattice point and the first lattice point of the target image are processed. Considering the potential energy based on the correlation between the second grid points of the reference image corresponding to the grid points one-to-one, the gradient of the potential energy based on the positions of the first and second grid points By the image energy force received by the second grid point, the first elastic energy force received from the elastic energy between the grid points adjacent to the second grid point, and the frictional force generated at the second grid point Using the first equation of motion related to the second lattice point as an image matching state transition model, a solution process is performed on this equation of motion, the equilibrium state of the second lattice point is obtained, and the local motion is determined. Mel.

  Furthermore, the image matching apparatus 100 according to the present embodiment introduces a parametric motion model that defines a mapping relationship by a plurality of parameters (motion parameters). That is, assuming that the movement of all grid points in the reference image screen (mapping relationship with the first grid point) is determined by a finite number of parametric motion models, the reference image screen is divided into a plurality of regions. Clustering (dividing) is performed, a motion parameter that is a parameter defining a motion model is estimated, and an optimal motion parameter is assigned to a second lattice point in the region. Then, the second motion equation introduced into the first motion equation as a correction term is used as a solution to correct the deviation between the local motion that is the solution processing result of the first motion equation and the motion based on the motion model. The motion is updated locally by processing. By repeating such motion model motion parameter estimation processing, motion parameter assignment processing, and local motion update (second motion equation solving processing), the accurate motion of the image is estimated, Image matching is performed by obtaining a mapping relationship between the first grid point and the second grid point.

  The first solution processing unit 101 between each of the plurality of first grid points on the target image and each of the plurality of second grid points corresponding to the first grid points on the reference image on a one-to-one basis. An image energy force received by the second lattice point due to a gradient of potential energy based on the correlation of the images, a first elastic energy force received from elastic energy between lattice points adjacent to the second lattice point, This is a processing unit that solves the equation of motion related to the second lattice point by the frictional force generated at the two lattice points by numerical analysis processing using the Euler method as a discrete variable method.

  The image energy force calculation unit 103 is a processing unit that calculates the image energy force received by the second lattice point based on the gradient of the potential energy based on the positions of the first lattice point and the second lattice point.

  The elastic energy force calculation unit 104 is a processing unit that calculates a first elastic energy force received from elastic energy between lattice points adjacent to the second lattice point. The frictional force calculation unit 105 is a processing unit that calculates the frictional force generated at the second grid point. Further, the elastic energy force calculation unit 104 performs the second grid point after the movement of the second grid point determined by the optimal motion parameter and the second grid point in the region where the label described later is the same. The elastic energy force is calculated, and the second elastic energy force is calculated so that the second lattice energy point and the third lattice point in regions having different labels are smaller than the second elastic energy force in the same region. Calculate In the present embodiment, the second lattice point and the third lattice point in regions having different labels are subjected to a process of setting the friction coefficient to 0 and the second elastic energy force to 0.

  The motion clustering processing unit 109 clusters (divides) the reference image into a plurality of regions by the clustering method from the local motion obtained by the solution processing by the first solution processing unit 101, and identifies the region for each divided region. A processing unit for assigning labels to be performed.

  The motion parameter estimation unit 110 estimates a motion parameter, which is a parameter that defines a motion model indicating the mapping relationship between the second grid point and the first grid point, in a region to which a label is assigned by numerical analysis. In the present embodiment, the motion parameter (first motion parameter) is estimated by performing numerical analysis by the least square method.

  The motion parameter correction unit 113 corrects the motion parameter by the steepest descent method so that the potential energy based on the image correlation between the point on the target image and the point on the reference image determined by the estimated motion parameter is minimized. To do.

  The motion parameter assignment unit 111 determines an optimum label for each second grid point in the region corresponding to the label from the estimated motion parameter, and obtains an optimum motion parameter (second motion parameter). The processing unit assigns the optimal motion parameter to each second lattice point in the region corresponding to the label.

  The correction term calculation unit 112 calculates a correction term for correcting the first equation of motion from the difference between the motion vector based on the allocated optimal motion parameter and the motion vector that is the result of the solution processing by the first solution processing unit. It is a processing unit.

  The first solution processing unit 108 receives the second energy energy from the image energy force and the elastic energy between the second lattice point and the third lattice point after the movement of the second lattice point determined by the optimum motion parameter. The processing unit solves the second equation of motion that defines the movement of the second lattice point by the numerical analysis using the elastic energy force, the frictional force, and the correction term obtained by the correction term calculation unit 112.

  The mapping processing unit 107 obtains the equilibrium state of the second lattice point from the result of the solution processing by the second solution processing unit 108, thereby obtaining the first lattice point on the target image and the second lattice point on the reference image. It is a processing part which calculates | requires the mapping relationship which is a corresponding relationship.

  The frame memory 106 is a memory for storing the input target image and reference image. Detailed processing of each unit will be described later.

  Next, an image matching state transition model introduced by the image matching apparatus 100 according to the present embodiment will be described.

The pixel values of the target image and the reference image at the point x∈E ² on the two-dimensional Euclidean space are shown by the formula (1).

（１）式は実数空間における連続画像モデルであり、本実施の形態ではデジタル画像を扱うため、２次元格子空間上の点ｎ∈Ｎ²について（２）式に示す離散画像モデルを考える。

Equation (1) is a continuous image model in real number space. In this embodiment, in order to handle digital images, a discrete image model shown in Equation (2) is considered for a point n∈N ² on a two-dimensional lattice space.

従って、連続画像モデルの画素値はデータ補間関数ｇを用いて、（３）式に示すように定義される。ここで、データ補間関数ｇとしては、例えば双曲形補間関数を使用することができる。

Therefore, the pixel value of the continuous image model is defined as shown in the equation (3) using the data interpolation function g. Here, as the data interpolation function g, for example, a hyperbolic interpolation function can be used.

本実施の形態では、対象画像から参照画像への画像マッチング問題を、対象画像上の曲面Ｘ⊂Ｅ²から参照画像上の曲面Ｙ⊂Ｅ²への写像として求めることにする。図２は、対象画像上の曲面から参照画像上の曲面への写像を示す説明図である。図３は、対象画像上の曲面Ｘの点ｘから参照画像上の曲面Ｙの点ｙへの写像を示す説明図である。すなわち、図２および３に示すように、画像マッチングは、対象画像上の曲面Ｘの点が参照画像上の曲面Ｙのどの点に対応しているのかを表す写像を求める問題として扱う。

In this embodiment, to obtaining an image matching problem from the target image to the reference image, as a mapping from the curved X⊂E ² on the target image to the curved surface Y⊂E ² on the reference image. FIG. 2 is an explanatory diagram showing mapping from the curved surface on the target image to the curved surface on the reference image. FIG. 3 is an explanatory diagram showing a mapping from a point x on the curved surface X on the target image to a point y on the curved surface Y on the reference image. That is, as shown in FIGS. 2 and 3, image matching is handled as a problem of obtaining a mapping that indicates which point of the curved surface Y on the reference image corresponds to the point of the curved surface X on the target image.

  Deriving such a mapping condition is as follows. Using the parameter (u, v), the point on the curved surface X is represented as x (u, v) εX, the point on the curved surface Y is represented as y (u, v) εY, and the parameter (u, v) and the point are represented. When x has the same scale, the point x (u, v) on the curved surface X can be expressed by the equation (4) by the parameter (u, v).

ここで、曲面Ｘから曲面Ｙへの求めるべき写像を（５）式で定義する。写像の対応関係は曲面Ｘ上の点ｘが一意に曲面Ｙに対応づけることができればよいので、写像ｇはすべての点が一意に決定する全単射とする。よって同一パラメータ（ｕ，ｖ）を有する点ｘ（ｕ，ｖ）、ｙ（ｕ，ｖ）は，（６）式を満たすことになる。

Here, the mapping to be obtained from the curved surface X to the curved surface Y is defined by equation (5). Since the mapping relationship is sufficient if the point x on the curved surface X can be uniquely associated with the curved surface Y, the mapping g is bijection in which all points are uniquely determined. Therefore, the points x (u, v) and y (u, v) having the same parameter (u, v) satisfy the expression (6).

このことは、すなわち、写像ｇについて考えることは、曲面Ｙそのものを考えることと同じである。一つの対象物が対象画像と参照画像に射影されて映っていると考えると、射影された対象画像上の像と参照画像上の像は、画素値が同じであると考えることができる。従って、対応づけるべき同一点の画素値が等しいと仮定すると、写像ｇに関して（７）式の条件式が導出される。かかる（７）式をパラメータｕ＝（ｕ，ｖ）を用いて表示すると、（８）式のようになる。

This means that thinking about the mapping g is the same as considering the curved surface Y itself. If one object is projected on the target image and the reference image, it can be considered that the projected image on the target image and the image on the reference image have the same pixel value. Therefore, assuming that the pixel values of the same point to be matched are equal, the conditional expression (7) is derived for the mapping g. When the equation (7) is displayed using the parameter u = (u, v), the equation (8) is obtained.

従って、画像マッチング問題は、（８）式の条件式を満たす写像ｇを求めるという問題に帰着することになる。

Therefore, the image matching problem results in a problem of obtaining a mapping g that satisfies the conditional expression (8).

  Next, a solution process for the image matching problem will be described. Expression (8) can be modified like Expression (9).

しかしながら、（９）式は、画像に含まれるノイズなどにより厳密には成立しないため、かかる式から直接写像ｇを求めることは困難である。このため、本実施の形態では、最小二乗法を利用している。最小二乗法によれば、（１０）式に示すエネルギー関数を考え、かかるエネルギー関数の最小化を行うことによりノイズが含まれていても条件を満足する写像ｇを求めることができる。

However, since the expression (9) is not strictly established due to noise included in the image, it is difficult to directly obtain the mapping g from such an expression. For this reason, the least square method is used in the present embodiment. According to the least square method, the energy function shown in the equation (10) is considered, and by minimizing the energy function, a mapping g that satisfies the condition can be obtained even if noise is included.

大数の法則が成立する程度の標本数があれば、最小二乗法により解が安定に求まると考えられるが、本実施の形態では、標本数が一つであるため、やはりノイズの影響を受けてしまう。解を安定して求めるためには、ベイズの定理により、近傍では滑らかに変化するという事前知識により事後分布の精度を高めるアプローチ、例えば、写像ｇに対して構造を規定して力学的なアプローチを考える。

If the number of samples is sufficient to satisfy the law of large numbers, the solution can be obtained stably by the least squares method. However, in this embodiment, since there is only one sample, it is still affected by noise. End up. In order to obtain a stable solution, an approach that improves the accuracy of the posterior distribution based on Bayes' theorem based on prior knowledge that it changes smoothly in the vicinity, for example, a dynamic approach by defining the structure for the mapping g Think.

  Since the mapping g is bijective, thinking about the mapping g is synonymous with considering the curved surface Y itself. As the mechanical structure, for example, a structure like a rigid body is conceivable, but in the case of an image, deformation or the like is conceivable, so that the structure of the rigid body is not appropriate. Therefore, in the present embodiment, an elastic body structure that can flexibly cope with image deformation or the like is introduced into the curved surface Y.

  The energy of the elastic body is expressed by equation (11).

また、（１０）式で示される画像エネルギーを曲面Ｙ上の点を用いて表すと、（１２）式で表される。

Further, when the image energy represented by the equation (10) is expressed using points on the curved surface Y, it is represented by the equation (12).

ポテンシャルエネルギーは、（１２）式の画像エネルギーと（１１）式の弾性エネルギーの線形結合により（１３）式のように定義される。

The potential energy is defined as in equation (13) by linear combination of image energy in equation (12) and elastic energy in equation (11).

このポテンシャルエネルギーを最小にすることは、画像マッチングの条件式である（７）式をできるだけ満たすように、かつ弾性体としての安定形態が取れるような曲面Ｙを探索することを意味する。このことをエネルギー最小化問題として定式化すると（１４）式で表される。

Minimizing this potential energy means searching for a curved surface Y that satisfies Equation (7), which is a conditional expression for image matching, as much as possible and that can take a stable form as an elastic body. When this is formulated as an energy minimization problem, it is expressed by equation (14).

従って、条件式（７）式を満たす写像ｇを求める画像マッチング問題は、エネルギー最小化問題の（１４）式を解法することに帰着する。

Therefore, the image matching problem for obtaining the mapping g that satisfies the conditional expression (7) results in solving the expression (14) of the energy minimization problem.

  Further, the necessary condition for optimality can be obtained from the equation (15) from the equation (14).

従来の手法では、この最適性必要条件の（１５）式を満たすような曲面Ｙ上の点ｙを求めることにより解法している。これらは静的なエネルギーの最小化であり、静的な釣り合いの位置を求める静的な最適化である。図４は、ポテンシャルエネルギーの最小化の点を探索する場合の静的な最適点の説明をするための模式図である。図４に示すように、静的な最適化では初期値に最も近い位置で質点のポテンシャルエネルギーの局所最適点（釣り合い位置）が求まる。従って、この局所最適点が大域的な最適点になるかどうかは初期値に依存することになる。言い換えれば、初期値の選定が好ましくない場合には、大域的に最適な点が得られない場合がある。

In the conventional method, the solution is obtained by obtaining a point y on the curved surface Y that satisfies the equation (15) of the optimality requirement. These are static energy minimizations and static optimizations that determine the position of the static balance. FIG. 4 is a schematic diagram for explaining a static optimum point when searching for a potential energy minimizing point. As shown in FIG. 4, in the static optimization, the local optimum point (balance position) of the potential energy of the mass point is obtained at the position closest to the initial value. Therefore, whether or not this local optimum point becomes a global optimum point depends on the initial value. In other words, when the selection of the initial value is not preferable, a globally optimum point may not be obtained.

  Therefore, in order to find a global optimum point and solve the equation (14), in this embodiment, each lattice point obtained by dividing the target image and the curved surface on the reference image into a lattice shape is considered as a mass point. The kinetic energy at the lattice point (mass point) of the curved surface Y on the reference image is introduced to perform global optimization. In other words, a time axis is introduced to the curved surface Y so that the curved surface Y can be deformed according to time, and a structure of motion is applied to the curved surface. FIG. 5 is a schematic diagram for explaining a global optimum point when searching for a potential energy minimizing point. As shown in FIG. 5, since the initial potential energy is converted into kinetic energy by the introduction of kinetic energy, the global optimality satisfying Eq. (15) in a wider range from the local optimal point to the lattice point (mass point). You can search for points. As a result, a wider range can be searched than in the static optimization, so that the dependency on the initial value is lowered, and there is a possibility that it is more robust against the influence of noise and the like.

Next, introduction of an equation of motion for global optimization of equation (15) will be described. First, the point y on the curved surface Y is expanded as shown in Equation (16) using the time τ introduced into the equation of motion. The total differentiation of the point y on the curved surface Y expressed by the equation (16) is the equation (17), but assuming that the parameter (u, v) ^T is independent of the time τ, the equation (18) is Since this holds, equation (17) is eventually expressed by equation (19). Further, the second derivative of the point y on the curved surface Y from the equation (19) is expressed by the equation (20). The total differentiation of the point y on the curved surface Y shown by the equation (19) and the second order differentiation shown by the equation (20) are expressed by the equation (21), respectively.

（１９）式は、曲面Ｙ上の点ｙの速度を表し、（２０）式は、曲面Ｙ上の点ｙの加速度を表している。従って、曲面Ｙ上の点ｙの運動エネルギーは、（２２）式で定義することができる。

Equation (19) represents the speed of the point y on the curved surface Y, and Equation (20) represents the acceleration of the point y on the curved surface Y. Therefore, the kinetic energy of the point y on the curved surface Y can be defined by equation (22).

このような運動エネルギーを導入したことは、曲面Ｙ上の点ｙを質点と考えて、曲面上の運動を定義したことになり、曲面Ｙ上の点ｙの運動という力学的な構造を規定したことを意味する。かかる曲面Ｙ上の点ｙの運動は運動方程式によって記述される。

The introduction of such kinetic energy means that the point y on the curved surface Y is considered as a mass point and the motion on the curved surface is defined, and the dynamic structure of the movement of the point y on the curved surface Y is defined. Means that. The motion of the point y on the curved surface Y is described by a motion equation.

  Next, the equation of motion of the point y on the curved surface Y is derived as follows according to the Lagrange method. Since the curved surface Y follows the equation (14) for energy minimization, the curved surface Y needs to move in a direction that minimizes the potential energy U. Lagrangian is defined by equation (23).

  Equation (23) means that the object falls so as to minimize the gravitational field, similar to the free fall of the object in the gravitational field that is the potential energy of the position. Since the Lagrangian equation is expressed by the equation (24), the equation of motion of the point y on the curved surface Y is expressed by the equation (25) from the equations (23) and (24).

（２５）式の運動方程式に従った運動は保存系であり、エネルギー保存則により解は収束しない。このため、系を非保存系にするために（２６）式で示される摩擦エネルギーを導入する。

The motion according to the equation of motion of Equation (25) is a conservation system, and the solution does not converge due to the energy conservation law. For this reason, in order to make the system non-conservative, the frictional energy represented by the equation (26) is introduced.

（２６）式で示される摩擦エネルギーを導入した場合のラグランジュ方程式は、（２７）式で表される。

The Lagrangian equation when the friction energy represented by the equation (26) is introduced is represented by the equation (27).

この（２７）式を解法することによって、摩擦エネルギーを導入した運動方程式が（２８）式で得られる。

By solving the equation (27), an equation of motion in which friction energy is introduced is obtained by the equation (28).

次に、動的な探索、（２８）式の運動方程式に従った解がエネルギー最小化の（１４）式に対する最適性必要条件である（１５）式を満足することについて説明する。

Next, it will be described that the dynamic search and the solution according to the equation of motion of the equation (28) satisfy the equation (15) which is the optimum requirement for the energy minimization (14).

  First, consider the case where the search is completed. In this case, since the point y stops completely, both the velocity (total derivative of y) and acceleration (second derivative of y) of the point y are both 0, so that the equation (29) is obtained from the equation of motion (28). Thus, equation (30) is established.

（３０）式が成立することは、最適性必要条件の（１５）式を満足することを意味している。従って、探索が完了した場合、そのときの点ｙは最適点であることがわかる。

The fact that the expression (30) is satisfied means that the expression (15) of the optimality requirement is satisfied. Therefore, when the search is completed, it can be seen that the point y at that time is the optimum point.

  Next, consider a case where the point y on the curved surface Y is stationary. In this case, the velocity at the point y (the total derivative of y) is 0, but the acceleration (the second derivative of y) is not 0, so that equation (31) is obtained from the equation of motion of equation (28).

ここで、点ｙが最適点であると仮定すると、最適性条件の（１５）式により、点ｙの加速度（ｙの２階微分）は０となり、最適点の探索が完了する。一方、点ｙが最適点でない場合には、点ｙの加速度（ｙの２階微分）が０とならないため、最適点の探索は完了していないことになる。すなわち、点ｙが静止した位置が最適点の位置である場合に（３１）式により最適性必要条件の（１５）式が満足することを意味している。

Assuming that the point y is the optimum point, the acceleration of the point y (second derivative of y) becomes 0 according to the equation (15) of the optimum condition, and the search for the optimum point is completed. On the other hand, when the point y is not the optimum point, the acceleration of the point y (second-order derivative of y) does not become 0, and the optimum point search is not completed. That is, when the position where the point y is stationary is the position of the optimum point, it means that the expression (15) of the optimality requirement is satisfied by the expression (31).

  Next, consider the case where frictional energy is acting on the point y. Equation (32) is derived from the nature of the frictional energy, which means that the movement of the point y stops and the search for the optimum point is completed over time.

このように、最適点の探索が完了した場合、点ｙが静止した場合、点ｙに摩擦エネルギーが作用している場合の上記説明により、（２８）式の運動方程式に従った点ｙの運動が、（１５）式の最適性必要条件を満足することがわかる。

As described above, when the search for the optimum point is completed, when the point y is stationary, or when the friction energy is acting on the point y, the motion of the point y according to the equation of motion of the equation (28) However, it is understood that the optimality requirement of the equation (15) is satisfied.

  On the other hand, equation (15) is in the form of first derivative, but the optimality requirement for second derivative is satisfied because the search may be completed at the peak of the potential energy peak shown in FIG. May not be. In this case, the optimality requirement for the second floor can be satisfied by assuming the equation (33) as a condition for the search for the optimum point not to be completed at the apex of the peak of potential energy. Specifically, the assumption of equation (33) can be achieved by realizing the condition of equation (34) with a search algorithm.

  From the above, it can be seen that the search for the optimum point according to the equation of motion of Equation (28) gives the optimum solution for the optimization problem shown by Equation (14). Such a search is more robust than a static search.

Here, in the equation of motion shown in the equation (28), since the partial differential with respect to the point y and the parameters (u, v) is included in the elastic energy E _k , the partial differential equation is converted into the differential equation by the finite difference method. Replace. When the right side of the elastic energy E _k is discretized by the finite difference method, the equation (35) is obtained.

ここで、ｎは離散化されたパラメータ（ｕ，ｖ）の格子点、ｙ_nはｎに対応する離散化されたＹ上の点、ｙ_n+(-1,0)，ｙ_n+(1,0)，ｙ_n+(0,-1)，ｙ_n+(0,1)は点ｙ_nの偏微分に関連する上下左右の格子点を示している。

Here, n is a lattice point of the discretized parameter (u, v), y _n is a point on Y that is discretized corresponding to _n , y _{n + (− 1,0)} , y _{n + (1,0 )} , Y _{n + (0, -1)} , y _{n + (0,1)} indicate the upper, lower, left, and right lattice points related to the partial differentiation of the point y _n .

The equation (35) when partially differentiated by y _n (36) below is obtained and such expression shows elastic energy force is the force generated by the elastic energy.

この弾性エネルギー力を（３７）式で示し、離散パラメータΔｕ，Δｖをそれぞれ１とすると、弾性エネルギー力は（３８）式で表されることになる。

When this elastic energy force is expressed by the equation (37) and the discrete parameters Δu and Δv are respectively 1, the elastic energy force is expressed by the equation (38).

ここで、画像エネルギーＥ_iは、パラメータ（ｕ，ｖ）に関して独立と考え、（２８）式の運動方程式もそのまま離散可能となり、（３９）式で表される。

Here, the image energy E _i is considered to be independent with respect to the parameters (u, v), and the equation of motion of the equation (28) can be discrete as it is, and is expressed by the equation (39).

この（３９）式の右辺第１項は、点ｙ_nにおける画像エネルギーの勾配の逆向きの力である画像エネルギー力であり、（４０）式のようにおく。

The first term on the right side of the equation (39) is an image energy force that is the force opposite to the gradient of the image energy at the point y _n , and is set as in the equation (40).

次に、この画像エネルギー力Ｆ_i（ｎ）、すなわち画像エネルギーの勾配の計算について説明する。本実施の形態では、画像エネルギーの勾配を、局所最適化の手法によって求めている。図６および図７は、画像エネルギー力の概念を示す説明図である。ここで、画像モデルＳ_cは連続画像モデルであり、実際にはサンプリングされた画像モデルＳ_pを利用しているため、局所最適化においてもサンプリングされた画像モデルに基づいておこなっている。

Next, calculation of the image energy force F _i (n), that is, the gradient of the image energy will be described. In this embodiment, the gradient of image energy is obtained by a local optimization technique. 6 and 7 are explanatory diagrams showing the concept of image energy force. Here, the image model S _c is a continuous image model, in practice because it uses the image model S _p sampled, is performed also based on the sampled image model in local optimization.

The sampling point closest to y _n (τ) is set as the local space center y _c , and y _c is obtained by equation (41).

そして、隣接空間Ｌを（４２）式で定義する。このとき、局所探索集合Ωは、（４３）式のように定義することができる。

And the adjacent space L is defined by (42) Formula. At this time, the local search set Ω can be defined as in equation (43).

そして、局所最適化の後、その方向へのベクトルを求め、求めたベクトルを正規化し勾配の大きさを乗算すると、画像エネルギー力として（４４）〜（４６）式が得られる。本実施の形態では、画像エネルギー力として（４４）〜（４６）式を使用している。

Then, after local optimization, a vector in that direction is obtained, and when the obtained vector is normalized and multiplied by the magnitude of the gradient, equations (44) to (46) are obtained as image energy forces. In the present embodiment, the equations (44) to (46) are used as the image energy force.

なお、実装上の扱いやすさ等を考慮してエネルギー関数を（４７）式のように定義して、（４７）式により導かれる画像エネルギー力を用いることもできる。

It is also possible to define an energy function as shown in Equation (47) in consideration of ease of handling in mounting and the like, and use the image energy force derived from Equation (47).

以上説明した画像エネルギーＦ_i（ｎ）、弾性エネルギー力Ｆ_k（ｎ）、摩擦力を使用して、格子点ｙ_n（τ）に対する運動方程式を画像マッチング状態遷移モデルとして生成すると、運動方程式は、（４８）式で表される。また、ｙ_n（０）が参照画像に等しいとし、運動エネルギーは０、すなわちｙ_n（０）の速度が０として考えると、運動方程式の初期条件は、（４９）式となる。

Using the image energy F _i (n), elastic energy force F _k (n), and frictional force described above to generate the equation of motion for the lattice point y _n (τ) as an image matching state transition model, the equation of motion is , (48). Further, _assuming that y _n (0) is equal to the reference image and the kinetic energy is 0, that is, the velocity of y _n (0) is 0, the initial condition of the equation of motion is the equation (49).

本実施の形態では、画像エネルギー力計算部１０３によって、（４４）〜（４６）式で示される画像エネルギー力Ｆ_iを計算し、弾性エネルギー力計算部１０４によって（３８）式で示される弾性エネルギー力Ｆ_kを計算し、摩擦力計算部１０５によって（４８）式の右辺第３項で示される摩擦力を計算する。そして、第１解法処理部１０１によって（４８）、（４９）式で示される運動方程式を数値解析処理によって解法している。

In the present embodiment, the image energy force calculation unit 103 calculates the image energy force F _i expressed by the equations (44) to (46), and the elastic energy force calculation unit 104 calculates the elastic energy expressed by the equation (38). The force _Fk is calculated, and the frictional force calculation unit 105 calculates the frictional force indicated by the third term on the right side of the equation (48). The first equation processing unit 101 solves the equations of motion represented by equations (48) and (49) by numerical analysis processing.

  Next, the solution of the equation of motion by the first solution processing unit 101 will be described. Since the equation of motion represented by the equations (48) and (49), that is, the ordinary differential equation, cannot generally be solved analytically, the time T sufficient for the system of equations of motion to converge is considered, The convergence state of the equation of motion is estimated by calculating the t = (0, T) interval by numerical analysis.

  In the present embodiment, numerical analysis processing by the Euler method is performed using the fact that the ordinary differential equation is uniquely determined by the discrete variable method when the initial value is determined.

  In the present embodiment, numerical analysis processing by the Euler method is performed, but the present invention is not limited to this. In addition to the Euler method, the discrete variable method includes various methods such as the Runge-Kutta method, the Brillesh-Store method, the predictor / corrector method, and the hidden Runge-Kutta method. A technique may be used. In the following, a numerical analysis process using the Euler method will be described as an example.

  Since the Euler method is a numerical solution to the first-order ordinary differential equation, the equations (48) and (49) are converted by applying the variable transformation of the equation (50) to the equations (48) and (49). Convert to first-order ordinary differential equations. Thereby, formulas (51) and (52) are obtained.

一方、常微分方程式（５３）式をオイラー法で数値解析処理を行って解法する場合には、（５４）式のようにｔ^(m)からｔ^(m+1)=ｔ^(m)+ｈへ解を進展させる。

On the other hand, when solving the ordinary differential equation (53) by performing numerical analysis processing using the Euler method, t ^(m) to t ^{(m + 1)} = t ^(m) + h as shown in equation (54). Develop the solution.

ここで、ｘ^(m)はｍステップであることを示し、ｈはステップ幅である。オイラー法のスキームを（５１）、（５２）式に適用すると（５５）、（５６）式の更新式が得られる。

Here, x ^(m) indicates m steps, and h is a step width. When the Euler scheme is applied to the equations (51) and (52), the updated equations (55) and (56) are obtained.

本実施の形態では、第１解法処理部１０１によって、（５１）、（５２）式の更新式を繰り返して解法することによって第１の運動方程式（４８）、（４９）式の系の収束状態を求めている。そして、マッピング処理部１０７では、この収束状態から、曲面Ｙをｙ_n＝ｙ_n ^(T)のように求め、さらに、（６）式よりｙ_n＝ｇ（ｎ）の関係から、写像ｇを求めている。

In the present embodiment, the first solution processing unit 101 repeatedly solves the update equations (51) and (52), thereby converging the first equations of motion (48) and (49). Seeking. Then, the mapping processing unit 107 obtains the curved surface Y from this convergence state as y _n = y _n ^(T) , and further calculates the mapping g from the relationship of y _n = g (n) from the equation (6). Looking for.

  In the update formulas (55) and (56), the motion state of each lattice point is determined by a local relationship, and the motion is estimated while the local state is in an equilibrium state. Introduce more global structure to improve accuracy.

  The image matching apparatus 100 according to the present embodiment introduces a parametric motion model that defines a mapping relationship by a plurality of parameters (motion parameters). The outline of motion estimation using a parametric motion model will be described below.

  It is assumed that the movement of all grid points in the screen of the reference image (mapping relationship with the first grid point) is determined by a finite number of parametric motion models. Assuming that the motion parameters are known, the problem of determining the motion of each grid point in the screen is to determine which parametric motion model each grid point belongs to, that is, to select the optimal motion parameter. Become. This is to optimize a plurality of motion parameters, which is an extremely easy method compared to solving an optimization problem with respect to the motion that can be taken on the entire screen.

  However, since the motion parameter is not known, the result of motion estimation according to the update formulas of the local dynamic matching processes (55) and (56) is used, and the result of the motion estimation is obtained using a clustering method. Then, it is divided into a plurality of clusters (groups). In the present embodiment, by performing numerical analysis by the k-means method by the motion clustering processing unit 109, the second lattice point is obtained from the result of motion estimation according to the update formulas (55) and (56). Clustered into multiple regions. As a clustering method, a method other than the k-means method may be used.

  Then, the motion parameter estimation unit 110 estimates the motion parameter (first motion parameter) of the parametric motion model from the local motion in each cluster (in the region) using the least square method. Note that a method other than the method of least squares can be used as a method of estimating the motion parameter.

  After estimating the motion parameter, the motion parameter correcting unit 113 corrects the estimated motion parameter. Thereafter, the motion parameter assignment unit 111 assigns the motion parameters of the motion model to which each lattice point belongs to each lattice point. Specifically, the motion parameter assignment unit 111 introduces a label indicating which motion parameter each lattice point belongs to, and assigns a label, thereby assigning the optimum motion parameter (second motion parameter). Make an assignment. Therefore, the motion parameter assignment unit 111 sets an evaluation function for the motion model, and sets a label indicated by the motion parameter that minimizes the evaluation function at each lattice point.

  For each grid point, the first solution processing unit 101 obtains local motion by the update formulas (55) and (56). At the same time, each grid point also has motion from the parametric motion model. Desired. For this reason, a correction term indicating the force for correcting the shift of the local movement of each grid point and the motion model is introduced, the correction term is obtained by the correction term calculation unit 112, and the second solution processing unit 108 A local update is performed by solving the second equation of motion obtained by adding the obtained correction term to the first equation of motion by numerical analysis. FIG. 8 is an explanatory diagram showing a state of motion estimation using a motion model on an image.

  As described above, in the image matching apparatus 100 according to the present embodiment, the motion parameter estimation processing by the motion parameter processing unit 110 described above, the motion parameter allocation processing by the motion parameter allocation unit 111, the correction term calculation unit 112, and By using the motion division that repeats local motion update (second motion equation solution processing) by the second solution processing unit 108, the accurate motion of the image is estimated, and the first grid point and Image matching is performed by obtaining a mapping relationship with the second grid point. FIG. 10 is an explanatory diagram showing a conceptual flow of image matching using motion division by introducing the parametric motion model according to the first embodiment.

  For example, as shown in FIG. 9, let us consider a state in which a part of an object on an image is image-matched with an incorrect motion. If there is an erroneous movement as shown in FIG. 9, it will be visually recognized as a boundary distortion or the like. In such a case, correct image matching can be performed on the entire object by assigning accurate motion information to an erroneous motion region. The optimal motion parameter assignment in the present embodiment, that is, the determination of the parametric motion model to which the lattice points belong corresponds to assigning accurate motion information to such erroneous motion regions. As a result, it is possible to prevent a problem that a part of the characters is broken on the image.

  Next, details of motion estimation using such a parametric motion model will be described. First, a parametric motion model is defined. Each parametric motion model is labeled with an arbitrary number of labels αεL⊂Z. Assume that the label of the lattice point n is z (n) εL, and a set of regions (each motion region) having the label α is Nα = {n | z (n) = α}. A parametric motion model for each region of the label α is defined by equation (57).

動きパラメータａα^(m)は種々の形式を採用することができる。例えば、動きパラメータａα^(m)として、アファインモデルを用いた場合には（５８）式、２次形式モデルを用いた場合には（５９）式で示される。ここで、ｍは、イテレーションの変数である。

The movement parameter aα ^(m) can take various forms. For example, when the affine model is used as the motion parameter aα ^(m) , the equation (58) is used, and when the quadratic form model is used, the equation (59) is used. Here, m is an iteration variable.

この他、動きパラメータとして、回転変換モデル、ユークリッド変換モデル、相似変換モデル、射影変換モデル、Ｌｉｅ変換モデルなどを用いてもよい。

In addition, a rotation transformation model, an Euclidean transformation model, a similarity transformation model, a projective transformation model, a Lie transformation model, or the like may be used as the motion parameter.

  Next, motion parameter estimation will be described. The motion parameter estimation processing is performed by the motion parameter estimation unit 110. Since the motion vector determined by the parametric motion model is expressed by Expression (60), using this and Expression (57), the motion vector determined by the parametric motion model is expressed by Expression (61).

一方、対応点により求まる動きベクトルは、（６２）式で示される。

On the other hand, the motion vector obtained from the corresponding point is expressed by equation (62).

パラメトリック動きモデルによって定まる動きベクトルと対応点により求まる動きベクトルとの誤差を（６３）式で定義すると、誤差エネルギーは（６４）式で表される。従って、各動きベクトルの定義式（６１）、（６２）式により、誤差エネルギーは（６５）式のようになり、結局（６６）式で示されることになる。

When the error between the motion vector determined by the parametric motion model and the motion vector obtained from the corresponding point is defined by Equation (63), the error energy is expressed by Equation (64). Therefore, the error energy is expressed by the equation (65) according to the definition equations (61) and (62) of each motion vector, and is eventually expressed by the equation (66).

一方、ラベルα領域内の誤差は、（６７）式で示され、このラベルα領域内の誤差が（６８）式に示すように最小になるように動きパラメータが推定される。

On the other hand, the error in the label α region is expressed by equation (67), and the motion parameter is estimated so that the error in the label α region is minimized as shown in equation (68).

本実施の形態では、動きパラメータと得られた動きベクトルの誤差がガウス分布に従っていると仮定して、最適な動きパラメータは最小二乗法を使用した数値解析によって求めている。すなわち、最小二乗法における正規方程式は、（６９−１）、（６９−２）、（６９−３）式に示されるため、動きパラメータは（７０−１）式で示される。

In the present embodiment, assuming that the error between the motion parameter and the obtained motion vector follows a Gaussian distribution, the optimum motion parameter is obtained by numerical analysis using the least square method. That is, since the normal equation in the least square method is represented by equations (69-1), (69-2), and (69-3), the motion parameter is represented by equation (70-1).

かかる（７０−１）式は連立一次方程式であり、このため、動きパラメータは、（７０）式により特異値分解やＬＵ分解などの手法で求めることができる。

Since the equation (70-1) is a simultaneous linear equation, the motion parameter can be obtained by a method such as singular value decomposition or LU decomposition according to the equation (70).

  Next, parameter correction processing by the parameter correction unit 113 will be described. The motion parameter estimation by the motion parameter estimation unit 110 estimates the motion parameter from the local motion vector, but it is not known whether the motion model defined by the estimated motion parameter is suitable for the image. That is, the likelihood of the motion parameter for the image is not considered.

  This becomes obvious as the following problem. An estimation error is included in the local motion estimation, and an estimation error is also included in the motion parameter estimation. For this reason, since the obtained motion model does not necessarily match the image, the accuracy of label allocation may be reduced. As a result, the correct label assignment is not performed, and the next local motion estimation is performed using the information, so that the error further increases. Since this can be regarded as an error divergence effect due to feedback, it is necessary to cancel the error feedback at any stage. By optimizing the likelihood of the motion parameter image, the error of the motion parameter image can be reduced, and the feedback can be canceled sufficiently.

The likelihood of the motion parameter for the image is expressed by equations (70-2) and (70-3).

Expressions (70-2) and (70-3) are equivalent to the likelihood of labeling. However, since the labeling is optimized by the label, the motion parameter is not optimized. As a motion parameter optimization problem, the equation (70-4) may be solved.

  In equation (70-4), the standard deviation term of noise is omitted because it does not affect the result of maximum likelihood estimation. It is known that this problem can be solved stably by the Lucas-Kanade algorithm. Here, since it is desired to perform maximum likelihood estimation individually for each label α, the Lucas-Kanade algorithm is applied for each label.

When a set of regions having labels, that is, each motion region is represented by the equation (70-5), the motion parameter optimization problem is defined for each region and represented by the equation (70-6).

The Lucas-Kanade algorithm is approximated by iterating the solution iteratively as shown in equation (70-9), assuming that the nearest neighbor (70-8) equation of the true solution (70-7) is known. This is a method for obtaining a solution, and it is sufficient to perform minimization according to equation (70-10).

When the equation (70-10) is Taylor-expanded around aα, the equation (70-11) is obtained. Differentiating this equation by Δaα yields the equation (70-12).

When the Newton method is applied to such an equation, the update of the solution is as shown in equation (70-13). Here, H is Hessian expressed by the equation (70-14).

Therefore, the solution update equation is expressed by equation (70-15).

  Here, for example, for a sufficiently small constant ε (> 0), an end condition is defined such as end when the expression (70-16) is satisfied. It may be configured to end after repeating a predetermined number of iterations.

Since the expression (70-10) treats each pixel equivalently, if an outlier is included in the label, the solution (average) is shifted due to the influence. Therefore, weighting is applied to the equation (70-10), and the problem indicated by the equation (70-17) is considered.

By setting the weight w (n) to be low in the case of an outlier, the weight w (n) is robust against the influence of the outlier. If the update equation is derived in the same manner as the solution update (70-13), the equation (70-18) is obtained. The weight w (n) can be determined using an M estimation framework or the like.

  Therefore, in the present embodiment, the motion parameter estimated by the motion parameter estimation unit 110 by the motion parameter correction unit 113 may be corrected by the equations (70-18), (70-14), and (70-15).

Next, calculation of the elastic force corresponding to discontinuity by the elastic energy force calculation unit 104 will be described. Different parametric motion models generally have different motion parameters. For this reason, distortion is less when the elastic constant at the boundary between the two is set to zero. The boundary portion is a portion having a different label. That is, the elastic constant between the lattice points n ₁ and n ₂ is expressed by the equation (71). The elastic force corresponding to the discontinuity (hereinafter referred to as “discontinuous corresponding elastic force”) can be defined by the equation (72) by the equations (71) and (38).

このため、弾性エネルギー力計算部１０４では、（７１）式に基づき弾性係数を決定し、後述する第２の運動方程式の項となる不連続対応弾性力を（７２）式により計算している。

For this reason, the elastic energy force calculation unit 104 determines an elastic coefficient based on the equation (71), and calculates the discontinuous corresponding elastic force, which is a term of the second equation of motion described later, using the equation (72).

Next, calculation of the correction term by the correction term calculation unit 112 will be described. When the motion determined by the parametric motion model is different from the motion of each grid point, it is necessary to perform correction to approximate each grid point to the parametric motion model. For this reason, a correction term for approximation to a second equation of motion described later is added, and the correction term calculation unit 112 calculates the correction term. Assuming that such approximate correction is based on linear elastic energy, the correction term is a linear function having a difference between the motion vector determined by the parametric motion model and the motion vector of each lattice point as input, as shown in equation (73). It is defined as the correction force indicated by the output. Here, k ₂ represents the strength of correction.

また、パラメトリック動きモデルによって定まる動きと各格子点の動きとの差が所定値以上となる場合には、（７３）式に示す補正項では補正力が必要以上に大きくなってしまう。このような場合には、パラメトリック動きモデルによって定まる動きベクトルと各格子点の動きベクトルの差を入力とした非線形関数の出力で示される補正項、例えば、（７４）式のような対数をとった補正項を採用することにより、補正項の値が膨大になることを防止することができる。

Further, when the difference between the motion determined by the parametric motion model and the motion of each grid point is a predetermined value or more, the correction force in the correction term shown in the equation (73) becomes larger than necessary. In such a case, a correction term indicated by the output of a nonlinear function having the difference between the motion vector determined by the parametric motion model and the motion vector of each grid point as an input, for example, a logarithm as shown in equation (74) is taken. By adopting the correction term, it is possible to prevent the value of the correction term from becoming enormous.

次に、動きクラスタリング処理部１０９による動きクラスタリング処理について説明する。更新式（５５）、（５６）式による局所的な動的マッチング処理が完了した時点では、ラベルの割り当て処理が行われていないため、動きパラメータの推定を行うことができない。このため、更新式（５５）、（５６）式による局所的な動的マッチング処理が完了して最初に動きパラメータの推定処理を行う場合には、動きパラメータの推定処理の実行前に、動きクラスタリング処理部１０９によって、更新式（５５）、（５６）式によって求められた動き推定をクラスタリングすることによりラベルの初期値として使用している。ここで、更新式（５５）、（５６）式による動き推定のクラスタリング処理としては、輝度等の画素値に基づいて画像をｋ−Ｍｅａｎｓ法によって複数個にクラスタリングを行って画像を分割し、そのクラスタリング結果によって複数個の動きパラメータを求め、求めた動きパラメータを初期値としてｋ−Ｍｅａｎｓ法による数値解析処理により動きパラメータのラベルを求める。

Next, motion clustering processing by the motion clustering processing unit 109 will be described. When the local dynamic matching process according to the update expressions (55) and (56) is completed, the label allocation process is not performed, so that the motion parameter cannot be estimated. Therefore, when the motion parameter estimation process is performed first after the local dynamic matching process according to the update expressions (55) and (56) is completed, the motion clustering process is performed before the motion parameter estimation process is executed. The motion estimation obtained by the update equations (55) and (56) is clustered by the processing unit 109 and used as an initial value of the label. Here, as a clustering process of motion estimation by the update formulas (55) and (56), an image is divided into a plurality of images by the k-Means method based on pixel values such as luminance, and the image is divided. A plurality of motion parameters are obtained from the clustering result, and the motion parameter labels are obtained by numerical analysis processing by the k-Means method using the obtained motion parameters as initial values.

  Next, a motion parameter assignment process by the motion parameter assignment unit 111, that is, a label assignment process will be described. In the motion parameter assignment process by the motion parameter assignment unit 111, the optimization problem shown in the equation (75) is released as a motion parameter assignment problem for the lattice point n.

ここで、Ｇはコスト関数であり、本実施の形態では、非特許文献「A. M. Tekalp, “digital video processing”, Prentice Hall, 1995，pp. 103」に示されている、（７６）式のようなＭＰＣ（Maximum Pel Count）を用いている。

Here, G is a cost function, and in this embodiment, as shown in the non-patent document “AM Tekalp,“ digital video processing ”, Prentice Hall, 1995, pp. 103, MPC (Maximum Pel Count) is used.

ここで、Ｂ⊂Λ²は、格子点ｎを中心とした局所領域(ブロック)を表す集合であり、Ｎ（Ｂ）は集合の要素数である。また、非特許文献「A. M. Tekalp, “digital video processing”, Prentice Hall, 1995，pp. 102-103」に示されているＭＳＥ(Mean Squared Error), ＭＡＤ(Mean Absolute Error)などを用いても良い。

Here, B⊂Λ ² is a set representing a local region (block) centered on the lattice point n, and N (B) is the number of elements in the set. Further, MSE (Mean Squared Error), MAD (Mean Absolute Error) and the like shown in non-patent literature “AM Tekalp,“ digital video processing ”, Prentice Hall, 1995, pp. 102-103” may be used. .

  Here, it is inferred that the accuracy of the motion parameter needs to be high in order to assign the motion parameter accurately. As shown in FIG. 11, the larger the area where the image matching is correct (the white area in FIG. 11) with respect to the matching target area, the higher the estimation accuracy of the motion parameter for that area. This is a characteristic of the least square method (the estimation accuracy improves as the number of samples increases due to the central limit theorem). For this reason, if this motion parameter is assigned to an area where image matching has not been performed (gray area in the figure), accurate image matching is likely to occur. Conversely, when the matching area is narrow, the motion parameter estimation accuracy decreases. In this case, even when the low-precision motion parameter is expanded, the deviation from the true motion parameter becomes large, and it is difficult to perform accurate image matching.

  Here, it is possible to optimize the entire screen. Focusing on the localization of image processing, the localization of the image is formulated by the Gibbs distribution equivalent to the Markov property, and the optimization problem that we want to find as a prior distribution is the posterior distribution maximization problem It can be formulated. In addition, the posterior distribution maximization problem is basically NP-hard, and it is difficult to obtain a global optimal solution, and is usually solved by simulated annealing or the like. In order to solve the problem more efficiently, the non-patent document “Y. Boykov, O. Veksler and R. Zabih,“ Fast Approximate Energy Minimization via Graph Cuts ”, IEEE Trans. On Pattern Analysis and Machine Int., Vol. 23, No. 11, pp. 1222-1239, 2001 ", and can be formulated using graph theory graph cuts.

  Next, the second equation of motion for which the second solution processing unit 108 performs solution processing will be described. As shown in the equation (77), the second equation of motion is obtained by adding a correction term of the equation (73) or (74) to the first equation of motion of the equations (55) and (56) and further elasticity. As energy force, it has the term of the elastic force corresponding to discontinuity shown in (71) and (72) types. The second solution processing unit 108 solves the second equation of motion by a numerical analysis process using a discrete variable method such as the Euler method, like the first equation of motion.

次に、以上のように構成された実施の形態１にかかる画像マッチング装置１００による画像マッチング処理について説明する。図１２は、実施の形態１にかかる画像マッチング装置１００による画像マッチングの全体処理の手順を示すフローチャートである。まず、カウント回数Ｌを１に初期化し（ステップＳ１２０１）、第１解法処理部１０１、画像エネルギー力計算部１０３、弾性エネルギー力計算部１０４、摩擦力計算部１０５による（５５）、（５６）式による局所的な動的マッチング処理を行う（ステップＳ１２０２）。

Next, an image matching process performed by the image matching apparatus 100 according to the first embodiment configured as described above will be described. FIG. 12 is a flowchart of an overall image matching process performed by the image matching apparatus 100 according to the first embodiment. First, the count number L is initialized to 1 (step S1201), and the first solution processing unit 101, the image energy force calculation unit 103, the elastic energy force calculation unit 104, and the frictional force calculation unit 105 perform the expressions (55) and (56). A local dynamic matching process is performed by (step S1202).

  When the dynamic matching processing is completed, motion clustering processing is performed by the motion clustering processing 109 (step S1203), and the reference image screen is divided into a plurality of frames to obtain initial values of motion parameter labels. Next, motion parameter estimation processing is performed by the motion parameter estimation unit 110 (step S1204), and motion parameters for each label area are estimated. Then, the motion parameter correcting unit 113 corrects the motion parameter estimated by the motion parameter estimating unit 110 using the equations (70-18), (70-14), and (70-15) (step S1205). Then, the motion parameter allocating unit 111 determines an optimal motion parameter from the corrected motion parameter, and the determined optimal motion parameter is allocated to each second grid point in the region corresponding to the label. Parameter assignment processing is performed (step S1206).

Next, the dynamic matching process with the correction term added, that is, the solution process of the second equation of motion of the equation (77) by the correction term calculation unit 112 and the second solution processing unit 108 is performed to obtain y _n ^(LT2) ). Obtained (step S1207) and stored in the frame memory 106. Then, the count number L is incremented by 1 (step S1208), and it is checked whether the count number L exceeds a predetermined number of times (step S1209).

  When the count number L does not exceed the predetermined number of times (step S1209: No), the processing from step S1204 to S1208 is repeated. That is, the motion parameter estimation processing by the motion parameter processing unit 110, the motion parameter correction processing by the motion parameter correction unit 113, the motion parameter allocation processing by the motion parameter allocation unit 111, the correction term calculation unit 112, and the second solution method The local motion update (second motion equation solving process) by the processing unit 108 is repeated.

On the other hand, if the count number L exceeds a predetermined number in step S1209 (step S1209: Yes), the second solution processing unit 108 sets y _n to y n for all grid points n. _n ^(T) is set (step S1210). Then, the mapping processing unit 107 obtains a correspondence relationship between the target image and the reference image, that is, a mapping (step S1211). Image matching is performed by estimating the exact motion of the image by the processing using such motion division and obtaining the mapping relationship between the first lattice point and the second lattice point.

  Next, the local dynamic matching process in step S1202 will be described. FIG. 13 is a flowchart of a local dynamic matching process according to the first embodiment. Specifically, the image matching processing is realized by solving the equations (55) and (56) shown above by numerical analysis.

First, the first solution processing unit 101 sets time τ ⁽⁰⁾ = 0 (step S1301), and sets initial values y _n ⁽⁰⁾ = n and v _n ⁽⁰⁾ = 0 (step S1302). Thereby, the equation (56) is executed.

Next, the image energy force calculation unit 103 calculates the image correlation potential force F _i ^(m) (n) in m steps for all the lattice points n (step S1303). The calculation processing of the image energy force F _i ^(m) (n) will be described later.

Next, the elastic energy force calculation unit 104 calculates the elastic energy force F _k ^(m) (n) in m steps with respect to all the lattice points n by the equation (38) (step S1304). Then, the frictional force calculation unit 105 calculates the frictional force [−μv _n ^(m) ] in m steps for all the lattice points n (step S1305).

Next, the first solution processing unit 101 uses the image energy force F _i ^(m) (n) and the elastic energy force F _k ^(m) (n) obtained in Steps S1303 to S1305 as the update equations of Expression (55). The frictional force [−μv _n ^(m) ] is updated (step S1306).

Next, the first solution processing unit 101 stores the value of y _n ^{(m) in} the frame memory 106 (step S1307). Then, the first solution processing unit 101 updates τ ^{(m + 1)} = τ ^(m) + h (step S1308), and increases m by 1 (step S1309). Then, it is determined whether or not τ ^{(m + 1)} has exceeded a predetermined time T (step S1310). If not, step S1303 to S1309 are repeatedly executed.

On the other hand, tau ^{(m + 1)} is the case beyond the T is the first solution processing unit 101, for all grid points n, sets the y _n ^(T) to y _n (step S1311).

  Next, the image energy force calculation process in step S1303 will be described. FIG. 14 is a flowchart showing a procedure of image energy force calculation processing by the image energy force calculation unit 103.

First, the image energy force calculation unit 103 calculates the sampling point closest to y _n (τ) as the local space center y _{c using} equation (41) (step S1401). Next, the image energy force calculation unit 103 sets the adjacent space L defined by the equation (42) (step S1402), and calculates the local search set Ω by the equation (43) (step S1403).

Next, the image energy force calculation unit 103 calculates y _min as a local optimization calculation using the equation (46) (step S1404), and uses the local optimization y _min obtained in step S1404 for normalization. Then, d = y _min −y _n ^(m) is calculated by the equation (45) (step S1405). Then, the image energy force calculation unit 103 calculates the image energy force F _i ^(m) (n) by using equation (44) using y _min , d, and the like (step S1406).

  Next, the motion clustering processing by the motion clustering processing unit 109 in step S1203 will be described. FIG. 15 is a flowchart of a motion clustering process performed by the motion clustering processing unit 109 according to the first embodiment.

  First, based on pixel values such as luminance as a result of local dynamic matching processing, the motion clustering processing unit 109 performs clustering on a plurality of reference images by the k-Means method to divide the images and label them. Obtained (step S1501).

  Then, based on the label based on the clustering result, an initial motion parameter is estimated to obtain a plurality of initial motion parameters (step S1502). Next, the motion parameter label α is obtained by numerical analysis processing by the k-Means method using the obtained initial motion parameter as an initial value (step S1503). Here, a plurality of labels α are obtained corresponding to a plurality of initial motion parameters. Such motion clustering processing is executed only once immediately after the local dynamic matching processing is completed, as shown in FIG.

  Next, the motion parameter estimation processing by the motion parameter estimation unit 110 in step S1204 will be described. FIG. 16 is a flowchart of a motion parameter estimation process performed by the motion parameter estimation unit 110 according to the first embodiment.

  First, the motion parameter estimation unit 110 applies A to all lattice points n with respect to the label α obtained by the motion clustering processing unit 109, and b using the equation (69-2). Is calculated (step S1601). Then, the solution process is performed by numerical analysis by the singular value decomposition method (or LU method) for the equation (69-1) (step S1602). Then, the processes in steps S1601 and S1602 are repeated for all labels α (step S1603). The motion parameter estimated in this way is corrected by equations (70-18), (70-14), and (70-15) by the motion parameter correcting unit 113 in step S1205.

  Next, the motion parameter assignment processing by the motion parameter assignment unit 111 in step S1206 will be described. FIG. 17 is a flowchart of a motion parameter assignment process performed by the motion parameter assignment unit 111 according to the first embodiment.

First, the motion parameter assignment unit 111 obtains α ⁰ that minimizes the cost function G ⁽ⁿ⁾ expressed by the equation (76) for all labels α (step S1701). This processing is to execute equation (75). Then, the motion parameter assignment unit 111 sets α ⁰ as a label for the lattice point n (step S1702). With this process, the optimum motion parameter is assigned to the lattice point n. Then, the motion parameter assignment unit 111 repeatedly executes the processes of steps S1701 and S1702 for all the lattice points n (step S1703). Thereby, an optimal motion parameter is assigned to all the lattice points n.

  Next, the dynamic matching process with the correction term added in step S1207 will be described. FIG. 18 is a flowchart of a dynamic matching process procedure to which correction terms according to the first embodiment are added. Specifically, the image matching process is realized by solving the equation (77) shown above by numerical analysis.

First, the second solution processing unit 108 sets time τ ⁽⁰⁾ = 0 (step S1801), and sets initial values y _n ⁽⁰⁾ = n and v _n ⁽⁰⁾ = 0 (step S1802). Thereby, the second expression of the expression (77) is executed.

Next, the image energy force calculation unit 103 calculates the image correlation potential force F _i ^(m) (n) in m steps for all grid points n (step S1803). The calculation processing of the image energy force F _i ^(m) (n) is performed in the same manner as the processing described with reference to FIG.

Next, the elastic energy force calculation unit 104 calculates the discontinuity-corresponding elastic force F _dk ^(m) (n) in m steps for all the lattice points n by the equations (71) and (72) (step S1804). . The calculation processing of the discontinuity corresponding elastic force F _dk ^(m) (n) will be described later.

Next, the correction term calculation unit 112 calculates the correction term F _d ^(m) (n) in m steps for all lattice points n by the equation (73) (step S1805). Then, the frictional force calculation unit 105 calculates the frictional force [−μ (τ ^(m) ) v _n ^(m) ] in m steps for all lattice points n (step S1806).

Next, the second solution processing unit 108 changes the equation (77) to the image energy force F _i ^(m) (n) and the discontinuity corresponding elastic force F _dk ^(m) ( n), the correction term F _d ^(m) (n), and the frictional force [−μ (τ ^(m) ) v _n ^(m) ] are updated (step S1807).

Next, the value of y _n ^(m) is stored in the frame memory 106 by the second solution processing unit 108 (step S1808). Then, the second solution processing unit 108 updates τ ^{(m + 1)} = τ ^(m) + h (step S1809), and increases m by 1 (step S1810). Then, it is determined whether or not τ ^{(m + 1)} has exceeded a predetermined time T ₁ (step S1811). If not, step S1803 to S1810 are repeatedly executed.

Next, the calculation processing of the elastic force F _dk ^(m) (n) corresponding to the discontinuity by the elastic energy force calculation unit 104 in step S1804 will be described. FIG. 19 is a flowchart illustrating a procedure of calculation processing of the discontinuous corresponding elastic force F _dk ^(m) (n) by the elastic energy force calculation unit 104 according to the first embodiment.

First, elastic energy force calculating unit 104 checks whether the label z grid point n ₁ in m steps ^(m) (n ₁₎ and the grid point n ₂ label z ^(m) (n ₂₎ is equal to ( Step S1901). This is to check whether the lattice point n ₁ and the lattice point n ₂ belong to a motion model having the same motion parameter.

When the grid point n ₁ label z ^(m) (n ₁₎ and the grid point n ₂ label z ^(m) (n ₂₎ are equal (step S1901: Yes), i.e., grid points n ₁ and the lattice point If n ₂ belongs to the motion model having the same motion parameter, the elastic coefficient k ^(m) (n ₁ , n ₂ ) between the lattice point n ₁ and the lattice point n ₂ is set to k ₁ , (step S1902). Then, the discontinuous corresponding elastic force F _dk ^(m) (n) is calculated by the equation (72) (step S1904).

On the other hand, in step S1901, if the grid point n ₁ label z ^(m) (n ₁₎ and the grid point n ₂ label z ^(m) (n ₂₎ are not equal (Step S1901: No), i.e., grid points If n ₁ and lattice point n ₂ belong to motion models having different motion parameters, the elastic coefficient k ^(m) (n ₁ , n ₂ ) is set to 0 (step S1903). As a result, the discontinuity-corresponding elastic force F _dk ^(m) (n) becomes zero according to equation (72). Note that the processing from step S1901 to S1903 is processing for executing equation (71).

  The discontinuity-corresponding elastic force calculated in this way is input to the term of the update equation (77) in step S1807.

  As described above, in the image matching apparatus 100 according to the first embodiment, after the local dynamic matching process according to the update equations (55) and (56) is finished, the reference image is divided into a plurality of regions and each region is divided. For each, the optimal motion parameter of the motion parameter for classifying the motion model defining the motion of the second grid point is obtained, and the solution processing by numerical analysis of the second motion equation (77) indicating the motion model is performed. Since the correspondence between the target image and the reference image is obtained, the discontinuity of the region boundary and the uniformity within the region can be expressed robustly, and more accurate image matching can be performed. As a result, an interpolation frame with higher image quality can be generated.

(Embodiment 2)
The image matching apparatus according to the second embodiment uses M estimation, which is a robust statistical technique, in motion parameter estimation processing.

  In the first embodiment, in the motion parameter estimation process, the motion parameter is estimated using the least square estimation method on the assumption that the error between the motion parameter and the obtained motion vector follows a Gaussian distribution. However, when it cannot be said that the error between the motion parameter and the motion vector simply follows a Gaussian distribution, it is preferable to use M estimation, which is a robust statistical technique. Such a case is, for example, a case where a large outlier is included in a part of an image, and in dynamic matching, the boundary between motion and motion tends to be in such a situation due to distortion.

  The configuration of the image matching apparatus according to the second embodiment is the same as that of the first embodiment. The overall image matching process, local dynamic matching process, motion clustering process, motion parameter assignment process, and dynamic matching process with correction terms in the second embodiment are the same as in the first embodiment. Is called.

ここで、ｃ（＞０）は定数であり、正規分布の場合にはｃ＝６．０となる。動きパラメータ推定部１１０は、重み関数ｗ(ｎ)を用いて標本平均ｚ⁽ⁱ⁾、標本偏差ｓ⁽ⁱ⁾を（８０）式のように計算する。 A standard normal variable z (n) in the M estimation is expressed by equation (78). According to Tukey's Biweight method, the weight function w (n) is expressed by equation (79).

Here, c (> 0) is a constant, and c = 6.0 in the case of normal distribution. The motion parameter estimation unit 110 calculates the sample average z ⁽ⁱ⁾ and the sample deviation s ⁽ⁱ⁾ using the weight function w (n) as shown in Equation (80).

ここで上付の(i)は、イテレーションを示している。また、動きパラメータ推定部１１０では、母分散は、標本偏差標本偏差ｓ⁽ⁱ⁾を（８１−１）、（８１−２）式で示すようにプラグインして推定することとする。

Here, the superscript (i) indicates an iteration. In the motion parameter estimation unit 110, the population variance is estimated by plugging in the sample deviation sample deviation s ⁽ⁱ⁾ as shown by the equations (81-1) and (81-2).

そして、動きパラメータ推定部１１０によって、母標準偏差の計算と重み関数の計算を交互にＩ回繰り返すことによってＭ推定による動きパラメータの推定が行われる。

Then, the motion parameter estimation unit 110 performs motion parameter estimation by M estimation by alternately repeating the calculation of the mother standard deviation and the weighting function I times.

  FIG. 20 is a flowchart of a motion parameter estimation process performed by the motion parameter estimation unit 110 according to the second embodiment. First, the motion parameter estimation unit 110 initializes the weight w (n) to 1 (step S2001). Then, with respect to the label α obtained by the motion clustering processing unit 109, A is calculated from the equation (81-2) and b is calculated from the equation (81-3) for all the lattice points n (step S2002). . Then, the equation (81-1) is solved by numerical analysis using the singular value decomposition method (or LU method) (step S2003).

Next, the motion parameter estimation unit 110 calculates the standard normal variable z (n) according to the equation (78) (step S2004). Then, using the weight function w (n), the sample average z ⁽ⁱ⁾ and the sample deviation s ⁽ⁱ⁾ are calculated according to the equation (80) (step S2005). Next, the motion parameter estimation unit 110 updates the weight w (n) (step S2006) and increases the iteration i by 1 (step S2007). Then, it is checked whether i has exceeded a predetermined number of times I (step S2008). If not, the processes from step S2002 to S2007 are repeatedly executed.

  On the other hand, if i exceeds the predetermined number I in step S2008, the processes in steps S1601 and S1602 are repeated for all labels α (step S2009).

  As described above, the image matching apparatus according to the second embodiment uses the M estimation, which is a robust statistical technique, in the motion parameter estimation process, and therefore, when the error between the motion parameter and the motion vector does not follow the Gaussian distribution. However, it becomes possible to robustly express the discontinuity of the region boundary and the uniformity within the region, and more accurate image matching can be performed. As a result, a higher quality interpolated frame can be generated. it can.

(Embodiment 3)
In the image matching device according to the third embodiment, when the second grid point is included in the occlusion area, which is an area that exists on the target image but does not exist on the reference image, the motion parameter is determined at the grid point. Is assigned, and motion parameters are assigned only to lattice points not included in the occlusion area.

  An area that exists on the target image but does not exist on the reference image is called a covered occlusion area, and is an area that cannot inherently perform image matching. In the image matching apparatus 100 according to the first embodiment, since motion parameters are also assigned to the covered occlusion area, the result of motion parameter estimation may be distorted.

  In the image matching apparatus according to the third embodiment, a covered occlusion area is detected in advance, and motion parameters are not assigned to lattice points in the occlusion area.

  FIG. 21 is a block diagram of a configuration of an image matching apparatus 2100 according to the third embodiment. The image matching apparatus 2100 according to the third embodiment is different from the first embodiment in that an occlusion detection unit 2101 is provided, and the other configuration is the same as that of the image matching apparatus according to the first embodiment. The overall image matching process, local dynamic matching process, motion clustering process, motion parameter estimation process, and dynamic matching process with correction terms in the third embodiment are performed in the same manner as in the first embodiment. Is called.

  The occlusion detection unit 2101 is a processing unit that detects a covered occlusion area. 22 and 23 are explanatory diagrams showing the state of the covered occlusion area. The covered occlusion area typically appears in front of a moving object, as shown in FIG. In the state shown in FIG. 22, the area where the spring is contracted is the covered occlusion area.

  Such a covered occlusion area can be formulated as follows. As shown in FIG. 23, 8 lattices around lattice points are considered. Then, the degree of deformation from the original square lattice is examined by the sum of the areas of the triangles formed from the eight lattices. Such an area S (n) is defined by an outer product of vectors as shown in equation (82). The area ratio r (n) is expressed by the equation (83), where the area of the square lattice is 4. In the equation (83), x represents the vector cross product calculation.

そして、この面積比が閾値Ｔよりも小さければ、図２２に示すような縮んでいる領域である、すなわちcoveredオクルージョン領域の可能性が高い。本実施の形態にかかるオクルージョン検出部２１０１では、（８３）式に示す面積比を閾値Ｔと比較することにより、オクルージョン領域か否かを判定している。そして、動きパラメータ割り当て部１１１では、オクルージョン検出部２１０１によりcoveredオクルージョン領域であると判定された格子点に関しては、動きパラメータの割り当てをおこなわないように構成している。

If the area ratio is smaller than the threshold value T, there is a high possibility that the area is shrunk as shown in FIG. 22, that is, a covered occlusion area. The occlusion detection unit 2101 according to the present embodiment determines whether or not it is an occlusion region by comparing the area ratio shown in the equation (83) with a threshold value T. The motion parameter assignment unit 111 is configured not to assign a motion parameter to a lattice point determined to be a covered occlusion region by the occlusion detection unit 2101.

  FIG. 24 is a flowchart of a motion parameter assignment process performed by the image matching apparatus 2100 according to the second embodiment.

  First, the occlusion detection unit 2101 calculates the total area S (n) of the triangles formed from the eight lattices around the lattice points using the equation (82) (step S2401). Next, the occlusion detection unit 2101 uses S (n) to calculate the area ratio r (n) using the equation (83) (step S2402).

Next, the occlusion detection unit 2101 determines whether the area ratio r (n) is smaller than a predetermined threshold T, that is, whether the area is an occlusion region (step S2403). When the area ratio r (n) is not smaller than the predetermined threshold T (step S2403: No), that is, when the area ratio r (n) is not an occlusion area, the motion parameter assignment unit 111 sets (76 for all labels α. ) Α ⁰ that minimizes the cost function G ⁽ⁿ⁾ shown in equation ⁽²⁾ is obtained (step S2404). Then, the motion parameter assignment unit 111 sets α ⁰ as a label for the lattice point n (step S2405). With this process, the optimum motion parameter is assigned to the lattice point n.

  On the other hand, in step S2403, if the area ratio r (n) is smaller than the predetermined threshold T (step S2403: Yes), that is, if it is an occlusion region, the optimal motion parameters in steps S2404 and S2405 are assigned. Not performed.

  Then, the motion parameter assignment unit 111 repeatedly executes the processing of steps S2401 to S2405 for all the lattice points n (step S2406). As a result, an optimal motion parameter is assigned to all lattice points n that are not in the occlusion region.

  As described above, in the image matching apparatus according to the third embodiment, the covered occlusion area is detected in advance, and no motion parameter is assigned to the lattice point of the occlusion area, so that distortion due to the covered occlusion is minimized. Therefore, more accurate image matching can be performed, and as a result, an interpolation frame with higher image quality can be generated.

(Embodiment 4)
In the image matching apparatus according to the fourth embodiment, motion parameters are not assigned to regions with low reliability.

  In the image matching apparatus according to the first embodiment, all optimization of the motion parameter assignment process is trusted. However, there is actually a region where the reliability of optimization of the motion parameter assignment process is low. For example, if there are two or more values that approximate the optimal value of the motion parameter, it can be said that the reliability is low. That is, when there are two or more values that approximate the optimal value for different motion parameters, the cost function surface is considered to be multimodal, and the reliability of the optimal value of the motion parameter is It is considered low. In the present embodiment, the motion parameter assignment process is not performed for such a low-reliability region.

  The configuration of the image matching apparatus according to the fourth embodiment is the same as that of the first embodiment. The overall image matching process, local dynamic matching process, motion clustering process, motion parameter estimation process, and dynamic matching process with correction terms in the fourth embodiment are performed in the same manner as in the first embodiment. Is called.

  The motion parameter assignment unit 111 according to the present embodiment determines whether or not a lattice point is included in the low reliability region, and performs a process that does not assign a motion parameter when the lattice point is included in the low reliability region. ing.

  The determination of the low reliability region can be formulated by the equation (84).

ここで、α⁰は、求めた最適なラベルであり、α¹は次に最適なラベルである。また、ε（＞０）は定数である。ただし、動きパラメータの値が近似する場合には、（８４）式に該当する場合があるため、（８５）式を満たすほど動きパラメータが近似する場合には、（８４）式は適用しないこととする。

Here, α ⁰ is the optimum label obtained, and α ¹ is the next optimum label. Also, ε (> 0) is a constant. However, when the motion parameter value approximates, it may correspond to the equation (84). Therefore, when the motion parameter approximates to satisfy the equation (85), the equation (84) is not applied. To do.

図２５は、実施の形態４にかかる動きパラメータ割り当て部１１１による動きパラメータ割り当て処理の手順を示すフローチャートである。まず、動きパラメータ割り当て部１１１は、全てのラベルαに対して、（７６）式に示すコスト関数Ｇ⁽ⁿ⁾が最小になるα⁰を求める（ステップＳ２５０１）。次に、動きパラメータ割り当て部１１１は、α⁰より一つ大きいα¹、すなわち次に最適なラベルを求める（ステップＳ２５０２）。

FIG. 25 is a flowchart of a motion parameter assignment process performed by the motion parameter assignment unit 111 according to the fourth embodiment. First, the motion parameter assignment unit 111 obtains α ⁰ that minimizes the cost function G ⁽ⁿ⁾ expressed by the equation (76) for all labels α (step S2501). Next, the motion parameter assignment unit 111 obtains α ^{1 which} is one greater than α ⁰ , that is, the next optimal label (step S2502).

  Next, it is examined whether or not the motion parameter is approximate, that is, whether or not the equation (85) is satisfied (step S2503). If the formula (85) is not satisfied, it is checked whether or not the low reliability region is satisfied based on whether the formula (86) is satisfied (step S2504).

If it is determined that equation (86) is not satisfied and the region is not a low-reliability region (step S2504: No), α ⁰ obtained in step S2501 is set as an optimum label for the lattice point n (step S2505). ). On the other hand, if it is determined in step S2504 that the expression (86) is satisfied and the region is a low-reliability region (step S2504: Yes), the process of step S2505 is not performed, and optimal label allocation, that is, optimal No assignment of motion parameters is performed.

If the motion parameter is approximated in step S2503, that is, if equation (85) is satisfied (step S2503: Yes), the low-reliability region determination process in step S2404 is not performed, and the determination is performed in step S2501. α ⁰ is set as an optimum label for the lattice point n (step S2505).

  With this process, the optimum motion parameter is assigned to the lattice point n. Then, the motion parameter assignment unit 111 repeatedly executes the processing from step S2501 to S2505 for all the lattice points n (step S2506).

  As described above, in the image matching apparatus according to the fourth embodiment, since motion parameters are not assigned to regions with low reliability, discontinuity of region boundaries and uniformity within regions can be expressed robustly. As a result, more accurate image matching can be performed, and as a result, an interpolation frame with higher image quality can be generated.

  The image matching devices according to the first to fourth embodiments include a control device such as a CPU, a storage device such as a ROM (Read Only Memory) and a RAM, an external storage device such as an HDD and a CD drive device, and a display device such as a display device. It has an input device such as a keyboard and a mouse, and has a hardware configuration using a normal computer.

  The image matching program executed by the image matching apparatus according to the first to fourth embodiments is a file in an installable or executable format, and is a CD-ROM, flexible disk (FD), CD-R, DVD (Digital Versatile Disk). And the like recorded on a computer-readable recording medium.

  Further, the image matching program executed by the image matching apparatus according to the first to fourth embodiments may be provided by being stored on a computer connected to a network such as the Internet and downloaded via the network. good. Further, the image matching program executed by the image matching apparatus of the present embodiment may be provided or distributed via a network such as the Internet.

  Further, the image matching program of the present embodiment may be configured to be provided by being incorporated in advance in a ROM or the like.

  The image matching program executed by the image matching apparatus according to the first to fourth embodiments includes the above-described units (first solution processing unit 101, image energy force calculation unit 103, elastic energy force calculation unit 104, and friction force calculation). Unit 105, second solution processing unit 108, motion clustering processing unit 109, motion parameter estimation unit 110, motion parameter allocation unit 111, correction term calculation unit 112, mapping processing unit 107, occlusion detection unit 2101) As the actual hardware, a CPU (processor) reads out and executes an image matching program from the storage medium, and the respective units are loaded onto the main storage device. Solution processing unit 101, image energy force calculation unit 103, elastic energy force Arithmetic unit 104, frictional force calculation unit 105, second solution processing unit 108, motion clustering processing unit 109, motion parameter estimation unit 110, motion parameter allocation unit 111, correction term calculation unit 112, and mapping process Unit 107 and occlusion detection unit 2101) are generated on the main memory.

  It should be noted that the present invention is not limited to the above-described embodiment as it is, and can be embodied by modifying the constituent elements without departing from the scope of the invention in the implementation stage. In addition, various inventions can be formed by appropriately combining a plurality of constituent elements disclosed in the above embodiments. For example, some components may be deleted from all the components shown in the embodiment. Furthermore, constituent elements over different embodiments may be appropriately combined.

1 is a block diagram illustrating a configuration of an image matching apparatus according to a first embodiment. It is explanatory drawing which shows the mapping from the curved surface on a target image to the curved surface on a reference image. It is explanatory drawing which shows the mapping from the point x of the curved surface X on a target image to the point y of the curved surface Y on a reference image. It is a schematic diagram for demonstrating the static optimal point in the case of searching the point of potential energy minimization. It is a schematic diagram for demonstrating the global optimal point in the case of searching for the point of potential energy minimization. It is explanatory drawing which shows the concept of image energy force. It is explanatory drawing which shows the concept of image energy force. It is explanatory drawing which shows the state on the image of the motion estimation using a motion model. It is explanatory drawing which shows the state by which the part of the object on an image is image-matched by the incorrect motion. It is explanatory drawing which shows the conceptual flow of the image matching which introduce | transduced the parametric motion model concerning Embodiment 1 and used motion division. It is a schematic diagram for demonstrating the relationship between a matching object area | region and the estimation accuracy of a motion parameter. 3 is a flowchart illustrating a procedure of overall image matching processing by the image matching apparatus 100 according to the first embodiment; 3 is a flowchart showing a procedure of local dynamic matching processing according to the first exemplary embodiment; 5 is a flowchart illustrating a procedure of image energy force calculation processing by an image energy force calculation unit 103. 10 is a flowchart showing a procedure of motion clustering processing by the motion clustering processing unit 109 according to the first embodiment. 3 is a flowchart showing a procedure of motion parameter estimation processing by a motion parameter estimation unit 110 according to the first exemplary embodiment. 4 is a flowchart illustrating a procedure of motion parameter assignment processing by a motion parameter assignment unit 111 according to the first embodiment. 6 is a flowchart showing a procedure of dynamic matching processing to which a correction term according to the first embodiment is added. 4 is a flowchart illustrating a procedure of calculation processing of a discontinuity-corresponding elastic force F _dk ^(m) (n) by an elastic energy force calculation unit 104 according to the first embodiment. 10 is a flowchart showing a procedure of motion parameter estimation processing by a motion parameter estimation unit 110 according to the second exemplary embodiment. It is a block diagram which shows the structure of the image matching apparatus 2100 concerning Embodiment 3. FIG. It is a schematic diagram for demonstrating a covered occlusion area | region. It is a schematic diagram for demonstrating a covered occlusion area | region. 12 is a flowchart illustrating a procedure of motion parameter assignment processing by the image matching apparatus 2100 according to the second embodiment. 14 is a flowchart illustrating a procedure of motion parameter assignment processing by a motion parameter assignment unit 111 according to the fourth embodiment.

Explanation of symbols

100, 2100 Image matching device 101 First solution processor 103 Image energy force calculator 104 Elastic energy force calculator 105 Friction force calculator 106 Frame memory 107 Mapping processor 108 Second solution processor 109 Motion clustering processor 110 Motion parameters Estimation unit 111 Motion parameter assignment unit 112 Correction term calculation unit 2101 Occlusion detection unit

Claims (17)

An image matching device for obtaining a correspondence between a target image and a reference image,
There is an image correlation between each of the plurality of first grid points on the target image and each of the plurality of second grid points corresponding to the first grid points on the reference image. Based on the potential energy gradient based on the position of each second grid point and the position of each first grid point corresponding to the second grid point. An image energy force calculator for calculating the image energy force received by the grid points;
An elastic energy force calculator that calculates a first elastic energy force received from elastic energy between each second lattice point and another second lattice point adjacent to the second lattice point;
A friction force calculator for calculating a friction force acting on each of the second lattice points;
A first solution processing unit that solves the first equation of motion related to each second lattice point by the image energy force, the first elastic energy force, and the friction force by numerical analysis;
A motion clustering processing unit that divides the reference image into a plurality of regions and assigns a label for identifying the region to each of the divided regions, as a result of the solution processing by the first solution processing unit;
A motion parameter estimator that estimates, by numerical analysis, a first motion parameter that is a parameter that defines a motion model indicating a mapping relationship between the second lattice point and the first lattice point in the region to which the label is assigned; ,
A motion parameter correction unit that corrects the first motion parameter so that the potential energy between the point on the target image and the point on the reference image determined by the estimated first motion parameter is minimized; ,
From the corrected first motion parameter, a second motion parameter that is the first motion parameter optimal for each second grid point in the region corresponding to the label is determined and determined. A motion parameter assigning unit that assigns the second motion parameter to each second grid point in the region corresponding to the label;
A correction term calculation unit that calculates a correction term for correcting the first equation of motion based on the assigned second motion parameter and the result of the solution processing by the first solution processing unit;
A second solution processing unit that solves the second equation of motion relating to each of the second lattice points obtained by adding the correction term to the first equation of motion by numerical analysis; and a solution by the second solution processing unit A mapping processing unit for obtaining a correspondence relationship between the target image and the reference image from a processing result;
An image matching apparatus comprising:   The image matching apparatus according to claim 1, wherein the motion parameter estimation unit estimates the motion parameter by a numerical analysis process using a least square method based on the label.   The image matching apparatus according to claim 1, wherein the motion parameter estimation unit estimates the motion parameter by numerical analysis based on M estimation based on the label.   The image matching apparatus according to claim 1, wherein the motion parameter correction unit corrects the motion parameter by obtaining a point at which the potential energy is minimum using a steepest descent method. The elastic energy force calculation unit further receives a first energy received from an elastic energy between the second lattice point and the third lattice point after the movement of the second lattice point determined by the second motion parameter. Calculate the elastic energy force of 2,
The second solution processing unit performs solution processing by numerical analysis on the second equation of motion based on the image energy force, the second elastic energy force, the friction force, and the correction term. The image matching apparatus according to claim 1.   The elastic energy force calculation unit calculates the second elastic energy force for the second lattice point and the third lattice point in the region where the label is the same, and the second elastic energy force in the region where the label is different. 6. The image matching according to claim 5, wherein the second elastic energy force is calculated to be smaller than that in the case where the label is in the same region with respect to the grid point and the third grid point. apparatus.   The correction term calculation unit calculates the correction term based on a difference between a motion vector based on the second motion parameter and a motion vector of the second grid point. Image matching device.   The correction term calculation unit calculates, as the correction term, an output of a linear function having a difference between a motion vector based on the second motion parameter and a motion vector of the second grid point as an input. The image matching device according to claim 7.   When the difference between the motion vector based on the second motion parameter and the motion vector of the second lattice point is equal to or greater than a predetermined value, the correction term calculation unit calculates a nonlinear function using the difference as an input. The image matching apparatus according to claim 7, wherein an output is calculated as the correction term.   The motion parameter assignment unit defines the label by a cost function based on a difference value between pixel values of the second grid point and the first grid point, and sets the label having the minimum cost function as the label. The image matching apparatus according to claim 1, wherein the second motion parameter is assigned to the second grid point.   The motion parameter allocating unit defines the label by the cost function to which the continuity of the label is added, and uses the label having the minimum cost function as the second motion parameter as the second lattice. The image matching apparatus according to claim 10, wherein the image matching apparatus is assigned to a point.   The motion parameter allocating unit is configured such that the value of the cost function of the label determined to be optimal and the value of the cost function of the second label next to the label are within a first range. The image matching apparatus according to claim 10, wherein the second motion parameter is not assigned to each second grid point in an area corresponding to the label.   The motion parameter allocating unit further includes a case where the value of the cost function of the label determined to be optimal and the value of the cost function of the label optimal next to the label are within a second range. The image matching apparatus according to claim 12, wherein the second motion parameter is assigned to each second grid point in an area corresponding to the label. An occlusion detection unit that detects an occlusion region that is present in the target image but is not present in the reference image;
The motion parameter assignment unit assigns the second motion parameter only to the second lattice point not included in the occlusion region when the second lattice point is included in the occlusion region. The image matching apparatus according to claim 1.   The occlusion detection unit is formed from an area formed from grid points around the second grid point and grid points around the first grid point corresponding to the second grid point in the reference image. The image matching apparatus according to claim 14, wherein the occlusion region is detected based on a ratio to a certain area. An image matching method for obtaining a correspondence between a target image and a reference image,
There is an image correlation between each of the plurality of first grid points on the target image and each of the plurality of second grid points corresponding to the first grid points on the reference image. Based on the potential energy gradient based on the position of each second grid point and the position of each first grid point corresponding to the second grid point. An image energy force calculation step for calculating the image energy force received by the grid points;
An elastic energy force calculation step of calculating a first elastic energy force received from elastic energy between each of the second lattice points and another second lattice point adjacent to the second lattice point;
A friction force calculating step for calculating a friction force acting on each of the second grid points;
A first solution processing step of performing a solution process by numerical analysis of the first equation of motion related to each of the second lattice points by the image energy force, the first elastic energy force, and the friction force;
From the result of the solution processing by the first solution processing step, a motion clustering step of dividing the reference image into a plurality of regions and assigning a label for identifying the region for each divided region;
A motion parameter estimation step of estimating, by numerical analysis, a first motion parameter that is a parameter that defines a motion model indicating a mapping relationship between the second lattice point and the first lattice point in the region to which the label is assigned; ,
A motion parameter correcting step of correcting the first motion parameter so that the potential energy between the point on the target image and the point on the reference image determined by the estimated first motion parameter is minimized; ,
From the modified first motion parameter, a second motion parameter that is the first motion parameter optimum for each second grid point in the region corresponding to the label is determined and determined. A motion parameter assigning step for assigning the second motion parameter to each second grid point in the region corresponding to the label;
A correction term calculation step for calculating a correction term for correcting the first equation of motion based on the assigned second motion parameter and the result of the solution processing by the first solution processing step;
A second solution processing step for solving the second motion equation relating to each of the second lattice points obtained by adding the correction term to the first motion equation by numerical analysis; and a solution by the second solution processing step A mapping processing step for obtaining a correspondence relationship between the target image and the reference image from a processing result;
An image matching method comprising: An image matching program for obtaining a correspondence between a target image and a reference image,
There is an image correlation between each of the plurality of first grid points on the target image and each of the plurality of second grid points corresponding to the first grid points on the reference image. Based on the potential energy gradient based on the position of each second grid point and the position of each first grid point corresponding to the second grid point. An image energy force calculation step for calculating the image energy force received by the grid points;
An elastic energy force calculation step of calculating a first elastic energy force received from elastic energy between each of the second lattice points and another second lattice point adjacent to the second lattice point;
A friction force calculating step for calculating a friction force acting on each of the second grid points;
A first solution processing step of performing a solution process by numerical analysis of the first equation of motion related to each of the second lattice points by the image energy force, the first elastic energy force, and the friction force;
From the result of the solution processing by the first solution processing step, a motion clustering step of dividing the reference image into a plurality of regions and assigning a label for identifying the region for each divided region;
A motion parameter estimation step of estimating, by numerical analysis, a first motion parameter that is a parameter that defines a motion model indicating a mapping relationship between the second lattice point and the first lattice point in the region to which the label is assigned; ,
A motion parameter correcting step of correcting the first motion parameter so that the potential energy between the point on the target image and the point on the reference image determined by the estimated first motion parameter is minimized; ,
Determining, from the modified first motion parameter, the second motion parameter that is the first motion parameter optimal for each second grid point in the region corresponding to the label; Assigning said second motion parameter to each second grid point in the region corresponding to said label;
A correction term calculation step for calculating a correction term for correcting the first equation of motion based on the assigned second motion parameter and the result of the solution processing by the first solution processing step;
A second solution processing step for solving the second motion equation relating to each of the second lattice points obtained by adding the correction term to the first motion equation by numerical analysis; and a solution by the second solution processing step A mapping processing step for obtaining a correspondence relationship between the target image and the reference image from a processing result;
Matching program that makes computer execute.

Images