CN101197043B - Shape-maintaining fitting algorithm in image processing - Google Patents

Abstract

The invention proposes a shape-preserving fitting algorithm in image processing. This method first inverts the approximate value of special points on the original surface according to the discrete image, then constructs a bicubic Coons interpolation surface on the area surrounded by four adjacent data points for the new data, in which the first-order partial values of the four corner points The derivative is estimated by the shape-preserving interpolation method, and the second-order partial derivatives of the four corner points are obtained by averaging the partial derivatives of its adjacent pixels. On this surface, the Coons surface is resampled according to the scaling ratio, and we use area sampling to obtain the gray value of the resampling point. By increasing the sampling density, more pixels are obtained, so as to achieve the purpose of image enlargement. The new method can effectively overcome the mosaic phenomenon and preserve clearer image details.

Description

Shape-maintaining fitting method during image is processed

(1) technical field

The present invention relates to the image interpolation technology, belong to image processing field.

(2) background technology

Along with the development of infotech, digital picture occurs in social life in the increasing field, and having become computing machine needs extremely important class data to be processed.And the basic operation of processing as image, image amplifies the background that has a wide range of applications, most of with digital picture as the demand that all can have image to amplify in the application of data.In magic magiscan, the CT image is amplified, can bring great convenience for diagnosis.Good image magnification method can make the image that obtains have higher quality, thereby provides very large facility for further working.

It is exactly the pixel that makes new advances according to Pixel Information interpolation discrete in the original image that image amplifies.The structure interpolation curved surface has several different methods, and such as the arest neighbors interpolation, bilinear interpolation is based on the interpolation of B é zier curved surface etc.The simplest method is the zeroth order interpolation, namely the isolated ground of each pixel is amplified in proportion, and the image that produces like this has serious mosaic.The interpolation curved surface of arest neighbors interpolation is surface of discontinuity, and enlarged image has obvious mosaic; It is continuous that the bilinear interpolation curved surface can reach CO, is continuous in the pixel region inside of original image, but discontinuous on the border, transition is level and smooth not between pixel.For solving mosaic phenomenon and border continuous problem, usually for the continuous mathematical model of original image structure, according to the scaling ratio model is resampled again, to obtain the image of corresponding size.The curved surface that obtains such as the interpolation based on B é zier curved surface is that C1 is continuous in whole image space, and the enlarged image that therefore obtains is level and smooth on the whole, reasonable effect is arranged, but the outline line in the image also can be owing to smooth effect thickens.Also have a lot of similarly methods of high-order interpolation, these methods can overcome mosaic phenomenon effectively, but because the smooth effect of high-order interpolation, so that the image detail loss is serious, the objects in images soft edge is unclear.For this reason, the interpolation algorithm based on fitted dividing curve is suggested.This method is at first carried out match to the outline line in the image, for the pixel region that has outline line to pass through, for outline line the zones of different that this pixel segmentation becomes is carried out respectively interpolation, to keep the clear of outline line.It is continuous that but the boundary except outline line of the pixel region that the interpolation curved surface of the method structure passes through at outline line can not reach C1.Based on the inhomogeneous interpolation image amplification method of cubic spline, the burst continuation algorithm equal segments interpolation algorithm of image scaling is suggested in succession, also is used for the processing of outline line in the optimized image.But these piecewise interpolations have also destroyed whole interpolation curved surface in the continuity of pixel adjoiner guaranteeing that outline line place curved surface is discontinuous with when giving prominence to profile.

(3) summary of the invention

Mosaic or the serious phenomenon of loss of detail occur easily when utilizing existing method to carry out zooming, how discussion of the present invention effectively overcomes mosaic phenomenon, makes the image behind the scaling keep clearly image detail.Discrete picture is done scaling operation, and one of effective method is to obtain the former curved surface that generates discrete picture, and former curved surface is directly done the scaling operation.Structure generates the former curved surface of discrete picture, need to know the value of every bit on the former curved surface, and each known discrete picture is an area sampling value, therefore can according to the counter approximate value of asking every bit on the former curved surface of area sampling value, then construct the Proximal surface of former curved surface according to approximate value.Basic thought of the present invention can be summarized as follows.At first according to the counter approximate value of asking particular point on the former curved surface of discrete picture, the approximate value of particular point forms a quadrilateral mesh, then construct bicubic Coons interpolation curved surface on the quadrilateral mesh, all bicubic Coons interpolation curved surface amalgamations form the Proximal surface of former curved surface together.At the zone structure bicubic Coons interpolation curved surface that adjacent four pixels surround, wherein, the single order partial derivative of four angle points and second order local derviation are estimated to obtain by the shape preserving interpolation method.According to the scaling ratio Coons curved surface is resampled on this curved surface, we adopt area sampling to obtain the gray-scale value of resample points.By increasing sampling density, obtain more pixel, thereby reach the purpose that image amplifies.

The present invention can be divided into three parts, A, according to the counter approximate value of asking former Surface Data point of discrete picture; The structure of B, Coons curved surface; C, resample according to the scaling ratio.The step of then utilizing the shape preserving interpolation Coons surface to carry out the image amplification is described below:

Step1: according to the counter approximate value of asking particular point on the former curved surface of discrete picture;

Step2: according to the value of the data point of obtaining obtain data point at u to the derivative that makes progress with v;

Step3: the second order local derviation of obtaining the data point place;

Step4: the bicubic interpolation patch calculates fast.

Step5: according to the scaling ratio Coons patch is resampled, obtain more accurate pixel;

(4) description of drawings

Fig. 1, dull non-protruding situation is for P _{I-2, j}, P _{I-1, j}, P _{I, j}And P _{I+1, j}The broken line of four some formation is dull non-protruding;

Fig. 2, the situation of convex polygon is for P _{I-2, j}, P _{I-1, j}, P _{I, j}And P _{I+1, j}Four points couple together and form a plane convex polygon;

Fig. 3, each picture element differentiate and bicubic patch organigram, the present invention be its each by adjacent four pixels, such as P _{I, j}, P _{I+1, j}, P _{I, j+1}, P _{I+1, j+1}Structure bicubic Coons interpolation curved surface is put P with adjacent image point on the zone that surrounds _{I-1, j}, P _{I+1, j}And P _{I+1, j}Three point estimation P _{I, j}The u at place is to derivative (P _{I, j}') _u, use pixel P _{I, j-1}, P _{I, j}And P _{I, j+1}Three point estimation P _{I, j}The v at place is to derivative (P _{I, j}') _v

Fig. 4, the area sampling synoptic diagram, wherein, the pixel of solid circle expression original image, the point that hollow circle expression resamples.For resample points 1, corresponding sample area is A, 2,3 points that resample, and its sample area is respectively B, C, and the value that resamples is determined jointly by its two adjacent Coons curved surfaces.

To be two kinds of methods amplify 2 times comparison to the CT image for Fig. 5, Fig. 6, and Fig. 5 is that bilinear interpolation method obtains, and Fig. 6 is obtained by method of the present invention.The quality of Fig. 5 and Fig. 6 two width of cloth images does not visually have significant difference.For comparing mosaic phenomenon, we amplify high power to image, see Fig. 7, Fig. 8.

To be two kinds of methods amplify 5 times comparison to the CT image for Fig. 7, Fig. 8, and Fig. 7 is that bilinear interpolation method obtains, and Fig. 8 is obtained by method of the present invention.On display, the mosaic phenomenon of Fig. 7 is very obvious, and Fig. 8 does not just have mosaic phenomenon.

(5) embodiment

The below at first illustrates the principles of science of the important step institute foundation in the above-mentioned algorithm:

The two-dimentional match that A, based on data excavate

1) mapping relations of area sampling and point sampling

Pixel value on the CT image is the area sampling value, and supposing has n * n pixel on one deck CT image, and then these pixels are the area sampling values that obtain at n * n unit area from space curved surface.Show that for example, during the scaling conversion, ideal situation is to obtain the curved surface that approaches of luv space curved surface carrying out image, then carry out resampling approaching curved surface, obtain to approach curved surface in the area sampling value of each new region.Therefore, if can obtain according to the area sampling value of CT image the more accurate point sampling value of this regional center point, just can according to these point sampling values construct precision higher approach curved surface.So, need to carry out pre-service to the area sampling value of CT image, obtain more accurate point sampling value.We adopt data mining technology to be asked by given pixel and approach curved surface.

Set up now the mapping relations of area sampling and point sampling, the situation of one dimension be discussed first:

If f (x) was n some P _i(i=0,1,2 ..., curve n-1),

(i=0,1,2 ..., n-1) be from this curve obtain corresponding to the zone [x _i, x _I+1] n area sampling value, then have following relation to set up:

${\∫}_{x_i}^{x_{i+1}}f\left(x\right)\mathrm{dx}={\stackrel{\‾}{P}}_i\×(x_{i+1}-x_i),i=\mathrm{0,1,2},...,n-1---\left(1\right)$

By (1) formula as can be known, by setting up interpolation point P _iWith the area sampling value

(i=0,1,2 ..., relation n-1), can by

Obtain P _j(i=0,1,2 ..., n-1), namely obtain more accurate point sampling value by the area sampling value.

The situation of two dimension is discussed again:

If P (x, y) represents a curved surface, P _Ij(i, j=0,1,2 ..., n-1) be n * n point on this curved surface, obviously, P (x, y) can put approximate representation by this n * n.Again (i, j=0,1,2 .., n-1) be from curved surface P (x, y) obtain corresponding to regional A _IjN * n area sampling value, then have following relation to set up:

$\underset{A_{\mathrm{ij}}}{\∫\∫}P(x,y)\mathrm{dxdy}={\stackrel{\‾}{P}}_{\mathrm{ij}}\×S_{\mathrm{ij}},i=\mathrm{0,1,2},...n-1---\left(2\right)$

Wherein, S _IjRepresent regional A _IjArea.

By (2) formula as can be known, by setting up some P _IjWith the area sampling value (i, j=0,1,2 ..., relation n-1), can by

Obtain P _Ij(i, j=0,1,2 ..., n-1), namely obtained the approximate value of point sampling value by the area sampling value.

By the mapping relations of above-mentioned zone sampling and point sampling as can be known, by curve construction/curved surface, can be obtained by the area sampling value of known CT image comparatively accurately point sampling value.

2) calculating of point sampling value

This algorithm utilizes the mapping relations of above-mentioned area sampling and point sampling, on the basis of CT area sampling data, obtains the point sampling value of more accurate CT pixel, adopts secondary Lagrange's interpolation curved surface to approach original curved surface in the mapping process.In addition, the interior pixels point of CT image has been taked different disposal routes with the boundary pixel point, made the sampled value of the pixel that obtains more reasonable and accurate.

If

(i, j=0,1,2 ..., n-1), be n * n pixel value on the tomographic image, A _Ij(i, j=0,1,2 ..., n-1) expression is with (x _j, y _i) centered by unit square, we wish to try to achieve point (x _j, y _i) the more accurate sampled value P that locates _Ij(i, j=0,1,2 ..., n-1).

Calculation level (x ₁, y ₁) the sampled value P that locates ₁₁Algorithm idea as follows:

1, by summit P _Ij(i, j=0,1,2) obtains a secondary lagrangian fit curved surface P (x, y);

2, according to the mapping relations of (2) formula, with P (x, y) respectively at A _IjCarry out integration on (i, j=0,1,2) zone, integrated value equals the known pixel value on this zone

(i, j=0,1,2) multiply by regional A _IjThe area of (i, j=0,1,2);

3, obtain nine yuan of linear function groups by the 2nd step, find the solution and to get P _IjThe value of (i, j=0,1,2).We only keep P ₁₁Value because it is that relative accuracy is the highest in nine new values, the value of other points can adopt same process to obtain.

Specifically being calculated as follows of its each step:

Step 1: utilize summit P _IjSecondary lagrangian fit curved surface of (i, j=0,1,2) structure, surface equation is:

$P(x,y)=\underset{i=0}{\overset{2}{\mathrm{\Σ}}}\underset{j=0}{\overset{2}{\mathrm{\Σ}}}{P}_{\mathrm{ij}}{L}_i\left(x\right)L_j\left(y\right)$

$=\frac{(x-1)(x-2)(y-1)(y-2)}{4}{P}_{00}-\frac{(x-1)(x-2)y(y-2)}{2}{P}_{01}+\frac{(x-1)(x-2)y(y-1)}{4}{P}_{02}---\left(3\right)$

$-\frac{x(x-2)(y-1)(y-2)}{2}{P}_{10}+x(x-2)y(y-2)P_{11}-\frac{x(x-2)y(y-1)}{2}{P}_{12}$

$+\frac{x(x-1)(y-1)(y-2)}{4}{P}_{20}-\frac{x(x-1)y(y-2)}{2}{P}_{21}+\frac{x(x-1)y(y-1)}{4}{P}_{22}$

Step 2: the lagrangian fit curved surface P (x, y) that (3) formula is represented is at regional A _IjCarry out integration on (i, j=0,1,2), and integrated value equals

(i, j=0,1,2) and A _IjThe product of (i, j=0,1,2) area, again A _Ij(i, j=0,1,2) is unit square, so integration satisfies following relational expression:

${\∫}_{x_j}^{x_{j+1}}{\∫}_{y_i}^{y_{i+1}}P(x,y)\mathrm{dxdy}={\stackrel{\‾}{P}}_{\mathrm{ij}}---\left(4\right)$

Step 3: can obtain one with P by (3) and (4) _Ij(i, j=0,1,2) is nine yuan of linear function groups of unknown quantity, separates this system of equations and can obtain P _IjThe value of (i, j=0,1,2).Wherein, we only keep P ₁₁Value.

In the above-mentioned steps, the integral process of (3) formula is actually Lagrangian basis function is carried out integration, is easy to obtain a general expression formula, calculates simple, convenient.In addition, the matrix of coefficients of the system of equations in the step 3 all is identical in all cases, therefore, can obtain first the inverse matrix of matrix of coefficients, then this inverse matrix is used for the solving equation group, thereby has improved the travelling speed of system, has reduced the complexity of system.Inverse matrix is as follows:

$\left[\begin{array}{ccccccccc}3\frac{335}{913}& -2\frac{153}{911}& \frac{467}{779}& -2\frac{153}{911}& 1\frac{79}{156}& -\frac{43}{130}& \frac{467}{779}& -\frac{43}{130}& \frac{70}{521}\\ \frac{593}{977}& -\frac{67}{182}& \frac{8}{67}& 1\frac{243}{443}& -1\frac{19}{250}& \frac{43}{182}& -\frac{184}{619}& \frac{13}{85}& -\frac{7}{97}\\ -\frac{201}{649}& \frac{17}{79}& -\frac{43}{910}& 1\frac{243}{443}& -1\frac{19}{250}& \frac{43}{182}& \frac{402}{649}& -\frac{235}{546}& \frac{43}{455}\\ \frac{593}{977}& 1\frac{243}{443}& -\frac{184}{619}& -\frac{67}{182}& -1\frac{19}{250}& \frac{13}{85}& \frac{8}{67}& \frac{43}{182}& -\frac{7}{97}\\ \frac{47}{412}& \frac{66}{251}& -\frac{53}{862}& \frac{66}{251}& \frac{269}{350}& -\frac{13}{119}& \frac{53}{862}& -\frac{13}{119}& \frac{225}{555}\\ -\frac{1}{19}& -\frac{91}{592}& \frac{13}{595}& \frac{66}{251}& \frac{269}{350}& -\frac{13}{119}& \frac{67}{637}& \frac{91}{296}& -\frac{35}{801}\\ -\frac{201}{649}& 1\frac{243}{443}& \frac{402}{649}& \frac{17}{79}& -1\frac{19}{250}& -\frac{235}{546}& -\frac{43}{910}& \frac{43}{182}& \frac{43}{455}\\ -\frac{1}{19}& \frac{66}{251}& \frac{67}{637}& -\frac{91}{592}& \frac{269}{350}& \frac{91}{296}& \frac{13}{595}& -\frac{13}{119}& -\frac{35}{801}\\ \frac{14}{553}& -\frac{91}{592}& -\frac{53}{862}& -\frac{91}{592}& \frac{269}{350}& \frac{91}{296}& -\frac{53}{862}& \frac{91}{296}& \frac{53}{431}\end{array}\right]$

The structure of B, Coons curved surface

1) utilize the shape preserving interpolation algorithm carry out each picture element at u to the derivative that makes progress with v.Put P with adjacent image point _{I-1, j}, P _{I, j}And P _{I+1, j}Three point estimation P _{I, j}The place u to derivative (P ' _{I, j}) _u, use pixel P _{I, j-1}, P _{I, j}And P _{I, j+1}Three point estimation P _{I, j}The place v to derivative (P ' _{I, j}) _vThe below with (P ' _{I, j}) _uBe example, carry out finding the solution of each picture element derivative.Along u to getting 3 P _{I-1, j}, P _{I, j}And P _{I+1, j}The structure interpolation curve, owing to the step-length such as be between 3, the curve of establishing three point interpolation is P (t), then its first order derivative is calculated as follows:

(P′ _i，j) _u＝(P _i+1，j-P _i-1，j)/2 (5)

To border data point P _{1, j}And P _{N, j}, utilize following formula calculate (P ' _{1, j}) _u=2 (P _{2, j}-P _{1, j})-(P ' _{2, j}) _u, (P ' _{N, j}) _u=2 (P _{N, j}-P _{N-1, j})-(P ' _{N-1, j}) _u

Directly estimate that by (5) formula three Hermite curves of derivative structure do not have the shape of data point suggestion sometimes, this is unacceptable to clinical practice.Therefore need to apply to the derivative that (5) formula is estimated certain constraint, make the cubic curve of structure have the shape of data point suggestion.Theoretical analysis is known, constructs three Hermite curves by the derivative of (5) formula definition, if curve has identical monotonicity and convexity with data point, then (P ' _{I, j}) _iNeed satisfy following condition

${\left(P_{i,j}^{\&prime;}\right)}_u=\left\{\begin{array}{cc}\min [\max (0,{\left(P_{i,j}^{\&prime;}\right)}_u),3\min ({\mathrm{\&Delta;P}}_{i-1,j},{\mathrm{\&Delta;P}}_{i,j})]& 0<\min ({\mathrm{\&Delta;P}}_{i-1,j},{\mathrm{\&Delta;P}}_{i,j})\\ \max \left[\min (0,{\left(P_{i,j}^{\&prime;}\right)}_u),3\max ({\mathrm{\&Delta;P}}_{i-1,j},{\mathrm{\&Delta;P}}_{i,j})\right]& 0>\max ({\mathrm{\&Delta;P}}_{i-1,j},{\mathrm{\&Delta;P}}_{i,j})\\ 0& 0>={\mathrm{\&Delta;P}}_{i-1,J}*{\mathrm{\&Delta;P}}_{i,j}\end{array}\right.---\left(6\right)$

Wherein, Δ P _{I-1, j}=P _{I, j}-P _{I-1, j}, Δ P _{I, j}=P _{I+1, j}-P _{I, j}

(6) characteristics of the derivative of formula definition be the maximal value of three Hermite curves of constructing less than the maximal value of given data point, minimum value is greater than the minimum value of given data point.In order to improve interpolation precision and to obtain more rational curve shape, we do reduction to the constraint of (6) formula and process, and estimate derivative by (5) formula first, then it are according to circumstances adjusted.Selection mode is as follows:

If 1. for P _{I-2, j}, P _{I-1, j}, P _{I, j}And P _{I+1, j}The broken line of four some formation is dull non-protruding situations, and we can directly select the derivative of (6) definition to estimate, just can satisfy guarantor's type requirement of Coons interpolation.

If 2. for P _{I-2, j}, P _{I-1, j}, P _{I, j}And P _{I+1, j}Four points couple together and form a plane convex polygon, and we need to do reduction to the constraint of (6) formula and process, and process as follows:

To P _{I-1, j}, P _{I, j}Value and (P ' _{I, j}) _uCarrying out interpolation, to get curve as follows:

P＝a ₁u ²+b ₁u+c ₁ (7)

Wherein

$a_1=\frac{P_{i-1,j}-P_{i,j}-{\left({P^{\&prime;}}_{i,j}\right)}_u}{{2x}_{i,j}-1},$ b ₁＝(P′ _i，j) _u-2a ₁u _i，j， ${\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="1.png"\; state="\{\{state\}\}">c</figure-callout>}}_1=P_{i,j}-a_1{u}_{i,j}^2-b_1{u}_{i,j}.$

To P _{I-1, j}, P _{I, j}Value and (P ' _{I-1, j}) _uCarrying out interpolation, to get curve as follows:

Q＝a ₂u ²+b ₂u+c ₂ (8)

By point (P _{I-2, j}, u _{I-2, j}) and (P _{I-1, j}, u _{I-1, j}) straight line that consists of is:

P1＝(P _i-1，j-P _i-2，j)(u-u _i-1，j)+P _i-1，j＝A ₁u+B ₁ (9)

By (P _{I, j}, u _{I, j}) and (P _{I+1, j}, u _{I+1, j}) straight line that consists of is:

P2＝(P _i+1，j-P _i，j)(u-u _i，j)+P _i，j＝A ₂u+B ₂ (10)

Remember that it is (u that two straight lines intersect intersection point _{Imid, j}, P _{Imid, j}), following formula is found the solution

P-P1＝a ₁u ²+b ₁u+c ₁-(P _i，j-P _i-1，j)(u-u _i-1，j)-P _i-1，j＝0 (11)

Then formula (15) has a solution at least, and namely interpolation curve and straight line P1 intersect at the end points place at least, and the end points solution is (P _{I-1, j}, u _{I-1, j}).And if have non-end points solution u _p, and satisfy u _{I-1, j}＜u _p＜u _{Imid, j}, prove that then interpolation curve P is not at a P _{I-1, j}, P _{Imid, j}, P _{I, j}In the convex closure that surrounds, this just need to adjust (P ' _{I, j}) _uAdjust (P ' _{I, j}) _uAfter, as long as make interpolation curve and straight line P ₁Tangent, just can make the interpolation curve P of acquisition at three some P _{I-1, j}, P _{Imid, j}And P _{I, j}Convex closure within, namely satisfy following condition:

$\left\{\begin{array}{c}{A}_{1}u+{\mathrm{<figure-callout\; id="1"\; label="B"\; filenames="1.png"\; state="\{\{state\}\}">B</figure-callout>}}_1=a_1{u}^2+b_{1}u+{\mathrm{<figure-callout\; id="1"\; label="c"\; filenames="1.png"\; state="\{\{state\}\}">c</figure-callout>}}_1\\ {2a}_{1}u+{\mathrm{<figure-callout\; id="1"\; label="b"\; filenames="1.png"\; state="\{\{state\}\}">b</figure-callout>}}_1=A_1\end{array}\right.---\left(12\right)$

The solving equation group can get after (12) (P ' _{I, j}) _u, then (P ' _{I, j}) _uCan satisfy the requirement of guarantor's type.

In like manner, following formula is found the solution:

Q-P2＝a ₂u ²+b ₂u+c ₂-(P _i+1，j-P _i，j)(u-u _i，j)-P _i，j＝0 (13)

Then formula (13) is except end points solution (P _{I, j}, u _{I, j}), may there be non-end points solution u _qIf satisfy u _{Imid, j}＜u _q＜u _{I, j}, prove that then interpolation curve Q is not at a P _{I-1, j}, P _{Imid, j}, P _{I, j}In the convex closure that surrounds, this just need to adjust (P ' _{I-1, j}) _uAdjust (P ' _{I-1, j}) _uAfter, as long as so that curve Q and straight line P2 tangent, Q is at a P _{I-1, j}, P _{Imid, j}, P _{I, j}In the convex closure that surrounds, namely satisfy following condition:

$\left\{\begin{array}{c}{A}_{2}u+B_2=a_2{u}^2+b_{2}u+c_2\\ {2a}_{2}u+b_2=A_2\end{array}\right.---\left(14\right)$

The solving equation group can get after (14) (P ' _{I-1, j}) _u, then (P ' _{I-1, j}) _uCan satisfy the requirement of guarantor's type.

For (P ' _{I, j}) _vFind the solution and adjustment and above-mentioned (P ' _{I, j}) _uFind the solution and adjust similar, do not repeat them here.

2) the second order local derviation is estimated

The calculating of second order local derviation is usually very complicated, and in the present invention, we have selected a kind of fairly simple method of estimation.By finding the solution u to the single order local derviation, find the solution v to the single order local derviation, get it and on average then get P _{I, j}The second order local derviation at place is calculated as follows:

(P″ _i，j) _uv＝(P _i+1，j-P _i-1，j+P _i，j+1+P _i，j-1)/4

3) structure P _{I, j}, P _{I+1, j}, P _{I, j+1}, P _{I+1, j+1}The bicubic Coons patch P that surrounds _Ij(u, v).

Bicubic Coons patch is the curved surface of being vowed, being led resultant second order local derviation information definition by the point of adjacent 4 corner points, is used for carrying out interpolation between one group of data point.For image, the Coons curved surface is applied to the interpolation of pixel gray-scale value.Get the image of a width of cloth m * n, for its each by adjacent four pixels, such as P _{I, j}, P _{I+1, j}, P _{I, j+1}, P _{I+1, j+1}Structure bicubic Coons interpolation curved surface is designated as P on the zone that surrounds _Ij(u, v).

Above-mentioned four parameter values corresponding to angle point are respectively (0,0), (1,0), (0,1), (1,1), then to above-mentioned four summit pixel value P (0,0), P (1,0), P (0,1), P (1,1), and four summits cut resultant second order local derviation totally 16 interpolation information carry out the bicubic Coons curved surface P (u, v) of interpolation, u, v ∈ [0,1] * [0,1] is defined as follows:

$P_{\mathrm{ij}}(u,v)=\left[\begin{array}{cccc}{F}_0\left(u\right)& F_1\left(u\right)& G_0\left(u\right)& G_1\left(u\right)\end{array}\right]C\left[\begin{array}{c}{F}_0\left(v\right)\\ F_1\left(v\right)\\ G_0\left(v\right)\\ G_1\left(v\right)\end{array}\right],$ u，v∈[0，1]×[0，1](15)

Wherein

$C=\left[\begin{array}{cccc}P\left(\mathrm{0,0}\right)& P\left(\mathrm{0,1}\right)& {P^{\&prime;}}_v\left(\mathrm{0,0}\right)& {p^{\&prime;}}_v\left(\mathrm{0,1}\right)\\ P\left(\mathrm{1,0}\right)& P\left(\mathrm{1,1}\right)& {P^{\&prime;}}_v\left(\mathrm{1,0}\right)& {P^{\&prime;}}_v\left(\mathrm{1,1}\right)\\ {P^{\&prime;}}_u\left(\mathrm{0,0}\right)& {P^{\&prime;}}_u\left(\mathrm{0,1}\right)& {P^{\&prime;\&prime;}}_{\mathrm{uv}}\left(\mathrm{0,0}\right)& {P^{\&prime;\&prime;}}_{\mathrm{uv}}\left(\mathrm{0,1}\right)\\ {P^{\&prime;}}_u\left(\mathrm{1,0}\right)& {P^{\&prime;}}_u\left(\mathrm{1,1}\right)& {P^{\&prime;\&prime;}}_{\mathrm{uv}}\left(\mathrm{1,0}\right)& {P^{\&prime;\&prime;}}_{\mathrm{uv}}\left(\mathrm{1,1}\right)\end{array}\right]$

P is the gray-scale value of given picture element, P _u, P _vRepresent that respectively u is to vowing P with guide v _UvThen represent the second order local derviation.In the practical application, only can obtain the gray-scale value of each pixel of piece image, by 1) and 2) part as can be known, utilize the method for shape preserving interpolation algorithm and simple method of estimation second order local derviation can obtain each picture element at u to the derivative and the second order local derviation that make progress with v, substitution formula (1) bicubic Coons patch P _Ij(u, v).

C, resample according to the scaling ratio

1) according to the scaling ratio Coons patch is resampled

If image needs amplifieroperation, we need to resample to original image according to the scaling ratio.The bicubic Coons interpolation curved surface of above-mentioned structure is applied to the gradation of image interpolation, and new sampled point is taken from this bicubic Coons interpolation curved surface sheet.

Suppose that the image enlargement factor is s, then original image I (u, v) is at u, and the v direction resamples according to the interval of l/s, namely obtains amplifying the image after s times.Digital picture I (u, the v) interpolation that is about to m * n is the image I of m ' * n ' ' (u ', v ').Wherein

m′＝m*s；

n′＝n*s；

Then resample points I ' (u ', v ') corresponding to the position of the pixel in the original image is:

u＝u′/s；

v＝v′/s；

The label of this pixel bicubic Coons patch at place in original image is i=int (u), j=int (v), and then according to this interpolation curved surface sheet (1) of hole of structure, the gray-scale value of resample points I ' (u ', v ') is P _Ij(u-i, v-j).

By the process that obtains the CT data as can be known, the CT data value of each pixel in fact is the mean value of the density of its place voxel, is an area sampling value, rather than this an accurate point sampling value.So during resampling, we also adopt area sampling to obtain the gray-scale value of resample points.Here we also adopt area sampling to obtain the gray-scale value of resample points.With the gray-scale value of the mean value on certain zonule, resample points place as resample points:

$I^{\&prime;}\left(u^{\&prime;}{v}^{\&prime;}\right)=\left(\underset{A}{\&Integral;\&Integral;}{P}_{\mathrm{ij}}(u,v)\mathrm{ds}\right)/S_A---\left(16\right)$

Wherein A is resample area.If on the border of former Coons patch, then the value of its area sampling should be determined by its two adjacent Coons curved surfaces resample points jointly just.That is:

$I^{\&prime;}\left(u^{\&prime;}{v}^{\&prime;}\right)=(\underset{B_{j-1}}{\&Integral;\&Integral;}{P}_{i,j-1}(u,v)\mathrm{ds}+\underset{B_j}{\&Integral;\&Integral;}{P}_{i,j}(u,v)\mathrm{ds})/S_B$

$I^{\&prime;}\left(u^{\&prime;}{v}^{\&prime;}\right)=(\underset{C_{i-1}}{\&Integral;\&Integral;}{P}_{i-1,j}(u,v)\mathrm{ds}+\underset{C_i}{\&Integral;\&Integral;}{P}_{i,j}(u,v)\mathrm{ds})/S_C$

By choosing different sampling densities, calculate the color value of sampled point by above-mentioned formula, just can obtain the image of any scaling ratio.

2) the quick calculating of image zooming

If integral domain is (x1, x2) ﹠amp; (y1, y2) then according to the expression-form (15) of Coons curved surface, finds the solution (16), gets (17) formula

I＝(IntegralF ₀(x2)-IntegralF ₀(x1))*(P _i，j*(IntegralF ₀(y2)-IntegralF ₀(y1))+P _i，j+1*(IntegralF1(y2)-IntegralF1(y1))+(P′ _i，j) _v*(IntegralG0(y2)-IntegralG0(y1))+(P′ _i，j+1) _v*(IntegralG1(y2)-IntegralG1(y1)))+(IntegralF1(x2)-IntegralF1(x1))*(P _i+1，j*(IntegralF0(y2)-IntegralF0(y1))+P _i+1，j+1*(IntegralF1(y2)-IntegralF1(y1))+(P′ _i+1，j) _v*(IntegralG0(y2)-IntegralG0(y1))+(P′ _i+1，j+1) _v*(IntegralG1(y2)-IntegralG1(y1)))+(IntegralG0(x2)-IntegralG0(x1))*((P′ _i，j) _u*(IntegralF0(y2)-IntegralF0(y1))+(P′ _i，j+1) _u*(IntegraIF1(y2)-IntegralF1(y1))+(P″ _i，j) _uv*(IntegralG0(y2)-IntegralG0(y1))+P″ _i，j+1) _uv*(IntegralG1(y2)-IntegralG1(y1)))+(IntegralG1(x2)-IntegralG1(x1))*((P′ _i+1，i) _u*(IntegralF0(y2)-IntegralF0(y1))+(P′ _i+1，j+1) _u*(IntegralF1(y2)-IntegralF1(y1)+(P″ _i+1，j) _uv*(IntegralG0(y2)-IntegralG0(y1))+(P″ _i+1，j+1) _uv*(IntegralG1(y2)-IntegralG1(y1)))(17)

Wherein, Integral represents the integration to the function of its back, for F ₀(), F ₁(), G ₀(), G ₁The integrated form of (), IntegralF ₀(t)=t ⁴/ 2-t ³+ t, IntegralF1 (t)=-t ⁴/ 2+t ³, IntegralG0 (t)=t ⁴/ 4-2t ³/ 3+t ²/ 2, IntegralG1 (t)=t ⁴/ 4-t ³/ 3

Above-mentioned expression formula is 4 order polynomial forms, and the process of each integration all will be calculated its 4 order polynomial form, and its calculated amount is too large.But according to the form that resamples as can be known, as long as enlargement factor determines, being assumed to be s, is the same for situation about resampling on all Coons patchs that are made of adjacent four pixels.That is to say, for the resample points of parameter identical (being assumed to be (uu, vv)) on the different Coons patchs, its integral domain is identical, and namely u is limited to (uu-1/2s, uu+1/2s) up and down, v is limited to (vv-1/2s, vv+1/2s) up and down.In conjunction with the form of (17) formula, if to the F in (15) formula ₀(t), F ₁(t), G ₀(t), G ₁(t) at (uu-1/2s, uu+1/2s) ﹠amp; Integration on (vv-1/2s, vv+1/2s) calculates in advance and stores, then in the process of the area sampling that calculates each resample points, and just need not be at every turn all to F ₀(), F ₁(), G ₀(), G ₁() recomputates its integrated form, thereby reduces its time complexity, accelerates computation process.For example image amplifies s doubly, for F ₀(), F ₁(), G ₀(), G ₁The integration of () is in the interval

Integration on (s is the scaling multiple) $\mathrm{IntegralF}0\left(t\right){|}_{i/s-1/2s}^{i/s+1/2s}1<=i<=s-1,$ $\mathrm{IntegralF}1\left(t\right){|}_{i/s-1/2s}^{i/s+1/2s}1<=i<=s-1,$ $\mathrm{IntegralG}0\left(t\right){|}_{i/s-1/2s}^{i/s+1/2s}1<=i<=s-1,$ $\mathrm{IntegralG}1\left(t\right){|}_{i/s-1/2s}^{i/s+1/2s}1<=i<=s-1,$ Calculate in advance and store,

Then can search for corresponding integrated value in record when (17) formula of calculating, (17) formula of bringing into gets final product to get integral result, and experimental result shows that its computing velocity can improve 40%.

We are applied in the present invention in the CT image processing system.The CT image is amplified, can bring great convenience for diagnosis.Its realization can realize by software programming.The false code of given first the present invention programming.

In the CT image processing system, if certain width of cloth CT image is amplified, the multiple that we will need enlarged image and needs to amplify passes to the ZoomBasedCoons method.From the process of amplifying, can see that the present invention selects shape preserving interpolation to calculate the single order local derviation and the second order local derviation is constructed Coons interpolation curved surface sheet, naturally guaranteed continuity and the shape-retaining ability of the image after the amplification, also just guaranteed the flatness of image.According to the ratio of amplifying, resample at the Coons interpolation curved surface sheet of structure, with the mean value of the picture element on the zonule, the resample points place gray-scale value as resample points.Bicubic Coons interpolation curved surface sheet reaches the secondary precision, so can obtain more image detail.So this invention can make the image after the amplification that obtains have higher quality, can effectively overcome mosaic phenomenon, make the image of generation that more clearly image detail be arranged.

As an example, need in the width of cloth diagnosis process to select the image used, respectively it is amplified 2 times and 5 times.Design sketch is seen accompanying drawing Fig. 7, Fig. 8.Can obtain reasonable effect from illustrating as can be known the present invention.

Claims (6)

1. A type-preserving fitting method in image processing, characterized in that: All pixels on the original image are taken from an original surface, and the approximate gray value of the corresponding point on the original surface is calculated according to the gray value of each pixel on the original image to form a series of new data. For new data, a bicubic Coons interpolation surface is constructed on the area surrounded by four adjacent data points, where, The first-order partial derivatives of the four corner points are estimated by the shape-preserving interpolation method, The second-order partial derivatives of the four corner points are obtained by averaging the partial derivatives of their adjacent pixels, Resampling is performed on the bicubic Coons surface according to the scaling, Use area sampling to obtain the gray value of the resampling point, and increase the sampling density to obtain more pixels. Among them, the first-order partial derivatives of the four corner points have quadratic approximation accuracy, and the solution method is as follows: Use the adjacent pixel points P _{i-1, j} , P _{i, j} and P _{i+1, j} to estimate the u derivative (P′ _{i, j} ) _u at P _i, j, and use pixel P _{i , j-1} , P _{i, j} and P _{i, j+1} three-point estimate of the v-director (P′ _{i, j} ) _v at P _i, j, for (P′ _{i, j} ) _u , along the u direction Take three points P _{i-1, j} , P _{i, j} and P _{i+1, j} to construct an interpolation curve. Since the three points are of equal step length, let the curve for interpolation of three points be P(t), then one of The order derivative is calculated as follows: (P' _{i, j} ) _u = (P _{i+1, j} - P _{i-1, j} )/2. 2. The shape-preserving fitting method in image processing according to claim 1, characterized in that: the step of forming new data comprises, utilizing the mapping relationship between area sampling and point sampling, adopting a quadratic Lagrangian interpolation surface To approximate the original surface to obtain more accurate point sampling values of pixel points, where the quadratic Lagrangian interpolation surface equation is:

                                                  

Integrate the Lagrangian interpolation surface P(x, y) represented by formula (1) on the area A _ij (i, j=0, 1, 2), and the integral value is equal to

Product with the area of A _ij (i, j=0, 1, 2):

From formula (1) and formula (2), a system of linear equations with P _ij (i, j=0, 1, 2) as the unknown quantity can be obtained. By solving this system of equations, P _ij (i, j =0, 1, 2), where only the value of P _ij is reserved, The integration process is actually to integrate the Lagrangian basis function, and the inverse matrix of the coefficient matrix of the equation system is as follows:

Use this inverse matrix to solve the system of equations. 3. the shape-preserving fitting method in the image processing according to claim 1, is characterized in that: in order to guarantee the shape-preserving property of interpolation curved surface, adopt following method to P _i, the u derivative (P' _i at j place _{, j} ) _u to limit, If for P _{i-2, j} , P _{i-1, j} , P _{i, j} and P _{i+1, j} the polyline segment formed by the four points is monotone and non-convex, for (P′ _{i, j} ) _u The restrictions are as follows:

If the four points of P _{i-2, j} , P _{i-1, j} , P _{i, j} and P _{i+1, j} are connected to form a planar convex polygon, the constraint of formula (3) needs to be weakened, That is, adjust (P′ _{i , j ) u} so that the interpolated curve and straight _line P _{i -2, j} P _{i-1, j} are tangent, so that the obtained interpolation curve is within the convex hull of the three points P _{i-1, j} , P _{imid, j} and P _{i, j} , then (P′ _{i, j} ) _u It can meet the requirement of type protection, where, Δu _i =u _i+1,j -u _i,j ,

Δu _i-1 =u _i,j -u _i-1,j ,

P _{imid, j} is the intersection point of straight line P _{i-2, j} P _{i-1, j} and straight line P _{i, j} P _{i+1, j} . 4. the type-preserving fitting method in the image processing according to claim 1 is characterized in that: get u to the average of first-order partial derivation and v to first-order partial derivation then get Pi _{, the second-order partial derivation at j} place guide: (P″ _{i, j} ) _uv = (P _{i+1, j} - P _{i-1, j} + P _{i, j+1} + P _{i, j-1} )/4. 5. the type-preserving fitting method in the image processing according to claim 1, is characterized in that: according to scaling ratio Coons curved surface sheet is carried out resampling, adopts area sampling to obtain the gray value of resampling point: The image magnification factor is s, then the original image I(u, v) is resampled at the interval of 1/s in the direction of u, v, that is, the image enlarged by s times is obtained, that is, the digital image I(u, v) of m×n ) is an image I'(u', v') interpolated as m'×n', where the position of the resampling point I'(u', v') corresponding to the pixel in the original image is: u=u'/s; v=v'/s; The label of the bicubic Coons surface sheet where the pixel is located in the original image is i=int (u), j=int (v), and the gray value of the resampling point is obtained by using area sampling, as follows:

If the resampling point is exactly on the boundary of the original Coons surface patch, the value of its area sampling should be jointly determined by its adjacent two bi-cubic Coons surfaces, namely:

6. The type-preserving fitting method in image processing according to claim 5, characterized in that: wherein, adopting area sampling to obtain the gray value of the resampling point comprises: utilizing the following method to accelerate the solution of the integral: Solve (4) to get formula (5) I=(IntegralF ₀ (x2)-IntegralF ₀ (x1))*(P _i,j *(IntegralF ₀ (y2)-IntegralF ₀ (y1))+P _i,j+1 *(IntegralF1(y2)-IntegralF1 (y1))+(P′ _i,j)v *(IntegralG0(y2)-IntegralG0(y1))+(P′i _,j+1 ) _v *(IntegralG1(y2)-IntegralG1(y1)))+ (IntegralF1(x2)-IntegralF1(x1))*(P _{i+1, j} *(IntegralF0(y2)-IntegralF0(y1))+P _{i+1, j+1} *(IntegralF1(y2)-IntegralF1(y1 ))+(P′ _{i+1, j} ) _v *(IntegralG0(y2)-IntegralG0(y1))+(P′ _{i+1, j+1} ) _v *(IntegralG1(y2)-IntegralG1(y1)) )+(IntegralG0(x2)-IntegralG0(x1))*((P′ _{i, j} ) _u *(IntegralF0(y2)-IntegralF0(y1))+(P′ _{i, j+1} ) _u *(IntegralF1( y2)-IntegralF1(y1))+(P″ _{i, j} ) _uv *(IntegralG0(y2)-IntegralG0(y1))+(P″ _{i, j+1} ) _uv *(IntegralG1(y2)-IntegralG1(y1 )))+(IntegralG1(x2)-IntegralG1(x1))*((P′ _{i+1, j} ) _u *(IntegralF0(y2)-IntegralF0(y1))+(P′ _{i+1, j+1} ) _u *(IntegralF1(y2)-IntegralF1(y1))+(P″ _{i+1, j} ) _uv *(IntegralG0(y2)-IntegralG0(y1))+(P″ _{i+1, j+1} ) _uv *(IntegralG1(y2)-IntegralG1(y1)))(5) Among them, F ₀ (), F ₁ (), G ₀ (), G ₁ () are

Integral represents the integral of the function behind it, then the integral form of F ₀ (), F ₁ (), G ₀ (), G ₁ () is as follows: IntegralF0(t)=t ⁴ /2-t ³ +t, IntegralF1(t)=-t ⁴ /2+t ³ , IntegralG0(t)=t ⁴ /4-2t ³ /3+t ² /2, IntegralG1(t)=t ⁴ /4-t ³ /3, for The integral of F ₀ (), F ₁ (), G ₀ (), G ₁ () is in the interval

Integral over 1<=i<=s-1

1<=i<=s-1,

1<=i<=s-1,

1<=i<=s-1,

1<=i<=s-1, calculated and stored in advance, where s is the zoom multiple, then when calculating the area sampling value, search for the corresponding integral value in the record, and bring it into formula (5) to get the integral result , among them, P _{i, j} , P _{i+1, j} , P _{i, j+1} , P _{i+1, j+1} are four adjacent points, and the u derivative at P _{i, j} is (P′ _{i, j} ) _u , the derivative of v is (P′ _{i, j} ) _v , the derivative of u at P _{i, j+1} is (P′ _{i, j+1} ) _u , and the derivative of v is (P′ _{i , j+1} ) _v , P _{i+1, the u derivative at j} is (P′ _{i+1, j} ) _u , and the v derivative is (P′ _{i+1, j} ) _v , P _i+1, The u derivative at _j+1 is (P′ _{i+1, j+1} ) _u , the v derivative is (P′ _{i+1, j+1} )v, x1, y1, x2, y2 are the upper and lower limits of the integral .

Images