CN104537611A - Optimized light intensity transmission phase recovery method based on Savitzky-Golay differential filters - Google Patents

Abstract

The invention discloses an optimized light intensity transmission phase recovery method based on Savitzky-Golay differential filters. The method comprises the following steps: acquiring the distribution of an object to be evaluated on a plurality of measuring planes through experiments, and calculating the light intensity axial differential through Savitzky-Golay differential filters of different orders; solving a light intensity transmission equation through the obtained light intensity axial differential and light intensity distribution on a focusing surface, and solving a group of phase distributions corresponding to the Savitzky-Golay differential filters of different orders; and decomposing the obtained phase distributions corresponding to the Savitzky-Golay differential filters of different orders by adopting a complementary filter group to obtain phase distributions corresponding to the Savitzky-Golay differential filters of different orders after filtering, and summing the phase distributions to obtain a final phase distribution. Through adoption of the optimized light intensity transmission phase recovery method, the contradiction of difficulty in giving consideration to both low and middle-frequency clouding noise and high-order nonlinearity in a conventional method is solved effectively; the noise resistance of a light intensity transmission equation method is enhanced greatly; and the accuracy of phase reconstruction is increased greatly.

Description

Based on the optimization light intensity transmission phase recovery method of Savitzky-Golay difference filter

Technical field

The invention belongs to phase recovery and the quantitative phase imaging technique of optical measurement, particularly a kind of optimization shaven head based on Savitzky-Golay difference filter transmits phase recovery method by force.

Background technology

Phase recovery is an important topic of optical measurement and imaging technique, and no matter in biomedical or field of industry detection, phase imaging technology is all playing an important role.Make a general survey of the progress of optical measurement nearly half a century, the most classical Method for Phase Difference Measurement should not belong to by non-interfering mensuration.But the shortcoming of interferometry is also fairly obvious: interferometry generally needs the light source (as laser) of high coherence, thus need comparatively complicated interference device; The introducing of extra reference path causes the requirement for measurement environment to become very harsh; The speckle coherent noise of the light source introducing of high coherence limits spatial resolution and the measuring accuracy of imaging system.

Difference and interferometry, another kind of very important phase measurement does not need by interference, and they are referred to as phase recovery.Because the PHASE DISTRIBUTION directly measuring light wave fields is very difficult, and the amplitude/intensity measuring light wave fields is very easy.Therefore, can be thought of as being recovered by intensity distributions this process of (estimation) phase place one mathematical " inverse problem ", i.e. phase retrieval problem.Phase recovery method also can be subdivided into process of iteration and direct method.Light intensity transmission equation method is the typical direct method of one in phase recovery method.Light intensity transmission equation is a Some Second Order Elliptic partial differential equation, that illustrates the quantitative relationship of the phase place of light wave in the variable quantity of light intensity is vertical with optical axis on optical axis direction plane.When the axial differential of light intensity and light distribution known, directly can obtain phase information by numerical solution light intensity transmission equation.Compare and interferometric method and iterative phase restoring method, its major advantage comprises: (1) non-interfering, only by measurement object plane light intensity direct solution phase information, does not need to introduce additional reference light; (2) non-iterative, obtains phase place by the direct solution differential equation; (3) can well white-light illuminating be applied to, as in traditional light field microscope kohler's illumination ( ); (4) without the need to Phase-un-wrapping, directly obtain the absolute profile of phase place, there are not 2 π phase place parcel problems in general interferometry; (5) optical system that need not be complicated, does not have harsh requirement for experimental situation, vibrates insensitive.

The committed step solved in light intensity transmission equation is the axial differential obtaining light intensity, and this amount is not directly measured, and needs to be obtained by numerical finite difference.Classic method is by the slight out-of-focus image of collection two width, and the out of focus distance of this two width image is equal relative to centre focus image and direction contrary, then utilizes centered finite difference methods to solve the estimation of the axial differential obtaining light intensity.But the Gradient estimates precision obtained like this, even if in the absence of noise, also only can reach the second order order of magnitude of out of focus distance.In this external actual measurement, the existence of noise is inevitable.Research shows, the practical manifestation of light intensity transmission equation phase recovery depends on the selection of out of focus distance in centered finite difference methods strongly.When out of focus distance was selected not at that time, result there will be very strong nebulous low-frequency noise.In order to improve the signal to noise ratio (S/N ratio) of measurement result, out of focus distance must be increased, so just will inevitably break the hypothesis of method of finite difference local linear, namely exchange the raising of signal to noise ratio (S/N ratio) with the cost of sacrificing estimated accuracy for.In order to address this problem, many researchers propose to adopt the ionization meter of multiaspect (>2) to estimate axial differential, so that correction of Nonlinear error or reduce the impact of noise more accurately.Although these methods have good effect under some particular condition, their performance depends on the grade of noise and the spatial frequency characteristic of object being measured very much.When given one group of light intensity view data, choose a kind of optimal method very difficult often.

Comparatively detailed background introduction is carried out below by for this problem.

Consider that a monochrome propagated along z-axis is concerned with paraxial light wave fields, its complex amplitude U (r) is wherein j is imaginary unit, the PHASE DISTRIBUTION of φ (r) for recovering.Light intensity transmission equation can be expressed as

$-k\frac{\∂I\left(r\right)}{\∂z}=\▿\·[I\left(r\right)\▿\mathrm{\φ}\left(r\right)]---\left(1\right)$

Wherein k is wave number 2 π/λ, λ is adopted optical wavelength, and r is lateral attitude coordinate (x, y).I (r) is the light distribution on focusing surface, thinks that it is positioned at the plane of z=0 without loss of generality.▽ is the Hamiltonian operator acting on r plane, for the axial differential of light intensity.Under normal circumstances, remove based on biplanar finite-difference formula axial differential ([1] M.Beleggia estimating light intensity by following, M.A.Schofield, V.V.Volkov, and Y.Zhu, " On the transport of intensity technique for phase retrieval, " Ultramicroscopy 102, 37-49 (2004). [2] D.Paganin, A.Barty, P.J.McMahon, and K.A.Nugent, " Quantitative phase-amplitude microscopy.III.The effects of noise, " Journal of Microscopy 214, 51-61 (2004) .)

$\frac{\∂I\left(r\right)}{\∂z}\≈\frac{I(r,\mathrm{\Δz})-I(r,-\mathrm{\Δz})}{2\mathrm{\Δz}}---\left(2\right)$

Note function light intensity function explicitly being expressed as out of focus distance, delta z here.Difference is adopted to carry out approximate differential based on biplanar finite-difference formula (2), so there is a problem: namely choosing difference result that much out of focus distance, delta z obtain, to carry out approximate differential be rational.Ideally muting, the answer of this problem is apparent.Here first qualitative explanation is carried out with some simple examples.As shown in Fig. 1 (a) to Fig. 1 (c), when defocusing amount Δ z more hour, the precision of this Difference formula is higher, and the spatial resolution of the phase place of correspondingly rebuilding is also higher.But when Δ z becomes greatly, the accuracy of the Difference formula of formula (2) will decline, and shows the decline of the spatial resolution of rebuilding phase place, has namely occurred " phase ambiguity " phenomenon, as shown in Fig. 1 (c) thereupon." phase ambiguity " is sometimes otherwise known as " nonlinearity erron ", this is because formula (2) itself utilizes local linear to be similar to approach the axial differential of light intensity, and the increase of Δ z directly results in the rising of nonlinear terms generation error in actual signal.So say from mathematical angle, Δ z should be reduced as much as possible, thus improve the precision of Difference formula.

But the factors such as the quantization effect of the noise existed in actual measurement, detector make this light intensity differential estimation problem complicated, as shown in Fig. 2 (a) to Fig. 2 (b).When Noise, it is too little that Δ z can not get, otherwise light intensity differential estimate will flooded by noise, as shown in Fig. 2 (a).Serious cloud low-frequency noise is just there will be in the phase place of now rebuilding.In order to improve signal to noise ratio (S/N ratio), have to increase Δ z, as shown in Fig. 2 (b).But Δ z get excessive time, phase ambiguity effect is also just obvious all the more.So in the presence of noise, must make compromise between noise and nonlinearity erron.The out of focus distance of this optimum must be closely-related with the feature of noise level and object itself.Can imagine, under certain conditions nonlinearity erron and low frequency aberration all less (now measuring error is minimum), distance is at this moment best out of focus distance.But the priori of the spatial frequency these two aspects of noise level and object is often all difficult to predict before measuring in actual conditions.In addition when object comparatively complicated (comprising High-frequency and low-frequency information) simultaneously; and noise comparatively serious time, often there will be a kind of implacable situation: the lower limit of the out of focus distance namely determined by the size of luminous intensity measurement noise can exceed the out of focus that determined by the most high spatial frequency of object apart from the upper limit.Briefly, phase low-frequency composition is affected by noise comparatively large, but insensitive to nonlinearity erron, therefore tends to choose larger out of focus distance.On the contrary, phase place high-frequency information is comparatively large by the impact of nonlinear effect, but has higher signal to noise ratio (S/N ratio), therefore tends to choose less out of focus distance.The contradiction of the two is often irreconcilable, needs to have to make certain " compromise ", allows error reduce as much as possible.But at this kind in particular cases, in any case select defocusing amount, all can not carry out accurately recovering to object phase based on biplanar axial differential method of estimation.

Consider the intensity signal only adopting two planes in biplanar axial differential method of estimation, wherein unique controllable parameter is exactly out of focus distance.In order to solve the problem, many researchers propose to adopt the ionization meter of multiaspect (>2) to estimate axial differential, are exactly in order to correction of Nonlinear error more neatly or the impact reducing noise.High-Order Finite-Difference Method (higher order finite difference) is by ([3] L.Waller such as Waller, L.Tian, and G.Barbastathis, " Transport of Intensity phase-amplitude imaging with higher order intensity derivatives; " Opt.Express 18,12552-12561 (2010) .) propose.For the situation of central difference, i.e. the light intensity I (r, i Δ z) of given multiple measurement plane, i=-n ..., 0 ..., n, the difference formula of High-Order Finite-Difference Method can be expressed as

$\frac{\∂I\left(r\right)}{\∂z}\≈\underset{i=-n}{\overset{n}{\mathrm{\Σ}}}\frac{a_{i}I(r,\mathrm{i\Δz})}{\mathrm{\Δz}}---\left(3\right)$

Wherein represent that coefficient is

$a_i=\frac{{(-1)}^{i+1}{(n!)}^2}{i(n+i)!(n-i)!}---\left(4\right)$

The core concept of High-Order Finite-Difference Method goes to approach the high-order Taylor expansion item in the axial differential of light intensity better by increasing luminous intensity measurement, for the luminous intensity measurement adopting 2n plane, if do not consider the impact of noise, the evaluated error of formula (3) is O (Δ z ²ⁿ) order.But when noise is deposited in case, High-Order Finite-Difference Method is often very responsive to noise.To disinthibite the impact of noise to more effectively utilize the luminous intensity measurement of multiple plane, Soto and Acosta ([4] M.Soto and E.Acosta, " Improved phase imaging from intensity measurements in multiple planes; " Appl.Opt.46,7978-7981 (2007) .) propose squelch method of finite difference (noise-reduction finite difference), its difference formula is still formula (3), but represents that coefficient is different:

$a_i=\frac{3i}{n(n+1)(2n+1)}---\left(5\right)$

This coefficient is the higher order term by ignoring in all Taylor expansion items, and minimum noise effect is derived and obtained.Contrast (4) can find with the coefficient of (5), and in High-Order Finite-Difference Method, the coefficient the closer to the light intensity of central plane will be large as much as possible, and this is to go estimation to approach higher order term in the axial differential of light intensity better.And the object of squelch method of finite difference is the quadratic sum minimizing coefficient, namely wish the light intensity the closer to central plane, coefficient is less.So there is contradiction in some sense in these two kinds of methods.In order to be in harmonious proportion this contradiction, ([5] R.Bie such as Bie, X.-H.Yuan, M.Zhao, and L.Zhang, " Method for estimating the axial intensity derivative in the TIE with higher order intensity derivatives and noise suppression; " Opt.Express 20,8186-8191 (2012) .) above-mentioned two kinds of methods are combined, propose squelch High-Order Finite-Difference Method, consider the impact of higher order term and noise in the Taylor expansion of noise simultaneously.For the luminous intensity measurement adopting 2n plane, higher order term is only considered m rank (m < 2n+1) by squelch High-Order Finite-Difference Method, so just can supply to optimize for squelch (reduction noise suppression factor) leaves " degree of freedom ".Except these are based on the method for finite difference, ([3] L.Waller such as Waller, L.Tian, and G.Barbastathis, " Transport of Intensity phase-amplitude imaging with higher order intensity derivatives; " Opt.Express 18,12552-12561 (2010) .) also advise adopting least square fitting method to process this problem, regard a discrete point range as by the light intensity measurement of each pixel in Different Plane, then by least square method, it is carried out curve fitting.By adjusting the exponent number of least square fitting, the change curve in light intensity axial area can be simulated as far as possible exactly, so just can suppress high-order error and noise effect simultaneously.Owing to finally can be obtained the analytical expression of the light intensity out of focus distance of each position by matching, light intensity also just can calculate easily at the axial differential value of center.

Above-mentioned four kinds of methods are the typical methods estimated based on the axial differential of multilevel light intensity.In general, increase compared to X-plane in high-order error and squelch based on multilevel method, under some particular condition, also obtain good effect, but some researchist has been found that the practical manifestation of above-mentioned many planar approach also depends on the spatial frequency characteristic of noise level and object being measured very much.When given one group of intensity data, choose a kind of optimal method very difficult often.For squelch High-Order Finite-Difference Method and least square fitting method, also need the order determining to approach (matching), this seems exactly the same with the problem of how best defocusing amount in biplane.Therefore there is the contradiction that low frequency cloud noise and high-order nonlinear are difficult to take into account simultaneously in classic method, the poor accuracy of its Phase Build Out.

Summary of the invention

A kind of optimization light intensity based on Savitzky-Golay difference filter is the object of the present invention is to provide to transmit phase recovery method, the contradiction that low frequency cloud noise and high-order nonlinear are difficult to take into account simultaneously can be efficiently solved, greatly improve the noise immunity of light intensity transmission equation method and the accuracy of Phase Build Out.

The technical solution realizing the object of the invention is: a kind of transmission of the optimization light intensity based on Savitzky-Golay difference filter phase recovery method, and step is as follows:

The first step, gathers the light distribution I of object to be asked on multiple measurement plane (r, i Δ z), i=-n ..., 0 ..., n, wherein Δ z represents the direct interval of these measurement planes, and n is data half-breadth;

Second step, utilizes the Savitzky-Golay difference filter of different order, calculates the axial differential of light intensity, namely press formula (6), utilize m=1 respectively, 3, ... the Savitzky-Golay difference filter of 2n-1 order, calculates the axial differential of respective light intensity

$\frac{\&PartialD;I\left(r\right)}{\&PartialD;z}\&ap;\underset{i=-n}{\overset{n}{\mathrm{\&Sigma;}}}\frac{a_{i}I(r,\mathrm{i\&Delta;z})}{\mathrm{\&Delta;z}}---\left(6\right)$

Wherein a _ifor filter coefficient, for the Savitzky-Golay difference filter on m rank, its coefficient a _ibe expressed as

$a_i=\underset{k=0}{\overset{m}{\mathrm{\&Sigma;}}}\frac{(2k+1){\left(2n\right)}^{\left(k\right)}}{{(2n+k+1)}^{(k+1)}}{P}_k^n\left(i\right)P_k^{n,1}\left(0\right)---\left(7\right)$

Wherein (a) ^(b)generalized factorial function, it is defined as (a) (a-1) ... (a-b+1), and (a) ⁽⁰⁾=1; be Ge Lan polynomial expression, it is defined as

$P_k^n\left(i\right)=\underset{j=0}{\overset{k}{\mathrm{\&Sigma;}}}\frac{{(-1)}^{j+k}{(j+k)}^{\left(2j\right)}{(n+i)}^{\left(j\right)}}{{(j!)}^2{\left(2n\right)}^{\left(j\right)}}---\left(8\right)$

Its s rank differential is defined as

$P_k^{n,s}\left(t\right)={\left(\frac{d^s}{d{x}^s}{P}_k^n\left(t\right)\right)}_{x=t}---\left(9\right)$

Obtain the axial differential of one group of light intensity like this m=1,3 ... 2n-1, they are by respectively

M=1,3 ... the Savitzky-Golay difference filter of 2n-1 order to estimate gained;

3rd step, the axial differential of the light intensity obtained by second step and the light distribution on focusing surface, solve light intensity transmission equation, finally solves and obtain one group of PHASE DISTRIBUTION corresponding to the Savitzky-Golay difference filter of different order;

4th step, the PHASE DISTRIBUTION corresponding to the Savitzky-Golay difference filter of different order 3rd step obtained adopts complementary filter group to decompose, the filtered PHASE DISTRIBUTION φ ' of Savitzky-Golay difference filter corresponding to different order obtained _m(r), m=1,3 ... 2n-1;

5th step, the filtered PHASE DISTRIBUTION φ ' of Savitzky-Golay difference filter corresponding to different order that the 4th step is obtained _m(r), m=1,3 ... 2n-1 summation finally PHASE DISTRIBUTION φ (r) to be asked

φ(r)＝φ′ ₁(r)+φ′ ₃(r)+...+φ′ _2n-1(r)。(19)

The present invention compared with prior art, its remarkable advantage: the present invention efficiently solves the contradiction that classic method medium and low frequency cloud noise and high-order nonlinear are difficult to take into account simultaneously, avoid the problem being difficult to determine filter order in Fixed-order Savitzky-Golay wave filter, the method carries out multistage subdifferential estimation based on Savitzky-Golay wave filter, the mode that frequency domain optimal filtering, spatial domain are recombinated again, improves the noise immunity of light intensity transmission equation method and the accuracy of Phase Build Out greatly.In addition, the present invention itself has, accuracy advantages of higher little without the need to manual adjustments parameter, calculated amount, the fields such as quantitative phase imaging, phase place be micro-that make it be suitable for very much.

Below in conjunction with accompanying drawing, the present invention is described in further detail.

Accompanying drawing explanation

Fig. 1 (a) be under noise-free case little out of focus apart under original light intensity image (left side), by two width image differences estimate the axial differential of the light intensity that obtains (in), and the PHASE DISTRIBUTION (right side) recovered.

Fig. 1 (b) be under noise-free case medium out of focus apart under original light intensity image (left side), by two width image differences estimate the axial differential of the light intensity that obtains (in), and the PHASE DISTRIBUTION (right side) recovered.

Fig. 1 (c) be under noise-free case large out of focus apart under original light intensity image (left side), by two width image differences estimate the axial differential of the light intensity that obtains (in), and the PHASE DISTRIBUTION (right side) recovered.

Fig. 2 (a) be in Noise situation little out of focus apart under original light intensity image (left side), by two width image differences estimate the axial differential of the light intensity that obtains (in), and the PHASE DISTRIBUTION (right side) recovered.

Fig. 2 (b) be in Noise situation medium out of focus apart under original light intensity image (left side), by two width image differences estimate the axial differential of the light intensity that obtains (in), and the PHASE DISTRIBUTION (right side) recovered.

Fig. 3 is based on the optimization light intensity transmission phase recovery method process flow diagram of Savitzky-Golay difference filter.

Fig. 4 (a) is experiment light path schematic diagram.

Fig. 4 (b) is the plot of light intensity collected.

Fig. 5 (a) be traditional double planar process recover the PHASE DISTRIBUTION that obtains, two interplanar spacings are 100 μm.

Fig. 5 (b) is the PHASE DISTRIBUTION that the inventive method is recovered to obtain.

Fig. 6 (a) is the hologram of identical object under test.

Fig. 6 (b) is the Fourier spectrum of hologram.

Fig. 6 (c) is that digital hologram rebuilds the PHASE DISTRIBUTION (in frame, part corresponds to the reconstructed results of method proposed by the invention, i.e. the corresponding part of Fig. 5 (b)) obtained.

Fig. 6 (d) phase place xsect contrasts, the line mark region respectively in corresponding diagram 6 (c).

Embodiment

Composition graphs 3, the present invention is based on the optimization light intensity transmission phase recovery method of Savitzky-Golay difference filter, concrete implementation step is as follows:

The first step, gathers the light distribution I of object to be asked on multiple measurement plane (r, i Δ z), i=-n ..., 0 ..., n, wherein Δ z represents the direct interval of these measurement planes, and n is data half-breadth.Harvester can adopt (but being not limited to) as Suo Shi Fig. 4 (a) 4f imaging system experimental provision, light source is by irradiating object under test after beam-expanding collimation.Object passes through a 4f system imaging in CCD plane.The lens that this 4f system is f by two focal lengths form, and are spaced apart 2f between two lens, and the distance of object distance first lens is f, and the distance of second lens distance CCD is also f.CCD camera is positioned in a horizontal moving stage, by mobile CCD, just can gather object to be asked on multiple measurement plane light distribution.

Second step, utilizes the Savitzky-Golay difference filter of different order, calculates the axial differential of light intensity, namely press formula (6), utilize m=1 respectively, 3, ... the Savitzky-Golay difference filter of 2n-1 order, calculates the axial differential of respective light intensity

$\frac{\&PartialD;I\left(r\right)}{\&PartialD;z}\&ap;\underset{i=-n}{\overset{n}{\mathrm{\&Sigma;}}}\frac{a_{i}I(r,\mathrm{i\&Delta;z})}{\mathrm{\&Delta;z}}---\left(6\right)$

Wherein a _ifor filter coefficient, for the Savitzky-Golay difference filter on m rank, its coefficient a _ican be expressed as

$a_i=\underset{k=0}{\overset{m}{\mathrm{\&Sigma;}}}\frac{(2k+1){\left(2n\right)}^{\left(k\right)}}{{(2n+k+1)}^{(k+1)}}{P}_k^n\left(i\right)P_k^{n,1}\left(0\right)---\left(7\right)$

Wherein (a) ^(b)generalized factorial function, it is defined as (a) (a-1) ... (a-b+1), and (a) ⁽⁰⁾=1. be Ge Lan polynomial expression, it is defined as

$P_k^n\left(i\right)=\underset{j=0}{\overset{k}{\mathrm{\&Sigma;}}}\frac{{(-1)}^{j+k}{(j+k)}^{\left(2j\right)}{(n+i)}^{\left(j\right)}}{{(j!)}^2{\left(2n\right)}^{\left(j\right)}}---\left(8\right)$

Its s rank differential is defined as

$P_k^{n,s}\left(t\right)={\left(\frac{d^s}{d{x}^s}{P}_k^n\left(t\right)\right)}_{x=t}---\left(9\right)$

Obtain the axial differential of one group of light intensity like this m=1,3 ... 2n-1, they are by m=1 respectively, 3 ... the Savitzky-Golay difference filter of 2n-1 order to estimate gained.

3rd step, obtains second step, by different order m=1, and 3 ... the Savitzky-Golay difference filter of 2n-1 estimates the axial differential of each light intensity obtained m=1,3 ... 2n-1, with light distribution I (r) on focusing surface=I (r, 0), utilizes formula (10) respectively, adopts Fast Fourier Transform (FFT) to solve light intensity transmission equation

$\mathrm{\&phi;}(x,y)=-k{\&dtri;}^{-2}\&dtri;\&CenterDot;\left[I^{-1}(x,y){\&dtri;\&dtri;}^{-2}\frac{\&PartialD;I(x,y)}{\&PartialD;z}\right]---\left(10\right)$

In formula inverse Laplace's operation symbol, for gradient operator, be vector dot, k is wave number, with operational symbol is all realized by Fourier transform, namely

${\&dtri;}^{-2}\{\&CenterDot;\}{F}^{-1}\left\{F\right\{\&CenterDot;\left\}\frac{1}{-4{\mathrm{\&pi;}}^2(u^2+v^2)}\right\}---\left(11\right)$

$\&dtri;\{\&CenterDot;\}=F^{-1}\left\{j2\mathrm{\&pi;uF}\right\{\&CenterDot;\},j2\mathrm{\&pi;vF}\{\&CenterDot;\left\}\right\}---\left(12\right)$

Wherein F represents Fourier transform, and (u, v) is the frequency domain coordinates corresponding with volume coordinate (x, y), and j is imaginary unit.As shown in Figure 3, finally solving what obtain is one group of PHASE DISTRIBUTION φ _m(x), m=1,3 ... 2n-1, it corresponds respectively to different order m=1, and 3 ... the Savitzky-Golay difference filter of 2n-1.

4th step, the PHASE DISTRIBUTION φ corresponding to the Savitzky-Golay difference filter of different order that the 3rd step is obtained _m(r), m=1,3 ... 2n-1 utilizes corresponding complementary filter group H respectively _m(e ^{j ω}), m=1,3 ... 2n-1 carries out filtering (wherein ω=Δ z λ (u ²+ v ²) represent normalization polar coordinates in frequency domain), obtain the filtered PHASE DISTRIBUTION φ ' of Savitzky-Golay difference filter that a group corresponds to different order _m(r), m=1,3 ... 2n-1.First PHASE DISTRIBUTION, such as formula shown in (13), is namely carried out Fourier transform, is then multiplied by H in a frequency domain by detailed process _m(e ^{j ω}), m=1,3 ... 2n-1, then inverse transformation is made the return trip empty territory:

φ′ _m(r)＝F ^-1{F{φ _m(r)}H _m(e ^jω)} (13)

Here the form of complementary filter group is that (form of complementary filter is not limited only to desirable 0-1 wave filter to desirable 0-1 wave filter here, also Butterworth filter can be adopted similarly, Gaussian filters etc., only adopt desirable 0-1 wave filter to be illustrated here), their frequency response H _m(e ^{j ω}), m=1,3 ... the sum total of 2n-1 mould is 1 (passing through original signal completely).For order m=3,5 ..., 2n-3, filters H _m(e ^{j ω}) be bandpass filter, its free transmission range is f _c ^m-2to f _c ^mbetween, namely

For order m=1, filters H ₁(e ^{j ω}) be low-pass filter, its cutoff frequency is f _c ¹, namely

For order m=2n-1, filters H _2n-1(e ^{j ω}) be a Hi-pass filter, its cutoff frequency is f _c ^2n-3, namely

F above in formula _c ^mrepresent the cutoff frequency of the 0.3dB point of m rank Savitzky-Golay difference filter.The concrete numerical value of this cutoff frequency realizes in the following way:

(1) if data half-breadth n is between 1-10, then f _c ^mcan directly find by reference table 1, such as: as data half-breadth n=5, for the filter freguency response H of order m=3 ₃(e ^{j ω}) cutoff frequency f _c ³tabling look-up can be 0.126.

(2) for larger n>=25 of data half-breadth, f _c ^mcan be calculated by formula (17).

${f_c}^m=\frac{m-0.3}{3.5n-1},\&GreaterEqual;25---\left(17\right)$

(3) time for data half-breadth 10 < n < 25, f _c ^mcan be calculated by formula (18).

${f_c}^m=\frac{m-0.3}{3.5n-2.5},10<n<25---\left(18\right)$

After this step, obtain the filtered PHASE DISTRIBUTION φ ' of Savitzky-Golay difference filter corresponding to different order _m(r), m=1,3 ... 2n-1.

5th step, the filtered PHASE DISTRIBUTION φ ' of Savitzky-Golay difference filter corresponding to different order that the 4th step is obtained _m(r), m=1,3 ... 2n-1 summation finally PHASE DISTRIBUTION φ (r) to be asked

φ(r)＝φ′ ₁(r)+φ′ ₃(r)+...+φ′ _2n-1(r) (19)

Comprehensive above-mentioned steps is known, and the present invention carries out frequency domain decomposition by the reconstruction phase place of the Savitzky-Golay difference filter to different order, optimal filtering, the mode of recombinating again carry out optimum angle recovery.Owing to combining the optimal frequency composition of each order Savitzky-Golay difference filter, it can efficiently solve the contradiction that classic method medium and low frequency cloud noise and high-order nonlinear are difficult to take into account simultaneously, avoid classic method, namely be difficult to the problem determining filter order in Fixed-order Savitzky-Golay wave filter, thus greatly can improve the noise immunity of light intensity transmission equation and the accuracy of Phase Build Out.

The 0.3dB normalized frequency (n=1 ~ 10) of the Savitzky-Golay difference filter under table 1 different order m and data half-breadth n

Embodiment

The validity of the transmission of the optimization light intensity based on Savitzky-Golay difference filter phase recovery method proposed by the invention is described by embodiment below.Experimental provision schematic diagram is as shown in Fig. 4 (a), and He-Ne laser light source (λ=632.8nm) is by irradiating object under test after beam-expanding collimation.Object passes through a 4f system imaging in CCD plane.The lens that this 4f system is f=25mm by two focal lengths form, and are spaced apart 2f between two lens, and the distance of object distance first lens is f.CCD camera is positioned in a horizontal moving stage, to be used for changing different out of focus distances.Sample in experiment is etched in a geometric scheme on transparent organic glass sheet and experiment light path, as shown in Fig. 4 (a).Acquire 51 width out-of-focus images by mobile CCD, out of focus scope is-1.250mm to 1.250mm, and the distance of two width adjacent images is spaced apart 50 μm.The present invention have chosen wherein some plot of light intensity pictures as shown in Fig. 4 (b).These images are collected (The Imaging Source DMK 41AU02, pixel dimension 4.65 μm) by a monochromatic CCD.Experimental result contrasts, and determinand is a phase object: Fig. 5 (a) is traditional double planar process, and two interplanar spacings are 100 μm; Fig. 5 (b) is the inventive method.Wherein every width figure bottom-right graph picture is the amplification display of homologous thread boxed area.As can be seen from experimental result, by the present invention the PHASE DISTRIBUTION recovering to obtain compared to traditional double planar process, no matter be obtained for very large lifting in signal to noise ratio (S/N ratio) or in the resolution of rebuilding phase place.It has not only accurately recovered the minor detail on sample, and ground unrest have also been obtained and suppresses preferably.

In order to quantitative test method proposed by the invention better recover the precision of phase place, adopt the digital holographic microscope (4x object lens, NA=0.13) of a Michelson structure to go to measure identical sample.Fig. 6 (a) and Fig. 6 (b) sets forth the hologram and Fourier spectrum thereof that record and obtain.Fig. 6 (c) finally rebuilds the PHASE DISTRIBUTION (frame portion corresponds to the reconstructed results of method proposed by the invention, i.e. the corresponding part of Fig. 5 (b)) obtained.Result with optimal frequency back-and-forth method holographic in order to more clearly comparative figures, we have chosen the phase place xsect being positioned at object edge place and contrast, as shown in Fig. 6 (d), the line mark region respectively in corresponding diagram 6 (c).Can find out that the two is coincidently fine: the inventive method has not only reconstructed larger phase hit exactly, also recovered object detail of the high frequency (as due to etch process imperfection the undulatory vertical groove that causes).

Claims (5)

1., based on an optimization light intensity transmission phase recovery method for Savitzky-Golay difference filter, it is characterized in that step is as follows:

The first step, gathers the light distribution I of object to be asked on multiple measurement plane (r, i Δ z), i=-n ..., 0 ..., n, wherein Δ z represents the direct interval of these measurement planes, and n is data half-breadth;

Second step, utilizes the Savitzky-Golay difference filter of different order, calculates the axial differential of light intensity, namely press formula (6), utilize m=1 respectively, 3, ... the Savitzky-Golay difference filter of 2n-1 order, calculates the axial differential of respective light intensity

$\frac{\&PartialD;I\left(r\right)}{\&PartialD;z}\&ap;\underset{i=-n}{\overset{n}{\mathrm{\&Sigma;}}}\frac{a_{i}I(r,\mathrm{i\&Delta;z})}{\mathrm{\&Delta;z}}---\left(6\right)$

Wherein a _ifor filter coefficient, for the Savitzky-Golay difference filter on m rank, its coefficient a _ibe expressed as

$a_i=\underset{k=0}{\overset{m}{\mathrm{\&Sigma;}}}\frac{(2k+1){\left(2n\right)}^{\left(k\right)}}{{(2n+k+1)}^{(k+1)}}{P}_k^n\left(i\right)P_k^{n,1}\left(0\right)---\left(7\right)$

Wherein (a) ^(b)generalized factorial function, it is defined as (a) (a-1) ... (a-b+1), and (a) ⁽⁰⁾=1; be Ge Lan polynomial expression, it is defined as

$P_k^n\left(i\right)=\underset{j=0}{\overset{k}{\mathrm{\&Sigma;}}}\frac{{(-1)}^{j+k}{(j+k)}^{\left(2j\right)}{(n+i)}^{\left(j\right)}}{{(j!)}^2{\left(2n\right)}^{\left(j\right)}}---\left(8\right)$

Its s rank differential is defined as

$P_k^{n,s}\left(t\right)={\left(\frac{d^s}{{\mathrm{dx}}^s}{P}_k^n\left(t\right)\right)}_{x=t}---\left(9\right)$

Obtain the axial differential of one group of light intensity like this m=1,3 ... 2n-1, they are by m=1 respectively, 3 ... the Savitzky-Golay difference filter of 2n-1 order to estimate gained;

3rd step, the axial differential of the light intensity obtained by second step and the light distribution on focusing surface, solve light intensity transmission equation, finally solves and obtain one group of PHASE DISTRIBUTION corresponding to the Savitzky-Golay difference filter of different order;

4th step, the PHASE DISTRIBUTION corresponding to the Savitzky-Golay difference filter of different order 3rd step obtained adopts complementary filter group to decompose, the filtered PHASE DISTRIBUTION φ ' of Savitzky-Golay difference filter corresponding to different order obtained _m(r), m=1,3 ... 2n-1;

5th step, the filtered PHASE DISTRIBUTION φ ' of Savitzky-Golay difference filter corresponding to different order that the 4th step is obtained _m(r), m=1,3 ... 2n-1 summation finally PHASE DISTRIBUTION φ (r) to be asked

φ(r)＝φ′ ₁(r)+φ′ ₃(r)+...+φ′ _2n-1(r)。(19)

2. the transmission of the optimization light intensity based on Savitzky-Golay difference filter phase recovery method according to claim 1 and 2, is characterized in that in the third step, the axial differential of each light intensity obtained by second step with light distribution I (r) on focusing surface=I (r, 0), formula (10) is utilized to adopt Fast Fourier Transform (FFT) to solve light intensity transmission equation respectively

$\mathrm{\&phi;}(x,y)=-k{\&dtri;}^{-2}\&dtri;\&CenterDot;\left[I^{-1}(x,y){\&dtri;\&dtri;}^2\frac{\&PartialD;I(x,y)}{\&PartialD;z}\right]---\left(10\right)$

▽ in formula ^-2be inverse Laplace's operation symbol, ▽ is gradient operator, and be vector dot, k is wave number, ▽ and ▽ ^-2operational symbol is all realized by Fourier transform, namely

${\&dtri;}^{-2}\{\&CenterDot;\}=F^{-1}\left\{\begin{array}{cc}F& \{\&CenterDot;\}\frac{1}{-4{\mathrm{\&pi;}}^2(u^2+v^2)}\end{array}\right\}---\left(11\right)$

$\&dtri;\{\&CenterDot;\}{=F}^{-1}\left\{\begin{array}{ccc}j2\mathrm{\&pi;uF}& \{\&CenterDot;\},j2\mathrm{\&pi;vF}& \{\&CenterDot;\}\end{array}\right\}---\left(12\right)$

Finally solving what obtain is one group of PHASE DISTRIBUTION, φ _m(x), m=1,3 ... 2n-1, it corresponds respectively to different order m=1, and 3 ... the Savitzky-Golay difference filter of 2n-1, wherein F represents Fourier transform, (u, v) is the frequency domain coordinates corresponding with volume coordinate (x, y), and j is imaginary unit.

3. the transmission of the optimization light intensity based on Savitzky-Golay difference filter phase recovery method according to claim 1 and 2, it is characterized in that in the 4th step, the PHASE DISTRIBUTION φ corresponding to the Savitzky-Golay difference filter of different order that the 3rd step is obtained _m(r), m=1,3 ... 2n-1 utilizes corresponding complementary filter group H respectively _m(e ^{j ω}), m=1,3 ... 2n-1 carries out filtering (wherein ω=Δ z λ (u ²+ v ²) represent normalization polar coordinates in frequency domain), obtain the filtered PHASE DISTRIBUTION φ ' of Savitzky-Golay difference filter that a group corresponds to different order _m(r), m=1,3 ... 2n-1, first PHASE DISTRIBUTION, such as formula shown in (13), is namely carried out Fourier transform, is then multiplied by H in a frequency domain by detailed process _m(e ^{j ω}), m=1,3 ... 2n-1, then inverse transformation is made the return trip empty territory:

φ′ _m(r)＝F ^-1{F {φ _m(r)}H _m(e ^jω)} (13)

4. the transmission of the optimization light intensity based on Savitzky-Golay difference filter phase recovery method according to claim 3, is characterized in that the form of complementary filter group adopts desirable 0-1 wave filter, their frequency response H _m(e ^{j ω}), m=1,3 ... the sum total of 2n-1 mould is 1, passes through original signal completely; For order m=3,5 ..., 2n-3, filters H _m(e ^{j ω}) be bandpass filter, its free transmission range is f _c ^m-2to f _c ^mbetween, namely

For order m=1, filters H ₁(e ^{j ω}) be low-pass filter, its cutoff frequency is f _c ¹, namely

For order m=2n-1, filters H _2n-1(e ^{j ω}) be a Hi-pass filter, its cutoff frequency is f _c ^2n-3, namely

F above in formula _c ^mrepresent the cutoff frequency of the 0.3dB point of m rank Savitzky-Golay difference filter.

5. the transmission of the optimization light intensity based on Savitzky-Golay difference filter phase recovery method according to claim 4, is characterized in that the cutoff frequency f of complementary filter group _c ^mconcrete Numerical Implementation step is as follows:

(1) if data half-breadth n is between 1-10, then f _c ^mreference table 1 is directly found;

(2) for larger n>=25 of data half-breadth, f _c ^mcalculated by formula (17),

${f_c}^m=\frac{m-0.3}{3.5n-1},n\&GreaterEqual;25---\left(17\right)$

(3) time for data half-breadth 10 < n < 25, f _c ^mcalculated by formula (18),

${f_c}^m=\frac{m-0.3}{3.5n-2.5},10<n<25---\left(18\right)$

The 0.3dB normalized frequency (n=1 ~ 10) of the Savitzky-Golay difference filter under table 1 different order m and data half-breadth n