KR102031620B1 - Image Matching Method - Google Patents

Abstract

An image registration method is disclosed. In an image matching method according to an embodiment of the present invention, in a method of image matching by an image processing apparatus, defining a second image after a predetermined time elapses with respect to a first image, warping annealing is performed on the second image. And applying a first image and a second image applied to the warping annealing to calculate a velocity field for minimizing an energy function consisting of a data fidelity function and a normalization function.

Description

Image Matching Method

BACKGROUND OF THE INVENTION 1. Field of the Invention The present invention relates to an image registration method, and more particularly, to an image registration method capable of obtaining a velocity field describing motion of an image by minimizing an energy function composed of data fidelity and normalization.

In general, image segmentation techniques include medical imaging, such as tumor and disease detection, vascular volume measurement, computer-assisted surgery, diagnostics, treatment planning, and study of anatomical structures, and satellite, road, forest, etc. It is widely used in various fields, such as geography such as object detection, biometrics such as face and fingerprint, forensic science, forensic science, precision measurement, and precision repair of bridges / buildings.

Normalization, on the other hand, is generally imposed to reduce the allowable space of the solution in some image analysis tasks formulated as inverse problems. When applying normalization to motion estimation, the normalization parameters were inferred from the observed data in such a way that the combined probability of the gradient field and the velocity field is maximum. Another global approach uses bidirectional filtering and integrated noise estimation to normalize the estimated motion. As a weighting factor for shaping image segmentation, there are many spatial adaptive shaping techniques incorporating image gradients in the form of edge indicator functions. Local variations in image intensity within a fixed size window were also used to modulate normalization, and normalization parameters were chosen based on variance.

Non-global normalization has also been proposed for optical flow, image segmentation, and image reconstruction based on static image statistics. Unlike static normalization, there is a dynamic adaptation method that estimates normalization parameters through a dynamic system where spatial normalization is applied globally. Another method was developed based on Morozov's principle of inconsistency, with the residuals bounded by the estimated noise. Anisotropic structural tensors were used based on either total variation or generalized total variation. Most adaptive normalization algorithms consider constant spatial statistics during optimization iterations regardless of the residuals.

However, this method is unstable when the computational cost is high and the noise variance is small. In addition, the multi-label model suffers from tasks involving user interactions such as bounding boxes, outlines, graffiti, or points, and is difficult due to inaccurate or redundant slices when using a large number of labels to force such user interaction. There was this.

Therefore, there is a demand for developing a technology capable of segmenting images without user interaction or prior knowledge.

In this regard, there is a Korean Patent Publication No. 10-2017-0081969 entitled "Texture Image Processing Apparatus and Method".

The technical problem to be solved by the present invention is to provide an image registration method that can match the image without user interaction or prior knowledge.

The technical problems of the present invention are not limited to the technical problems mentioned above, and other technical problems not mentioned will be clearly understood by those skilled in the art from the following description.

According to an aspect of the present invention, there is provided a method of matching an image, the method of matching an image by an image processing apparatus, the method comprising: defining a second image after a predetermined time elapses with respect to the first image; Applying a warping annealing to a second image, and calculating a velocity field that minimizes an energy function composed of a data fidelity function and a normalization function using the first image and a second image applied to the warping annealing.

Preferably, the second image f ₂ (x) is defined by the equation described below. [Equation]

Here, f ₁ (x) is the first image, u (x) is a very small deformation of the image area, Is a noise process.

Preferably, f ₁ (x + u (x)) may be defined by the equation described below.

[Equation]

Here, ∇f ₁ represents a spatial gradient of an image f _1.

Preferably, the second image to which the warping annealing is applied may be defined by the equation described below.

[Equation]

Where τ is a control parameter taking into account the degree of distortion between the forward and reverse directions.

Preferably, the data fidelity function may be defined by the equation described below.

[Equation]

Where λ (x) is the relative weight function ego, Ρ (u) is the residual, β is a control parameter related to the change of the residual ρ (u) greater than '0', and α is 0 <α <1, which controls the degree of leanness of the relative weighting function λ. Constant parameters, Is , Is the hoover loss function.

Preferably, the normalization function may be defined by the equation described below.

[Equation]

here, , silver , The Huber function ( ), U (x) = (u ₁ (x), u ₂ (x)) is a component of the velocity field.

Preferably, the energy function can be defined by the equation described below.

[Equation]

Where θ is a scalar magnification factor with θ> 0, w = (w ₁ , w ₂ ) is a double variable for equivalent constraint u = v.

Advantageously, the energy function can be simplified to an unconstrained enlarged Lagrangian, such as the equation described below.

[Equation]

Where r and z = (z ₁ , z ₂ ) are auxiliary variables to be minimized.

Advantageously, the final energy function, simplified to Lagrangian, can be optimized using a gradient based alternating directional number technique (ADMM).

According to the present invention, the balance between data fidelity and normalization can be dynamically determined by using a relative weighting function based on the residual.

In addition, images can be matched without user interaction or prior knowledge.

Effects of the present invention are not limited to the above-mentioned effects, and other effects not mentioned will be clearly understood by those skilled in the art from the following description.

1 is a flowchart illustrating an image registration method according to an embodiment of the present invention.
2 is an exemplary diagram for describing image registration according to an embodiment of the present invention.
3 is a graph for explaining the annealing effect in the warfil according to an embodiment of the present invention.
4 is a diagram visually comparing a motion estimation method of a Middlebury data set using another algorithm according to an embodiment of the present invention.

As the invention allows for various changes and numerous embodiments, particular embodiments will be illustrated in the drawings and described in detail in the written description. However, this is not intended to limit the present invention to specific embodiments, it should be understood to include all modifications, equivalents, and substitutes included in the spirit and scope of the present invention. In describing the drawings, similar reference numerals are used for similar elements.

Terms such as first, second, A, and B may be used to describe various components, but the components should not be limited by the terms. The terms are used only for the purpose of distinguishing one component from another. For example, without departing from the scope of the present invention, the first component may be referred to as the second component, and similarly, the second component may also be referred to as the first component. The term and / or includes a combination of a plurality of related items or any item of a plurality of related items.

When a component is referred to as being "connected" or "connected" to another component, it may be directly connected to or connected to that other component, but it may be understood that other components may be present in between. Should be. On the other hand, when a component is said to be "directly connected" or "directly connected" to another component, it should be understood that there is no other component in between.

The terminology used herein is for the purpose of describing particular example embodiments only and is not intended to be limiting of the present invention. Singular expressions include plural expressions unless the context clearly indicates otherwise. In this application, the terms "comprise" or "have" are intended to indicate that there is a feature, number, step, operation, component, part, or combination thereof described in the specification, and one or more other features. It is to be understood that the present invention does not exclude the possibility of the presence or the addition of numbers, steps, operations, components, components, or a combination thereof.

Unless defined otherwise, all terms used herein, including technical or scientific terms, have the same meaning as commonly understood by one of ordinary skill in the art. Terms such as those defined in the commonly used dictionaries should be construed as having meanings consistent with the meanings in the context of the related art and shall not be construed in ideal or excessively formal meanings unless expressly defined in this application. Do not.

Hereinafter, exemplary embodiments of the present invention will be described in detail with reference to the accompanying drawings.

1 is a flowchart illustrating an image registration method according to an embodiment of the present invention, and FIG. 2 is an exemplary diagram for describing an image registration according to an embodiment of the present invention.

Referring to FIG. 1, the image processing apparatus defines a second image after a predetermined time elapses with respect to the first image (S110).

Space x∈Ω and time A series of images f (x; t) taken from will be described.

Video registration problem is a pair of images Wow Velocity field describing motion between Aim to calculate. The desired velocity field u is obtained by minimizing the energy function consisting of data fidelity and normalization.

In general, the energy function is a data fidelity function ( ) And normalization function ( Is defined using That is, the energy function is defined by Equation 1 described below.

[Equation 1]

In general, the tradeoff between data fidelity and normalization is assumed to be constant throughout the entire process of space (ie, the entire image area) and time (ie, generally iterative) optimization. However, neither is preferable. The spatially adaptable regularity will be described with reference to FIG. 2. A residual between the division result (b) in which the input image (a) and (a) is divided is displayed as (c), and the dispersion of the residual is displayed as (d). Looking at (c) and (d), the residuals or variances (gray levels: big light, small dark) are not constant in space. This may also apply to other imaging issues such as motion estimation (optical flow) and image reconstruction, where local dispersion of residuals often differs in space and optimization. Therefore, spatially adapted normalization beyond static image features is needed as studied in local contrast changes. The normalization of these tasks varies with space, but the transformations are constant through iteration because they are linked to image statistics. Thus, the present invention can implement spatial adaptive normalization using residuals that change during iteration.

That is, the present invention uses a normalization technique that can spatially adaptively and semi-automatically determine the weights between the data fidelity term and the normalization term. This method determines the relative weighting function without prior information based on the residual between the input data and the estimated solution.

That is, the relative weight function lambda is defined by equation (2) described below.

[Equation 2]

here, Where ρ (u) is the residual, β is a control parameter related to the change of the residual ρ (u) greater than '0', and α is 0 <α <1, which controls the degree of leanness of the relative weighting function λ It can mean a constant parameter. The relative weight function λ determines the balance between data fidelity and normalization. For example, as λ becomes larger, data fidelity is considered more important, so the boundaries remain clear, and as λ becomes smaller, more normalization is imposed, resulting in an over smoothing effect.

The relative weighting function λ between data fidelity and normalization is adaptively applied at each point x according to the residual ρ (u (x)) determined by the local suitability of the data for the model. Adaptation rules based on relative weighting functions (λ) have stronger equalization when the residuals are large, Normalization is strong when is small and weak and equivalent when residual is small. Is designed to be designed during the energy optimization process.

With positive Lagrange multiplier (α) The range of normalizes the weight 1-λ with the data fidelity ρ (u) so that normalization is imposed everywhere In the definition of a robust hoover loss function with a threshold parameter μ> 0 ( ).

The hoover loss function can be defined by equation (3) described below.

[Equation 3]

The use of the hoover loss function compared to L ₂ norm has the advantage of having a continuous derived value, unlike the undifferentiated L ₁ norm leading to stair artifacts while maintaining geometric features such as edges better. Further, Huber loss enables an efficient convex optimization due to the approximate equivalence of the operator L ₁ norm (proximal operator).

On the other hand, for data fidelity, an additional noise process ( Consider an image matching model based on the brightness coherence hypothesis.

Therefore, f ₂ (x) may be defined as in Equation 4 below.

[Equation 4]

Where u is an extremely small variant of the image area.

Previous velocity field solution to linearize f ₁ (x + u (x)) in equation (4) The first Taylor series expansion can be applied to. Then, f ₁ (x + u (x)) may be expressed as Equation 5 described below.

[Equation 5]

Here, ∇f ₁ represents a spatial gradient of the image f ₁ , and for the sake of simplicity, the superscript notation for transposition of the slope is omitted in ∇f ₁ . Then, linearization of the brightness coherence condition of Equations 2 and 3 leads to an image matching equation such as Equation 6 described below.

[Equation 6]

Here, f _t = f ₂ -f ₁ represents the time derivative of f.

Noise process ( Hoover Loss Function with Threshold μ> 0 Assuming a binomial distribution leading to, the data fidelity can be defined by equation (7) described below.

[Equation 7]

Here, λ (x) is determined based on the residual ρ (u) as a relative weighting function.

ρ (u) may be defined by Equation 8 described below.

[Equation 8]

Huber loss function for normalization Use the standard smoothness term with the threshold η> 0.

Then, the normalization function may be defined by Equation 9 described below.

[Equation 9]

here, May be defined by Equation 10 described below.

[Equation 10]

here, silver , The Huber function ( , U (x) = (u ₁ (x), u ₂ (x)) is a component of the velocity field. Normalization is required in areas where caliber problems are obvious (e.g., u is not unique because it is not unique in homogeneous areas) or in occlusal areas where u is not defined. However, normalization does not affect solutions where u is well defined.

On the other hand, to better cope with large displacements, an annealing method between forward and reverse warping is used in the image registration calculation.

That is, when calculating u, the control parameter Consider forward and backward deformation of the furnace domain. Then, f ₂ (x) may be defined as in Equation 11 described below.

[Equation 11]

Here, τ considers the degree of distortion between the forward and reverse directions. We introduce a simple annealing process for the control parameter τ whose value gradually changes from 0.5 to 1 during the optimization process. Considering the annealing in the twisting direction, the data fidelity is modified by Equation 12 described below.

[Equation 12]

Is defined by Equation 13 described below.

[Equation 13]

Where λ is the residual Is determined based on, and the initial velocity field u ₀ is omitted for ease of representation.

The energy function is then expressed as Equation 14 described below.

 [Equation 14]

Here, the initial value is τ = 0.5, which represents a symmetrical form of energy, and τ → 1 increases sequentially with the annealing process. A simple example of the annealing process is based on optimization iterations.

The energy function for the image registration of equation (14) is minimized for the velocity field u. The intermediate solution of u is used repeatedly as the initial solution u _o , and image warping is applied accordingly.

When the partition variable u = v and the variable v = (v ₁ , v ₂ ) are applied to the energy function of Equation 14, Equation 14 may be defined as Equation 15 described below.

 [Equation 15]

Where θ is a scalar magnification factor with θ> 0, w = (w ₁ , w ₂ ) is a double variable for the equivalent constraint u = v. Data Fidelity of Equation (34) And normalization of equation (31) Is a non-smooth function Moreau-Yosida can be optimized efficiently by normalization. That is, the data fidelity of Equation 13 And normalization of equation (10) Is a normalized form And Is replaced by.

Then, Equation 13 is defined as Equation 16 described below, and Equation 10 is defined as Equation 17 described below.

[Equation 16]

[Equation 17]

Where r and z = (z ₁ , z ₂ ) are auxiliary variables to be minimized.

The energy function then becomes Lagrangian (Eq. ).

Equation 18

Energy function By minimizing, we can obtain the desired velocity field u = (u ₁ , u ₂ ). To optimize the energy function, we can use gradient based alternating number technique (ADMM) algorithm. For the control parameter τ for warping annealing, 0.5 to 1 is given the step size at each iteration Divided by.

The ADDM algorithm is shown in Table 1.

TABLE 1

According to the ADDM algorithm shown in Table 1 Proceed to minimize.

Hereinafter, the robustness and effects of the normalization technique for image registration according to the present invention will be described using experimental results. In the experiment, qualitative and quantitative evaluations were performed based on the Middlebury optical flow data set. The mean endpoint error (AEE) and mean angle error (AAE) are used for quantitative evaluation.

3 is a graph for explaining the annealing effect in the warfil according to an embodiment of the present invention. (a) is a graph comparing the energy of motion estimation with or without annealing, and (b) is a graph comparing the average endpoint error of motion estimation with or without annealing. to be.

The effect of the annealing method on the degree of distortion between the forward and backward directions is described. To emphasize the role of the annealing parameters, we apply the image matching algorithm according to the present invention without using adaptive normalization and use the image pairs with the largest discrepancy (Grove3 sequence) on average in the data set. The result is performed by a motion estimation algorithm with or without warping annealing parameter τ. It can be seen that the algorithm using annealing in warping provides faster convergence and better accuracy compared to the baseline without the warping annealing scheme. Changes in the warping annealing parameters in the optimization iteration result in changes in the energy, and as a result It causes an energy curve fluctuation such as (a) with = 0.005.

4 is a diagram visually comparing a motion estimation method of a Middlebury data set using another algorithm according to an embodiment of the present invention.

The algorithm according to the invention is compared with a Huber variant with the traditional Horn-Schunck model (HS) [51], the TV-L1 model (TV) [63] and the L1 data fidelity (HTV) [24]. (a) the input image, (b) the ground truth, (c) the HS, (d) the TV, (e) the HTV, and (f) the method of the present invention. This visual comparison indicates that the method according to the invention is more accurate than the other algorithms quantitatively evaluated in Table 2. Where AEE and AAE are calculated for each case.

TABLE 2

Occlusion is not specially considered in the calculation of optical flows to emphasize the role of adaptive normalization, which implicitly deals with occlusions that result in higher residuals due to mismatches. The parameters of each method are optimally chosen in terms of errors, μ = 0.01, η = 0.3, α = 0.01, β = 10, = 0.005, θ = 0.1 was used.

So far I looked at the center of the preferred embodiment for the present invention. Those skilled in the art will appreciate that the present invention can be implemented in a modified form without departing from the essential features of the present invention. Therefore, the disclosed embodiments should be considered in descriptive sense only and not for purposes of limitation. The scope of the present invention is shown in the claims rather than the foregoing description, and all differences within the scope will be construed as being included in the present invention.

Claims (9)

In the image processing apparatus to match the image,
Defining a second image after a predetermined time elapses with respect to the first image;
Applying annealing in warping to the second image; And
Calculating a velocity field for minimizing an energy function consisting of a data fidelity function and a normalization function using the first image and a second image applied to the warping annealing;
The second image ( f ₂ (x)) is an image registration method, characterized in that defined by the equation
[Equation]

Here, f ₁ (x) is the first image, u (x) is a very small deformation of the image area, Is a noise process. delete The method of claim 1,
F ₁ (x + u (x)) is defined by the equation described below.
[Equation]

Here, ∇f ₁ represents a spatial gradient of an image f _1. The method of claim 1,
The second image to which the warping annealing is applied is defined by the equation described below.
[Equation]

Where τ is a control parameter taking into account the degree of distortion between the forward and reverse directions. The method of claim 1,
The data fidelity function is defined by the equation described below
[Equation]

Where λ (x) is the relative weight function ego, Ρ (u) is the residual, β is a control parameter related to the change of the residual ρ (u) greater than '0', and α is 0 <α <1, which controls the degree of leanness of the relative weighting function λ. Constant parameters, Is , Is the hoover loss function. The method of claim 1,
The normalization function is an image matching method, characterized in that defined by the equation described below
[Equation]

here, , silver , The Huber function ( ), U (x) = (u ₁ (x), u ₂ (x)) is a component of the velocity field. The method of claim 1,
The energy function is defined by the equation described below
[Equation]

Where θ is a scalar magnification factor with θ> 0, w = (w ₁ , w ₂ ) is a double variable for equivalent constraint u = v. The method of claim 7, wherein
And said energy function is simplified to an unconstrained enlarged Lagrangian, as shown below.
[Equation]

Where r and z = (z ₁ , z ₂ ) are auxiliary variables to be minimized. The method of claim 8,
The final energy function, simplified to Lagrangian, is optimized using a gradient based alternating direction number technique (ADMM).