Optical diffraction tomography resolving method, device and equipment

Abstract

The invention provides an optical diffraction tomography resolving method, device and equipment, and belongs to the technical field of optical diffraction tomography. The method comprises the following steps: acquiring an off-axis interference image of a target sample; preprocessing the off-axis interference image to obtain a target refractive index distribution diagram; and taking the target refractive index distribution diagram as input, carrying out iterative reconstruction on the refractive index distribution condition of the target sample through a preset algorithm until the difference between a theoretical scattered field calculated based on the refractive index distribution condition and an actual scattered field measured by an experiment meets a preset condition, and obtaining a resolving result diagram, the preset algorithm is generated by combining total variation regularization and a beam propagation method, and is realized by using a GPU (Graphics Processing Unit) under CUDA (Compute Unified Device Architecture). High-precision and high-efficiency recovery of the three-dimensional refractive index of a sample is realized, and the resolving efficiency and accuracy of optical diffraction tomography are improved.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of optical diffraction tomography, and particularly relates to an optical diffraction tomography method, device and equipment. BACKGROUND

[0002] Optical diffraction tomography (ODT) is an emerging three-dimensional label-free quantitative phase imaging technology, which has attracted extensive attention and rapid development in the fields of optical microscopic imaging and biomedical research in recent years. As a cutting-edge technology combining computational imaging and multi-angle illumination, ODT can accurately reconstruct the three-dimensional refractive index distribution inside the sample by adjusting the incident angle and collecting scattered light field data under multiple viewing angles, combining Rytov or Born approximation inversion algorithm, and realize high-sensitivity and high-resolution imaging of cell structure, organelle and even tissue microstructure.

[0003] In traditional bright-field microscopic imaging, transparent or weakly absorbing samples (such as living cells, tissue sections, etc.) are often difficult to image due to lack of sufficient absorption contrast, resulting in low image contrast and difficulty in observing internal structures. Although phase contrast microscopes and differential interference microscopes can enhance phase contrast to some extent, they are still limited to two-dimensional projection imaging and cannot provide three-dimensional structural information. ODT can achieve three-dimensional reconstruction of the internal structure of transparent samples by obtaining diffraction patterns through multi-angle illumination and combining inversion algorithms, breaking through the limitations of traditional imaging methods.

[0004] Compared with imaging techniques relying on dyes or fluorescent probes, ODT has the true advantage of label-free, avoiding the interference of exogenous labels on cell state, and circumventing problems such as fluorescent bleaching and phototoxicity, making the system more stable and reliable in applications such as live cell imaging and long-term dynamic observation. In addition, ODT is compatible with various optical architectures such as off-axis holography, common-path interference, and digital holographic microscopy, and has good system integration capability and engineering scalability. With the increasing demand for high-resolution, large-field-of-view, high-throughput three-dimensional imaging in scientific research and clinical applications, ODT is gradually becoming an important tool in the fields of biophysics, tissue engineering, pathological diagnosis, and drug development.

[0005] However, due to the physical limitations of the numerical aperture (NA) of the optical system, there is an inherent gap in the angular spectrum data collected by the imaging system, resulting in a serious lack of axial spectrum in the recovered three-dimensional refractive index image, which manifests as image compression or stretching in the axial direction. The classic inversion algorithm based on Rytov or Born approximation cannot fully suppress the artifacts caused by this spectral deficiency. Therefore, in order to improve image quality, researchers have introduced regularization algorithms from image processing in recent years.

[0006] The regularization idea was first proposed by Tikhonov for solving the ill-posed inverse problem. The Tikhonov regularization model is widely used in image denoising and restoration, but it tends to produce smooth solutions, which is not conducive to the preservation of image edges and texture details. To balance image denoising and structure detail restoration, various regularization strategies that integrate image prior models have been developed, allowing key structural information to be preserved while effectively suppressing noise.

[0007] However, the regularization strategy that integrates image prior models requires multiple iterations, which is computationally intensive and time-consuming for image calculation. Current algorithms are mostly implemented on a central processing unit (CPU) platform based on scripting languages such as Matlab or Python, which face limitations in parallel computing capacity, inflexible memory scheduling, and low computational efficiency. When faced with the reconstruction task of dozens of multi-angle large-size images, a single calculation often takes several hours or even longer, which seriously restricts the application efficiency of ODT in practical scenarios. SUMMARY

[0008] The present application provides an optical diffraction tomography calculation method, device and equipment, which realizes high-precision and high-efficiency recovery of three-dimensional refractive index of a sample, and improves the calculation efficiency and accuracy of optical diffraction tomography.

[0009] The present application provides an optical diffraction tomography calculation method, comprising: obtaining an off-axis interference image of a target sample; preprocessing the off-axis interference image to obtain a target refractive index distribution map; inputting the target refractive index distribution map, iteratively reconstructing the refractive index distribution of the target sample by a preset algorithm until the difference between the theoretical scattering field calculated based on the refractive index distribution and the actual scattering field measured in the experiment meets a preset condition, obtaining a calculation result map, and the preset algorithm is generated by combining total variation regularization and beam propagation method and is implemented using GPU under CUDA architecture.

[0010] According to the optical diffraction tomography solving method provided in the application, the off-axis interference image is preprocessed to obtain a target refractive index distribution map, including: acquiring first complex amplitude information of the target sample according to the off-axis interference image; acquiring an initial refractive index distribution map of the target sample according to the first complex amplitude information; extracting a minimum effective region containing the target sample from the initial refractive index distribution map through a preset image segmentation network to obtain a reference image; acquiring a two-dimensional cross-sectional image in the XZ direction in the reference image; performing Z-axis direction compression and clipping operations on the two-dimensional cross-sectional image according to a spatial transformation network to obtain a target refractive index distribution map, wherein the target refractive index distribution map includes an effective region of the target sample in the Z-axis direction in the two-dimensional cross-sectional image.

[0011] According to the optical diffraction tomography solving method provided in the application, the Z-axis direction compression and clipping operations on the two-dimensional cross-sectional image according to the spatial transformation network to obtain a target refractive index distribution map, including: inputting the two-dimensional cross-sectional image into the spatial transformation network to obtain an affine transformation matrix predicted by the spatial transformation network, wherein the affine transformation matrix contains a scaling factor and an offset amount in the Z-axis direction; performing compression and clipping operations on the two-dimensional cross-sectional image according to the affine transformation matrix to obtain a reference two-dimensional cross-sectional image, wherein the scaling factor indicates the scaling of the two-dimensional cross-sectional image in the Z-axis direction, and the offset amount is used to indicate the center position of the clipping region; acquiring a sampling coordinate of the reference two-dimensional cross-sectional image; performing sampling processing on the reference two-dimensional cross-sectional image based on the sampling coordinate through a bilinear interpolation strategy to obtain a target refractive index distribution map.

[0012] According to the optical diffraction tomography solving method provided in the application, the acquiring of the first complex amplitude information of the target sample according to the off-axis interference image, including: performing Fourier transform on a first off-axis interference image corresponding to a current frame to obtain a target spectrum map; acquiring an incident wave vector according to the target spectrum map; acquiring a second off-axis interference image corresponding to a subsequent frame; after processing the second off-axis interference image according to the incident wave vector, filtering through a low-pass spatial filter to obtain the first complex amplitude information of the target sample.

[0013] According to the optical diffraction tomography solving method provided in the application, the initial refractive index distribution of the target sample is obtained according to the first complex amplitude information, and the method comprises the following steps: obtaining second complex amplitude information corresponding to an off-axis interference image in the case of no sample; determining a Rytov phase according to the first complex amplitude information and the second complex amplitude information; obtaining a component of the incident wave vector in multiple directions and a coordinate of the target sample in a frequency domain according to the target spectrum; and obtaining the initial refractive index distribution of the target sample according to the Rytov phase, the coordinate of the target sample in the frequency domain and the component of the incident wave vector in multiple directions.

[0014] According to the optical diffraction tomography solving method provided in the application, the initial refractive index distribution is taken as input, and the refractive index distribution of the target sample is iteratively reconstructed through a preset algorithm, and the method comprises the following steps: obtaining a phase propagation kernel, the phase propagation kernel being generated based on a light beam propagation theory; calculating a target theoretical scattering field corresponding to the current refractive index distribution according to the phase propagation kernel; calculating a target gradient value according to the target theoretical scattering field; updating the current refractive index distribution of the target sample through the target gradient value; and repeating the above steps until the difference between the target theoretical scattering field calculated based on the refractive index distribution and an actual scattering field measured in an experiment is less than a preset value.

[0015] According to the optical diffraction tomography solving method provided in the application, the phase propagation kernel is obtained, and the method comprises the following steps: obtaining an interval of each propagation layer in a light field propagation process; obtaining a surrounding refractive index of the target sample and a frequency domain coordinate of the target sample; and obtaining a phase propagation kernel according to the interval of each propagation layer, the frequency domain coordinate and the surrounding refractive index.

[0016] According to the optical diffraction tomography solving method provided in the application, the target theoretical scattering field corresponding to the current refractive index distribution is calculated according to the phase propagation kernel, and the method comprises the following steps: obtaining a video memory index position corresponding to each pixel in the initial refractive index distribution; calculating a theoretical scattering field of each illumination angle propagated to a current propagation layer through forward and reverse two-dimensional Fourier transform according to the video memory index position, the phase propagation kernel and a theoretical scattering field of a previous propagation layer, respectively; determining a next propagation layer as the current propagation layer, and repeating the above steps until a theoretical scattering field of each propagation layer corresponding to each illumination angle is obtained, respectively; and determining the target theoretical scattering field according to the theoretical scattering field of each propagation layer.

[0017] According to the optical diffraction tomography solution method provided by the present application, calculating the target gradient value based on the theoretical scattering field includes: obtaining the theoretical scattering field of the last propagation layer corresponding to each illumination angle; performing the following operations for each illumination angle: obtaining the residual between the theoretical scattering field of the last propagation layer corresponding to the current illumination angle and the actual scattering field; reversely propagating the residual layer by layer to obtain the residual of each propagation layer; sequentially calculating the gradient value based on the residual of each propagation layer and the theoretical scattering field of each propagation layer; repeating the above steps until the gradient value corresponding to each illumination angle is obtained; and obtaining the target gradient value based on the gradient value corresponding to each illumination angle.

[0018] According to the optical diffraction tomography solution method provided in the present application, the current refractive index distribution of the target sample is updated through the target gradient value based on the FISTA algorithm, including: obtaining a regularization coefficient, an update step size and a current refractive index distribution; obtaining a proximity operator based on the current refractive index distribution, the regularization coefficient, the update step size and the target gradient value; and updating the current refractive index distribution of the target sample based on the proximity operator.

[0019] According to the optical diffraction tomography solution method provided by the present application, determining the theoretical scattering field based on the light field of the last propagation layer includes: obtaining the value of a fidelity term at each illumination angle based on the theoretical scattering field of each propagation layer at each illumination angle, wherein the fidelity term is used to indicate the difference between the theoretical scattering field and the actual scattering field; determining the value of a final fidelity term based on the value of the fidelity term at each illumination angle; obtaining the value of a total variation regularization term based on the value of the final fidelity term; and determining the theoretical scattering field based on the value of the final fidelity term and the value of the total variation regularization term.

[0020] The present application also provides an imaging device, comprising a laser scanning module, an interference imaging module and a sample carrying displacement stage; The sample carrying displacement stage is used to control the position of the carried target sample in the XYZ three-dimensional directions; The laser scanning module is used to provide an illumination beam with multiple illumination angles according to the position of the target sample; The interference imaging module is used to collect the scattered light field of the target sample under the illumination light beams at the multiple illumination angles and form an off-axis interference image, and the off-axis interference image is used to indicate the three-dimensional refractive index distribution information of the target sample.

[0021] According to the imaging device provided in the present application, the laser scanning module comprises a single longitudinal mode laser, a 1:2 optical fiber beam splitter, a first optical fiber flange, a first collimating lens, a two-dimensional scanning galvanometer, a first tube lens and an immersion water objective arranged in sequence; the single longitudinal mode laser emits laser light, and the laser light is divided into object light and reference light through the 1:2 optical fiber beam splitter, the object light is emitted through the first optical fiber flange and then collimated into the two-dimensional scanning galvanometer through the first collimating lens; the rotation center of the two-dimensional scanning galvanometer is located at the front focal plane of the first tube lens, and is used to provide a plurality of illumination beams with different angles of illumination according to the position of the target sample; the first tube lens and the immersion water objective form a 4-f system, and are used to image the scanning surface of the two-dimensional scanning galvanometer to the sample plane.

[0022] According to the imaging device provided in the present application, the interference imaging module comprises an optical fiber channel, and an immersion oil objective, a second tube lens, a non-polarized cube beam splitter and a CMOS camera arranged in sequence; the immersion oil objective collects scattered light of the object light scattered by the target sample, the scattered light is adjusted through the second tube lens and then enters the non-polarized cube beam splitter; the non-polarized cube beam splitter combines the adjusted scattered light with reference light from the optical fiber channel, so that an off-axis interference image is formed on the CMOS camera.

[0023] According to the imaging device provided in the present application, the optical fiber channel comprises a second optical fiber flange, a second collimating lens and a mirror arranged in sequence; the reference light is emitted through the second optical fiber flange, collimated through the second collimating lens, and then enters the non-polarized cube beam splitter through the mirror.

[0024] The present application also provides an optical diffraction tomography solving device, comprising: An acquisition unit is configured to acquire an off-axis interference image of a target sample. A preprocessing unit is configured to pre-process the off-axis interference image to obtain a target refractive index distribution map. A reconstruction unit is configured to take the target refractive index distribution map as input, iteratively reconstruct the refractive index distribution of the target sample by using a preset algorithm until the difference between a theoretical scattering field calculated based on the refractive index distribution and an actual scattering field measured in an experiment meets a preset condition, and obtain a solving result map, wherein the preset algorithm is generated by combining total variation regularization and beam propagation method and is implemented by using a GPU under a CUDA architecture.

[0025] The present application also provides an electronic device comprising a memory, a processor and a computer program stored in the memory and executable on the processor, wherein the processor implements the optical diffraction tomography solving method as described above when executing the computer program.

[0026] The application further provides a non-transitory computer-readable storage medium, which stores a computer program, and the computer program is executed by a processor to implement the optical diffraction tomography solving method.

[0027] The application further provides a computer program product, which comprises a computer program, and the computer program is executed by a processor to implement the optical diffraction tomography solving method.

[0028] The optical diffraction tomography solving method, device and equipment provided by the application first acquire an off-axis interference image of a target sample, then pre-process the off-axis interference image to obtain a target refractive index distribution map, finally take the target refractive index distribution map as input, and iteratively reconstruct the refractive index distribution of the target sample by using a preset algorithm until the difference between a theoretical scattering field calculated based on the refractive index distribution and an actual scattering field measured in an experiment meets a preset condition, and a solving result map is obtained. The preset algorithm is generated by combining a total variation regularization and a beam propagation method, and is implemented by using a GPU under a CUDA architecture. The three-dimensional refractive index of the sample can be recovered with high precision and high efficiency, and the optical diffraction tomography solving efficiency and accuracy are improved. BRIEF DESCRIPTION OF DRAWINGS

[0029] In order to more clearly illustrate the technical solutions of the present application or the prior art, the following will briefly introduce the drawings needed to be used in the embodiments or prior art description. Obviously, the drawings in the following description are some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor.

[0030] Figure 1 is one of the flowcharts of the optical diffraction tomography solving method provided by the application.

[0031] Figure 2 is the optical path schematic diagram of the imaging device provided by the application.

[0032] Figure 3 is the sample detection schematic diagram provided by the application.

[0033] Figure 4 is the second flowchart of the optical diffraction tomography solving method provided by the application.

[0034] Figure 5 is the third flowchart of the optical diffraction tomography solving method provided by the application.

[0035] Figure 6 is the principle schematic diagram of the beam propagation model provided by the application.

[0036] Figure 7 is a time contrast chart provided by the present application.

[0037] Figure 8 is a result contrast chart provided by the present application.

[0038] Figure 9 is a functional unit composition block diagram of an optical diffraction tomography solving device provided by the present application.

[0039] Figure 10 is a structural schematic diagram of an electronic device provided by the present application. DETAILED DESCRIPTION

[0040] In order to make the objects, technical solutions and advantages of the present application clearer, the technical solutions in the present application will be described clearly and completely below in combination with the drawings in the present application. Obviously, the described embodiments are part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of protection of the present application.

[0041] The terms "first", "second", and the like in the specification and claims of the present application and the above-described drawings are used to distinguish different objects, rather than to describe a specific order. In addition, the terms "include" and "have" and any variations thereof are intended to cover non-exclusive inclusion. For example, a process, method, system, product or device including a series of steps or units is not limited to the listed steps or units, but can optionally include steps or units not listed, or can optionally include other steps or units inherent to the process, method, product or device.

[0042] Reference herein to "an embodiment" means that a particular feature, structure, or characteristic described in connection with the embodiment can be included in at least one embodiment of the present application. The appearance of the phrase in various places in the specification does not necessarily all refer to the same embodiment, nor is it necessarily independent or alternative embodiments to each other. Those skilled in the art explicitly and implicitly understand that the embodiments described herein can be combined with other embodiments.

[0043] The current regularization strategy needs to be iterated for multiple rounds, and the calculation amount is large, and the image solving takes a long time. The current algorithm is mostly implemented on a central processing unit (CPU) platform based on a script language such as Matlab or Python, and faces problems such as limited parallelism, inflexible memory scheduling, and low calculation efficiency. When facing the reconstruction task of dozens of multi-angle large-size images, a single solving often takes several hours or even longer, which seriously restricts the application efficiency of ODT in actual scenarios.

[0044] To solve the above problems, the application provides an optical diffraction tomography solving method, device and equipment, and the embodiments of the application are described in detail below with reference to the drawings.

[0045] Please refer to Figure 1 The optical diffraction tomography solving method comprises the following steps.

[0046] S101, an off-axis interference image of a target sample is acquired.

[0047] The off-axis interference image can be acquired by the imaging device provided by the application.

[0048] S102, the off-axis interference image is preprocessed to obtain a target refractive index distribution map.

[0049] S103, the target refractive index distribution map is taken as input, and the refractive index distribution of the target sample is iteratively reconstructed by a preset algorithm until the difference between the theoretical scattering field calculated based on the refractive index distribution and the actual scattering field measured in the experiment meets a preset condition.

[0050] The preset algorithm is generated based on a total variation (TV) regularization constraint beam propagation method (BPM) and is implemented using a GPU under a compute unified device architecture (CUDA). That is, the preset algorithm is a TV-BPM algorithm. In a specific implementation, the related data of the TV-BPM algorithm can be stored on the GPU memory, and the CUDA language is used to perform operations on the GPU. The difference between the theoretical scattering field and the actual scattering field measured in the experiment can be represented by an objective function, and the preset condition can refer to the minimum value of the objective function.

[0051] After obtaining the refractive index distribution after the final iterative reconstruction, the data in the GPU can be transmitted to the memory, and a TIFFSTACK format picture is output as a result.

[0052] The scheme realizes the fast processing of large-scale light intensity image data acquired under multi-angle illumination conditions through parallel computing under the GPU architecture, effectively overcomes the problems of low axial resolution, serious spectrum loss and other problems of the traditional classical inversion algorithm based on Rytov or Born approximation, and the technical bottlenecks of the conventional regularization method, such as many iteration times, high computational complexity, slow running speed and the like. The method can realize high-quality and high-efficiency three-dimensional quantitative label-free imaging, has good system integration and engineering application prospect, and is suitable for multiple fields such as biological imaging, medical diagnosis and high-throughput microscopic imaging.

[0053] In a specific implementation, refer to Figure 2 The present application provides an imaging device, comprising a laser scanning module 1, an interference imaging module 2 and a sample carrying displacement stage 3; the sample carrying displacement stage 3 is used to control the position of the carried target sample in XYZ three-dimensional directions; the laser scanning module 1 is used to provide a plurality of illumination beams with different illumination angles according to the position of the target sample; and the interference imaging module 2 is used to collect the scattered light field of the target sample under the plurality of illumination beams with different illumination angles and form an off-axis interference image, which is used to indicate the three-dimensional refractive index distribution information of the target sample.

[0054] In one possible embodiment, the laser scanning module 1 comprises a single longitudinal mode laser 1.1, a 1:2 optical fiber beam splitter 1.2, a first optical fiber flange 1.3, a first collimating lens 1.4, a two-dimensional scanning galvanometer 1.5, a first tube lens 1.6 and an immersion water objective 1.7 arranged in sequence; the single longitudinal mode laser 1.1 emits laser light, which is divided into object light and reference light through the 1:2 optical fiber beam splitter 1.2, and the object light is emitted from the first optical fiber flange 1.3 and collimated into the two-dimensional scanning galvanometer 1.5 through the first collimating lens 1.4; the rotation center of the two-dimensional scanning galvanometer 1.5 is located at the front focal plane of the first tube lens 1.6, and is used to provide a plurality of illumination beams with different illumination angles according to the position of the target sample; and the first tube lens 1.6 and the immersion water objective 1.7 constitute a 4-f system, which is used to image the scanning surface of the two-dimensional scanning galvanometer 1.5 to the sample plane.

[0055] In one possible embodiment, the laser scanning module 1 comprises a single longitudinal mode laser 1.1, a 1:2 optical fiber beam splitter 1.2, a first optical fiber flange 1.3, a first collimating lens 1.4, a two-dimensional scanning galvanometer 1.5, a first tube lens 1.6 and an immersion water objective 1.7 arranged in sequence; the single longitudinal mode laser 1.1 emits laser light, which is divided into object light and reference light through the 1:2 optical fiber beam splitter 1.2, and the object light is emitted from the first optical fiber flange 1.3 and collimated into the two-dimensional scanning galvanometer 1.5 through the first collimating lens 1.4; the rotation center of the two-dimensional scanning galvanometer 1.5 is located at the front focal plane of the first tube lens 1.6, and is used to provide a plurality of illumination beams with different illumination angles according to the position of the target sample; and the first tube lens 1.6 and the immersion water objective 1.7 constitute a 4-f system, which is used to image the scanning surface of the two-dimensional scanning galvanometer 1.5 to the sample plane.

[0056] In one possible embodiment, the interference imaging module 2 comprises an optical fiber channel, and an immersion oil objective 2.1, a second tube lens 2.2, a non-polarized cube beam splitter 2.3 and a CMOS camera 2.7 arranged in sequence; the immersion oil objective 2.1 collects the scattered light of the object light scattered by the target sample, the scattered light is adjusted by the second tube lens 2.2 and then enters the non-polarized cube beam splitter 2.3; and the non-polarized cube beam splitter 2.3 combines the adjusted scattered light with the reference light from the optical fiber channel, so as to form an off-axis interference image on the CMOS camera 2.7.

[0057] In a possible embodiment, the optical fiber channel comprises a second optical fiber flange 2.4, a second collimating lens 2.5 and a mirror 2.6 arranged in sequence; the reference light is emitted through the second optical fiber flange 2.4, collimated by the second collimating lens 2.5, and then emitted into the non-polarization cube beam splitter through the mirror 2.6.

[0058] In specific implementations, the imaging device can further comprise a synchronization control module for coordinating the timing matching between the galvanometer scanning angle and the camera exposure, so as to realize accurate synchronization of multi-angle illumination and image acquisition. Through the control mechanism, the system can realize fast, efficient and consistent angle scanning and image acquisition, and ensure the accuracy and quality of subsequent three-dimensional image reconstruction.

[0059] In the present application, the above imaging device can realize high-speed acquisition of multi-angle and full-field complex amplitude information, and provide a high-quality and stable data source for subsequent regularization image reconstruction algorithm based on GPU acceleration, thereby effectively improving the speed and three-dimensional image quality of optical diffraction tomography.

[0060] As can be seen, in the embodiment, the off-axis interference image of the target sample is first acquired, then the off-axis interference image is preprocessed to obtain a target refractive index distribution map, and finally the target refractive index distribution map is taken as input to iteratively reconstruct the refractive index distribution of the target sample by a preset algorithm until the difference between the theoretical scattering field calculated based on the refractive index distribution and the actual scattering field measured in the experiment meets a preset condition, to obtain a calculation result map. The preset algorithm is generated in combination with total variation regularization and beam propagation method and is realized using GPU under CUDA architecture. The three-dimensional refractive index of the sample can be recovered with high precision and high efficiency, and the calculation efficiency and accuracy of optical diffraction tomography are improved.

[0061] In a possible embodiment, the preprocessing of the off-axis interference image to obtain the target refractive index distribution map comprises: acquiring first complex amplitude information of the target sample according to the off-axis interference image; acquiring an initial refractive index distribution map of the target sample according to the first complex amplitude information; extracting a minimum effective region containing the target sample from the initial refractive index distribution map through a preset image segmentation network to obtain a reference image; acquiring a two-dimensional cross-sectional image in the XZ direction in the reference image; performing Z-axis direction compression and clipping operation on the two-dimensional cross-sectional image according to a spatial transformation network to obtain a target refractive index distribution map, wherein the target refractive index distribution map includes the effective region of the target sample in the Z-axis direction in the two-dimensional cross-sectional image.

[0062] Among them, the off-axis interference image collected by the off-axis interference mode can be processed, and the complex amplitude information of the sample is demodulated by using Fourier and spatial domain filtering methods, that is, the phase and amplitude data of the sample scattering light field are obtained. When obtaining the initial refractive index distribution map, Rytov approximation can be used for three-dimensional diffraction tomographic reconstruction to obtain the three-dimensional refractive index distribution map of the sample, that is, the initial refractive index distribution map.

[0063] In extracting the reference image, the preset image segmentation network can be a pre-trained YOLOV11 (You Only Look Once version 11) image segmentation network. Based on the network, the obtained focal plane refractive index map can be identified and analyzed, and the minimum effective region (ROI) containing the target sample can be accurately located and extracted, so as to effectively reduce the data redundancy and computational resource overhead in subsequent processing. The focal plane refractive index map can be a two-dimensional cross-sectional image in the XZ direction.

[0064] In a specific implementation, YOLOV11 uses an improved CSPDarknet as a backbone network, which has stronger feature extraction capability. In the feature fusion part, an efficient path aggregation structure (such as PANet or BiFPN) is introduced, which significantly enhances the detection performance of multi-scale and small-size targets. At the same time, YOLOV11 introduces the Anchor-Free mechanism, discards the traditional anchor frame design, and directly regresses the target boundary box through the center point and its width and height, effectively reducing the dependence on hyperparameters, improving the training efficiency and generalization ability of the network. The model also supports multi-task expansion such as image segmentation and pose estimation, and is flexible in deployment, especially suitable for embedded systems and edge computing scenarios.

[0065] In terms of training data labeling, the two-dimensional optical diffraction map obtained by Rytov approximation reconstruction in the foregoing can be used as the segmentation basis. Since the sample cells in the map have higher phase contrast and clear boundary features, the recognition accuracy and robustness of the neural network can be effectively enhanced. As shown in Figure 3 It is shown that the YOLOV11 neural network is used to realize automatic identification and high-precision cutting of multiple cell regions in the image, which greatly improves the efficiency and feasibility of subsequent regularization iteration.

[0066] In a specific implementation, in a regularization iterative reconstruction process based on a beam propagation method (BPM), a three-dimensional voxel grid is usually preset for simulating light field propagation, and the grid has equal interval sampling in X, Y and Z directions to ensure the accuracy of calculation. To ensure that the entire sample is completely contained in the calculation area, the traditional method usually sets a large pixel depth in the Z axis (i.e. optical axis) direction. However, this way increases the calculation overhead and storage pressure while containing boundary redundant voxels, which is not conducive to high-throughput processing and real-time imaging requirements. Therefore, the present scheme performs compression and clipping operations in the Z axis direction, automatically identifies and intercepts the effective area of the sample in the Z direction, so as to achieve the purpose of simplifying data input and improving calculation efficiency.

[0067] That is, referring to Figure 4 It can be seen that the preprocessing step of the present scheme includes: first obtaining the light field complex amplitude through off-axis interference, then preliminarily reconstructing based on Rytov approximation to obtain an initial refractive index distribution map, then performing image segmentation and effective area extraction on the initial refractive index distribution map to obtain a reference image, and then performing Z direction compression and clipping on the reference image to obtain a preprocessed target refractive index distribution map. So that the regularization iterative reconstruction is finally performed based on the target refractive index distribution map.

[0068] As can be seen, in the present embodiment, based on the above preprocessing step, the redundant calculation burden of the non-target area is significantly reduced, the calculation efficiency of the regularization reconstruction and the overall processing speed of the system are improved, and at the same time, under the premise of ensuring the image quality, the image misjudgment and resource waste caused by stretching artifacts are avoided.

[0069] In one possible embodiment, the Z-axis direction compression and clipping operation on the two-dimensional cross-sectional image according to the spatial transformation network to obtain the target refractive index distribution map comprises: inputting the two-dimensional cross-sectional image into the spatial transformation network to obtain an affine transformation matrix predicted by the spatial transformation network, the affine transformation matrix containing a scaling factor and an offset in the Z-axis direction; performing compression and clipping operation on the two-dimensional cross-sectional image according to the affine transformation matrix to obtain a reference two-dimensional cross-sectional image, the scaling factor indicating the scaling of the two-dimensional cross-sectional image in the Z-axis direction, and the offset indicating the center position of the clipping area; obtaining the sampling coordinates of the reference two-dimensional cross-sectional image; performing sampling processing on the reference two-dimensional cross-sectional image based on the sampling coordinates through a bilinear interpolation strategy to obtain the target refractive index distribution map.

[0070] In a specific implementation, in the structure of the spatial transformation network (STN), the LocalizationNet is composed of three convolution modules, but its structure is not limited to the traditional standard convolution operation. Each module is composed of different two-dimensional convolution layers, ReLU activation functions, and maximum pooling layers, which are used to extract multi-scale spatial features of the image. First, the global features of the input image are extracted through the Localization Net, and an affine transformation matrix is predicted to describe the axial geometric deformation of the image. The affine matrix can be limited to only contain the scaling factor s and the offset ty in the Z direction, which is expressed as: wherein the parameter s is used to compensate for the Z-axis stretching artifact commonly existing in the initial refractive index distribution map, so that the compressed image is closer to the true sample thickness; and ty is used to align the center of the cropped region to the center of the image, which is particularly suitable for sample images with uneven stretching (e.g., more serious Z negative direction stretching).

[0071] According to the above affine matrix, a sampling coordinate grid (Grid Generator) is generated to determine the sampling position of each output pixel in the input image, and the corresponding transformation formula is: wherein is the transformed coordinate, is the original coordinate. Since the sampling coordinates are usually non-integer positions, a bilinear interpolation strategy (Bilinear Sampler) is used to calculate the pixel value, and the calculation formula is: wherein is the output image after sampling, is the image after transformation. The size of the image is .

[0072] The two-dimensional image compressed by the STN module can also be sent to the YOLOV11 network for target region detection to identify the image block containing the effective cell structure. Then the corresponding three-dimensional data body (along the Z direction) is cropped synchronously to extract the sample block with accurate spatial structure and minimum data volume. In a specific implementation, the Rytov reconstruction result can be used as the initial training input of the YOLOV11 network, and the image obtained after iterative reconstruction by the TV-BPM algorithm can be used as a supervision label (ground truth) to train the parameter regression module in the STN.

[0073] i.e., for example Figure 5As shown, when performing Z-axis direction compression clipping, first, input the picture to the STN module, the picture is a reference image. Then extract the features through the Localization Net, and then perform vectorization and full connection layer regression operation, and then construct an affine matrix and perform image transformation, and then sample the transformed image by bilinear interpolation, and finally use YOLOV11 to detect the processed image, and finally output the effective area, that is, the target refractive index distribution map.

[0074] It can be seen that, in the embodiment, the Z-axis direction compression clipping significantly reduces the required computing resources in the regularization reconstruction, effectively improves the convergence speed of the algorithm iteration, optimizes the image quality in the Z-axis direction while maintaining the imaging accuracy, and provides key support for large field of view, high resolution and label-free three-dimensional imaging.

[0075] In one possible embodiment, the first complex amplitude information of the target sample is obtained according to the off-axis interference image, including: performing Fourier transform on the first off-axis interference image corresponding to the current frame to obtain a target spectrum; obtaining an incident wave vector according to the target spectrum; obtaining a second off-axis interference image corresponding to a subsequent frame; after processing the second off-axis interference image according to the incident wave vector, filtering through a low-pass spatial filter to obtain the first complex amplitude information of the target sample.

[0076] Among them, the intensity map of the off-axis hologram can be first Fourier transformed. The incident wave vector at this time is found by finding the maximum value in the spectrum, that is, the wave vector of the reference light . Subsequently, other off-axis holograms are directly processed, and the direct current term and the conjugate term are filtered out through a low-pass spatial filter, and finally the digital wavefront is obtained by normalizing the reference light image , that is, the complex amplitude information: Among them, IDFT is inverse Fourier transform, L is a low-pass filter, and DFT is Fourier transform.

[0077] It can be seen that, in the embodiment, although the off-axis hologram has low space-bandwidth product utilization, the complex amplitude of the light field can be obtained only once, which can improve the efficiency of acquiring the phase and amplitude data of the sample scattered light field.

[0078] In a possible embodiment, the obtaining the initial refractive index distribution of the target sample according to the first complex amplitude information comprises: obtaining second complex amplitude information corresponding to an off-axis interference image without the sample; determining a Rytov phase according to the first complex amplitude information and the second complex amplitude information; obtaining components of the incident wave vector in multiple directions and coordinates of the target sample in a frequency domain according to the target spectrum; and obtaining the initial refractive index distribution of the target sample according to the Rytov phase, the coordinates of the target sample in the frequency domain, and the components of the incident wave vector in the multiple directions.

[0079] Where, although the diffraction tomography three-dimensional reconstruction based on Rytov approximation will cause the image axial stretching problem and the poor axial imaging quality due to the problem of spectrum missing. But the picture recovered based on Rytov approximation has good imaging quality in two-dimensional direction, which can facilitate subsequent accurate cell recognition. Moreover, the speed of the result recovered by the method is very fast, and as an initial solution substituted into subsequent regularization iteration, the method can also greatly accelerate the convergence speed of the regularization iterative algorithm. The complex amplitude of the scattering field of the sample can be obtained from step one Similarly, we can measure the complex amplitude of the background field without the sample Then the Rytov phase can be obtained : The Fourier diffraction theorem under Rytov approximation is: Where, is the information of the sample in the frequency domain, is the component of the incident wave vector along different directions, is the wave number of light in the medium. is the coordinate in the frequency domain. According to the above formula, the initial three-dimensional refractive index distribution of the sample can be quickly solved.

[0080] It can be seen that in the embodiment, after the complex amplitude information is obtained, the Rytov approximation is used for three-dimensional diffraction tomography reconstruction to obtain the initial three-dimensional refractive index distribution of the sample. And the distribution is used as the initial estimation of the subsequent regularization iterative algorithm, which helps to improve the convergence speed and reconstruction accuracy of the algorithm.

[0081] In one possible embodiment, the initial refractive index distribution map is used as input, and the refractive index distribution of the target sample is iteratively reconstructed using a preset algorithm, including: obtaining a phase propagation kernel, which is generated based on beam propagation theory; calculating a target theoretical scattering field corresponding to the current refractive index distribution according to the phase propagation kernel; calculating a target gradient value according to the target theoretical scattering field; updating the current refractive index distribution of the target sample using the target gradient value; and repeating the above steps until the difference between the target theoretical scattering field calculated based on the refractive index distribution and the actual scattering field measured experimentally is less than a preset value.

[0082] Among them, the idea of ​​beam propagation theory is to process the propagation process of the light field layer by layer, such as Figure 6 As shown, the interval between each layer is , the light field distribution of each layer (i.e. the theoretical scattered field in this scheme) can be obtained by calculating the previous layer, and the refractive index distribution is also processed layer by layer, so that the light field distribution of each layer can be obtained by calculating the incident light field layer by layer. The distribution of the light field through each layer is for: in, is the wave number in vacuum, is the refractive index distribution, is the distribution in the two-dimensional angular frequency domain, is the refractive index of the medium, is the fast Fourier transform, is the fast inverse Fourier transform.

[0083] When we want to infer the refractive index distribution inside the sample from the measured scattered light field information, we can construct an inverse problem model. This model is based on the difference between the scattered light field and the background light field, and continuously corrects the refractive index estimation through iterative optimization, so that the theoretical scattered field it generates gradually approaches the experimentally measured scattered field. Specifically, let is the actual scattered field measured experimentally, is estimated from the current refractive index The calculated theoretical scattered field. The gap between the theoretical scattered field and the actual scattered field measured experimentally (i.e., the objective function) satisfies the preset conditions and can be expressed as minimizing the following error norm: It can be seen that in the embodiment, the target data after the three-dimensional recognition and compression processing is taken as input and substituted into the total variation regularization algorithm (TV-BPM) based on the beam propagation method (BPM) for iterative reconstruction. The algorithm is accelerated by using GPU parallel computing under the CUDA architecture, effectively improving the iteration efficiency. The introduction of the Rytov approximate initial solution greatly improves the convergence speed of the TV-BPM algorithm, and can shorten the reconstruction time from several hours to several minutes.

[0084] In one possible embodiment, the phase propagation kernel is obtained by: obtaining the interval of each propagation layer in the light field propagation process; obtaining the surrounding refractive index of the target sample and the frequency domain coordinates of the target sample; and obtaining the phase propagation kernel according to the interval of each propagation layer, the frequency domain coordinates and the surrounding refractive index.

[0085] The phase propagation kernel obtained according to the phase propagation theory is: The phase propagation kernel can be generated by writing a CUDA kernel function . Wherein is the surrounding refractive index of the sample, is the wavelength of the incident wave, is the two-dimensional frequency domain coordinates, is the interval between different layers.

[0086] According to the numerical aperture, magnification of the oil immersion objective 2.1 and the pixel size of the camera 2.7, the maximum frequency that can be received by the system can be calculated. The frequency interval corresponding to each pixel , the objective can collect the maximum frequency corresponding to the number of pixels in the image . When the lateral frequency is within the range that can be collected by the objective, the value is 1; when the lateral frequency exceeds the range that can be collected by the objective, the value is 0. In this way, the surrounding refractive index of the sample, the wavelength of the incident wave and other information can be determined.

[0087] It can be seen that in the embodiment, the TV-BPM algorithm is used for iterative reconstruction, and GPU parallel computing is used under the CUDA architecture to accelerate the implementation, effectively improving the iteration efficiency.

[0088] In a possible embodiment, the calculating the target theoretical scattering field corresponding to the current refractive index distribution includes: obtaining a video memory index position corresponding to each pixel in the initial refractive index distribution; calculating a theoretical scattering field propagated to a current propagation layer at each illumination angle according to the video memory index position, the phase propagation kernel and a theoretical scattering field of a previous propagation layer through forward and inverse two-dimensional Fourier transform; determining a next propagation layer as the current propagation layer, and repeating the above steps until a theoretical scattering field of each propagation layer corresponding to each illumination angle is obtained; and determining the target theoretical scattering field according to the theoretical scattering field of each propagation layer.

[0089] wherein the kernel function is assigned a block of size (32, 32), a grid, wherein a size of a single image collected by default is pixels, a total number of layer slices set. An index corresponding to each pixel in the GPU can be calculated by the following method: pixel column thread number can be calculated by pixel row thread number and pixel layer number Similarly, a corresponding video memory index position is calculated according to Other kernel functions mentioned later also use this method to establish the index of the pixel and the video memory position. The kernel function is written and the cufftExecC2C function is used to perform forward and inverse two-dimensional Fourier transform calculation on the GPU to calculate layer by layer. The current refractive index distribution The light field propagated to the Jth layer at the Lth angle is: According to the above method, the total exit field under the Lth angle is calculated layer by layer until the theoretical scattering field of the last propagation layer is obtained. In particular, the basic four arithmetic operations (addition, subtraction, multiplication and division) between the matrix and the corresponding pixels (elements) of the matrix are included in the kernel function, at this time the matrix data type is a type unique to CUDA, and there is no four arithmetic operations about this type itself or the four arithmetic operations between this type and the float type in the CUFFT library, so a related structure body needs to be written to perform operator overloading.

[0090] It can be seen that in this embodiment, the BPM regularization algorithm is implemented in the CUDA language, and three-dimensional parallel computing is performed through the GPU's thread block and thread grid structures; the algorithm adopts a video memory pre-allocation strategy, and uses shared memory and register management structures for key variables such as refractive index, propagation phase, complex amplitude, etc., which significantly improves the efficiency of multi-angle and multi-layer image reconstruction.

[0091] In a possible embodiment, calculating the target gradient value based on the theoretical scattering field includes: obtaining the theoretical scattering field of the last propagation layer corresponding to each illumination angle; performing the following operations for each illumination angle: obtaining the residual between the theoretical scattering field of the last propagation layer corresponding to the current illumination angle and the actual scattering field; backpropagating the residual layer by layer to obtain the residual of each propagation layer; sequentially calculating the gradient value based on the residual of each propagation layer and the theoretical scattering field of each propagation layer; repeating the above steps until the gradient value corresponding to each illumination angle is obtained; and obtaining the target gradient value based on the gradient value corresponding to each illumination angle.

[0092] Among them, when the gradient value corresponding to each illumination angle is calculated, the total theoretical scattered field corresponding to each illumination angle can be calculated first. and the actual measured light field The residual . Then the residual is passed back layer by layer: Then, based on the residual of each layer, the gradient is obtained: Then calculate sequentially until , where conj represents the complex conjugate operation, represents the gradient of the jth layer under L illumination angles.

[0093] Finally, the target gradient value is determined based on the gradient value corresponding to each illumination angle: in, is the total number of irradiation angles.

[0094] It can be seen that in this embodiment, since the calculation of gradients at different angles can be performed in parallel, the TV-BPM algorithm accelerated by the GPU based on CUDA programming can greatly improve the computing efficiency.

[0095] In a possible embodiment, the FISTA algorithm is used to update the current refractive index distribution of the target sample based on the target gradient value, including: obtaining a regularization coefficient, an update step and a current refractive index distribution; obtaining a proximity operator according to the current refractive index distribution, the regularization coefficient, the update step and the target gradient value; and updating the current refractive index distribution of the target sample according to the proximity operator.

[0096] The refractive index distribution can be updated by using the FISTA algorithm In a specific implementation, the intermediate variable can be set as The initial value of the intermediate variable is set as First, the following is calculated: wherein is an update step.

[0097] Then, the following is calculated: wherein is a proximity operator, is a regularization coefficient.

[0098] The following formula is used to calculate: Finally, the updated refractive index distribution is calculated as: The intermediate variable is then updated as for the next round of updating the refractive index distribution.

[0099] As can be seen, in this embodiment, the TV-BPM algorithm accelerated by the GPU based on CUDA programming can greatly improve the calculation efficiency.

[0100] In a possible embodiment, the theoretical scattering field is determined according to the light field of the last propagation layer, including: obtaining a value of a fidelity term for each illumination angle according to the theoretical scattering field of each propagation layer under each illumination angle, the fidelity term being used to indicate the difference between the theoretical scattering field and the actual scattering field; determining a value of a final fidelity term according to the value of the fidelity term for each illumination angle; obtaining a value of a total variation regularization term according to the value of the final fidelity term; and determining the theoretical scattering field according to the value of the final fidelity term and the value of the total variation regularization term.

[0101] The objective function is expressed as minimizing the following error norm. On this basis, in order to enhance the stability and robustness of the solution, the total variation (TV) regularization term is introduced, and the weighted regularization objective function is further constructed as follows: in: is the fidelity term, which represents the difference between the theoretical light field and the actual light field. Since the scattered light field values ​​at different angles are measured in the experiment, the fidelity term here can be in the form of an average value: That is, the values ​​of the fidelity terms at L angles are calculated each time, and the average is taken as the fidelity term of the target optimization function. Represents the total variation regularization term, which is used to adjust the refractive index distribution Apply global constraints. Here we use the second-order norm representation: Operator is a differential vector operator used to find the differential of the refractive index along the spatial direction.

[0102] It can be seen that in this embodiment, the introduction of the total variation regularization term can enhance the stability and robustness of the solution.

[0103] To further clarify the effect of this solution, please refer to Figure 7 ,Depend on Figure 7 It can be seen that for different demonstration data sizes, based on 150 BPM iterations and the computational time on different platforms, it can be seen that the computational time of this scheme using CUDA to perform iterative reconstruction is much shorter than that based on scripting languages ​​such as Matlab and Python on the CPU platform.

[0104] Also see Figure 8 ,Depend on Figure 8 It can be seen that the refractive index distribution map reconstructed based on the start regularization iteration method of this scheme is more accurate than the refractive index distribution map reconstructed based on the existing inversion algorithm such as Rytov approximation.

[0105] An optical diffraction tomography solution device provided by the present application is described below. The optical diffraction tomography solution device described below corresponds to the optical diffraction tomography solution method described above.

[0106] See also Figure 9The optical diffraction tomography solving device 900 comprises: an acquisition unit 901 configured to acquire an off-axis interference image of a target sample; a preprocessing unit 902 configured to preprocess the off-axis interference image to obtain a target refractive index distribution map; and a reconstruction unit 903 configured to take the target refractive index distribution map as input, iteratively reconstruct a refractive index distribution condition of the target sample by using a preset algorithm until a difference between a theoretical scattering field calculated based on the refractive index distribution condition and an actual scattering field measured in an experiment meets a preset condition, and obtain a solving result map, wherein the preset algorithm is generated by combining a total variation regularization and a beam propagation method and is implemented by using a GPU under a CUDA architecture.

[0107] In one possible implementation, in the preprocessing of the off-axis interference image to obtain the target refractive index distribution map, the preprocessing unit 902 is specifically configured to: acquire first complex amplitude information of the target sample according to the off-axis interference image; acquire an initial refractive index distribution map of the target sample according to the first complex amplitude information; extract a minimum effective region containing the target sample from the initial refractive index distribution map by using a preset image segmentation network to obtain a reference image; acquire a two-dimensional cross-sectional image in an XZ direction of the reference image; and perform Z-axis direction compression and clipping operations on the two-dimensional cross-sectional image according to a spatial transformation network to obtain the target refractive index distribution map, wherein the target refractive index distribution map includes an effective region of the target sample in the Z-axis direction in the two-dimensional cross-sectional image.

[0108] In one possible implementation, in the preprocessing of the off-axis interference image to obtain the target refractive index distribution map, the preprocessing unit 902 is specifically configured to: input the two-dimensional cross-sectional image into the spatial transformation network to obtain an affine transformation matrix predicted by the spatial transformation network, wherein the affine transformation matrix contains a scaling factor and an offset in the Z-axis direction; perform compression and clipping operations on the two-dimensional cross-sectional image according to the affine transformation matrix to obtain a reference two-dimensional cross-sectional image, wherein the scaling factor indicates scaling of the two-dimensional cross-sectional image in the Z-axis direction, and the offset is used to indicate a center position of a clipping region; acquire a sampling coordinate of the reference two-dimensional cross-sectional image; and perform sampling processing on the reference two-dimensional cross-sectional image based on the sampling coordinate by using a bilinear interpolation strategy to obtain the target refractive index distribution map.

[0109] In a possible embodiment, in the obtaining the first complex amplitude information of the target sample according to the off-axis interference image, the preprocessing unit 902 is specifically configured to: perform Fourier transform on a first off-axis interference image corresponding to a current frame to obtain a target spectrum; obtain an incident wave vector according to the target spectrum; obtain a second off-axis interference image corresponding to a subsequent frame; and perform filtering on the second off-axis interference image according to the incident wave vector, and then perform filtering through a low-pass spatial filter to obtain the first complex amplitude information of the target sample.

[0110] In a possible embodiment, in the obtaining the initial refractive index distribution map of the target sample according to the first complex amplitude information, the preprocessing unit 902 is specifically configured to: obtain second complex amplitude information, the second complex amplitude information corresponding to an off-axis interference image in a case without a sample; determine a Rytov phase according to the first complex amplitude information and the second complex amplitude information; obtain a component of the incident wave vector in multiple directions and a coordinate of the target sample in a frequency domain according to the target spectrum; and obtain the initial refractive index distribution map of the target sample according to the Rytov phase, the coordinate of the target sample in the frequency domain, and the component of the incident wave vector in the multiple directions.

[0111] In a possible embodiment, in the inputting the initial refractive index distribution map and iteratively reconstructing the refractive index distribution of the target sample through a preset algorithm, the reconstruction unit 903 is specifically configured to: obtain a phase propagation kernel, the phase propagation kernel being generated based on a light beam propagation theory; calculate a target theoretical scattering field corresponding to the current refractive index distribution according to the phase propagation kernel; calculate a target gradient value according to the target theoretical scattering field; update the current refractive index distribution of the target sample through the target gradient value; and repeat the above steps until a difference between a target theoretical scattering field calculated based on the refractive index distribution and an actual scattering field measured in an experiment is less than a preset value.

[0112] In a possible embodiment, in the obtaining the phase propagation kernel, the reconstruction unit 903 is specifically configured to: obtain an interval of each propagation layer in a light field propagation process; obtain a surrounding refractive index of the target sample and a coordinate of the target sample in a frequency domain; and obtain a phase propagation kernel according to the interval of each propagation layer, the coordinate in the frequency domain, and the surrounding refractive index.

[0113] In a possible implementation, in the calculation of the target theoretical scattering field corresponding to the current refractive index distribution according to the phase propagation kernel, the reconstruction unit 903 is specifically configured to: obtain a video memory index position corresponding to each pixel in the initial refractive index distribution map; calculate a theoretical scattering field propagated to a current propagation layer at each illumination angle according to the video memory index position, the phase propagation kernel and a theoretical scattering field of a previous propagation layer through forward and inverse two-dimensional Fourier transform respectively; determine a next propagation layer as the current propagation layer, and repeat the above steps until a theoretical scattering field of each propagation layer corresponding to the each illumination angle is obtained; and determine the target theoretical scattering field according to the theoretical scattering field of each propagation layer.

[0114] In a possible implementation, in the calculation of the target gradient value according to the theoretical scattering field, the reconstruction unit 903 is specifically configured to: obtain the theoretical scattering field of the last propagation layer corresponding to the each illumination angle; for each illumination angle, perform the following operations: obtain a residual error of the theoretical scattering field of the last propagation layer corresponding to the current illumination angle and the actual scattering field; perform inverse transmission on the residual error layer by layer to obtain a residual error of each propagation layer; and calculate a gradient value according to the residual error of each propagation layer and the theoretical scattering field of each propagation layer in sequence; repeat the above steps until a gradient value corresponding to each illumination angle is obtained; and obtain the target gradient value according to the gradient value corresponding to each illumination angle.

[0115] In a possible implementation, in the updating of the current refractive index distribution of the target sample based on the FISTA algorithm and the target gradient value, the reconstruction unit 903 is specifically configured to: obtain a regularization coefficient, an update step length and a current refractive index distribution; obtain a proximity operator according to the current refractive index distribution, the regularization coefficient, the update step length and the target gradient value; and update the current refractive index distribution of the target sample according to the proximity operator.

[0116] In a possible implementation, in the determination of the theoretical scattering field according to the light field of the last propagation layer, the reconstruction unit 903 is specifically configured to: obtain a value of a fidelity term at each illumination angle according to the theoretical scattering field of each propagation layer at each illumination angle, the fidelity term being used to indicate a difference between the theoretical scattering field and the actual scattering field; determine a value of a final fidelity term according to the value of the fidelity term at each illumination angle; obtain a value of a total variation regularization term according to the value of the final fidelity term; and determine the theoretical scattering field according to the value of the final fidelity term and the value of the total variation regularization term.

[0117] Please refer to Figure 10 , Figure 10 is a structural schematic diagram of an electronic device provided in the present application. As shown in Figure 10As shown, the electronic device can include a processor 1010, a communications interface 1020, a memory 1030, and a communications bus 1040, wherein the processor 1010, the communications interface 1020, and the memory 1030 complete mutual communication through the communications bus 1040. The processor 1010 can invoke a logical instruction in the memory 1030 to execute an optical diffraction tomography solving method, which includes: obtaining an off-axis interference image of a target sample; pre-processing the off-axis interference image to obtain a target refractive index distribution map; taking the target refractive index distribution map as input, iteratively reconstructing the refractive index distribution of the target sample through a preset algorithm until the difference between a theoretical scattering field calculated based on the refractive index distribution and an actual scattering field measured in an experiment meets a preset condition, and obtaining a solving result map, wherein the preset algorithm is generated in combination with total variation regularization and a beam propagation method and is implemented using a GPU under CUDA.

[0118] In addition, the logical instruction in the memory 1030 described above can be implemented in the form of a software functional unit and sold or used as an independent product, and can be stored in a computer readable storage medium. Based on such understanding, the technical solutions of the present application essentially or the part that contributes to the prior art or part of the technical solutions can be embodied in the form of a software product. The computer software product is stored in a storage medium and includes a plurality of instructions for causing a computer device (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present application. The aforementioned storage medium includes: a U disk, a mobile hard disk, a read-only memory (ROM, Read-Only Memory), a random access memory (RAM, Random Access Memory), a magnetic disk or an optical disk, and various media that can store program codes.

[0119] On the other hand, the present application also provides a non-transitory computer readable storage medium having a computer program stored thereon, wherein the computer program is executed by a processor to implement the optical diffraction tomography solving method provided by the above-mentioned methods, which includes: obtaining an off-axis interference image of a target sample; pre-processing the off-axis interference image to obtain a target refractive index distribution map; taking the target refractive index distribution map as input, iteratively reconstructing the refractive index distribution of the target sample through a preset algorithm until the difference between a theoretical scattering field calculated based on the refractive index distribution and an actual scattering field measured in an experiment meets a preset condition, and obtaining a solving result map, wherein the preset algorithm is generated in combination with total variation regularization and a beam propagation method and is implemented using a GPU under CUDA.

[0120] In yet another aspect, the present application also provides a computer program product comprising a computer program which, when executed by a processor, implements any of the above optical diffraction tomography solving methods, the method comprising: obtaining an off-axis interference image of a target sample; pre-processing the off-axis interference image to obtain a target refractive index distribution map; taking the target refractive index distribution map as input, iteratively reconstructing the refractive index distribution of the target sample by a preset algorithm until the difference between a theoretical scattering field calculated based on the refractive index distribution and an actual scattering field measured in experiment meets a preset condition, obtaining a solving result map, the preset algorithm being generated in combination with total variation regularization and beam propagation method and implemented using GPU under CUDA.

[0121] The apparatus embodiments described above are merely illustrative, wherein the units described as separate components may or may not be physically separate, and the components displayed as units may or may not be physical units, i.e., may be located in one place, or may be distributed on multiple network units. Part or all of the modules can be selected to achieve the purpose of the embodiment scheme according to actual needs. Those skilled in the art can understand and implement without creative labor.

[0122] Through the description of the above embodiments, those skilled in the art can clearly understand that each embodiment can be realized by means of software and necessary general hardware platform, and of course can also be realized by hardware. Based on such understanding, the above technical solutions can be embodied in the form of software product, which can be stored in a computer readable storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., and includes a plurality of instructions to make a computer device (which can be a personal computer, server, or network device, etc.) execute the method described in each embodiment or some parts of the embodiment.

[0123] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present application, and not to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that: it can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement for some technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present application.

Claims

1. An optical diffraction tomography solution method, characterized in that: include: Acquire an off-axis interferometric image of the target sample; Preprocessing the off-axis interference image to obtain a target refractive index distribution map; The target refractive index distribution map is used as input, and the refractive index distribution of the target sample is iteratively reconstructed using a preset algorithm until the gap between the theoretical scattering field calculated based on the refractive index distribution and the actual scattering field measured experimentally meets a preset condition, thereby obtaining a solution result map. The preset algorithm is generated by combining total variation regularization and the beam propagation method and is implemented using a GPU under CUDA.

2. The method according to claim 1, characterized in that The preprocessing of the off-axis interference image to obtain a target refractive index distribution map includes: Acquire first complex amplitude information of the target sample according to the off-axis interference image; acquiring an initial refractive index distribution map of the target sample according to the first complex amplitude information; Extracting the minimum effective area containing the target sample from the initial refractive index distribution map through a preset image segmentation network to obtain a reference image; Acquire a two-dimensional cross-sectional image in the XZ direction of the reference image; The two-dimensional section image is compressed and cropped in the Z-axis direction according to the spatial transformation network to obtain a target refractive index distribution map, which includes the effective area of ​​the target sample in the Z-axis direction in the two-dimensional section image.

3. The method according to claim 2, characterized in that The performing Z-axis compression and cropping operations on the two-dimensional section image according to the spatial transformation network to obtain a target refractive index distribution map includes: Inputting the two-dimensional section image into the spatial transformation network to obtain an affine transformation matrix predicted by the spatial transformation network, wherein the affine transformation matrix includes a scaling factor and an offset in the Z-axis direction; performing compression and cropping operations on the two-dimensional slice image according to the affine transformation matrix to obtain a reference two-dimensional slice image, wherein the scaling factor indicates a scaling condition of the two-dimensional slice image in the Z-axis direction, and the offset indicates a center position of a cropping area; Obtaining sampling coordinates of the reference two-dimensional section image; The reference two-dimensional section image is sampled and processed based on the sampling coordinates using a bilinear interpolation strategy to obtain a target refractive index distribution map.

4. The method according to claim 2, characterized in that The acquiring first complex amplitude information of the target sample according to the off-axis interference image includes: Performing Fourier transform on the first off-axis interference image corresponding to the current frame to obtain a target spectrum diagram; Obtaining an incident wave vector according to the target spectrum; acquiring a second off-axis interference image corresponding to a subsequent frame; After processing the second off-axis interference image according to the incident wave vector, filtering is performed through a low-pass spatial filter to obtain first complex amplitude information of the target sample.

5. The method according to claim 4, characterized in that The obtaining of the initial refractive index distribution map of the target sample according to the first complex amplitude information includes: Acquiring second complex amplitude information, where the second complex amplitude information corresponds to an off-axis interference image in the absence of a sample; determining a Rytov phase according to the first complex amplitude information and the second complex amplitude information; Acquire the components of the incident wave vector in multiple directions and the coordinates of the target sample in the frequency domain according to the target spectrum diagram; An initial refractive index distribution map of the target sample is obtained according to the Rytov phase, the coordinates of the target sample in the frequency domain, and the components of the incident wave vector in multiple directions.

6. The method according to any one of claims 1 to 5, characterized in that The step of taking the initial refractive index distribution map as input and iteratively reconstructing the refractive index distribution of the target sample using a preset algorithm includes: obtaining a phase propagation kernel, wherein the phase propagation kernel is generated based on a beam propagation theory; Calculating a target theoretical scattering field corresponding to the current refractive index distribution according to the phase propagation kernel; Calculating a target gradient value according to the target theoretical scattering field; Updating the current refractive index distribution of the target sample by using the target gradient value; The above steps are repeated until the difference between the target theoretical scattering field calculated based on the refractive index distribution and the actual scattering field measured experimentally is smaller than a preset value.

7. The method according to claim 6, characterized in that The obtaining of the phase propagation kernel comprises: Obtain the interval of each propagation layer during the light field propagation process; Acquiring the surrounding refractive index of the target sample and the frequency domain coordinates of the target sample; A phase propagation kernel is obtained according to the interval of each propagation layer, the frequency domain coordinates and the surrounding refractive index.

8. The method according to claim 6, characterized in that Calculating the target theoretical scattering field corresponding to the current refractive index distribution according to the phase propagation kernel includes: Obtaining a video memory index position corresponding to each pixel in the initial refractive index distribution map; Calculating the theoretical scattering field propagated to the current propagation layer at each illumination angle by forward and inverse two-dimensional Fourier transform according to the memory index position, the phase propagation kernel and the theoretical scattering field of the previous propagation layer; Determine the next propagation layer as the current propagation layer, and repeat the above steps until the theoretical scattered field of each propagation layer corresponding to each illumination angle is obtained respectively; The target theoretical scattering field is determined according to the theoretical scattering field of each propagation layer.

9. The method according to claim 6, characterized in that Calculating the target gradient value according to the theoretical scattering field includes: Obtaining a theoretical scattered field of the last propagation layer corresponding to each illumination angle; For each illumination angle, do the following: Obtaining a residual between a theoretical scattered field of the last propagation layer corresponding to a current illumination angle and the actual scattered field; The residual is transferred back layer by layer to obtain the residual of each propagation layer; calculating gradient values ​​according to the residual of each propagation layer and the theoretical scattering field of each propagation layer in sequence; Repeat the above steps until the gradient value corresponding to each illumination angle is obtained; A target gradient value is obtained according to the gradient value corresponding to each illumination angle.

10. The method according to claim 9, characterized in that Updating the current refractive index distribution of the target sample by using the target gradient value includes: Get the regularization coefficient, update step size and current refractive index distribution; Acquire a proximity operator according to the current refractive index distribution, the regularization coefficient, the update step size, and the target gradient value; The current refractive index distribution of the target sample is updated according to the proximity operator.

11. The method according to claim 8, characterized in that The determining the theoretical scattered field according to the light field of the last propagation layer includes: Obtaining a value of a fidelity term at each illumination angle according to the theoretical scattering field of each propagation layer at each illumination angle, wherein the fidelity term is used to indicate the difference between the theoretical scattering field and the actual scattering field; Determining the value of the final fidelity item according to the value of the fidelity item at each illumination angle; Obtaining a value of a total variation regularization term according to the value of the final fidelity term; The theoretical scattered field is determined according to the value of the final fidelity term and the value of the total variation regularization term.

12. An imaging device, characterized in that: It includes a laser scanning module, an interference imaging module and a sample carrying displacement stage; The sample carrying displacement stage is used to control the position of the carried target sample in the XYZ three-dimensional directions; The laser scanning module is used to provide an illumination beam with multiple illumination angles according to the position of the target sample; The interference imaging module is used to collect the scattered light field of the target sample under the illumination light beams at the multiple illumination angles and form an off-axis interference image, and the off-axis interference image is used to indicate the three-dimensional refractive index distribution information of the target sample.

13. The device according to claim 1, characterized in that The laser scanning module includes a single longitudinal mode laser, a one-to-two fiber beam splitter, a first fiber flange, a first collimating lens, a two-dimensional scanning galvanometer, a first tube lens and a water immersion objective lens, which are arranged in sequence; The single longitudinal mode laser emits laser light, which is divided into object light and reference light by the one-to-two optical fiber beam splitter. The object light is emitted by the first optical fiber flange and then collimated by the first collimating lens to enter the two-dimensional scanning galvanometer. The rotation center of the two-dimensional scanning galvanometer is located at the front focal plane of the first tube mirror, and is used to provide an illumination beam with multiple illumination angles according to the position of the target sample; The first tube lens and the water immersion objective lens form a 4-f system, which is used to conjugately image the scanning surface of the two-dimensional scanning galvanometer onto the sample plane.

14. The device according to claim 13, characterized in that The interference imaging module includes an optical fiber channel, and an immersion oil objective lens, a second tube lens, a non-polarizing cubic beam splitter and a CMOS camera arranged in sequence; The oil immersion objective lens collects scattered light after the object light is scattered by the target sample, and the scattered light is adjusted by the second tube lens and then incident on the non-polarizing cubic beam splitter; The non-polarizing cubic beam splitter combines the adjusted scattered light with the reference light from the fiber optic path, so that an off-axis interference image is formed on the CMOS camera.

15. The device according to claim 14, characterized in that The optical fiber channel includes a second optical fiber flange, a second collimating lens and a reflector arranged in sequence; The reference light is emitted from the second optical fiber flange and then collimated by the second collimating lens, and the collimated reference light is emitted into the non-polarizing cubic beam splitter through the reflector.

16. A device for calculating optical diffraction tomography, characterized in that: include: an acquisition unit, configured to acquire an off-axis interference image of a target sample; a preprocessing unit, configured to preprocess the off-axis interference image to obtain a target refractive index distribution map; A reconstruction unit is configured to take the target refractive index distribution map as input and iteratively reconstruct the refractive index distribution of the target sample using a preset algorithm until the difference between a theoretical scattering field calculated based on the refractive index distribution and an actual scattering field measured experimentally meets a preset condition, thereby obtaining a solution result map. The preset algorithm is generated by combining total variation regularization and beam propagation method and is implemented using a GPU under the CUDA architecture.

17. An electronic device comprising a memory, a processor, and a computer program stored in the memory and running on the processor, characterized in that: When the processor executes the computer program, the optical diffraction tomography solution method according to any one of claims 1 to 11 is implemented.

18. A non-transitory computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the optical diffraction tomography solution method according to any one of claims 1 to 11 is implemented.

19. A computer program product comprising a computer program, characterized in that When the computer program is executed by a processor, the optical diffraction tomography solution method according to any one of claims 1 to 11 is implemented.

Images