CN105046657A - Image stretching distortion adaptive correction method - Google Patents

Abstract

The invention provides an image stretching distortion adaptive correction method comprising: dividing an image into uniform square grids and establishing an input-output coordinate; determining a person area and a background area by using a face detection algorithm, establishing homographic constraint for each grid in order to correct the stretching distortion appearing in the person area, establishing homographic compatible constraint for adjacent grids in order to guarantee the continuity of the grids; determining a linear area by using a linear detection algorithm and establishing linear maintaining constraint for the linear area in order to protect a corrected line against bending; guaranteeing the uniform continuity of the whole image by means of global coordinate smoothness and uniformity constraint; using an energy function to express the constraint, adding weights and adding the energy equations, computing a coordinate value with the minimum energy by using a least square method so as to obtain a corrected image coordinate value; performing rendering by using bilinear mapping to obtain a final corrected image.

Description

A kind of image stretch distortion self-adapting correction method

Technical field

The present invention relates to a kind of digital image processing techniques, particularly relate to a kind of image stretch distortion self-adapting correction method.

Background technology

Camera wide-angle lens head can obtain the image compared with Large visual angle, but also easily introduces pattern distortion.Wherein a kind of distortion carries out imaging when pin-hole imaging model followed by wide-angle lens, and when having object near camera in imaging scene and be in field of view edge, in image after shooting, this object can seriously be stretched due to perspective effect, causes the distortion of this geometry of objects.Particularly this object be face or personage time, because the distortion of human eye to face and character image is comparatively responsive, this stretch distortion is more obvious.Therefore, needing the correct image to there is object stretch distortion, making image meet human eye vision.

But the straight line in image scene can be bent because correcting the stretching of object, therefore, need making other scenes or object remain unaffected at correction stretch distortion simultaneously.

In addition, panoramic picture is also a kind of wide angle picture in essence.When panoramic picture adopts linear perspective projection, this image also there will be stretch distortion.

In existing technology, there is the overall projection correction method by changing projection centre, having by image background to be separated the method corrected respectively again with figure picture, having projection correction's method of the maintenance picture material of optimization.

Summary of the invention

The object of the invention is to for the deficiencies in the prior art, provide a kind of image stretch to distort self-adapting correction method, the object especially stretch distortion that occurs of personage in correcting image while keeping background straight line in image straight.

The object of the invention is to be achieved through the following technical solutions: a kind of image stretch distortion self-adapting correction method, said method comprising the steps of:

Steps A: input the image with object stretch distortion.

Step B: set up grid and coordinate, is divided into square uniform grid by the image with object stretch distortion, and grid intersecting point coordinate is designated as x _i,j=(x _i,j, y _i,j), i, j are grid intersection point sequence number; If grid intersecting point coordinate x _i,jcoordinate after correction is u _i,j=(u _i,j, v _i,j).

Step C: structure u _i,jbound energy function.

Step C1: structure homography bound energy function.Calculate the angle rotation of people's object area being needed to field of view center to rotate computing formula is:

Wherein L is the horizontal-shift distance of people's object area center to field of view center, and f is focal length, L and f is all converted in units of pixel count.Then the homography conversion of people's object area can be written as:

$u_{i,j}=f\left(\frac{U}{W}\right)+L$

$v_{i,j}=f\left(\frac{V}{W}\right)$

The homography conversion described by above formula is designated as H (x _i,j), then can define the energy function of personage's area grid point homography constraint:

$E_{h1}=\underset{x_{i,j}\∈P_1}{\Σ}|u_{i,j}-H\left(x_{i,j}\right){|}^2$

Wherein, P ₁the set of behaviour object area net point.

The homography bound energy function of background area net point is defined as:

$E_{h2}=\underset{x_{i,j}\∈P_2}{\Σ}|u_{i,j}-x_{i,j}{|}^2$

Wherein, P ₂for the set of background area net point.

Total homography bound energy function is

$E_h=E_{h1}+E_{h2}=\underset{x_{i,j}\∈P_1}{\Σ}|u_{i,j}-H\left(x_{i,j}\right){|}^2+\underset{x_{i,j}\∈P_2}{\Σ}|u_{i,j}-x_{i,j}{|}^2$

Step C2: structure homography compatibility restriction energy function.

Definition P _hcthe net point on total border between behaviour object area and background area, H ₁, H ₂be respectively the homography conversion of background area and people's object area.Then the energy function of homography compatibility restriction is:

$E_{hc}=\underset{x_{i,j}\∈P_{hc}}{\Σ}\left(|u_{i,j}-H_1\left(x_{i,j}\right){|}^2+|u_{i,j}-H_2\left(x_{i,j}\right){|}^2\right)$

Step C3: structure straight line keeps bound energy function.First the straight-line segment except personage in background area in line segment detection algorithm detected image is used.The straight line set detected is defined as L.

Have in the grid of straight line process at each and define a virtual point

$u_{i,j}^L={\mathrm{au}}_{i,j}+{\mathrm{bu}}_{i+1,j}+{\mathrm{cu}}_{i,j+1}+{\mathrm{du}}_{i+1,j+1}$

Wherein coefficient a, b, c, d are calculated through bilinear interpolation by the x-y coordinate figure corresponding with u-v coordinate.Computing formula is:

$a=\frac{\left(x_{i+1,j}-x_{i,j}^L\right)\left(y_{i,j+1}-y_{i,j}^L\right)}{\left(x_{i+1,j}-x_{i,j}\right)\left(y_{i,j+1}-y_{i,j}\right)}$

$b=\frac{\left(x_{i,j}^L-x_{i,j}\right)\left(y_{i,j+1}-y_{i,j}^L\right)}{\left(x_{i+1,j}-x_{i,j}\right)\left(y_{i,j+1}-y_{i,j}\right)}$

$c=\frac{\left(x_{i,j}^L-x_{i,j}\right)\left(y_{i,j}^L-y_{i,j}\right)}{\left(x_{i+1,j}-x_{i,j}\right)\left(y_{i,j+1}-y_{i,j}\right)}$

$d=\frac{\left(x_{i+1,j}-x_{i,j}^L\right)\left(y_{i,j}^L-y_{i,j}\right)}{\left(x_{i+1,j}-x_{i,j}\right)\left(y_{i,j+1}-y_{i,j}\right)}$

Wherein, be taken as the mid point of line segment in gridblock.

For any straight line, define its starting point and terminal is respectively u _startand u _end.

Point u on calculated curve _i,jand this is to the distance D between subpoint on straight line, and D is constrained to zero.Can be obtained by geometric relationship, the calculating formula of this distance D is

D＝|(u _i,j-u _start) ^Tn|

$n=\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right]\frac{(u_{end}-u_{start})}{\left|\right|u_{end}-u_{start}\left|\right|}$

Wherein, n=(n ₁, n ₂) be the unit normal vector of straight line.Calculate distance D, then the energy function of definable straight line maintenance constraint is

$E_l=\underset{x_{i,j}\∈L}{\Σ}{\left({\left(u_{i,j}-u_{start}\right)}^{T}n\right)}^2=\underset{x_{i,j}\∈L}{\Σ}{\left(n_1\left(u_{i,j}-u_{start}\right)+n_2\left(v_{i,j}-v_{start}\right)\right)}^2;$

Step C4: structure world coordinates smoothness constraint energy function.The Jacobi matrix of definition u-v coordinate is

$J=\frac{\∂u}{\∂x}=\left[\begin{array}{cc}\frac{\∂u}{\∂x}& \frac{\∂u}{\∂y}\\ \frac{\∂v}{\∂x}& \frac{\∂v}{\∂y}\end{array}\right];$

When every in Jacobi matrix meets Cauchy-Riemann equations:

$\begin{array}{c}\frac{\∂u}{\∂x}=\frac{\∂v}{\∂y}\\ \frac{\∂u}{\∂y}=-\frac{\∂v}{\∂x}\end{array}$ Time,

The continuously smooth of coordinate can be realized.Above formula discretize is obtained:

u _i+1,j-u _i,j＝v _i,j+1-v _i,j

u _i,j+1-u _i,j＝-(v _i+1,j-v _i,j)；

Can coordinates computed smoothness constraint energy function by above formula:

E _s＝Σ((u _i+1,j-u _i,j)-(v _i,j+1-v _i,j)) ²+Σ((v _i+1,j-v _i,j)+(u _i,j+1-u _i,j)) ²；

Step C5: structure world coordinates uniform restraint energy function.By Laplace operator discretize, the energy function of definable world coordinates uniform restraint is

$E_u=\mathrm{\Σ}\left|\right|\left[\begin{array}{c}{u}_{i+1,j}\text{+}{u}_{i\text{-}1,j}\text{+}{u}_{i,j\text{+1}}\text{+}{u}_{i,j\text{-1}}{\text{-4u}}_{i,j}\\ v_{i+1,j}{\text{+v}}_{i\text{-}1,j}{\text{+v}}_{i,j\text{+1}}{\text{+v}}_{i,j\text{-1}}{\text{-4v}}_{i,j}\end{array}\right]|{|}^{\text{2}};$

Step D: linear least-squares solves.For each energy function adds weight coefficient and each energy function addition is obtained total energy function be:

$E=w_h^2{E}_h+w_{hc}^2{E}_{hc}+w_l^2{E}_l+w_s^2{E}_s+w_u^2{E}_u$

Make E=0, then E _h, E _hc, E _l, E _s, E _uall be required to be 0, can system of linear equations be constructed thus

Ax＝b

Wherein A is constraint matrix, comprises the u in bound energy function _i,jeach front coefficient.X is for comprising unknown quantity u to be asked _i,jand v _i,jcolumn vector, b comprises the constant term of separate equation.

Solve x with linear least square, solution formula is

x＝(A ^TA) ^-1A ^Tb

Step e: solve net point coordinate u _i,jand v _i,jafter, the coordinate adopting bilinear map to calculate except net point after arbitrfary point (x, y) correction is

(u,v)＝(1-x)(1-y)p1+x(1-y)p2+(1-x)yp3+xyp4。

Wherein, p1, p2, p3, p4 are respectively 4 net point coordinates of point (x, y) place grid, and p1 is upper left corner net point coordinate, and p2 is upper right corner net point coordinate, and p3 is lower left corner net point coordinate, and p4 is lower right corner net point coordinate.

Step F: after calculating calibration coordinate (u, v).By from input coordinate value (x, y) remapping and interpolation to the coordinate figure (u, v) calculated, obtain correcting image.

Method disclosed by the invention compares existing bearing calibration following feature: (1) can correct the stretch distortion of field of view edge place object, can keep again background straight line not by bending; (2) correct result images content consistent continuously, occur without the cutting imagination; (3) after correcting, the visual field loss of image is less; (4) bearing calibration completes automatically, without the need to manual intervention.

Accompanying drawing explanation

Fig. 1 is the process flow diagram of a kind of image stretch distortion correction method of the present invention;

Fig. 2 is the example image with stretch distortion described in the invention;

Fig. 3 is the foundation of grid and the definition schematic diagram of coordinate;

Fig. 4 is the division schematic diagram of image-region;

Fig. 5 is the definition schematic diagram of virtual point;

Fig. 6 is the calculating schematic diagram of distance D;

Fig. 7 is bilinear map schematic diagram;

Fig. 8 is the example image after the inventive method corrects;

Fig. 9 is that the example image after the inventive method corrects cuts out the image after into rectangle.

Embodiment

Image stretch distortion self-adapting correction method of the present invention comprises the image accepted with object stretch distortion; With uniform grid image divided and set up input and output coordinate system; Input picture is detected face location and size by Face datection algorithm, and determines rectangular area residing for personage and remaining background area by this information; Be respectively regional and set up homography constraint, to correct the stretch distortion occurred in people's object area, and set up homography compatibility restriction for adjacent area, the continuity of guarantee regional; Determine linearity region by line detection algorithm, set up straight line for linearity region and keep constraint, to ensure that bending does not appear in the straight line after correcting; Set up the level and smooth and homogeneity constraint guarantee entire image congruous continuity of world coordinates; Above-mentioned constraint is expressed as energy function, and adds weight, each energy equation is added; Energy function is minimized to construct system of linear equations, uses least square method to calculate this system of linear equations, obtain the image coordinate value after correcting; Use bilinear texture to map to carry out playing up obtaining final image.

For understanding implementation procedure of the present invention better, elaborate a kind of image stretch distortion self-adapting correction method of the present invention below in conjunction with accompanying drawing.The process flow diagram of concrete implementation step is see Fig. 1, and applying method of the present invention, to realize the concrete steps of image stretch distortion correction as follows:

Steps A: input the image with object stretch distortion.Fig. 2 is an example of the image with personage's stretch distortion, and distorted position double-head arrow is indicated.

Step B: set up grid and coordinate.For Fig. 2, as shown in Figure 3, image is divided into square uniform grid, grid intersecting point coordinate is designated as x _i,j=(x _i,j, y _i,j), i, j are grid intersection point sequence number.After scalloping, coordinate to be solved is designated as u _i,j=(u _i,j, v _i,j).

Step C: structure constraint energy function.

Step C1: structure homography bound energy function.Occur that the reason of stretch distortion is that object is in image f iotaeld-of-view edge, when rotation camera make it just to object and make object be in the shooting of visual field central authorities time, then object tensionless winkler foundation distortion.Therefore, the visual angle being changed image local area by homography conversion (perspective transform) is needed, to correct the stretch distortion of object in this region.First, use image Face datection algorithm to detect face location and size, determine rectangular area residing for personage and remaining background area by this information, the people's object area determined for Fig. 2 and background area are as shown in Figure 4.Calculate the angle rotation of people's object area being needed to field of view center to rotate computing formula is:

Wherein L is the horizontal-shift distance of people's object area center to field of view center, and f is focal length, L and f is all converted in units of pixel count.Then the homography conversion of people's object area can be written as:

$u_{i,j}=f\left(\frac{U}{W}\right)+L$

$v_{i,j}=f\left(\frac{V}{W}\right)$

The homography conversion described by above formula is designated as H (x _i,j), then can define the energy function of personage's area grid point homography constraint:

$E_{h1}=\underset{x_{i,j}\∈P_1}{\Σ}|u_{i,j}-H\left(x_{i,j}\right){|}^2$

Wherein P ₁the set of behaviour object area net point.Background area does not need homography to convert, and therefore the homography bound energy function of this regional network lattice point is defined as:

$E_{h2}=\underset{x_{i,j}\∈P_2}{\Σ}|u_{i,j}-x_{i,j}{|}^2$

Wherein P ₂for the set of background area net point.Total homography bound energy function is

$E_h=E_{h1}+E_{h2}=\underset{x_{i,j}\∈P_1}{\Σ}|u_{i,j}-H\left(x_{i,j}\right){|}^2+\underset{x_{i,j}\∈P_2}{\Σ}|u_{i,j}-x_{i,j}{|}^2$

Step C2: structure homography compatibility restriction energy function.There is public limit people's object area and background area.In order to prevent regional from causing region joining part to occur inconsistent phenomenon after homography conversion, homography compatibility restriction must be added at the adjacent part in region.Definition P _hcfor the net point on interregional total border, H ₁, H ₂be respectively the homography conversion in figure and ground region.Then the energy function of homography compatibility restriction is:

$E_{hc}=\underset{x_{i,j}\∈P_{hc}}{\Σ}\left(\text{|}{u}_{i,j}-H_1\left(x_{i,j}\right){\text{|}}^{\text{2}}\text{+}|u_{i,j}-H_2\left(x_{i,j}\right){|}^2\right)$

Described homography conversion H ₁for:

$u_{i,j}=f\left(\frac{U}{W}\right)+L$

$v_{i,j}=f\left(\frac{V}{W}\right)$

Described homography conversion H ₂it is the unit matrix of 3 × 3.

Step C3: structure straight line keeps bound energy function.Ensure that image cathetus still keeps straight after calibrated, need to add straight line and keep constraint.First the straight-line segment of background area except personage in line segment detection algorithm detected image is used.The straight line set detected is defined as L.

Due to be not any straight line all through net point, therefore as shown in Figure 5, each have in the grid of straight line process define a virtual point

$u_{i,j}^L={\mathrm{au}}_{i,j}+{\mathrm{bu}}_{i+1,j}+{\mathrm{cu}}_{i,j+1}+{\mathrm{du}}_{i+1,j+1}$

Wherein coefficient a, b, c, d are calculated through bilinear interpolation by the x-y coordinate figure corresponding with u-v coordinate.Computing formula is:

$a=\frac{\left(x_{i+1,j}-x_{i,j}^L\right)\left(y_{i,j+1}-y_{i,j}^L\right)}{\left(x_{i+1,j}-x_{i,j}\right)\left(y_{i,j+1}-y_{i,j}\right)}$

$b=\frac{\left(x_{i,j}^L-x_{i,j}\right)\left(y_{i,j+1}-y_{i,j}^L\right)}{\left(x_{i+1,j}-x_{i,j}\right)\left(y_{i,j+1}-y_{i,j}\right)}$

$c=\frac{\left(x_{i,j}^L-x_{i,j}\right)\left(y_{i,j}^L-y_{i,j}\right)}{\left(x_{i+1,j}-x_{i,j}\right)\left(y_{i,j+1}-y_{i,j}\right)}$

$d=\frac{\left(x_{i+1,j}-x_{i,j}^L\right)\left(y_{i,j}^L-y_{i,j}\right)}{\left(x_{i+1,j}-x_{i,j}\right)\left(y_{i,j+1}-y_{i,j}\right)}$

Wherein, be taken as the mid point of line segment in gridblock.

For any straight line, define its starting point and terminal is respectively u _startand u _end.As shown in Figure 6, in order to by stretching for curve be straight line, need the some u on calculated curve _i,jand this is to the distance D between subpoint on straight line, and D is constrained to zero.Can be obtained by geometric relationship, the calculating formula of this distance D is

D＝|(u _i,j-u _start) ^Tn|

$n=\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right]\frac{\left(u_{end}-u_{start}\right)}{\left|\right|u_{end}-u_{start}\left|\right|}$

Wherein, n=(n ₁, n ₂) be the unit normal vector of straight line.Calculate distance D, then the energy function of definable straight line maintenance constraint is

$E_l=\underset{x_{i,j}\∈L}{\Σ}{\left({\left(u_{i,j}-u_{stan}\right)}^{T}n\right)}^2={\underset{x_{i,j}\∈L}{\Σ}\left(n_1\left(u_{i,j}-u_{start}\right)+n_2\left({\mathrm{\ν}}_{i,j}-v_{start}\right)\right)}^2$

Step C4: structure world coordinates smoothness constraint energy function.Enable picture material congruous continuity, need define world coordinates smoothness constraint, the Jacobi matrix of definition u-v coordinate is

$J=\frac{\∂u}{\∂x}=\left[\begin{array}{cc}\frac{\∂u}{\∂x}& \frac{\∂u}{\∂y}\\ \frac{\∂v}{\∂x}& \frac{\∂v}{\∂y}\end{array}\right]$

When every in Jacobi matrix meets Cauchy-Riemann equations:

$\frac{\∂u}{\∂x}=\frac{\∂v}{\∂y}$

$\frac{\∂u}{\∂y}=-\frac{\∂v}{\∂x}$

Time, the continuously smooth of coordinate can be realized.Because grid is discrete, above formula discretize is obtained

u _i+1,j-u _i,j＝v _i,j+1-v _i,j

u _i,j+1-u _i,j＝-(v _i+1,j-v _i,j)

Can coordinates computed smoothness constraint energy function by above formula:

E _s＝Σ((u _i+1,j-u _i,j)-(v _i,j+1-v _i,j)) ²+Σ((v _i+1,j-v _i,j)+(u _i,j+1-u _i,j)) ²

Step C5: structure world coordinates uniform restraint energy function.Only consider that coordinate smoothness constraint can make correction of a final proof image peripheral meshes more much larger than grid in the middle part of image.Therefore, be necessary to add world coordinates uniform restraint, the grid ratio inside and outside image is consistent.The solution of Laplace's equation has even character, therefore adopts Laplace operator here.By Laplace operator discretize, the energy function of definable coordinate uniform restraint is

$E_u=\mathrm{\Σ}\left|\right|\left[\begin{array}{c}{u}_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}-4u_{i,j}\\ v_{i+1,j}+v_{i-1,j}+v_{i,j+1}+v_{i,j-1}-{\mathrm{vu}}_{i,j}\end{array}\right]|{|}^2$

Step D: linear least-squares solves.For each energy function adds weight coefficient and each energy function addition is obtained total energy function be:

$E=w_h^2{E}_h+w_{hc}^2{E}_{hc}+w_l^2{E}_l+w_s^2{E}_s+w_u^2{E}_u$

Make energy function minimum, need E=0 be made.E=0, then E _h, E _hc, E _l, E _s, E _uall be required to be 0, can system of linear equations be constructed thus

Ax＝b

Wherein A is constraint matrix, comprises the u in bound energy function _i,jeach front coefficient. _xfor comprising unknown quantity u to be asked _i,jand v _i,jcolumn vector, b comprises the constant term of separate equation.Because this system of linear equations is overdetermination, useable linear least square method pair _xsolve, solution formula is

x＝(A ^TA) ^-1A ^Tb

Step e: solve net point coordinate u _i,jand v _i,jafter, bilinear map need be adopted to calculate all the other coordinates except net point.As shown in Figure 7, known net point p1, p2, p3, p4 coordinate, in grid, the calculating formula of arbitrfary point (u, v) is

(u,v)＝(1-x)(1-y)p1+x(1-y)p2+(1-x)yp3+xyp4

Step F: after calculating calibration coordinate (u, v).Through from input coordinate value remapping and interpolation to the coordinate figure calculated, correction of a final proof image can be obtained, as shown in Figure 8.Correcting image irregularity boundary, need cut out as exporting after rectangle, as shown in Figure 9.

Claims (1)

1. an image stretch distortion self-adapting correction method, is characterized in that, said method comprising the steps of:

Steps A: input the image with object stretch distortion.

Step B: set up grid and coordinate, is divided into square uniform grid by the image with object stretch distortion, and grid intersecting point coordinate is designated as x _i,j=(x _i,j, y _i,j), i, j are grid intersection point sequence number; If grid intersecting point coordinate x _i,jcoordinate after correction is u _i,j=(u _i,j, v _i,j).

Step C: structure u _i,jbound energy function.

Step C1: structure homography bound energy function.Calculate the angle rotation of people's object area being needed to field of view center to rotate computing formula is:

Wherein L is the horizontal-shift distance of people's object area center to field of view center, and f is focal length, L and f is all converted in units of pixel count.Then the homography conversion of people's object area can be written as:

$u_{i,j}=f\left(\frac{U}{W}\right)+L$

$v_{i,j}=f\left(\frac{V}{W}\right)$

The homography conversion described by above formula is designated as H (x _i,j), then can define the energy function of personage's area grid point homography constraint:

$E_{h1}=\underset{x_{i,j}\∈P_1}{\Σ}{|u_{i,j}-H\left(x_{i,j}\right)|}^2$

Wherein, P ₁the set of behaviour object area net point.

The homography bound energy function of background area net point is defined as:

$E_{h2}=\underset{x_{i,j}\∈P_2}{\Σ}|u_{i,j}-x_{i,j}{|}^2$

Wherein, P ₂for the set of background area net point.

Total homography bound energy function is

$E_h=E_{h1}+E_{h2}=\underset{x_{i,j}\∈P_1}{\Σ}|u_{i,j}-H\left(x_{i,j}\right){|}^2+\underset{x_{i,j}\∈P_2}{\Σ}|u_{i,j}-x_{i,j}{|}^2$

Step C2: structure homography compatibility restriction energy function.

Definition P _hcthe net point on total border between behaviour object area and background area, H ₁, H ₂be respectively the homography conversion of background area and people's object area.Then the energy function of homography compatibility restriction is:

$E_{hc}=\underset{x_{i,j}\∈P_{hc}}{\mathrm{\Σ}}\left({\left|u_{i,j}-H_1\left(x_{i,j}\right)\right|}^2+{|u_{i,j}-H_2\left(x_{i,j}\right)|}^2\right)$

Step C3: structure straight line keeps bound energy function.First the straight-line segment except personage in background area in line segment detection algorithm detected image is used.The straight line set detected is defined as L.

Have in the grid of straight line process at each and define a virtual point

$u_{i,j}^L={\mathrm{au}}_{i,j}+{\mathrm{bu}}_{i+1,j}+{\mathrm{cu}}_{i,j+1}+{\mathrm{du}}_{i+1,j+1}$

Wherein coefficient a, b, c, d are calculated through bilinear interpolation by the x-y coordinate figure corresponding with u-v coordinate.Computing formula is:

$a=\frac{(x_{i+1,j}-x_{i,j}^L)(y_{i,j+1}-y_{i,j}^L)}{(x_{i+1,j}-x_{i,j})(y_{i,j+1}-y_{i,j})}$

$b=\frac{(x_{i,j}^L-x_{i,j})(y_{i,j+1}-y_{i,j}^L)}{(x_{i+1,j}-x_{i,j})(y_{i,j+1}-y_{i,j})}$

$c=\frac{(x_{i,j}^L-x_{i,j})(y_{i,j}^L-y_{i,j})}{(x_{i+1,j}-x_{i,j})(y_{i,j+1}-y_{i,j})}$

$d=\frac{(x_{i+1,j}-x_{i,j}^L)(y_{i,j}^L-y_{i,j})}{(x_{i+1,j}-x_{i,j})(y_{i,j+1}-y_{i,j})}$

Wherein, be taken as the mid point of line segment in gridblock.

For any straight line, define its starting point and terminal is respectively u _startand u _end.

Point u on calculated curve _i,jand this is to the distance D between subpoint on straight line, and D is constrained to zero.Can be obtained by geometric relationship, the calculating formula of this distance D is

D＝|(u _i,j-u _start) ^Tn|

$n=\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right]\frac{(u_{end}-u_{start})}{||u_{end}-u_{start}||}$

Wherein, n=(n ₁, n ₂) be the unit normal vector of straight line.Calculate distance D, then the energy function of definable straight line maintenance constraint is

$E_l=\underset{x_{i,j}\∈L}{\Σ}{\left({(u_{i,j}-u_{start})}^{T}n\right)}^2=\underset{x_{i,j}\∈L}{\mathrm{\Σ}}{\left(n_1\right(u_{i,j}-u_{start})+n_2(v_{i,j}-v_{start}))}^2;$

Step C4: structure world coordinates smoothness constraint energy function.The Jacobi matrix of definition u-v coordinate is

$J=\frac{\∂u}{\∂x}=\left[\begin{array}{cc}\frac{\∂u}{\∂x}& \frac{\∂u}{\∂y}\\ \frac{\∂v}{\∂x}& \frac{\∂v}{\∂y}\end{array}\right];$

When every in Jacobi matrix meets Cauchy-Riemann equations:

$\begin{array}{c}\frac{\∂u}{\∂x}=\frac{\∂v}{\∂y}\\ \frac{\∂u}{\∂y}=-\frac{\∂v}{\∂x}\end{array}$ Time,

The continuously smooth of coordinate can be realized.Above formula discretize is obtained:

$\begin{array}{c}{u}_{i+1,j}-u_{i,j}=v_{i,j+1}-v_{i,j}\\ u_{i,j+1}-u_{i,j}=-\left(v_{i+1,j}-v_{i,j}\right)\end{array};$

Can coordinates computed smoothness constraint energy function by above formula:

E _s＝Σ((u _i+1,j-u _i,j)-(v _i,j+1-v _i,j)) ²+Σ((v _i+1,j-v _i,j)+(u _i,j+1-u _i,j)) ²；

Step C5: structure world coordinates uniform restraint energy function.By Laplace operator discretize, the energy function of definable world coordinates uniform restraint is

$E_u=\mathrm{\Σ}{||\left[\begin{array}{c}{u}_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}-4u_{i,j}\\ v_{i+1,j}+v_{i-1,j}+v_{i,j+1}+v_{i,j-1}-4v_{i,j}\end{array}\right]||}^2;$

Step D: linear least-squares solves.For each energy function adds weight coefficient and each energy function addition is obtained total energy function be:

$E=w_h^2{E}_h+w_{hc}^2{E}_{hc}+w_l^2{E}_l+w_s^2{E}_s+w_u^2{E}_u$

Make E=0, then E _h, E _hc, E _l, E _s, E _uall be required to be 0, can system of linear equations be constructed thus

Ax＝b

Wherein A is constraint matrix, comprises the u in bound energy function _i,jeach front coefficient.X is for comprising unknown quantity u to be asked _i,jand v _i,jcolumn vector, b comprises the constant term of separate equation.

Solve x with linear least square, solution formula is

x＝(A ^TA) ^-1A ^Tb

Step e: solve net point coordinate u _i,jand v _i,jafter, the coordinate adopting bilinear map to calculate except net point after arbitrfary point (x, y) correction is

(u,v)＝(1-x)(1-y)p1+x(1-y)p2+(1-x)yp3+xyp4。

Wherein, p1, p2, p3, p4 are respectively 4 net point coordinates of point (x, y) place grid, and p1 is upper left corner net point coordinate, and p2 is upper right corner net point coordinate, and p3 is lower left corner net point coordinate, and p4 is lower right corner net point coordinate.

Step F: after calculating calibration coordinate (u, v).By from input coordinate value (x, y) remapping and interpolation to the coordinate figure (u, v) calculated, obtain correcting image.