CN104331616A - Solving light intensity transfer equation-based digital hologram demodulation method - Google Patents

Abstract

The invention discloses a solving light intensity transfer equation-based digital hologram demodulation method. The method comprises the following steps of at first, recreating a digital hologram, calculating to obtain complex amplitude distribution on an object plane, normalizing the complex amplitude distribution on the object plane by using amplitude, propagating the normalized complex amplitude for a small distance by using an angular spectrum numerical propagation method to obtain an out-of-focus light intensity image, taking the out-of-focus light intensity image as input data of a solving light intensity transfer equation, and solving the light intensity transfer equation by using a fast Fourier transform-based solving method to obtain phase distribution of an object. By the method, the continuous phase of the object can be directly demodulated from an abaxial hologram, the complex phase unwrapping procedure of the traditional digital hologram is effectively avoided, and tilt aberration introduced by an off-axis structure of an object parameter optical path and secondary phase distortion introduced by an microobjective in the object parameter optical path can also be removed.

Description

Based on the digital hologram demodulation method solving light intensity transmission equation

Technical field

The invention belongs to the hologram treatment technology in the digital hologram imaging in optical measurement, particularly a kind of digital hologram demodulation method based on solving light intensity transmission equation.

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.Particularly in the micro-field of biomedicine, optical microscope has become requisite instrument in biology and medical research, is the important foundation that cytology and cell biology are set up and developed.But, due to the low absorption to light wave of most of biological cell, adopt common light field optical microscopy to be difficult to realize the blur-free imaging of cell.On the other hand, biological cell PHASE DISTRIBUTION but for the extraction of cell micro observation, structural information and dynamic behavior research very important.So phase measurement or quantitative phase imaging solve this problem important technology.

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.Digital hologram, as the representational method of most in interferometric method, because its principle is simple, image taking speed is fast, record-playback is flexible, measuring accuracy is high, is widely used in the fields such as optical detection, material science, biomedical imaging.But, rebuild in digital hologram the phase place obtained be wrapped in (-π, π] between, corresponding to the principal value interval of arctan function.In order to obtain continuous phase distribution, Phase-un-wrapping must be carried out to it.The algorithm of many Phase-un-wrapping is had to be suggested ([1] D.C.Ghiglia and M.D.Pritt at present, Two-dimensional phase unwrapping:theory, algorithms, and software (Wiley New York:, 1998) .), but when there is noise and phase place " bad point " (phase residue) in phase diagram, as when running into the complexity, the non-continuous objects that there is large area phase place unconnected area, the performance of these algorithms is desirable all not to the utmost.Due to various interference inevitable in measuring process, not only regional area there will be phase place flaw sometimes, and error even also can along removing parcel direction Spreading and diffusion.Above method also also exists shortcomings in addition, as complicated in algorithm, and calculated amount is large and go parcel consuming time long etc., is difficult to be applied in the occasion of high speed real-time measurement.On the other hand, the quadratic phase introduced with microcobjective owing to tilting from axle distorts and is superimposed upon in the phase place of object.In order to accurately obtain the phase information of object itself, need to carry out physics or numerical compensation ([2] T.Colomb to these aberrations, F.Montfort, J.K ü hn, N.Aspert, E.Cuche, A.Marian, F.Charriere, S.Bourquin, P.Marquet, and C.Depeursinge, " Numerical parametric lens for shifting, magnification, and complete aberration compensation in digital holographic microscopy, " J.Opt.Soc.Am.A 23, 3177-3190 (2006) .), this be usually difficult to operation and very consuming time.

Summary of the invention

The object of the present invention is to provide a kind of digital hologram demodulation method based on solving light intensity transmission equation, the continuous phase of object can be recalled by Directly solution from off-axis hologram, not only efficiently avoid and numerous and diversely in conventional digital holography go encapsulation process, also effectively can remove the quadratic phase introduced due to microcobjective in the inclined aberration introduced from axle construction of thing ginseng light path and object light light path and distort.

The technical solution realizing the object of the invention is: a kind of digital hologram demodulation method based on solving light intensity transmission equation, and step is as follows:

The first step, reconstructing digital hologram calculates the COMPLEX AMPLITUDE U on object plane ₀(x, y);

Second step, by the COMPLEX AMPLITUDE U on this object plane ₀(x, y) is normalized by formula (7), obtains the complex amplitude U on the object plane after normalization ₀' (x, y),

$U_0^{\′}(x,y)=\frac{U_0(x,y)}{\left|U_0(x,y)\right|}---\left(7\right)$

In formula || absolute value of a complex number operation is got in representative, and object under test is changed into pure phase position thing by this step, and supposes that the light distribution of object plane is constant 1;

3rd step, by the complex amplitude U after normalization ₀' (x, y) adopts angular spectrum Numerical Propagation Method to propagate, and obtains out of focus plot of light intensity as I _{Δ z};

4th step, by out of focus plot of light intensity as I _{Δ z}as the input data solving light intensity transmission equation, through type (11) calculates the axial differential of light intensity

$\frac{\∂I(x,y)}{\∂z}\≈\frac{I_{\mathrm{\Δz}}(x,y)-1}{\mathrm{\Δz}}---\left(11\right)$

5th step, utilizes the solving method based on Fast Fourier Transform (FFT) to solve light intensity transmission equation, obtains the PHASE DISTRIBUTION of object, i.e. through type (12), adopts Fast Fourier Transform (FFT) to solve light intensity transmission equation, obtains final PHASE DISTRIBUTION

$\mathrm{\φ}(x,y)=F^{-1}\left\{\frac{k(u^2+v^2)F\{\∂I(x,y)/\∂z\}}{4{\mathrm{\π}}^2[{(u^2+v^2)}^2+\mathrm{\γ}]}\right\}---\left(12\right)$

Wherein F represents Fourier transform, and γ is regularization parameter, and regularization parameter γ has two kinds of selections, if need the quadratic phase distortion of removing in PHASE DISTRIBUTION, then its value is γ=(max (D _x, D _y)) ^-1, wherein D _xwith D _yit is the physical size of image; If do not need the quadratic phase distortion of removing in PHASE DISTRIBUTION, then the desirable dimensionless being greater than 0 of regularization parameter.

The present invention compared with prior art, its remarkable advantage: not only efficiently avoid and numerous and diversely in conventional digital holography go encapsulation process, the out of focus aberration also introduced for the inclined aberration introduced from axle construction and the microcobjective in object light light path removed due to thing ginseng light path provides simply a kind of and effective means.Because it only adopts four Fourier transforms to realize acquisition and the aberration compensation of absolute phase, it is simple many efficiently that it compares conventional phase solution packaging method, and there will not be the phase place flaw owing to going to wrap up improper generation.Have broad application prospects.

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

Accompanying drawing explanation

Fig. 1 (a) is the topological diagram that Laser Scanning Confocal Microscope measurement obtains.

Fig. 1 (b) is the line section display figure through convex center place as shown in Fig. 1 (a).

Fig. 1 (c) is the wrapped phase figure adopting the common road holographic mirror measurement of single beam splitter to obtain.

Fig. 1 (d) is the continuous phase bitmap that path dependent phase solution pack obtains.

Fig. 1 (e) is that Quality Map guidance method (quality-map guided) path has nothing to do the continuous phase bitmap that phase unwrapping pack obtains.

Fig. 1 (f) is the continuous phase bitmap obtained by the inventive method.

Fig. 1 (g) is surface topography 3-D display figure.

Fig. 1 (h) is the line section display figure corresponding to dotted line position in Fig. 1 (f).

Fig. 2 (a) is the wrapped phase distribution plan obtained by frequency spectrum translation method.

Fig. 2 (b) is through the continuous phase distribution plan obtained by Phase-un-wrapping in Fig. 2 (a).

Fig. 2 (c) is not through the wrapped phase distribution plan of slant correction.

Fig. 2 (d) is (numerical value inverse) wrapped phase distribution plan that the inventive method obtains.

Fig. 2 (e) is that the inventive method is in conjunction with the continuous phase distribution plan obtained after Tikhonov regularization.

Fig. 2 (f) is the 3-D display figure of Fig. 2 (e) PHASE DISTRIBUTION.

Embodiment

The present invention is based on and solve in the digital hologram demodulation method of light intensity transmission equation, simplified by the complex amplitude of object light during digital hologram record and reference light and be denoted as O and R, then the intensity distributions of hologram is expressed as

I _H(x,y)＝|O| ²+|R| ²+RO ^*+R ^*O (1)

In formula, * represents complex conjugate, and first two of formula (1) become zero-order image, are that object light is added with the direct of reference light intensity.Section 3 RO ^*real image, or-1 grade of picture, it contains the information of object light conjugate image, shows as a mirror image of object itself.Section 4 R ^*o is the virtual image, or+1 grade of picture, contains the raw information of object light itself, directly proportional with original objects, is to be beneficial to most directly to observe.Both real image and the virtual image are referred to as twin image.For the hologram from axle interrecord structure, hypothetical reference light wave is an ideal plane ripple, and has an angle theta relative to object light, namely

R (x, y)=| R|exp (-jk sin θ x) (2) then formula (1) can be expressed as

I _h(x, y)=| O| ²+ | R| ²+ | R|Oexp (-jk sin θ x)+| R|O ^*in exp (jk sin θ x) (3) formula, j is imaginary unit, and k is wave number, and θ is the angle of thing ginseng light.If consider the Fourier transform of hologram, so the effect of phase factor exp (-jk sin θ x) is namely that the Fourier spectrum of the virtual image and real image has been carried out mirror image translation at frequency domain.When angle theta is enough large, zero-order image be centrally located at frequency domain initial point, and two interference terms will be separated on frequency plane, and relative frequency domain initial point is mirror image distribution.

The concrete implementation step that the present invention is based on the digital hologram demodulation method solving light intensity transmission equation is described below in detail:

The first step, reconstructing digital hologram I _h(x, y) obtains the COMPLEX AMPLITUDE on object plane.Specific operation process is: first by hologram I _h(x, y) Fourier transform is carried out, then frequency domain bandpass filter is adopted to remove to filter out separately frequency spectrum corresponding to the virtual image (described in front, for off-axis hologram, the virtual image, real image are separated on frequency plane with zero-order term, so adopt frequency domain bandpass filter to be easy to be screened), then inverse Fourier transform territory of making the return trip empty obtains the COMPLEX AMPLITUDE U of its hologram plane _z(x, y).This process can be expressed as formula (4)

U _z(x, y)=| R|Oexp (-jk sin θ x)=F ^-1{ F{I _h(x, y) } × H (u, v) } Fourier transform and inverse transformation are referred to as F and F respectively in (4) formula ^-1, subscript z represents the distance between hologram plane and object plane; (u, v) representative is relative to the frequency domain coordinates of (x, y); H (u, v) represents the transition function of frequency domain bandpass filter, accurately can jig out the spectral regions at virtual image place.Certainly in order to avoid the loss of high frequency detail, its size should be large as much as possible, and utilize the transition function of relative smooth, as adopted Gauss (Gaussian) or peaceful (Hanning) window of the Chinese, can effectively avoid rebuilding gibbs (Gibbs) effect produced due to interrupted spectra in image.

Then by the COMPLEX AMPLITUDE U of hologram plane _z(x, y), adopts formula (5) by the distance of its " anti-spread "-z, thus has just got back on the object plane at itself place, obtain the COMPLEX AMPLITUDE U on object plane ₀(x, y)

U ₀(x, y)=F ^-1{ F{U _z(x, y) } H _-z(u, v) }, in (5) formula, Fourier transform and inverse transformation are referred to as F and F respectively ^-1, H _-z(u, v) is angular spectrum transition function, and its form is

$H_{-z}(u,v)=\exp [-j\frac{2\mathrm{\πz}}{\mathrm{\λ}}\sqrt{1-{\left(\mathrm{\λu}\right)}^2-{\left(\mathrm{\λv}\right)}^2}]---\left(6\right)$

In formula, λ is optical wavelength, and j is imaginary unit, and z represents the distance between hologram plane and object plane.

Second step, by the COMPLEX AMPLITUDE U on the object plane of first step gained ₀(x, y) is normalized by formula (7), obtains the complex amplitude U on the object plane after normalization ₀' (x, y),

$U_0^{\′}(x,y)=\frac{U_0(x,y)}{\left|U_0(x,y)\right|}---\left(7\right)$

In formula || absolute value of a complex number operation is got in representative, and object under test is changed into pure phase position thing by this step, and supposes that the light distribution of object plane is constant 1.

3rd step, by the complex amplitude U on the object plane after normalization ₀' (x, y) (Δ z here will obtain very little (being recommended as 1 μm) to adopt formula (8) and angular spectrum Numerical Propagation Method represented by formula (9) to propagate the distance of Δ z, to ensure the accuracy of finite difference approximation in the 4th step), obtain out of focus plot of light intensity as I _{Δ z}.

U _Δz(x,y)＝F ^-1{F{U ₀′(x,y))}H _Δz(u,v)}, (8)

I _Δz(x,y)＝|U ₀(x,y)| ²(9)

In formula, Fourier transform and inverse transformation are referred to as F and F respectively ^-1, H _-z(u, v) is angular spectrum transition function, and its form is

$H_{\mathrm{\Δz}}(u,v)=\exp \left[j\frac{2\mathrm{\π\Δz}}{\mathrm{\λ}}\sqrt{1-{\left(\mathrm{\λu}\right)}^2-{\left(\mathrm{\λv}\right)}^2}\right]---\left(10\right)$

4th step, by out of focus plot of light intensity as I _{Δ z}through type (11) calculates the axial differential of light intensity

$\frac{\∂I(x,y)}{\∂z}\≈\frac{I_{\mathrm{\Δz}}(x,y)-1}{\mathrm{\Δz}}---\left(11\right)$

5th step, through type (12), adopts Fast Fourier Transform (FFT) to solve light intensity transmission equation, obtains final PHASE DISTRIBUTION

$\mathrm{\φ}(x,y)=F^{-1}\left\{\frac{k(u^2+v^2)F\{\∂I(x,y)/\∂z\}}{4{\mathrm{\π}}^2[{(u^2+v^2)}^2+\mathrm{\γ}]}\right\}---\left(12\right)$

Wherein F represents Fourier transform, and γ is regularization parameter.Regularization parameter γ has two kinds of selections, if need the quadratic phase distortion of removing in PHASE DISTRIBUTION, then its value is γ=(max (D _x, D _y)) ^-1, wherein D _xwith D _yit is the physical size of image.If do not need the quadratic phase distortion of removing in PHASE DISTRIBUTION, then the desirable very little numerical value (being greater than the dimensionless of 0) of regularization parameter, as being less than 10 ^-10.

The PHASE DISTRIBUTION that through type (12) obtains is exactly directly continuous print, there is not phase place parcel, thus also just avoids complicated Phase-un-wrapping operation and separate a series of phase place flaw problems of wrapping up operation itself and may cause.In addition because inclined aberration itself can't cause the change of the axial differential signal of light intensity, the inventive method automatically can remove inclined aberration, and this is its another large advantage.After inclined aberration is accurately removed, remaining quadratic phase distortion is then effectively eliminated by arranging regularization parameter γ.In addition, because the present invention only adopts single width off-axis hologram, algorithm itself realizes by means of Fast Fourier Transform (FFT), and very simple and fast, makes it be very suitable for the measurement of dynamic sample.

Embodiment

By an embodiment, the validity of digital hologram demodulation method based on solving light intensity transmission equation proposed by the invention is described below.Testing sample is the reflection protrusion array of a rule.This array adopts on a silicon substrate in the pure oxygen environment of induced with laser non-intrusion type structure process technology under 0.75 atmospheric pressure to process.The object choosing this sample is to show the robustness of algorithm of the present invention under complicated phase place and noise situations.The accurate characterization of this microstructure appearance is most important in the application of integrated optics, MEMS (micro electro mechanical system) and area of solar cell for it.This microstructure can be obtained by Laser Scanning Confocal Microscope or atomic force microscope measurement, but these methods all need a large amount of sweep times.

Fig. 1 (a) is the topological diagram adopting Laser Scanning Confocal Microscope measurement to obtain; Fig. 1 (b) shows the sectional view through convex center position.The common road holographic mirror that this experiment have employed based on single beam splitter acquires off-axis hologram.Due to thing ginseng, altogether line structure itself can the quadratic phase distortion introduced of auto-compensation object lens, so only need to adopt frequency spectrum translation method to compensate droop error.The wrapped phase that Fig. 1 (c) obtains after giving and removing, can find out comprising a large amount of " bad point " (carrying out contour integration result not for 0 with free routing around it), cause classic method to be difficult to correctly remove parcel to this phase diagram.Fig. 1 (d) be adopt a kind of simple de-packaging method (line by line the method for development) relevant to path to obtain continuous phase, due to error accumulation and diffusion, in result, there is serious " bracing wire " effect.Adopt more complicated Quality Map guidance method can prevent propagation and the diffusion of error, but it still fundamentally can not avoid the appearance of mistake.As shown in Fig. 1 (e), some phase places change region (particularly in the region of minute protrusions close to top, being shown in magnification region) faster and have occurred significantly going to wrap up mistake.And adopt the hologram of inventive method demodulation to give intact continuous phase, there is not any defect areas and error propagation phenomenon, as shown in Fig. 1 (f).After being changed by phase height, the three-dimensional topology of sample can be obtained, be shown in Fig. 1 (g).Fig. 1 (h) gives the line section display corresponding to dotted line position in Fig. 1 (f), meets very well result with Laser Scanning Confocal Microscope.Due to the inventive method only adopt 4 Fast Fourier Transform (FFT)s (twice for numerical value propagate, other twice for solving light intensity transmission equation), its computational complexity is only 4N ²logN.In this embodiment, hologram is of a size of 1280 × 960, adopts the notebook computer of 2.5GHz dominant frequency to process under MATLAB environment.Method process altogether required time is 321.3ms, is only 1/6th (1.971s) that traditional Quality Map guides phase unwrapping pack.

The present invention adopts and measures the macrophage be fixed on microslide based on the holographic microscope of Michelson interferometer structure.The frequency spectrum of the quadratic phase distortion broadening that thing ginseng light phase curvature mismatch causes ± 1 grade, this allows traditional frequency spectrum translation method be difficult to correctly compensate inclined aberration.Wrapped phase distribution after Fig. 2 (a) is given and compensated by frequency spectrum translation method; The donut having occurred much bias can be found out in phase diagram, this represent droop error and do not remove completely.Give in Fig. 2 (b) and adopt Quality Map guidance method to go to wrap up the continuous phase distribution obtained.Obviously can see in figure that some go to wrap up error.In addition due to the existence of quadratic phase distortion, the phase information accurately obtaining object itself is difficult to.In order to show the calibration capability of algorithm of the present invention for aberration, we are directly applied to not in the filtered hologram of slope compensation (without any frequency spectrum translation).Now in wrapped phase by comprising the inclination carrier phase of a unusual high frequency, as shown in Fig. 2 (c).Because the inventive method directly can obtain continuous phase, conveniently show and compare, its inverse is returned wrapped phase, as shown in Fig. 2 (d).Compared to Fig. 2 (a), can find out that method of the present invention eliminates high frequency inclined aberration exactly automatically, only leave the donut being positioned at picture centre caused by quadratic phase distortion.If (regularization factors is γ=2 × 10 to introduce regularization in the process solving light intensity transmission equation ^-3μm ^-1), bending background will flatten, effectively as shown in Fig. 2 (e) Yu Fig. 2 (f).Result shows that the inventive method not only effectively eliminates phase distortion, does not also exist any owing to going to wrap up improper caused phase place flaw in result.

Claims (3)

1., based on the digital hologram demodulation method solving light intensity transmission equation, it is characterized in that step is as follows:

The first step, reconstructing digital hologram calculates the COMPLEX AMPLITUDE U on object plane ₀(x, y);

Second step, by the COMPLEX AMPLITUDE U on this object plane ₀(x, y) is normalized by formula (7), obtains the complex amplitude U ' on the object plane after normalization ₀(x, y),

$U_0^{\′}(x,y)=\frac{U_0(x,y)}{\left|U_0(x,y)\right|}---\left(7\right)$

In formula || absolute value of a complex number operation is got in representative, and object under test is changed into pure phase position thing by this step, and supposes that the light distribution of object plane is constant 1;

3rd step, by the complex amplitude U ' after normalization ₀(x, y) adopts angular spectrum Numerical Propagation Method to propagate, and obtains out of focus plot of light intensity as I _{Δ z};

4th step, by out of focus plot of light intensity as I _{Δ z}as the input data solving light intensity transmission equation, through type (11) calculates the axial differential of light intensity

$\frac{\∂I(x,y)}{\∂z}\≈\frac{I_{\mathrm{\Δz}}(x,y)-1}{\mathrm{\Δz}}---\left(11\right)$

5th step, utilizes the solving method based on Fast Fourier Transform (FFT) to solve light intensity transmission equation, obtains the PHASE DISTRIBUTION of object, i.e. through type (12), adopts Fast Fourier Transform (FFT) to solve light intensity transmission equation, obtains final PHASE DISTRIBUTION

$\mathrm{\φ}(x,y)=F^{-1}\left\{\frac{k(u^2+v^2)F\{\∂I(x,y)/\∂z\}}{4{\mathrm{\π}}^2[{(u^2+v^2)}^2+\mathrm{\γ}]}\right\}---\left(12\right)$

Wherein F represents Fourier transform, and γ is regularization parameter, and regularization parameter γ has two kinds of selections, if need the quadratic phase distortion of removing in PHASE DISTRIBUTION, then its value is γ=(max (D _x, D _y)) ^-1, wherein D _xwith D _yit is the physical size of image; If do not need the quadratic phase distortion of removing in PHASE DISTRIBUTION, then the desirable dimensionless being greater than 0 of regularization parameter.

2. the digital hologram demodulation method based on solving light intensity transmission equation according to claim 1, is characterized in that in a first step, reconstructing digital hologram I _h(x, y) obtains the COMPLEX AMPLITUDE U on object plane _zthe process of (x, y) is: first by hologram I _h(x, y) carries out Fourier transform, then adopts frequency domain bandpass filter to remove the frequency spectrum filtered out separately corresponding to the virtual image, then inverse Fourier transform territory of making the return trip empty obtains the COMPLEX AMPLITUDE U of its hologram plane _z(x, y), this process is expressed as formula (4)

U _z(x,y)＝|R|Oexp(-jksinθx)＝F ^-1{F{I _H(x,y)}×H(u,v)} (4)

In formula, Fourier transform and inverse transformation are referred to as F and F respectively ^-1, subscript z represents the distance between hologram plane and object plane; (u, v) representative is relative to the frequency domain coordinates of (x, y); H (u, v) represents the transition function of frequency domain bandpass filter;

Then by the COMPLEX AMPLITUDE U of hologram plane _z(x, y), adopts formula (5) by the distance of its " anti-spread "-z, thus has just got back on the object plane at itself place, obtain the COMPLEX AMPLITUDE U on object plane ₀(x, y)

U ₀(x,y)＝F ^-1{F{U _z(x,y)}H _-z(u,v)}, (5)

In formula, Fourier transform and inverse transformation are referred to as F and F respectively ^-1, H _-z(u, v) is angular spectrum transition function, and its form is

$H_{-z}(u,v)=\exp [-j\frac{2\mathrm{\πz}}{\mathrm{\λ}}\sqrt{1-{\left(\mathrm{\λu}\right)}^2-{\left(\mathrm{\λv}\right)}^2}]---\left(6\right)$

In formula, λ is optical wavelength, and j is imaginary unit, and z represents the distance between hologram plane and object plane.

3. the digital hologram demodulation method based on solving light intensity transmission equation according to claim 1, is characterized in that in the third step, by the complex amplitude U ' on the object plane after normalization ₀(x, y) adopts formula (8) to propagate the distance of Δ z with the angular spectrum Numerical Propagation Method represented by formula (9), obtains out of focus plot of light intensity as I _{Δ z},

U _Δz(x,y)＝F ^-1{F{U′ ₀(x,y))}H _Δz(u,v)}, (8)

I _Δz(x,y)＝|U ₀(x,y)| ²(9)

In formula, Fourier transform and inverse transformation are referred to as F and F respectively ^-1, H _-z(u, v) is angular spectrum transition function, and its form is

$H_{\mathrm{\Δz}}(u,v)=\exp [-j\frac{2\mathrm{\π\Δz}}{\mathrm{\λ}}\sqrt{1-{\left(\mathrm{\λu}\right)}^2-{\left(\mathrm{\λv}\right)}^2}]---\left(10\right)$