CN113393472B - 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
CN113393472B
CN113393472B CN202110595634.6A CN202110595634A CN113393472B CN 113393472 B CN113393472 B CN 113393472B CN 202110595634 A CN202110595634 A CN 202110595634A CN 113393472 B CN113393472 B CN 113393472B
Authority
CN
China
Prior art keywords
phase
quality
holographic
diagram
image
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.)
Active
Application number
CN202110595634.6A
Other languages
Chinese (zh)
Other versions
CN113393472A (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 method for unpacking is mainly to accurately judge the quality and the physical complexity of each point by using a similar quality map (QualityMap), start from a high-quality point and end at a low-quality starting point for unpacking 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 GDA0003682129170000021
in the formula I H (x, y) is the intensity information of the original image, I O (x, y) is the object light intensity information, I R (x, y) is the reference light intensity information, phi O (x, y) is object light phase information, phi R And (x, y) is reference light phase information.
And carrying out fast Fourier transform on the light intensity of the holographic image to obtain:
Figure GDA0003682129170000022
in the formula F [ ■]Denotes a fast Fourier transform, α (x, y) ═ I O (x，y)+I R (x, y) represents object light I O (x, y) and reference light I R (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 GDA0003682129170000023
in order to be a real image,
Figure GDA0003682129170000024
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 +1 Is the +1 order spectrum corresponding to the real image;
Figure GDA0003682129170000025
is a virtual image, F -1 Is 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 GDA0003682129170000031
wherein Im (■) represents f +1 The imaginary part of Re (■) represents f +1 Real part of (f) +1 Is 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) in the holographic phase pattern H (x, y) is expressed as a matrixX；
(2) Carrying out zero average value transformation on each row of the matrix X again;
(3) calculate the covariance matrix of matrix X
Figure GDA0003682129170000032
Wherein X T Is a transposed matrix of matrix X;
(4) computing eigenvectors e of the covariance matrix C 1 ，e 2 ……e t And a characteristic value λ, t ═ 1,2 … … N × M;
(5) sorting the eigenvalues from large 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 GDA0003682129170000033
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 GDA0003682129170000034
and
Figure GDA0003682129170000035
is the gradient of the phase difference between the pixel wrap points of the phase along two different directions x and y respectively,
Figure GDA0003682129170000036
and
Figure GDA0003682129170000037
respectively representing two averages of two phase-deflector function within the window. From this equation, a mass diagram Q (m,n), the quality diagram Q (M, N) is a two-dimensional data matrix diagram 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 GDA0003682129170000041
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 zoning composite phase unwrapping method based on digital holographic microscopic 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 a sample at the position of a focal plane of a 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 (b) of fig. 2, and a digital holographic image is collected on the platform and is marked as I H (x, y) as shown in (c) of FIG. 2.
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 I H "(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 phi H 。
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 GDA0003682129170000051
Figure GDA0003682129170000052
wherein A is O (x，y)、A R (x, y) are the light intensity information of the object light and the reference light respectively,
Figure GDA0003682129170000061
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:
sheet (x, y) ═ U O (x，y)+U R (x，y)
The intensity of the two coherent light beams can be expressed by the following equation:
Figure GDA0003682129170000062
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 I H (x, y) which converts this two-dimensional continuous light intensity distribution into an M x N matrix I by means of a digital imaging device H ' (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 digital holograms 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:
I H (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, which is called 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:
I H (x，y)＝|O| 2 +|R| 2 +|R|Oexp(-jksinθx)+|R|O * exp(jksinθx)
wherein | O | 2 Representing the luminous intensity information of the reference object, | R | 2 Representing 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 θ x) indicates the phase factor, and θ 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 obtained H ″(x，y)。
The light intensity recorded by the digital hologram is expressed as:
Figure GDA0003682129170000071
in the formula I H (x, y) is the intensity information of the original image, I O (x, y) is the object light intensity information, I R (x, y) is the reference light intensity information, phi O (x, y) is object light phase information, phi R And (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 GDA0003682129170000072
in the formula F [ ■]Denotes a fast Fourier transform, α (x, y) ═ I O (x，y)+I R (x, y) represents object light I O (x, y) and reference light I R (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 GDA0003682129170000073
in order to be a real image,
Figure GDA0003682129170000074
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 +1 Is the +1 order spectrum corresponding to the real image;
Figure GDA0003682129170000075
is a virtual image, F -1 The 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 I H ″(x，y)：
I H ″(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 +1 And (x, y) performing positive and negative tangent function transformation to extract phase information of the digital hologram:
Figure GDA0003682129170000081
wherein Im (■) represents f +1 The imaginary part of Re (■) represents f +1 The real part of (a); a spectrogram obtained by fourier-transforming phase information extracted from the HELA digital hologram of (c) in fig. 2 is shown in (a) in fig. 3, a spectrogram obtained by extracting + 1-order spectrum centered in (b) in fig. 3, and a phase spectrogram obtained by inverse fourier-transforming + 1-order spectrum into phase information is shown in (c) in fig. 3.
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 GDA0003682129170000082
the PCA principal component analysis method is used for carrying out specific work of phase aberration removal:
(1) holographic phase pattern phi H (X, y) is denoted as matrix X, which is the phase diagram φ H The 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) the covariance matrix of the matrix X is calculated
Figure GDA0003682129170000083
Wherein X T Is a transposed matrix of matrix X;
(4) finding eigenvectors e of the covariance matrix C 1 ，e 2 ……e t And 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 hologram mn And obtaining a quality diagram.
The phase derivative variance may be defined by:
Figure GDA0003682129170000091
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 GDA0003682129170000092
and
Figure GDA0003682129170000093
is the gradient of the phase difference between the pixel wrap points of the phase along two different directions x and y respectively,
Figure GDA0003682129170000094
and
Figure GDA0003682129170000095
which are used to represent two averages of two phase-deflector-rate functions within the window, respectively. According to the formula, a quality map Q (M, N) corresponding to the holographic phase map can be 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) is shown in (b) of fig. 5.
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 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, a unwrapping method based on path tracking is more recommended, and in high-quality areas, a phase unwrapping method based on a 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 diagram to be greater than or equal to 0.5 as follows:
Figure GDA0003682129170000101
and 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 when the mass ratio of the region is greater than or equal to the mass ratio x r of the quality map (x is a random index), judging the region as a low-quality region; 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 (a) in FIG. 5, and are divided according to the dividing mode in the phase diagram, as shown in (b) in FIG. 5, and are divided according to the positions, as shown in (c) in FIG. 5, and finally, the wrapping phases of different quality regions are unwrapped respectively, as shown in (d) in FIG. 5;
step 5.4: after unwrapping, the four unwrapped images are merged to form the final phase image, as 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 mass division schematic diagram of fig. 5, the conventional mass diagram has an obvious distribution of high and low mass regions, wherein the mass division position shown in (b) in fig. 5 obviously indicates that the bright regions are low mass regions (upper left, lower left, and lower right), and the dark regions are high mass regions (upper right), and then the wrapped phase diagram (for phase accuracy, phase de-aberration is required before) is divided according to the mass division schematic diagram to obtain the wrapped phase diagrams of the four sub-regions, as shown in (d) in fig. 6. Since the unwrapping method based on the minimum norm has the characteristic of locally extending parts, the most accurate phase information will not be obtained in a low-quality region if the method is adopted, as shown in (a) of fig. 6. 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 (b) of fig. 6, the phase diagram obtained by the branch-cut method retains more phase information but the central area has a more obvious line-pulling phenomenon, as shown in (c) of fig. 6, the phase diagram obtained by the mass-oriented diagram 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 regional composite unwrapping on the wrapped phase diagram corresponding to the low-quality region, carries out branch cutting unwrapping on the wrapped phase diagram corresponding to the low-quality region, carries out weighted least square unwrapping on the phase wrapped diagram corresponding to the high-quality region, and then splices and fuses the wrapped phase diagrams, as shown in (d) in fig. 6, the unwrapped phase accuracy is greatly optimized.
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 (7)

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 +1 level 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, wherein the Fourier transforming the original holographic image into the spectrogram, and the specific method for extracting the +1 level spectrum comprises the following steps:
the light intensity of the original holographic image is:
Figure FDA0003682129160000011
in the formula I H (x, y) is the intensity information of the original image, I O (x, y) is the object light intensity information, I R (x, y) is reference light intensity information, o O (x, y) is object light phase information, C R (x, y) is reference light phase information;
and carrying out fast Fourier transform on the light intensity of the holographic image to obtain:
Figure FDA0003682129160000012
in the formula F [ ■]Denotes a fast Fourier transform, α (x, y) I O (x,y)+I R (x, y) represents object light I O (x, y) and reference light I R (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 FDA0003682129160000013
in order to be a real image,
Figure FDA0003682129160000014
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 +1 Is the +1 level spectrum corresponding to the real image;
Figure FDA0003682129160000015
is a virtual image, F -1 Is the-1 order spectrum corresponding to the virtual image;
directly and independently extracting a +1 level spectrum from a Fourier spectrogram by intercepting an image;
the extracted holographic phase map is:
Figure FDA0003682129160000021
wherein Im (■) represents f +1 The imaginary part of Re (■) represents f +1 Real part of (f) +1 Is a spatial light intensity distribution
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 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 phi H (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 FDA0003682129160000022
Wherein X T Is a transposed matrix of matrix X;
(4) computing eigenvectors e of the covariance matrix C 1 ,e 2 ……e t And 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.
3. 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 FDA0003682129160000023
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 FDA0003682129160000024
and
Figure FDA0003682129160000025
is the gradient of the phase difference between the pixel wrap points of the phase along two different directions x and y respectively,
Figure FDA0003682129160000026
and
Figure FDA0003682129160000027
respectively representing two average values of two phase partial derivative rate functions in the window, and obtaining a quality diagram Q (M, N) corresponding to the holographic phase diagram according to the formula, wherein the quality diagram Q (M, N) is a data matrix diagram of two-dimensional data 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).
4. 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.
5. The partitioned composite phase unwrapping method based on digital holographic microscopy imaging as claimed in claim 4, 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.
6. The zoned composite phase unwrapping method based on digital holographic microscopy imaging as recited in claim 4, wherein a quality threshold is set to 0.5.
7. The zoned composite phase unwrapping method based on digital holographic microscopy imaging as claimed in claim 4, wherein the mass ratio is calculated by the formula:
Figure FDA0003682129160000031
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 CN113393472A (en) 2021-09-14
CN113393472B true 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)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115078413B (en) * 2022-06-30 2023-05-05 广西大学 Novel phase-shift electroholography based on transmission electron microscope
CN115452743A (en) * 2022-09-13 2022-12-09 南京理工大学 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

Also Published As

Publication number Publication date
CN113393472A (en) 2021-09-14

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
US6747771B2 (en) Off-axis illumination direct-to-digital holography
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
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
CN103322940A (en) Method for acquiring microscopic image in three-dimensional shape
CN110455834A (en) X-ray single exposure imaging device and method based on light 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
CN114241072A (en) Laminated imaging reconstruction method and system
Wu et al. Autofocusing algorithm for pixel-super-resolved lensfree on-chip microscopy
CN113298700A (en) High-resolution image reconstruction method in scattering scene
Zhang et al. Deep-learning-enhanced digital holographic autofocus imaging
KR100717414B1 (en) Off-axis illumination direct-to-digital holography
US20220020116A1 (en) Holographic ultra resolution imaging
Abdelsalam et al. Digital holographic shape measurement using Fizeau microscopy
CN110514137B (en) Phase unwrapping method, device, system, computer equipment and storage medium
CN108710205B (en) A kind of optical scanner holography self-focusing method based on edge gray difference function
Chen et al. High-resolution reconstruction of spectrum-overlapped off-axis holography by deflecting reference beam of Gaussian symmetry
Buitrago-Duque et al. Speckle Reduction in Digital Holographic Microscopy by Physical Manipulation of the Pupil Function
Qiu The feasibility of automatic focusing in digital holography by using Fresnel transform as numerical holographic reconstruction algorithm
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

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