CN113393472A - Zoning composite phase unwrapping method based on digital holographic microscopic imaging - Google Patents

Zoning composite phase unwrapping method based on digital holographic microscopic imaging Download PDF

Info

Publication number
CN113393472A
CN113393472A CN202110595634.6A CN202110595634A CN113393472A CN 113393472 A CN113393472 A CN 113393472A CN 202110595634 A CN202110595634 A CN 202110595634A CN 113393472 A CN113393472 A CN 113393472A
Authority
CN
China
Prior art keywords
phase
quality
holographic
image
map
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110595634.6A
Other languages
Chinese (zh)
Other versions
CN113393472B (en
Inventor
左超
王小凤
吴育民
李卓识
孙佳嵩
冯世杰
蔡泽伟
李加基
范瑶
张润南
张玉珍
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University of Science and Technology
Original Assignee
Nanjing University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University of Science and Technology filed Critical Nanjing University of Science and Technology
Priority to CN202110595634.6A priority Critical patent/CN113393472B/en
Publication of CN113393472A publication Critical patent/CN113393472A/en
Application granted granted Critical
Publication of CN113393472B publication Critical patent/CN113393472B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/11—Region-based segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/141—Discrete Fourier transforms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/136—Segmentation; Edge detection involving thresholding
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20048—Transform domain processing
    • G06T2207/20056—Discrete and fast Fourier transform, [DFT, FFT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20212—Image combination
    • G06T2207/20221—Image fusion; Image merging
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30168—Image quality inspection

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Computing Systems (AREA)
  • Discrete Mathematics (AREA)
  • Holo Graphy (AREA)

Abstract

The invention provides a zoning composite phase unwrapping method based on digital holographic microscopic imaging, which comprises the following steps of: collecting an original holographic image of a tested sample; extracting a holographic phase diagram; carrying out phase dephasing on the holographic phase diagram by adopting a PCA principal component analysis method to obtain the holographic phase diagram of the sample after dephasing; calculating phase derivative variance based on path transmission unwrapping to obtain a quality map; dividing the distribution of the quality map according to the quality factor of the quality map, unwrapping different quality areas by adopting different phase unwrapping algorithms, and fusing and splicing the quality areas to obtain unwrapped phases. The invention maintains the accuracy of the unwrapped phase obtained by resolving the high-quality area and strengthens the accuracy of the phase obtained by resolving the low-quality area.

Description

Zoning composite phase unwrapping method based on digital holographic microscopic imaging
Technical Field
The invention belongs to the optical microscopic imaging technology, and particularly relates to a zoning composite phase unwrapping method based on digital holographic microscopic imaging.
Background
The traditional phase unwrapping method based on digital holographic microscopic imaging mainly comprises two types of methods based on path tracking and minimum norm. Commonly used methods based on path tracking are the "cutch" method and the quality-oriented graph method proposed by Goldstein et al. The branch cutting method is mainly characterized in that a branch cutting line is established at a pixel point corresponding to two error phase values of the phase diagram to distinguish the phase diagram in an error area, so that the end points of the branch cutting line can be bypassed to perform unwrapping, and the accuracy is improved. Although the 'branch-and-cut' method can effectively avoid a 'pole' with a large error, a new 'pole' must be judged and found out first so as to establish an error segmentation line, so that the process is complex and time-consuming, and is not beneficial to unpacking in an image with overlarge pixels; the Quality-oriented unpacking method mainly utilizes a similar Quality Map (Quality Map) to accurately judge the Quality and the physical complexity of each point, starts from a high-Quality point and finishes at a low-Quality starting point, and unpacks the parcel based on path integration. The minimum unwrapping method based on norm is the least square unwrapping method which is a global unwrapping mode, error points cannot be eliminated in the unwrapping process, but the 'dead points' are bypassed, so that the dead points still exist in the unwrapped phase diagram, but tend to be smooth in numerical calculation, and the method has the main defect that a high-precision phase diagram cannot be unwrapped in a low-quality phase diagram.
Disclosure of Invention
The invention aims to provide a partitioned composite phase unwrapping method based on digital holographic microscopy imaging, which aims to solve the influence of 'dead spots' in a digital holographic unwrapping algorithm on unwrapping, maintain high-quality areas to obtain wrapping accuracy, simultaneously strengthen the wrapping accuracy in low-quality areas and effectively inhibit the propagation of errors.
The technical scheme for realizing the purpose of the invention is as follows: a zoning composite phase unwrapping method based on digital holographic microscopic imaging comprises the following specific steps:
step 1, collecting an original holographic image of a tested sample by adopting a microscopic imaging system built by a transmission type holographic optical path;
step 2, Fourier transforming the original holographic image into a spectrogram, extracting a +1 level spectrum, moving the spectrum to the center for inverse Fourier transformation to obtain a holographic image, focusing the holographic image by using an angular spectrum propagation principle, and extracting a holographic phase diagram;
step 3, performing phase dephasing on the holographic phase diagram by adopting a PCA principal component analysis method to obtain the holographic phase diagram of the sample after dephasing;
step 4, unwrapping based on path transmission, and calculating phase derivative variance according to the pixel size of the holographic phase diagram after phase distortion is eliminated to obtain a quality diagram;
and 5, dividing the distribution of the quality map according to the quality factor of the quality map, unwrapping different quality areas by adopting different phase unwrapping algorithms, and fusing and splicing the quality areas to obtain unwrapped phases.
Preferably, the original holographic image is fourier-transformed into a spectrogram, and the specific method for extracting the +1 level spectrum is as follows:
the light intensity of the original holographic image is:
Figure BDA0003090961620000021
in the formula IH(x, y) is the intensity information of the original image, IO(x, y) is the object light intensity information, IR(x, y) is the reference light intensity information, phiO(x, y) is object light phase information, phiRAnd (x, y) is reference light phase information.
And carrying out fast Fourier transform on the light intensity of the holographic image to obtain:
Figure BDA0003090961620000022
in the formula F [ ■]Denotes a fast Fourier transform, α (x, y) ═ IO(x,y)+IR(x, y) represents object light IO(x, y) and reference light IR(x, y) a direct current component after the correlation, a ═ F [ α (x, y)]The direct current component of the two phases of dry light is subjected to fast Fourier transform to obtain a result;
Figure BDA0003090961620000023
in order to be a real image,
Figure BDA0003090961620000024
indicating the intensity of the real image, exp [ i φ ]O(x,y)-iφR(x,y)]Is the argument distribution of the real image, i denotes the imaginary unit, F+1Is the +1 order spectrum corresponding to the real image;
Figure BDA0003090961620000025
is a virtual image, F-1Is the-1 order spectrum corresponding to the virtual image;
the +1 order spectrum is extracted separately from the fourier spectrogram directly by clipping the image.
Preferably, the extracted holographic phase map is:
Figure BDA0003090961620000031
wherein Im (■) represents f+1The imaginary part of Re (■) represents f+1Real part of (f)+1Is a spatial light intensity distribution.
Preferably, the specific method for performing phase dephasing on the holographic phase diagram by using the PCA principal component analysis method comprises the following steps:
(1) holographic phase pattern phiH(X, y) is denoted as matrix X;
(2) carrying out zero average value transformation on each row of the matrix X again;
(3) calculate the covariance matrix of matrix X
Figure BDA0003090961620000032
Wherein XTIs a transposed matrix of matrix X;
(4) computing eigenvectors e of the covariance matrix C1,e2……etAnd a characteristic value λ, t ═ 1,2 … … N × M;
(5) sorting the eigenvalues from big to small, and selecting the largest k of the eigenvalues; respectively taking the corresponding k eigenvectors as row vectors to form an eigenvector matrix P; y is PX which is the data matrix after the dimension of the matrix X is reduced to k dimension;
(6) and performing least square fitting on the obtained reduced-dimension data matrix Y to remove the phase difference.
Preferably, the phase derivative variance is calculated as:
Figure BDA0003090961620000033
wherein (m, n) is a pixel gradient range of the subscript value (i, j) of each pixel point in a square having a size represented by k × k pixels in the central range,
Figure BDA0003090961620000034
and
Figure BDA0003090961620000035
is the gradient of the phase difference between the pixel wrap points of the phase along two different directions x and y respectively,
Figure BDA0003090961620000036
and
Figure BDA0003090961620000037
respectively representing two averages of two phase-deflector function within the window. According to the formula, a quality map Q (M, N) corresponding to the holographic phase map is obtained, wherein the quality map Q (M, N) is a two-dimensional data matrix map with M rows and N columns, each element of the matrix is a quality factor, and the quality factor is a value of [0, 1]]The numerical value of (c).
Preferably, the specific method for partitioning the distribution of the quality map according to the quality factor of the quality map is as follows:
setting a quality threshold value, and determining the quality ratio of the quality factor in the quality map which is greater than or equal to the set threshold value;
in the quality map Q (M, N) of M multiplied by N, small rectangles of a multiplied by b are searched from four top corners of the quality map to the center, and a is more than or equal to 1 and less than or equal to M; b is more than or equal to 1 and less than or equal to N, when the mass ratio r 'of the small rectangle is more than or equal to x r, x is less than or equal to a random index, the searching is stopped, four small-mass diagram position coordinates determined after four vertexes are respectively searched are obtained, the areas of the four small-mass diagrams are compared, and one block with the smallest area is selected as a final zoning position to obtain a zoning mass diagram Q' (m, N); dividing the quality map Q (m, n) into four rectangular areas according to the position coordinates of the partitioned quality map Q' (m, n), reserving the overlapped part of L rows and L columns of each area, and determining the coordinate positions of the four areas;
respectively calculating the mass ratio of the four regions, and judging the region as a low-mass region when the mass ratio of the regions is more than or equal to the mass ratio x r of the mass diagram; and when the area quality ratio is smaller than the quality area quality ratio, judging the area as a high quality area.
Preferably, the different quality regions are unwrapped by using a different phase unwrapping algorithm, and the specific method comprises the following steps: the quality diagram Q (m, n) is divided into position coordinates of four areas and is brought into the phase diagram, the phase diagram is divided into four small phase diagrams with different qualities according to the dividing coordinates, the small phase diagrams are unwrapped in the low-quality areas by adopting a path tracking method, and the unwrapped phase diagrams in the high-quality areas are unwrapped by adopting a minimum norm method to obtain unwrapped phase diagrams of four different quality areas.
Preferably, the quality threshold is set to 0.5.
Preferably, the mass ratio is calculated by the formula:
Figure BDA0003090961620000041
compared with the prior art, the invention has the following remarkable advantages: (1) the invention maintains the accuracy of the unwrapped phase obtained by resolving the high-quality area and strengthens the accuracy of the phase obtained by resolving the low-quality area; (2) the invention effectively inhibits the propagation of errors and greatly optimizes the effect of the finally obtained unwrapped phase diagram.
The present invention is described in further detail below with reference to the attached drawing figures.
Drawings
FIG. 1 is a flow chart of a partitioned composite phase unwrapping method based on digital holographic microscopy imaging.
FIG. 2 is a system diagram of a partitioned composite phase unwrapping method based on digital holographic microscopy imaging.
Fig. 3 is a process diagram of phase information extracted from the digital hologram of HELA cell according to the present invention.
FIG. 4 is a comparison of HELA cell digital holograms before and after phase difference analysis by PCA principal component analysis.
FIG. 5 is a schematic sectional view of the HELA cell digital hologram sectional composite unwrapping method of the present invention.
FIG. 6 is a comparison of phase information obtained by different phase unwrapping methods.
FIG. 7 is a numerical processing diagram of cell segmentation using the phase map obtained by the present invention.
Fig. 8 is a diagram of final partitioning of the quality map.
Detailed Description
The invention mainly focuses on the research of the unwrapping method, provides zoning composite unwrapping and can obtain more accurate complete phase information.
As shown in fig. 1, a zonal compound phase unwrapping method based on digital holographic microscopy imaging is based on hologram acquisition by a microscope system built by a transmission type holographic optical path. Taking the light path based on transmission type holographic imaging as an example, the holographic microscope imaging system comprises a laser light source, a beam splitter, a reflector, a phase objective, a sample and a CCD camera. The sample selects living cells (HELA cells), the hologram of the living cells is collected, the phase information of the hologram is extracted according to the digital holographic imaging principle, the phase distortion is further corrected, and the hologram phase image after phase difference removal is used as the information of phase unwrapping for subsequent work.
The method comprises the following specific steps:
step 1: collecting a digital holographic microscopic image: a microscopic imaging system built by a transmission type holographic optical path is adopted, and a CCD camera is used for collecting a holographic image of a tested sample under a microscope. Two beams of interference light are introduced in the imaging process and are respectively marked as object light O and reference light R, and an image under the microscope lens is recorded by a CCD camera. As shown in fig. 2(a), after the microscopic imaging system is constructed, the following adjustment operations are performed on the optical path:
1. the CCD camera, the beam splitter prism and a 20cm collimating lens positioned in front of the CCD camera are connected with each other, and the camera is subjected to focusing adjustment.
2. Placing the sample at the position of the focal plane of the microscope objective;
3. adjusting a light path, wherein a laser emits a light beam, the light beam is divided into object light and reference light by a beam splitter, and the object light emits parallel light from a beam splitter prism through a reflector, a microscope objective and a sample;
4. the reference light emits parallel light from the beam splitter prism through the reflector;
5. after the assembly, a live cell digital holographic microscope is formed, as shown in fig. 2(b), and a digital holographic image is collected on the platform and marked as IH(x, y) as shown in FIG. 2 (c).
Step 2: extracting the wrapping phase of the holographic image: fourier transforming the digital holographic image into a spectrogram, extracting +1 level spectrum, moving to the center, and performing inverse Fourier transformation to obtain a holographic image, which is marked as IH"(x, y). And the holographic image is focused by using the angular spectrum propagation principle, and a holographic phase image is extracted and recorded as phiH。
Step 1 shows that object light O and reference light R are introduced in the holographic measurement, and according to the interference measurement theory, the two introduced coherent light beams must be coherent light beams with the same frequency and the same polarization state, and the complex amplitudes thereof can be respectively expressed as:
Figure BDA0003090961620000051
Figure BDA0003090961620000052
wherein A isO(x,y)、AR(x, y) are the light intensity information of the object light and the reference light respectively,
Figure BDA0003090961620000061
According to the principle of light wave superposition, the complex amplitude distribution of the phase information of the object light and the reference light respectively is the interference superposition of two coherent lights:
H(x,y)＝UO(x,y)+UR(x,y)
the intensity of the two coherent light beams can be expressed by the following formula:
Figure BDA0003090961620000062
the first two terms in the formula represent the phase intensity of two coherent lights, and the third term represents a hidden interference information element, which includes phase intensity information and amplitude distribution information of the optical wavefront with the same frequency.
Considering that the recording process of the digital hologram is a process of interference fringe digitization, the intensity distribution of the digital hologram is IH(x, y) which converts this two-dimensional continuous light intensity distribution into an M x N matrix I by means of a digital imaging deviceH' (x ', y ') (x ' ═ 0, 1, … …, M-1; y ' ═ 0, 1, … …, N-1), where the pixel size of the digital detector can be expressed as Δ x Δ y. Therefore, the recording of the digital hologram has a spatial discretization process. To further analyze the individual information recorded by the digital hologram, its light intensity is marked by the following formula:
IH(x,y)＝|O|2+|R|2+RO*+R*O
the first two terms of the above equation are the direct addition of the intensity of the object light and the reference light, referred to as the zero order image. Third part RO*Is a real image, or +1 level image, which represents information of the object light conjugate image. Fourth term R*O is a virtual image, or-1 order image, representing the information of the object light itself, and is most conveniently viewed directly. The twin image is composed of a real image and a virtual image, both of which can resolve information of object light, but the combination of the two can easily cause information loss, especially phase information. Therefore, in order to extract the phase information better, an off-axis light path recording method is adopted, and the method is common and can separate the required information well.
The light intensity distribution of the digital hologram is further expressed as:
IH(x,y)＝|O|2+|R|2+|R|Oexp(-jksinθx)+|R|O*exp(jksinθx)
wherein | O |2Representing the luminous intensity information of the reference object, | R |2Representing reference light intensity information, | R | Oexp (-jksin θ x) is a real image RO*Expanded representation, | R | O*exp (jksin θ x) is the virtual image R*The expansion of O indicates that exp (-jksin theta x) indicates the phase factor and theta is the angle of the reference light relative to the object light.
It is easy to see that the effective information to be acquired is hidden in the real image (+1 level image) or the virtual image (-1 level image), and the phase factor is just to perform frequency domain mirror image translation on the fourier spectrum of the real image or the virtual image with respect to the zero-level image, when the θ angle is large enough, the zero-level image is located at the origin of the frequency domain, and the real image and the virtual image (also called as two interference images) present mirror image distribution with respect to the origin, so that one of the real image or the virtual image is screened out and then reconstructed from the frequency domain back to the spatial domain, and then the complete spatial light intensity distribution I capable of extracting the amplitude information and the phase information can be obtainedH”(x，y)。
The light intensity recorded by the digital hologram is expressed as:
Figure BDA0003090961620000071
in the formula IH(x, y) is the intensity information of the original image, IO(x, y) is the object light intensity information, IR(x, y) is the reference light intensity information, phiO(x, y) is object light phase information, phiRAnd (x, y) is reference light phase information.
The fast Fourier transform is carried out on the light intensity, so that the following can be obtained:
Figure BDA0003090961620000072
in the formula F [ ■]Denotes a fast Fourier transform, α (x, y) ═ IO(x,y)+IR(x, y) represents object light IO(x, y) and reference light IR(x, y) a direct current component after the correlation, a ═ F [ α (x, y)]The direct current component of the two phases of dry light is subjected to fast Fourier transform to obtain a result;
Figure BDA0003090961620000073
in order to be a real image,
Figure BDA0003090961620000074
indicating the intensity of the real image, exp [ i φ ]O(x,y)-iφR(x,y)]Is the argument distribution of the real image, i denotes the imaginary unit, F+1Is the +1 order spectrum corresponding to the real image;
Figure BDA0003090961620000075
is a virtual image, F-1The result is obtained by carrying out fast Fourier transform on the virtual image, and is a-1 level spectrum corresponding to the virtual image;
the +1 level spectrum (or-1 level spectrum) can be extracted from the Fourier spectrogram directly by intercepting the image, and the real image or the virtual image of the digital hologram is separated. In the simulation of MATLAB, in order to separate a real image and a virtual image, a certain piece of spectrum information can be captured in a screenshot mode, but the information obtained in screenshot is a user-defined pixel size, is different from the original image and is not beneficial to the next calculation, so that the image obtained in screenshot needs to be subjected to pixel recovery through numerical value conversion. The most intense direct current component of a spectrogram is usually located at a dot (center) of an image, so that the + 1-level spectrogram after pixel recovery is moved to the dot, and a pixel point corresponding to the maximum value of the + 1-level spectrogram needs to be obtained.
Then, the spatial light intensity distribution I of the holographic image for extracting the amplitude information and the phase information can be obtained by carrying out inverse Fourier transform on the spatial light intensity distribution IH″(x，y)：
IH″(x，y)＝Φ-1[F+1]＝Φ-1[F(f+1(x,y))]＝f+1(x,y)
Wherein phi-1[■]Representing the inverse Fourier transform, and finally f+1And (x, y) performing positive and negative tangent function transformation to extract phase information of the digital hologram:
Figure BDA0003090961620000081
wherein Im (■) represents f+1The imaginary part of Re (■) represents f+1The real part of (a); the spectrogram obtained by fourier-transforming the phase information extracted from the HELA digital hologram of fig. 2(c) is shown in fig. 3(a), the spectrogram obtained by extracting the +1 level spectrum from the spectrogram at the center position is shown in fig. 3(b), and the spectrogram obtained by inverse fourier-transforming the +1 level spectrum to extract the phase information is shown in fig. 3 (c).
And step 3: generating a wrapped phase map: performing phase dephasing on the holographic phase diagram by adopting a PCA (principal component analysis) method to obtain a complete wrapped phase diagram, eliminating the influence of phase distortion, and obtaining the holographic phase diagram recorded as phi 'of the sample after the phase distortion is eliminated'H。
The dimensionality reduction operation is carried out by a Principal Component Analysis (PCA) method, and mainly the covariance matrix form is symmetrically triangulated. The assumed original data matrix is X, the corresponding covariance data matrix is C, the data matrix P is defined as a matrix formed by combining a group of original bases according to different step sequences, and the assumed matrix Y is the original matrix data after the matrix X is subjected to basis transformation. And the covariance matrix relationship of Y is represented by matrix C and matrix D, respectively, then this covariance matrix relationship can be represented as follows:
Figure BDA0003090961620000082
the PCA principal component analysis method is used for carrying out specific work of phase aberration removal:
(1) holographic phase pattern phiH(X, y) is denoted as matrix X, which is the phase diagram φHThe two-dimensional arrays of (X, y) are arranged into a matrix X with m rows and n columns;
(2) carrying out zero average value transformation on each row of the matrix X again;
(3) calculate the covariance matrix of matrix X
Figure BDA0003090961620000083
Wherein XTIs a transposed matrix of matrix X;
(4) finding eigenvectors e of the covariance matrix C1,e2……etAnd a characteristic value λ, t ═ 1,2 … … N × M;
(5) sorting the eigenvalues from large to small, and selecting the largest k of the eigenvalues. Then, respectively taking the k eigenvectors corresponding to the k eigenvectors as row vectors to form an eigenvector matrix P; y is PX which is the data matrix after the dimension of the matrix X is reduced to k dimension;
(6) and performing least square fitting on the obtained new dimension-reduced data matrix Y to remove the phase difference.
The phase difference results obtained are shown in fig. 4(a), and the HELA cell phase compensation results obtained by correcting the phase difference are shown in fig. 4 (b).
And 4, step 4: calculating wrapped phase quality map phase: unwrapping based on path propagation, calculating phase derivative variance Z according to pixel size of wrapped phase hologrammnAnd obtaining a quality diagram.
The phase derivative variance may be defined by:
Figure BDA0003090961620000091
this expression indicates that (m, n) is the pixel gradient range of the subscript value (i, j) of each pixel point in a square of the size represented by k × k pixels in the central range,
Figure BDA0003090961620000092
and
Figure BDA0003090961620000093
is the gradient of the phase difference between the pixel wrap points of the phase along two different directions x and y respectively,
Figure BDA0003090961620000094
and
Figure BDA0003090961620000095
which are used to represent two averages of two phase-deflector-rate functions within the window, respectively. According to the formula, a mass diagram Q (M, N) corresponding to the holographic phase diagram can be obtained, wherein the mass diagram Q (M, N) is a two-dimensional data matrix with M rows and N columnsIn the figure, each element of the matrix is a quality factor, and the quality factor is a value of [0, 1]The numerical value of (c) is shown in FIG. 5 (b).
And 5: partitioning and compositely unwrapping according to the quality diagram: and dividing the quality map into a low-quality area and a high-quality area according to the quality factor of the quality map, wherein different phase unwrapping algorithms are adopted for different quality areas. Obtaining an unwrapped phase diagram in a low-quality region by adopting an unwrapping method based on path tracking; and obtaining unwrapped phase diagrams of different quality areas by adopting a phase unwrapping method based on the minimum norm in the high quality area. And fusing and splicing the quality areas to obtain the unwrapped phase.
And judging whether the quality is good or bad according to the quality factors in the obtained quality map Q (m, n). A quality map Q (M, N) solved by the phase derivative variance is a matrix map with a quality factor value of M × N between [0, 1], when the quality factor of the quality map Q (M, N) is greater than 0.5, it indicates that the point tends to have poor quality, and when a region in a quality map has a point with a dense Q (M, N) element greater than 0.5, the region is a low-quality region, which is not favorable for the traditional phase unwrapping method based on the minimum norm, the traditional phase unwrapping method based on the minimum norm will polish a local region to be the same as the whole region during unwrapping, actively drive the unwrapped phase to be "smooth" in the region, and greatly lose phase accuracy; therefore, in low quality areas, it is more recommended to adopt a unwrapping method based on path tracking, while in high quality areas, a phase unwrapping method based on minimum norm can be considered.
The zoning composite unwrapping principle based on the quality map is as follows:
step 5.1: setting the quality threshold value to be 0.5, and determining the quality proportion of the quality factor in the quality map to be greater than or equal to 0.5 as follows:
Figure BDA0003090961620000101
step 5.2: determining a high-quality region and a low-quality region in a quality map by taking the mass ratio r as a basic reference, wherein the specific method comprises the following steps: firstly, in an M multiplied by N quality map Q (M, N), small rectangles of a multiplied by b are searched from four top corners of the quality map to the center, and a is more than or equal to 1 and less than or equal to M; b is more than or equal to 1 and less than or equal to N, when the mass ratio r 'of the small rectangle is more than or equal to x r (x is more than or equal to 0.5 and less than or equal to 1 and is a random index), the searching is stopped, the position coordinates of four small-mass graphs determined after the four vertexes are respectively searched are obtained, the areas of the four small-mass graphs are compared, and one block with the smallest area is selected as the final zoning position to obtain a zoning mass graph Q' (m, N). Then, based on the position coordinates of the partitioned quality map Q' (m, n), the quality map Q (m, n) is divided into four rectangular regions, and each region is reserved with an overlapping portion of L rows and L columns, which can be used for subsequent splicing and fusion, so as to determine the coordinate positions of the four regions, as shown in fig. 8.
Step 5.3: calculating a mass ratio r' of each region, and judging the region as a low-quality region when the mass ratio of the region is more than or equal to the mass ratio x r of the mass diagram (x is a random index); and when the area quality ratio is smaller than the quality area quality ratio, judging the area as a high quality area. Finally, the position coordinates of the four regions divided into the quality diagram Q (m, n) are brought into the phase diagram, the phase diagram is divided into four small phase diagrams with different qualities according to the dividing coordinates, namely, a high-quality region and a low-quality region are divided according to the mode shown in the figure 5(a), and are divided according to the dividing mode in the phase diagram, as shown in the figure 5(b), and are divided according to the positions, as shown in the figure 5(c), and finally, the wrapping phases of the different quality regions are respectively unwrapped, as shown in the figure 5 (d);
step 5.4: after unwrapping, the four unwrapped images are merged and fused to form the final phase image, which is shown in fig. 6 (d). The invention maintains the accuracy of the unwrapped phase obtained by solving the high-quality area, enhances the accuracy of the phase obtained by solving the low-quality area, effectively inhibits the propagation of errors and greatly optimizes the effect of the finally obtained unwrapped phase diagram.
In order to compare the unwrapping effect of the partitioned composite unwrapping method based on the digital holographic microscopy imaging system, fig. 6 shows the effect of the phase diagram obtained by partitioned composite unwrapping and the phase diagram obtained by the traditional unwrapping method. As can be seen from the schematic diagram of mass division of the partition in fig. 5, the conventional mass diagram has an obvious distribution of high-quality and low-quality regions, wherein it can be obviously seen from the mass division position shown in fig. 5(b), the bright regions are low-quality regions (upper left, lower left, and lower right), and the dark regions are high-quality regions (upper right), and then the wrapped phase diagram (for the accuracy of the phase, which needs to be phase-deblurred before) is divided according to the schematic diagram of mass division to obtain the wrapped phase diagrams of the four sub-regions, as shown in fig. 6 (d). Since the unwrapping method based on the minimum norm has the characteristic of locally extending parts, if the method is adopted, the most accurate phase information cannot be obtained in the low-quality region, as shown in fig. 6 (a). The unwrapping method based on path tracking can better cope with low-quality areas, but if the phase is too large, the error propagation is more obvious, as shown in fig. 6(b), the phase diagram obtained by the branch-cut method retains more phase information but the central area has a more obvious wire-pulling phenomenon, as shown in fig. 6(c), the phase diagram obtained by the mass-guided graph method has an obvious dark and bright area, and the error cannot be well suppressed. In order to inhibit the propagation of errors, the method carries out partitioned composite unwrapping on the wrapped phase diagram corresponding to the low-quality area, carries out branch-cut unwrapping on the wrapped phase diagram corresponding to the low-quality area, carries out weighted least square unwrapping on the phase wrapped diagram corresponding to the high-quality area, and then splices and fuses the phase wrapped phase diagram, as shown in fig. 6(d), the wrapped phase accuracy is greatly optimized and unwrapped.
According to the method, only the areas are divided according to the quality factors of the quality images, the high-quality areas and the low-quality areas are respectively unwrapped, the accuracy of the unwrapped phase is obviously improved, the influence of 'dead spots' on an unwrapping algorithm can be reduced, the holographic phase information is correct, and meanwhile the integrity of the holographic phase information is greatly improved.
In order to verify the application effect of the partitioned composite unwrapping method based on the digital holographic microscopic imaging system, the solved Hela cell phase information graph is further subjected to numerical processing, and the cell segmentation is performed by using a flood back-off algorithm, as shown in FIG. 7, the solved phase information of the digital holographic cell graph is accurate, the cell segmentation can be performed well, and the cell morphology is retained to the greatest extent.

Claims (9)

1. A zoning composite phase unwrapping method based on digital holographic microscopic imaging is characterized by comprising the following specific steps:
step 1, collecting an original holographic image of a tested sample by adopting a microscopic imaging system built by a transmission type holographic optical path;
step 2, Fourier transforming the original holographic image into a spectrogram, extracting a +1 level spectrum, moving the spectrum to the center for inverse Fourier transformation to obtain a holographic image, focusing the holographic image by using an angular spectrum propagation principle, and extracting a holographic phase diagram;
step 3, performing phase dephasing on the holographic phase diagram by adopting a PCA principal component analysis method to obtain the holographic phase diagram of the sample after dephasing;
step 4, unwrapping based on path transmission, and calculating phase derivative variance according to the pixel size of the holographic phase diagram after phase distortion is eliminated to obtain a quality diagram;
and 5, dividing the distribution of the quality map according to the quality factor of the quality map, unwrapping different quality areas by adopting different phase unwrapping algorithms, and fusing and splicing the quality areas to obtain unwrapped phases.
2. The zoned composite phase unwrapping method based on digital holographic microscopy imaging as claimed in claim 1, wherein fourier transform of the original holographic image into a spectrogram is performed, and the specific method for extracting +1 level spectrum is as follows:
the light intensity of the original holographic image is:
Figure FDA0003090961610000011
in the formula IH(x, y) is the intensity information of the original image, IO(x, y) is the object light intensity information, IR(x, y) is the reference light intensity information, phiO(x, y) is object light phase information, phiRAnd (x, y) is reference light phase information.
And carrying out fast Fourier transform on the light intensity of the holographic image to obtain:
Figure FDA0003090961610000012
in the formula F [ ■]Denotes a fast Fourier transform, α (x, y) ═ IO(x,y)+IR(x, y) represents object light IO(x, y) and reference light IR(x, y) a direct current component after the correlation, a ═ F [ α (x, y)]The direct current component of the two phases of dry light is subjected to fast Fourier transform to obtain a result;
Figure FDA0003090961610000021
in order to be a real image,
Figure FDA0003090961610000022
indicating the intensity of the real image, exp [ i φ ]O(x,y)-iφR(x,y)]Is the argument distribution of the real image, i denotes the imaginary unit, F+1Is the +1 order spectrum corresponding to the real image;
Figure FDA0003090961610000023
is a virtual image, F-1Is the-1 order spectrum corresponding to the virtual image;
the +1 order spectrum is extracted separately from the fourier spectrogram directly by clipping the image.
3. The zoned composite phase unwrapping method based on digital holographic microscopy imaging as claimed in claim 1, wherein the extracted holographic phase map is:
Figure FDA0003090961610000024
wherein Im (■) represents f+1The imaginary part of Re (■) represents f+1Real part of (f)+1Is a spatial light intensity distribution。
4. The partitioned composite phase unwrapping method based on digital holographic microscopy imaging as claimed in claim 1, wherein the specific method for performing phase dephasing on the holographic phase map by PCA principal component analysis comprises:
(1) holographic phase pattern phiH(X, y) is denoted as matrix X;
(2) carrying out zero average value transformation on each row of the matrix X again;
(3) calculate the covariance matrix of matrix X
Figure FDA0003090961610000025
Wherein XTIs a transposed matrix of matrix X;
(4) computing eigenvectors e of the covariance matrix C1,e2……etAnd a characteristic value λ, t ═ 1,2 … … N × M;
(5) sorting the eigenvalues from big to small, and selecting the largest k of the eigenvalues; respectively taking the corresponding k eigenvectors as row vectors to form an eigenvector matrix P; y is PX which is the data matrix after the dimension of the matrix X is reduced to k dimension;
(6) and performing least square fitting on the obtained reduced-dimension data matrix Y to remove the phase difference.
5. The zoned composite phase unwrapping method based on digital holographic microscopy imaging as claimed in claim 1, wherein the phase derivative variance is calculated by the formula:
Figure FDA0003090961610000026
wherein (m, n) is a pixel gradient range of the subscript value (i, j) of each pixel point in a square having a size represented by k × k pixels in the central range,
Figure FDA0003090961610000027
and
Figure FDA0003090961610000028
is the gradient of the phase difference between the pixel wrap points of the phase along two different directions x and y respectively,
Figure FDA0003090961610000029
and
Figure FDA00030909616100000210
respectively representing two averages of two phase-deflector function within the window. According to the formula, a quality map Q (M, N) corresponding to the holographic phase map is obtained, wherein the quality map Q (M, N) is a two-dimensional data matrix map with M rows and N columns, each element of the matrix is a quality factor, and the quality factor is a value of [0, 1]]The numerical value of (c).
6. The zoning composite phase unwrapping method based on digital holographic microscopy imaging as claimed in claim 1, wherein the specific method for zoning the distribution of the quality map according to the quality factor of the quality map is as follows:
setting a quality threshold value, and determining the quality ratio of the quality factor in the quality map which is greater than or equal to the set threshold value;
in the quality map Q (M, N) of M multiplied by N, small rectangles of a multiplied by b are searched from four top corners of the quality map to the center, and a is more than or equal to 1 and less than or equal to M; b is more than or equal to 1 and less than or equal to N, when the mass ratio r 'of the small rectangle is more than or equal to x r, x is less than or equal to a random index, the searching is stopped, four small-mass diagram position coordinates determined after four vertexes are respectively searched are obtained, the areas of the four small-mass diagrams are compared, and one block with the smallest area is selected as a final zoning position to obtain a zoning mass diagram Q' (m, N); dividing the quality map Q (m, n) into four rectangular areas according to the position coordinates of the partitioned quality map Q' (m, n), reserving the overlapped part of L rows and L columns of each area, and determining the coordinate positions of the four areas;
respectively calculating the mass ratio of the four regions, and judging the region as a low-mass region when the mass ratio of the regions is more than or equal to the mass ratio x r of the mass diagram; and when the area quality ratio is smaller than the quality area quality ratio, judging the area as a high quality area.
7. The partitioned composite phase unwrapping method based on digital holographic microscopy imaging as claimed in claim 6, wherein different quality regions are unwrapped by using different phase unwrapping algorithms, and the specific method is as follows: the quality diagram Q (m, n) is divided into position coordinates of four areas and is brought into the phase diagram, the phase diagram is divided into four small phase diagrams with different qualities according to the dividing coordinates, the small phase diagrams are unwrapped in the low-quality areas by adopting a path tracking method, and the unwrapped phase diagrams in the high-quality areas are unwrapped by adopting a minimum norm method to obtain unwrapped phase diagrams of four different quality areas.
8. The zoned composite phase unwrapping method based on digital holographic microscopy imaging as recited in claim 6, wherein a quality threshold is set to 0.5.
9. The zoned composite phase unwrapping method based on digital holographic microscopy imaging as claimed in claim 6, wherein the mass ratio is calculated by the formula:
Figure FDA0003090961610000031
CN202110595634.6A 2021-05-29 2021-05-29 Zoning composite phase unwrapping method based on digital holographic microscopic imaging Active CN113393472B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110595634.6A CN113393472B (en) 2021-05-29 2021-05-29 Zoning composite phase unwrapping method based on digital holographic microscopic imaging

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110595634.6A CN113393472B (en) 2021-05-29 2021-05-29 Zoning composite phase unwrapping method based on digital holographic microscopic imaging

Publications (2)

Publication Number Publication Date
CN113393472A true CN113393472A (en) 2021-09-14
CN113393472B CN113393472B (en) 2022-08-16

Family

ID=77619521

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110595634.6A Active CN113393472B (en) 2021-05-29 2021-05-29 Zoning composite phase unwrapping method based on digital holographic microscopic imaging

Country Status (1)

Country Link
CN (1) CN113393472B (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115078413A (en) * 2022-06-30 2022-09-20 广西大学 Novel phase shift electronic holography based on transmission electron microscope
WO2024055602A1 (en) * 2022-09-13 2024-03-21 南京理工大学 Lens-free single-frame phase recovery method based on partially coherent light-emitting diode illumination

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103063155A (en) * 2012-12-12 2013-04-24 浙江师范大学 Rapid package removal method of digital microscopic holographic phase diagram
EP3163251A1 (en) * 2015-10-30 2017-05-03 KLA-Tencor Corporation Method and system for regional phase unwrapping with pattern-assisted correction
CN111561877A (en) * 2020-04-24 2020-08-21 西安交通大学 Variable resolution phase unwrapping method based on point diffraction interferometer

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103063155A (en) * 2012-12-12 2013-04-24 浙江师范大学 Rapid package removal method of digital microscopic holographic phase diagram
EP3163251A1 (en) * 2015-10-30 2017-05-03 KLA-Tencor Corporation Method and system for regional phase unwrapping with pattern-assisted correction
CN111561877A (en) * 2020-04-24 2020-08-21 西安交通大学 Variable resolution phase unwrapping method based on point diffraction interferometer

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115078413A (en) * 2022-06-30 2022-09-20 广西大学 Novel phase shift electronic holography based on transmission electron microscope
WO2024055602A1 (en) * 2022-09-13 2024-03-21 南京理工大学 Lens-free single-frame phase recovery method based on partially coherent light-emitting diode illumination

Also Published As

Publication number Publication date
CN113393472B (en) 2022-08-16

Similar Documents

Publication Publication Date Title
CN113393472B (en) Zoning composite phase unwrapping method based on digital holographic microscopic imaging
EP3065001B1 (en) Holographic microscope and data processing method for high-resolution hologram image
Palacios et al. 3D image reconstruction of transparent microscopic objects using digital holography
CN109916522B (en) Aberration compensation method based on hologram continuation and implementation device thereof
CN105184295B (en) A kind of holoscan space length extracting method based on wavelet transformation and connected domain
CN113418469B (en) Spectrum confocal scanning common-path digital holographic measurement system and measurement method
CN106054570A (en) Method for realizing large-phase reconstruction of single digital hologram by adopting intensity transmission equation
CN114660060A (en) Macroscopic Fourier laminated super-resolution imaging method based on matrix scanning
He et al. Robust phase aberration compensation in digital holographic microscopy by self-extension of holograms
Wu et al. Autofocusing algorithm for pixel-super-resolved lensfree on-chip microscopy
CN114241072A (en) Laminated imaging reconstruction method and system
Chen et al. High-resolution reconstruction of spectrum-overlapped off-axis holography by deflecting reference beam of Gaussian symmetry
WO2021003380A1 (en) Calibration-free phase shifting procedure for self-interference holography
CN117053716A (en) Automatic detection method for outline parameters of circular aperture interferogram
Zhang et al. Deep-learning-enhanced digital holographic autofocus imaging
US12112452B2 (en) Holographic ultra resolution imaging
KR100717414B1 (en) Off-axis illumination direct-to-digital holography
Abdelsalam et al. Digital holographic shape measurement using Fizeau microscopy
CN110514137B (en) Phase unwrapping method, device, system, computer equipment and storage medium
Buitrago-Duque et al. Speckle Reduction in Digital Holographic Microscopy by Physical Manipulation of the Pupil Function
Jin et al. Spatial consistency calibration based on phase difference minimization for parallel slightly off-axis digital holographic microscopy
CN110715931B (en) Automatic detection method and detection device for defects of transparent sample
Ilhan et al. Autofocusing in digital holography
Sato et al. Lens-less holographic microscope with high resolving power and no-distortion
Guan et al. Fast autofocusing in off-axis digital holography based on search region segmentation and dichotomy

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant