CN106097291A - A kind of ART flame based on gradient total variation section restructing algorithm - Google Patents

Abstract

The invention discloses a kind of ART flame based on gradient total variation section restructing algorithm, flame image is carried out pretreatment；Generate projection picture；Iteration weight matrix W is sought according to ART criterion；Image is reconstructed according to ART Algebraic Iterative Method；Solve gradient and the gradient total variation of each iterative image；Image is adjusted according to radial direction TV method；Output reconstruct image.Compared with prior art, the present invention can be effectively improved section reconstruction quality and the reconstructed velocity of flame image；Especially in the case of projection angle is less, the decorrelation of this algorithm is more preferable, and gradient image is more sparse, and flame restructuring procedure good stability, and algorithm robustness is strong.

Description

A kind of ART flame based on gradient total variation section restructing algorithm

Technical field

The present invention relates to flame image reconstruction field, particularly relate to a kind of combustion in IC engine flame section restructing algorithm.

Background technology

The air pollution that motor vehicle exhaust emission in recent years causes is on the rise, it has also become the major reason that haze weather is formed One of, the energy-saving and emission-reduction of internal combustion engine have caused global concern.Combustion flame as the important sign of IC engine cylinder combustion process, Significant to research flame combustion mechanism, reduction burning blowdown and raising engine performance.Contemporary optics method and meter Calculate the development of machine technology so that utilizing optic visualization technology to be directly observed internal combustion engine flame combustion process becoming can Energy.At present internal combustion engine combustion flame is studied, mostly utilize the cylinder flame single width two dimension obtained from piston base or top to throw Shadow image carries out two dimensional character analysis, but owing to flame has translucence, the flame image utilizing the method to shoot is whole In-cylinder combustion flame is in the two-dimensional overlay result in specific projection face, it is impossible to characterize the three dimensions feature of flame.It is therefore desirable to Flame is carried out three-dimensional configuration reconstruct, and the quantitative analysis for flame provides enough information.

Scientific research personnel has carried out corresponding research to the reconstruct of flame three-dimensional configuration both at home and abroad, mainly utilizes high-speed synchronous video camera Shoot the flame two-dimensional projection image in multiple directions, and based on computed tomography (Computed Tomography, letter Claiming CT) flame two dimension slicing is reconstructed to form three-dimensional image by technology.The Samuel of University of Toronto uses for reference chromatography Imaging technique and the thought of CT technology image reconstruction, use the method for layering manufacture that nature flame is carried out three-dimensionalreconstruction. Flame is shot by J.Floyd etc. based on 5 CCD camera and mirror-image system, and utilizes the CT technology 10 width flames to obtaining View is reconstructed, and its reformatted slices effectiveness comparison is coarse, awaits the most perfect.The Tadashi Ito of Japan etc. introduce The switch of one foreign radiation sources controls the acquisition of flame thermal radiation projection and uses CT technology to divide to the temperature reconstructing flame Cloth image.Burner hearth flame three-dimensionalreconstruction aspect can law court scholar etc. be entered by the Zhou Huaichun of Chinese University of Science and Technology and the Cen of Zhejiang University Go numerous studies, used digitized camera head from hearth combustion spatial extraction burning two dimension flame radiation image, and utilize CT technology carries out the reconstruct of flame three-dimensional temperature field.Therefore, to the reconstruct of flame three-dimensional configuration it is crucial that to flame two-dimensional projection Image carries out section reconstruct, and its section reconstruction accuracy determines the effect of flame three-dimensionalreconstruction, but the studies above method is the most sharp Directly flame section is reconstructed by CT technology, the most effectively utilizes the sparse characteristic of flame, cause it at projection angle relatively In the case of little, quality reconstruction is poor.

Summary of the invention

For problem above, the present invention proposes a kind of ART flame based on radial direction TV section restructing algorithm, and the method is Optimization to flame two dimension slicing restructing algorithm, according to the radial diffusion feature of combustion flame, it is contemplated that flame image becomes complete Differ from the openness of (Total Variation, TV) territory, tradition TV expanded to polar coordinate system model from cartesian coordinate system, It is incorporated into algebraic reconstruction method (Algebraic reconstruction technique is called for short ART) flame two dimension slicing reconstruct In.

The present invention proposes ART flame based on radial direction TV section restructing algorithm, and the method comprises the following steps:

Step 1, flame sectioning image is carried out pretreatment: include greyscale transformation and smoothing denoising operation；

Step 2, the iteration weight matrix W trying to achieve the reconstruction of flame sectioning image acquisition corresponding projection image array P, tool Body step is as follows:

By the two dimension discrete volume elements independent for N=n × n of flame section, each volume elements has and specifically represents this list Value x of unit's brightness size_j(j=1,2 ..., N), corresponding one of each pixel in two-dimensional projection image passes this section The ray L of volume elements_k(k=1,2 ..., M), M is the ray sum of all projections, equal to projection angle number and certain Angles Projections The product of number of rays；And value P of this pixel_k(k=1,2 ..., M) be on this ray all volume elements to this ray brightness tribute The superposition offered；If ω_jkFor volume elements j to ray L_kWeight factor, represent volume elements and interradial dependency, volume elements j correlation Line L_kContribution be:

P_jk=ω_jk·x_j

In section, all volume elements are to L_kBrightness contribution and being expressed as:

$P_k=\underset{j=1}{\overset{N}{\Σ}}{P}_{jk}=\underset{j=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_{jk}\·x_j,k=1,2,...,M$

It is expressed as follows by the form of matrix:

Wx=P

P is that M ties up projection matrix, represents the pixel value of two-dimensional projection；X is that N-dimensional treats reformatted slices data vector；W is M × N Dimension projection weight matrix；

Step (3), setting primary iteration number of times i=1 and the maximum Maxcount=15 of iterations, and iteration is set Initial vector f⁽⁰⁾=0,0 ..., 0}；

Step 4, calculating reconstruct image, equation group formula is as follows:

$f^{\left(m\right)}=f^{(m-1)}+{\mathrm{\ω}}_{jk}(P_k-\underset{j=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_{jk}{f}^{(m-1)})/\underset{j=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_{jk}^2$

Wherein m=1,2 ..., M represents the equation number of the equation group of complete iteration each time；Public according to above-mentioned equation group Formula, calculates M × N equation group；By initial vector f⁽⁰⁾=0,0 ..., 0} substitutes into first equation of equation group and obtains f⁽¹⁾, by Second equation obtains f⁽²⁾, the like, finally tried to achieve f by m-th equation^(M), this completes a complete iteration；

Step 5, solve f^(M)GradientFormula is as follows:

$\▿f^{\left(M\right)}=(D_{\mathrm{\ρ}}{x}_{\mathrm{\ρ},\mathrm{\θ}},D_{\mathrm{\θ}}{x}_{\mathrm{\ρ},\mathrm{\θ}})$

Wherein D_ρAnd D_θIt is the discrete differential operator along footpath, pole and polar angle direction respectively, has equation below to calculate:

$D_{\mathrm{\ρ}}{x}_{\mathrm{\ρ},\mathrm{\θ}}=x_{\mathrm{\ρ},\mathrm{\θ}}-x_{{\mathrm{\ρ}}_{-},\mathrm{\θ}}$

$D_{\mathrm{\θ}}{x}_{\mathrm{\ρ},\mathrm{\θ}}=x_{\mathrm{\ρ},\mathrm{\θ}}-x_{{\mathrm{\ρ},\mathrm{\θ}}_{-}}$

Wherein ρ_-And ρ₊Represent respectively and reduce along footpath, pole and the change in increase direction, footpath, pole, θ_-And θ₊Represent the most respectively along pole Angle is clockwise and the anticlockwise change of polar angle；

Step 6, solve f^(M)The partial derivative of gradient total variationAnd to f^(M)Carry out gradient Total variation optimizes, such as following formula:

$f^{\left(M\right)}\′=f^{\left(M\right)}-q\·\▿{||f^{\left(M\right)}||}_{R-TV}$

Wherein, q is for declining step factor, and in the present invention, value is 0.1/k, and uses f^(M)' replace original f^(M)；

If step 7 i≤Maxcount, by f^(M)' carry out complete iteration next time according to the equation group of step 4, otherwise, Stop iteration, output reconstruct image.

ART flame based on the radial direction TV section restructing algorithm of the present invention, it is possible to be effectively improved the section weight of flame image Structure quality and reconstructed velocity；Especially in the case of projection angle is less, the decorrelation of this algorithm is more preferable, gradient image More sparse, and flame restructuring procedure good stability, algorithm robustness is strong.

Accompanying drawing explanation

Fig. 1 is x_ρ,θThe searching process schematic of neighbor pixel；Reference is: (1-1), searching(1- 2), x is found_ρ+,θ；(1-3) find(1-4), find

Fig. 2 is ART flame based on the radial direction TV section restructing algorithm flow chart of the present invention；

Fig. 3 is the gray level image of two width simulation flame sections, and reference is: (3-1), image 1；(3-2), image 2；

Fig. 4 is the PSNR change curve along with iterations reconstructing image under two image difference projection angles, accompanying drawing It is labeled as: (4-1), image 1 are under 18 jiaos；(4-2), image 1 is under 36 jiaos；(4-3), image 1 is under 72 jiaos；(4-4), image 2 under 18 jiaos；(4-5), image 2 is under 36 jiaos；(4-6), image 2 is under 72 jiaos；

Fig. 5 is reconstruct result schematic diagram.

Detailed description of the invention

For making the object, technical solutions and advantages of the present invention clearer, below in conjunction with accompanying drawing to embodiment party of the present invention Formula is described in further detail.

The sparse reconstruct of compressed sensing and algebraically iterative reconstruction method all have and go out original complete graph from less is more data reconstruction As the feature of data, and petrol engine combustion flame image presents the sparsity structure form of brightness aggregation block, and it presents relatively in TV territory Good sparse characteristic, therefore can be incorporated into total variation optimization in ART algorithm and gasoline engine flame is carried out sectioning image weight Structure, to improve the quality reconstruction of the section of flame.

First, set up gradient total variation model: the center of flame image is set to initial point, set up polar coordinate system, meter Calculate gradient and the gradient total variation of flame image；At definition initial point, gradient is 0, for pixel x_ρ,θSit in pole Position in mark system represents with M (ρ, θ), and its partial gradient is expressed as:

$\▿x_{\mathrm{\ρ},\mathrm{\θ}}=(D_{\mathrm{\ρ}}{x}_{\mathrm{\ρ},\mathrm{\θ}},D_{\mathrm{\θ}}{x}_{\mathrm{\ρ},\mathrm{\θ}})---\left(1\right)$

Wherein D_ρAnd D_θIt is the discrete differential operator along footpath, pole and polar angle direction respectively, has equation below to calculate:

$\begin{array}{c}{D}_{\mathrm{\ρ}}{x}_{\mathrm{\ρ},\mathrm{\θ}}=x_{\mathrm{\ρ},\mathrm{\θ}}-x_{{\mathrm{\ρ}}_{-},\mathrm{\θ}}\\ D_{\mathrm{\θ}}{x}_{\mathrm{\ρ},\mathrm{\θ}}=x_{\mathrm{\ρ},\mathrm{\θ}}-x_{\mathrm{\ρ},{\mathrm{\θ}}_{-}}\end{array}---\left(2\right)$

Solved gradient total variation by gradient, need footpath, pole ρ and polar angle θ both direction image slices vegetarian refreshments gradient Vector field homoemorphism is added, i.e. the mould of gradient vector is image pixel along component flat of footpath, pole ρ and polar angle θ both direction gradient The evolution summation of side's sum, computing formula is as follows:

$\begin{array}{c}{||x||}_{R-TV}=\underset{\mathrm{\ρ},\mathrm{\θ}}{\Σ}\sqrt{D_{\mathrm{\ρ}}{x_{\mathrm{\ρ},\mathrm{\θ}}}^2+D_{\mathrm{\θ}}{x_{\mathrm{\ρ},\mathrm{\θ}}}^2}\\ =\underset{\mathrm{\ρ},\mathrm{\θ}}{\Σ}\sqrt{{(x_{\mathrm{\ρ},\mathrm{\θ}}-x_{{\mathrm{\ρ}}_{-},\mathrm{\θ}})}^2+{(x_{\mathrm{\ρ},\mathrm{\θ}}-x_{\mathrm{\ρ},{\mathrm{\θ}}_{-}})}^2}\end{array}---\left(3\right)$

During program realizes, the total variation optimized algorithm of the derivative approximation of Sidky as reference, is then incited somebody to action radially by this algorithm Gradient TV adjustment derivative approximation thought is calculated as follows formula and is realized, and so far establishes gradient total variation mould Type.

$\begin{array}{c}\▿{||x||}_{R-TV}=\frac{\∂{||x||}_{R-TV}}{\∂x_{\mathrm{\ρ},\mathrm{\θ}}}\≈\frac{2x_{\mathrm{\ρ},\mathrm{\θ}}-x_{{\mathrm{\ρ}}_{\_},\mathrm{\θ}}-x_{\mathrm{\ρ},{\mathrm{\θ}}_{\_}}}{\sqrt{{(x_{\mathrm{\ρ},\mathrm{\θ}}-x_{{\mathrm{\ρ}}_{\_},\mathrm{\θ}})}^2+{(x_{\mathrm{\ρ},\mathrm{\θ}}-x_{\mathrm{\ρ},{\mathrm{\θ}}_{\_}})}^2+\mathrm{\ϵ}}}\\ +\frac{x_{\mathrm{\ρ},\mathrm{\θ}}-x_{{\mathrm{\ρ}}_{+},\mathrm{\θ}}}{\sqrt{{(x_{{\mathrm{\ρ}}_{+},\mathrm{\θ}}-x_{\mathrm{\ρ},\mathrm{\θ}})}^2+{(x_{{\mathrm{\ρ}}_{+},\mathrm{\θ}}-x_{\mathrm{\ρ},{\mathrm{\θ}}_{\_}})}^2+\mathrm{\ϵ}}}\\ +\frac{x_{\mathrm{\ρ},\mathrm{\θ}}-x_{\mathrm{\ρ},{\mathrm{\θ}}_{+}}}{\sqrt{{(x_{\mathrm{\ρ},{\mathrm{\θ}}_{+}}-x_{\mathrm{\ρ},\mathrm{\θ}})}^2+{(x_{\mathrm{\ρ},{\mathrm{\θ}}_{+}}-x_{{\mathrm{\ρ}}_{-},{\mathrm{\θ}}_{+}})}^2+\mathrm{\ϵ}}}\end{array}---\left(4\right)$

During the gradient of gradient total variation model calculates, for pixel x_ρ,θThe seeking of neighbor pixel Look for process as follows: for the discrete differential in direction, footpath, pole, first look for footpath, pole and reduce the pixel set of 1 unitPoint in set is positioned at zero as the center of circle, and radius is on the circle of ρ-1.At these pixels In, find the gradient component that the pixel vectors pixel minimum with preimage vegetarian refreshments vector angle is the minimizing direction, footpath, pole of needs Neighbor pixel, is expressed asThe gradient component neighbor pixel in increase direction, footpath, pole is found by above-mentioned similar method, It is expressed asFor the discrete differential in polar angle direction, first look for the pixel set constant in footpath, polePoint in set is positioned at zero as the center of circle, and radius is on the circle of ρ.In these pixels, With preimage vegetarian refreshments x_ρ,θFor starting point, find the pixel that pixel vectors is minimum with preimage vegetarian refreshments vector angle clockwise, be needs Polar angle reduce direction gradient component neighbor pixel, be expressed asPolar angle is found counterclockwise by above-mentioned similar method Increase the gradient component neighbor pixel in direction, be expressed asAs procedure described above and formula (1) can solve required critical path To gradient.

The concrete implementation step of ART flame based on the gradient total variation section restructing algorithm of the present invention is as follows:

Step 1: flame sectioning image is carried out pretreatment: include greyscale transformation and smoothing denoising operation；

Step 2: try to achieve iteration weight matrix W acquisition corresponding projection image array P that flame sectioning image is rebuild；

Step 3: set primary iteration number of times i=1 and the maximum Maxcount=15 of iterations, and iteration is set Initial vector f⁽⁰⁾=0,0 ..., 0}；

Step 4: calculate reconstruct image according to ART restructing algorithm, the method for relaxation that the general Kaczmarz of employing proposes solves, Choosing relaxation factor is constant 1,

$f^{\left(m\right)}=f^{(m-1)}+{\mathrm{\ω}}_{jk}(P_k-\underset{j=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_{jk}{f}^{(m-1)})/\underset{j=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_{jk}^2---\left(5\right)$

Wherein m=1,2 ..., M represents the equation number of the equation group of complete iteration each time；Equation according to formula (5) Group, calculates M × N equation group；By initial vector f⁽⁰⁾=0,0 ..., 0} substitutes into first equation of formula (4) equation group and obtains f⁽¹⁾, second equation obtain f⁽²⁾, the like, finally tried to achieve f by m-th equation^(M), this completes and the most completely change Generation；

Step 5: according to above-mentioned gradient total variation model, solve f^(M)Gradient

Step 6: according to derivative approximation thought, seek the total differential of gradient total variationAnd to f^(M)Enter Row gradient total variation optimizes, and i.e. uses gradient descent method to be adjusted such as following formula, and wherein q is for declining step factor, this Bright middle value is 0.1/k, and uses f^(M)' replace original f^(M)

$f^{\left(M\right)}\′=f^{\left(M\right)}-q\·\▿{||f^{\left(M\right)}||}_{R-TV}---\left(6\right)$

Step 7: if i≤Maxcount, by f^(M)' carry out complete iteration next time according to formula (5), otherwise stop iteration Output reconstruct image.

In gradient total variation model during radial direction gradient calculation, for pixel x_ρ,θThe searching of neighbor pixel Journey is as shown in Figure 1.

The specific embodiment of the present invention describes in detail as follows:

1) the different simulation flame sectioning image of selection two width, is intended as shown in Figure 3:

2), intend choosing the maximum Maxcount=15 of primary iteration number of times i=1 and iterations, decline step factor q For 0.1/k；

3) then according to the algorithm flow of technique scheme calculates gradient and the footpath of flame sectioning image respectively To gradient total variation, and the reconstruct image reconstructing out to ART is optimized adjustment, completes flame section restructuring procedure；

4) performance test

In test, the ART-R-TV algorithm of the present invention and the restructuring procedure of traditional ART-TV algorithm are all at 18,36 With carry out under 72 projection angles.The knot of the meansigma methods of 15 iteration of MSE and PSNR of the reconstruct image obtained in test Fruit is the most as shown in table 1.And the PSNR of two kinds of algorithms change curve in an iterative process is compared as shown in Figure 4.

MSE and the PSNR contrast of table 1, image 1 and the image 2 reconstruction result under three projection angles

As can be seen from Table 1, two width images utilize the reconstruction quality of algorithm of the present invention to be above passing in restructuring procedure System algorithm.Under 18,36 and 72 projection angles during iteration 15 times, the average raising of PSNR value of the reconstruct image of image 1 2.93db, 3.36db and 3.68db, the PSNR value of the reconstruct image of image 2 the most averagely improve 3.77db, 4.14db and 4.34db.Meanwhile, the MSE contrast of the reconstruct image of two width images obtains the algorithm of the present invention and puts down compared to the MSE of traditional algorithm Average is less and reconstruction accuracy is higher.

Contrast by the PSNR situation of change of Fig. 4, it can be seen that along with increasing of projection angle, the reconstruct figure of two kinds of algorithms The PSNR value of picture all can increase.And the PSNR of the reconstruct image of tradition ART-TV algorithm improves along with the increase of iterations Fluctuating unstable, growth curve mostly also is broken line.The ART-R-TV algorithm of the present invention reconstruct image PSNR with It is more stable that iterations increases smooth growth, it is ensured that restructuring procedure is stable at an average level, makes algorithm robustness more By force.Meanwhile, when same projection angle number is issued to identical PSNR, the iterations required for the algorithm of the present invention will be less than passing The ART-TV algorithm of system.

Table 2, two images reconstruct image under 18,36 and 72 projection angles respectively

Claims (1)

1. ART flame based on a gradient total variation section restructing algorithm, it is characterised in that this algorithm includes following step Rapid:

Step (1), flame sectioning image is carried out pretreatment: include greyscale transformation and smoothing denoising operation；

Step (2), the iteration weight matrix W trying to achieve the reconstruction of flame sectioning image acquisition corresponding projection image array P, specifically Step is as follows:

By the two dimension discrete volume elements independent for N=n × n of flame section, each volume elements has that specifically to represent this unit bright Value x of degree size_j(j=1,2 ..., N), corresponding one of each pixel in two-dimensional projection image passes this section volume elements Ray L_k(k=1,2 ..., M), M is the ray sum of all projections, equal to projection angle number and certain Angles Projections ray The product of number；And value P of this pixel_k(k=1,2 ..., M) it is that on this ray, this ray brightness is contributed by all volume elements Superposition；If ω_jkFor volume elements j to ray L_kWeight factor, represent volume elements and interradial dependency, volume elements j is to ray L_k Contribution be:

P_jk=ω_jk·x_j

In section, all volume elements are to L_kBrightness contribution and being expressed as:

$P_k=\underset{j=1}{\overset{N}{\Σ}}{P}_{jk}=\underset{j=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_{jk}\·x_j,k=1,2,...,M$

It is expressed as follows by the form of matrix:

Wx=P

P is that M ties up projection matrix, represents the pixel value of two-dimensional projection；X is that N-dimensional treats reformatted slices data vector；W is that M × N-dimensional is thrown Shadow weight matrix；

Step (3), setting primary iteration number of times i=1 and the maximum Maxcount=15 of iterations, and arrange at the beginning of iteration Beginning vector f⁽⁰⁾=0,0 ..., 0}；

Step (4), calculating reconstruct image, equation group formula is as follows:

$f^{\left(m\right)}=f^{(m-1)}+{\mathrm{\ω}}_{jk}(P_k-\underset{j=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_{jk}{f}^{(m-1)})/\underset{j=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_{jk}^2$

Wherein m=1,2 ..., M represents the equation number of the equation group of complete iteration each time；According to above-mentioned equation group formula, meter Calculate M × N equation group；By initial vector f⁽⁰⁾=0,0 ..., 0} substitutes into first equation of equation group and obtains f⁽¹⁾, by second Equation obtains f⁽²⁾, the like, finally tried to achieve f by m-th equation^(M), this completes a complete iteration；

Step (5), solve f^(M)GradientFormula is as follows:

$\▿f^{\left(M\right)}=(D_{\mathrm{\ρ}}{x}_{\mathrm{\ρ},\mathrm{\θ}},D_{\mathrm{\θ}}{x}_{\mathrm{\ρ},\mathrm{\θ}})$

Wherein D_ρAnd D_θIt is the discrete differential operator along footpath, pole and polar angle direction respectively, has equation below to calculate:

D_ρx_ρ,θ=x_ρ,θ-x_ρ-,θ

D_θx_ρ,θ=x_ρ,θ-x_ρ,θ-

Wherein ρ_-And ρ₊Represent respectively and reduce along footpath, pole and the change in increase direction, footpath, pole, θ_-And θ₊Represent the most respectively along polar angle up time Pin and the anticlockwise change of polar angle；

Step (6), solve f^(M)The partial derivative of gradient total variationAnd to f^(M)Carry out gradient complete It is deteriorated and optimizes, such as following formula:

$f^{\left(M\right)}\′=f^{\left(M\right)}-q\·\▿{||f^{\left(M\right)}||}_{R-TV}$

Wherein, q is for declining step factor, and in the present invention, value is 0.1/k, and uses f^(M)' replace original f^(M)；

Step (7) if i≤Maxcount, by f^(M)' carry out complete iteration next time according to the equation group of step (4), otherwise, Stop iteration, output reconstruct image.