Three-dimensional reconstruction method of cryoelectron microscope image sequence

Abstract

The invention discloses a three-dimensional reconstruction method of a cryoelectron microscope image sequence, and belongs to the technical field of electron microscope image processing. The method comprises the following steps: acquiring a cryoelectron microscope image sequence, and preprocessing; performing three-dimensional reconstruction on the preprocessed data, firstly designing a kernel function, processing a core compute-intensive task in tomographic reconstruction by utilizing a compute unified device architecture CUDA, then optimizing data reading efficiency through a GPU memory type, and customizing a CUDA kernel function based on a voxel scattering strategy; after calculation is completed, the reconstructed 3D volume data is subjected to data merging through MPI, and the calculation efficiency of a large-scale reconstruction task is improved; a memory management strategy based on slices is adopted, and MPI multi-thread parallel writing-in is carried out on the merged reconstructed daughter data; and outputting the reconstructed data. Compared with an existing reconstruction method, the method has the advantages that the reconstruction efficiency can be remarkably improved, and the reconstruction precision is not lost.

Description

Technical Field

[0001] The present invention belongs to the technical field of electron microscope image processing, and in particular, relates to a three-dimensional reconstruction method for a cryo-electron microscope image sequence. Background Art

[0002] With the in-depth application of cryo-electron tomography (cryo-ET) in the field of life science research, the amount of data generated has increased exponentially. At the same time, with the continuous advancement of sampling technology, the scale of projection data obtained in the data acquisition link is also becoming increasingly large. For example, the amount of image data collected by a single scan of a high-resolution cryo-electron microscope has reached TB level. However, these massive amounts of high-resolution projection data are extremely time-consuming when used for three-dimensional reconstruction, and the reconstruction speed has gradually become a bottleneck in the entire research process, especially when faced with large-scale data sets, the problem of low computational efficiency has become more and more prominent.

[0003] In the current computing technology system, the application of highly parallel architecture has brought significant changes to the field of biological computing. Many biological computing tasks can be efficiently performed by making full use of the multi-core parallel computing capabilities of the central processing unit (CPU) through the message passing interface (MPI), and by leveraging the powerful parallel computing advantages of the graphics processing unit (GPU) with the help of the compute unified device architecture (CUDA). For example, in the tomographic reconstruction algorithm, TOMO3D has integrated CPU acceleration technology, which has improved the reconstruction speed to a certain extent by optimizing the parallel execution of the algorithm on a multi-core CPU. At the same time, like RELION [2] 、EMAN [3] 、M [4] and IMOD [5] Computational frameworks such as , use the powerful computing power of GPU to support reconstruction tasks. They achieve relatively efficient 3D reconstruction by rationally scheduling GPU computing resources. In addition, many reconstruction algorithms optimized for GPU have been developed to provide high-performance fault processing solutions. These algorithms have achieved good reconstruction results in specific scenarios by deeply mining the hardware characteristics of GPU.

[0004] Although certain progress has been made in the application of the above-mentioned technologies, existing tools still have many shortcomings. On the one hand, the hardware architecture on which many existing tools rely is relatively old, and they fail to fully utilize the full computing potential of modern multi-core CPUs and advanced GPUs. For example, some reconstruction tools designed based on early hardware architectures cannot effectively utilize multiple cores of the CPU for parallel computing when facing new multi-core CPUs due to insufficient algorithm parallelization, resulting in a waste of computing resources. The present invention deeply optimizes the algorithm, which can perfectly adapt to modern multi-core CPUs and advanced GPUs, fully tap the potential of the hardware, and greatly improve computing efficiency. On the other hand, existing tools have limitations in terms of user interaction, which usually restrict users from flexibly and effectively adjusting 3D volume data parameters, making it difficult for researchers to finely control the reconstruction process according to actual needs when facing different experimental conditions and sample characteristics, greatly weakening the flexibility and adaptability of the tools. Summary of the invention

[0005] The purpose of the present invention is to provide a three-dimensional reconstruction method for a cryo-electron microscopy image sequence to accelerate the three-dimensional reconstruction process of cryo-electron tomography and make up for the shortcomings of the prior art.

[0006] To achieve the above object, the present invention is achieved through the following technical solutions:

[0007] A three-dimensional reconstruction method for a cryo-electron microscopy image sequence comprises the following steps:

[0008] S1: Acquire cryo-EM image sequences and perform preprocessing;

[0009] S2: Perform 3D reconstruction on the preprocessed data. First, design the kernel function and use the Compute Unified Device Architecture (CUDA) to handle the core computationally intensive tasks in tomographic reconstruction. Then optimize the data reading efficiency through the GPU memory type and customize the CUDA kernel function based on the voxel splatting strategy.

[0010] S3: After the calculation is completed, the reconstructed 3D volume data is merged using MPI to improve the computational efficiency of large-scale reconstruction tasks; a slice-based memory management strategy is adopted, and MPI multi-threading is used to write the merged reconstructed sub-volume data in parallel; and the reconstructed data is output.

[0011] Furthermore, in S1: the image sequence includes a contrast transfer function correction (CTF) and an aligned tilt sequence and its corresponding tilt angles; MPI technology is used for multi-threaded parallel file operations to read the image sequence, effectively alleviating the I / O bottleneck; at the same time, MPI is also used for preprocessing tasks, including filtering calculations and geometric adjustments of three-dimensional volume data, to improve data quality and optimize subsequent reconstruction processes.

[0012] Furthermore, the preprocessing includes:

[0013] (1) Filter calculation

[0014] In the weighted back projection method (Weighted Back Projection, WBP) and the filtered back projection method (Filtered Back Projection, FBP), filtering calculation is involved; WBP performs weighted processing on the projection image before back projection, thereby effectively reducing the blurring effect in the BPT process, and finally obtaining a clearer and more accurate 3D reconstruction result. The filtering formula is as follows:

[0015]

[0016] Where H(·) is a weighting function used to assign weights to each projection; σ Nx is the standard deviation of the Gaussian decay part; weight is a predefined threshold (the default value is 0.05);

[0017] Similar to the WBP algorithm, FBP relies on a Fourier transform-based filtering method to act on projection data, aiming to eliminate artifacts and enhance the structural details of the 3D reconstruction results. The present invention provides two filters, RamLak filter and Shepp-Logan filter:

[0018] The RamLak filter is mainly used to suppress high-frequency noise while retaining low-frequency information, thereby improving the clarity of the reconstructed image.

[0019] The Shepp-Logan filter is designed to reduce the generation of artifacts and make the reconstruction results more precise and accurate.

[0020] The filtering formulas of the two methods are:

[0021]

[0022] Where H(·) is the filter function, where H RL (·) represents RamLak filter, H SL(·) represents the Shepp-Logan filter, ω represents the frequency, ω c is the cut-off frequency.

[0023] (2) Geometry adjustment

[0024] Applying the geometry model, the spatial positioning of the sample in 3D space is precisely adjusted by adjusting the tilt angle θ, pitch angle φ, and Z shift t to ensure efficient and accurate 3D reconstruction. The adjustment formula is as follows:

[0025]

[0026] Furthermore, in the custom CUDA kernel function in S2, each voxel is "distributed" to the detector and the pixel values ​​of four adjacent points are updated to generate a projection image. Subsequently, these four points are used in the projection data to update the corresponding voxels in the reconstructed volume. Both projection and back-projection use transposition operations to enhance the accuracy of the results. This strategy simplifies the calculation and speeds up the processing compared to the traditional ray-casting method.

[0027] Voxel scattering is a method of mapping voxel values ​​onto an image. In this process, each voxel is "scattered" onto a plane (or detector) and weighted according to its range of influence in the image. The contribution of each voxel is usually calculated by interpolation. Assuming that the size of a voxel is Δx×Δy×Δz, and the center of the voxel is at (x,y,z), and the value of the voxel is V(x,y,z), then the contribution of the voxel to a pixel in the image plane is calculated by the interpolation function W(r):

[0028]

[0029] Where I(p) is the projection value of pixel p on the detector, V(r i ) is the voxel at position r i The value at W(r i ) is the voxel pixel r on the image plane i The interpolation function of the present invention is bilinear difference, that is,

[0030] W(r i )=(1-α)(1-β)·V(x1,y1)+α(1-β)·V(x2,y1)+(1-α)β·V(x1,y2)+αβ·V(x2,y2)

[0031] in is the normalized offset of the voxel in the x and y directions, V(x* ,y * ) is the midpoint of the slice image of the three-dimensional body on the z axis (x * ,y * )’s pixel value.

[0032] Furthermore, the reconstruction algorithms in S2 include Back Projection Technique (BPT), Weighted Back Projection (WBP), Filtered Back Projection (FBP), Simultaneous Iterative Reconstructive Technique (SIRT), Simultaneous Algebraic Reconstruction Technique (SART) and Alternating Direction Method of Multipliers (ADMM). The present invention designs differentiated implementation methods for different reconstruction algorithms, and the core goal is to efficiently complete the three-dimensional reconstruction task.

[0033] Furthermore, in S2, efficient function execution requires a reasonable grid configuration. Selecting the appropriate grid and block size for kernel function grid configuration is crucial to ensure the correctness of the calculation results and the efficient use of streaming multiprocessors (SMs) on the GPU. The optimal grid and block configuration can enhance parallelism, increase resource utilization, and thus improve overall performance.

[0034] Furthermore, the kernel function grid configuration of S2 specifically includes:

[0035] S2-1: Forward projection: A two-dimensional (2D) grid is used to generate the projection image. The size of the grid is determined by two parameters: the first dimension is calculated by dividing the number of voxels on the side of the volume by the number of threads per thread block, and the other dimension is the number of projections (if this value is 1, it is a 1D grid). To improve computational efficiency, each thread in the grid is responsible for calculating the projection of a specific part of the volume. These threads iterate along the y direction to calculate the projection value, thereby reducing the time required to calculate the projection coordinates.

[0036] S2-2: Difference calculation: A two-dimensional grid (2D) is used to calculate the difference between the calculated value and the actual projection. The dimension of the grid is determined by two parameters: the number of pixels in one dimension of the projected image and the number of projections in another dimension. Each thread is responsible for calculating the difference between a pixel point and its corresponding actual projection.

[0037] S2-3: Backprojection: A three-dimensional (3D) grid is used to perform the backprojection task. The dimensions of the grid are calculated based on three parameters: the width, thickness, and length of the volume. This ensures that the grid accurately represents the entire volume to be reconstructed. To optimize processing and avoid computational conflicts, each thread in the grid is assigned to a specific voxel and is responsible for calculating the backprojection value of that voxel at all projection angles. This approach improves efficiency by evenly distributing the computational load across multiple threads.

[0038] S2-4: Update Volume: Use a one-dimensional grid (1D) to update the three-dimensional volume. The number of threads in the grid is equal to the total number of voxels in the volume. Each thread is responsible for updating one voxel.

[0039] Furthermore, CUDA-based parallel computing: the highly parallel architecture of the GPU is used to efficiently calculate the input data and execute the core reconstruction algorithm; the reconstruction parameters can be customized: key parameters such as the number of iterations, relaxation factors, and geometric settings can be adjusted to optimize the reconstruction effect.

[0040] Compared with the prior art, the present invention has the following beneficial effects:

[0041] The present invention significantly improves the efficiency of the 3D reconstruction process while maintaining accuracy by adopting a parallel acceleration solution on the CPU and GPU sides. The acceleration solution mainly adopted by the present invention is divided into two parts: on the one hand, data preprocessing and file operations are optimized through parallel processing on the CPU side, thereby improving overall efficiency; on the other hand, on the GPU side, the core reconstruction tasks are accelerated through parallel computing, significantly improving processing speed and accuracy. Experimental results show that the present invention can significantly improve computing efficiency in most cases, especially in terms of computing time, showing superior performance compared to traditional methods.

[0042] The present invention improves the reconstruction speed through innovative design to cope with the challenge of increasing data volume; at the same time, it provides users with more flexible parameter adjustment space, so that they can finely control the reconstruction process according to different experimental conditions and sample characteristics to meet the needs of diverse experimental scenarios. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 It is a data processing workflow diagram of the present invention.

[0044] Figure 2 The following is a comparative analysis of the reconstruction results of SIRT technology using Fourier ring correlation (FRC) curves; wherein (A) real data. (B) reprojection generated by the present invention (reconstruction along the y axis). (C) reprojection generated by the present invention (reconstruction along the Z axis). (D) reprojection of a three-dimensional volume generated by IMOD. (E) reprojection of a three-dimensional volume generated by TOMO3D. (F) FRC curve. DETAILED DESCRIPTION

[0045] The technical solution of the present invention is further described below in conjunction with embodiments and drawings.

[0046] Example 1

[0047] Use SIRT to perform efficient 3D cryo-EM reconstruction. The experimental data sets include: BBb, EMPAIR-10453, EMPAIR-10045, and EMPAIR-10164. The specific process is as follows: Figure 1 As shown in the figure, first, distributed file reading is performed on the CPU through MPI, and then the geometric structure of the sample is adjusted. Subsequently, the computationally intensive reconstruction task is performed on the GPU to fully utilize its parallel computing capabilities, and two types of GPU memory management methods are adopted. Finally, the reconstructed volume data is distributed and output through the CPU. The method includes the following steps:

[0048] Step 1: Data input and preprocessing

[0049] Step 1.1 Input data: including the tilt sequence and its corresponding tilt angle after CTF correction and alignment.

[0050] Step 1.2 MPI acceleration mechanism:

[0051] (1) In this embodiment, MPI technology is used to perform multi-threaded parallel file operations, specifically for reading tilted sequences and writing merged reconstructed sub-volume data, effectively alleviating the I / O bottleneck.

[0052] (2) At the same time, MPI is also used for preprocessing tasks. In the SIRT method, three-dimensional volume data is aligned according to a customized angle to improve data quality and optimize the subsequent reconstruction process.

[0053] Step 2: Parallel reconstruction on GPU

[0054] Step 2.1 Data transfer: Transfer the preprocessed data from the CPU to the GPU.

[0055] Step 2.2 Customization of reconstruction parameters: Users can adjust key parameters such as the number of iterations, relaxation factors, and geometric settings to optimize the reconstruction effect.

[0056] Step 2.3 CUDA-based parallel computing: Utilize the highly parallel architecture of GPU to efficiently calculate the input data and execute the core reconstruction algorithm according to the user-defined parameters.

[0057] The formula of the SIRT reconstruction method used in this embodiment is:

[0058]

[0059] in, represents the intensity value of the i-th projected pixel at the k-th iteration, is the value of the j-th voxel at the k-th iteration, is the updated value of the j-th voxel at the k+1th iteration, and denotes the weight of the jth voxel’s contribution to the i-th projection pixel at the kth and k+1th iterations, respectively. i represents the measured intensity value of the i-th projection pixel, λ is the regularization parameter used to control the magnitude of the update at each iteration. n is the number of voxels used for projection, which is a subset of all voxels, N is the total number of voxels, and M is the number of projection pixels.

[0060] Based on the above formula and the grid configuration in S2 above, the specific function configuration of the SIRT method is expressed as:

[0061] S2-1: Forward Projection

[0062] For each projection i (computed in parallel), use the current volume estimate x j Perform forward projection to obtain simulated projection data

[0063]

[0064] This step maps the volume data to the projection space and obtains the predicted projection results.

[0065] S2-2: Difference calculation: Next, calculate the error of each projection i, that is, the true projection value p i With forward projection The result is normalized by the projection weight:

[0066]

[0067] This operation calculates the deviation of each projection and normalizes it.

[0068] S2-3: Back Projection

[0069] In the back-projection stage, for each voxel j (computed in parallel), the total projection error r i and weight w ij, Calculate the back-projection result:

[0070]

[0071] This process transfers the error information of the projection space back to the volume space, providing corrections for subsequent updates.

[0072] S2-4: Update 3D volume

[0073] Finally, using the back-projection result y j Update volume data. For each voxel j (parallel calculation), update according to the following formula:

[0074]

[0075] After multiple iterations, the reconstructed volume data x is finally obtained j .

[0076] Step 3: Data return and result combination

[0077] After the calculation of step 3.1 is completed, the reconstructed 3D volume data is transferred from the GPU back to the CPU.

[0078] Step 3.2 uses MPI to merge data to improve the computational efficiency of large-scale reconstruction tasks.

[0079] Step 4: Result output and storage management

[0080] In order to complete the reconstruction more efficiently, a slice-based memory management strategy is adopted, and MPI multi-threaded parallel writing is used.

[0081] Example 2

[0082] Use BPT for efficient 3D cryo-EM reconstruction. The experimental data sets include: BBb, EMPAIR-10453, EMPAIR-10045, and EMPAIR-10164. BPT is the simplest 3D volume reconstruction method. The specific steps are as follows:

[0083] Step 1: Data input and preprocessing

[0084] Step 1.1 Input data: including the tilt sequence and its corresponding tilt angle after CTF correction and alignment.

[0085] Step 1.2 MPI acceleration mechanism:

[0086] (1) In this embodiment, MPI technology is used to perform multi-threaded parallel file operations, specifically for reading tilted sequences and writing merged reconstructed sub-volume data, effectively alleviating the I / O bottleneck.

[0087] (2) At the same time, MPI is also used for preprocessing tasks. In BPT, the three-dimensional volume data is aligned according to the data input by the user. This process can enhance the high-frequency information in the projection data and help improve the edge quality of the reconstructed image.

[0088] Step 2: Parallel reconstruction on GPU

[0089] Step 2.1 Data transfer: Transfer the preprocessed data from the CPU to the GPU.

[0090] Step 2.2 Customization of reconstruction parameters: Users can adjust key parameters such as the number of iterations, relaxation factors, and geometric settings to optimize the reconstruction effect.

[0091] Step 2.3 CUDA-based parallel computing: Utilize the highly parallel architecture of GPU to efficiently calculate the input data and execute the core reconstruction algorithm according to the user-defined parameters.

[0092] According to the BPT algorithm, the reconstruction steps are:

[0093] (1) Back projection

[0094] For each voxel j in the volume (computed in parallel), the contributions are summed for all projection angles i:

[0095]

[0096] denominator+=w[i][j]

[0097] (2) Update the 3D body

[0098] After the accumulation is completed, each voxel j is normalized and the three-dimensional volume is updated:

[0099]

[0100] This step ensures that the value of each voxel is properly normalized, and finally the reconstructed 3D volume x is obtained.

[0101] Step 3: Data return and result combination

[0102] After the calculation of step 3.1 is completed, the reconstructed 3D volume data is transferred from the GPU back to the CPU.

[0103] Step 3.2 uses MPI to merge data to improve the computational efficiency of large-scale reconstruction tasks.

[0104] Step 4: Result output and storage management

[0105] In order to complete the reconstruction more efficiently, a slice-based memory management strategy is adopted, and MPI multi-threaded parallel writing is used.

[0106] Example 3

[0107] Use FBP for efficient 3D cryo-EM reconstruction. The experimental data sets include: BBb, EMPAIR-10453, EMPAIR-10045, and EMPAIR-10164. The FBP algorithm is a variant of the BPT algorithm. The specific steps are as follows:

[0108] Step 1: Data input and preprocessing

[0109] Step 1.1 Input data: including the tilt sequence and its corresponding tilt angle after CTF correction and alignment.

[0110] Step 1.2 MPI acceleration mechanism:

[0111] (1) In this embodiment, MPI technology is used to perform multi-threaded parallel file operations, specifically for reading tilted sequences and writing merged reconstructed sub-volume data, effectively alleviating the I / O bottleneck.

[0112] (2) At the same time, MPI is also used for preprocessing tasks, including RamLak filtering or Shepp-Logan filtering and alignment of three-dimensional volume data in FBP. This process can enhance the high-frequency information in the projection data and help improve the edge quality of the reconstructed image.

[0113] Step 2: Parallel reconstruction on GPU

[0114] Step 2.1 Data transfer: Transfer the preprocessed data from the CPU to the GPU.

[0115] Step 2.2 Customization of reconstruction parameters: Users can adjust key parameters such as the number of iterations, relaxation factors, and geometric settings to optimize the reconstruction effect.

[0116] Step 2.3 CUDA-based parallel computing: Utilize the highly parallel architecture of GPU to efficiently calculate the input data and execute the core reconstruction algorithm according to the user-defined parameters.

[0117] According to the FBP algorithm, the reconstruction steps are:

[0118] (1) Back projection

[0119] In this step, each voxel j in the volume is processed (using parallel computing) and the information of all projection angles is accumulated:

[0120]

[0121] denominator+=w[i][j]

[0122] The purpose of this operation is to transform the filtered projection data Reverse mapping to volume space, while counting the sum of weights received by each voxel, in preparation for subsequent normalization.

[0123] (2) Update the 3D body

[0124] Each voxel j is normalized to update the 3D volume, which is specifically implemented as follows:

[0125]

[0126] This step ensures that the value of each voxel is properly normalized, and finally the reconstructed 3D volume x is obtained.

[0127] Step 3: Data return and result combination

[0128] After the calculation of step 3.1 is completed, the reconstructed 3D volume data is transferred from the GPU back to the CPU.

[0129] Step 3.2 uses MPI to merge data to improve the computational efficiency of large-scale reconstruction tasks.

[0130] Step 4: Result output and storage management

[0131] In order to complete the reconstruction more efficiently, a slice-based memory management strategy is adopted, and MPI multi-threaded parallel writing is used.

[0132] Example 4

[0133] WBP is used for efficient 3D cryo-EM reconstruction. The experimental data sets include: BBb, EMPAIR-10453, EMPAIR-10045, and EMPAIR-10164. The WBP algorithm flow and specific implementation are the same as the FBP algorithm except for the filter selection and FBP algorithm, so they are not repeated here.

[0134] Example 5

[0135] SART is used for efficient 3D cryo-EM reconstruction. The experimental data sets include: BBb, EMPAIR-10453, EMPAIR-10045, and EMPAIR-10164. The SART algorithm uses an iterative update strategy to reconstruct volume data. The method includes the following steps:

[0136] Step 1: Data input and preprocessing

[0137] Step 1.1 Input data: including the tilt sequence and its corresponding tilt angle after CTF correction and alignment.

[0138] Step 1.2 MPI acceleration mechanism:

[0139] (1) In this embodiment, MPI technology is used to perform multi-threaded parallel file operations, specifically for reading tilted sequences and writing merged reconstructed sub-volume data, effectively alleviating the I / O bottleneck.

[0140] (2) At the same time, MPI is also used for preprocessing tasks. In the SART method, three-dimensional volume data is aligned according to a customized angle to improve data quality and optimize the subsequent reconstruction process.

[0141] Step 2: Parallel reconstruction on GPU

[0142] Step 2.1 Data transfer: Transfer the preprocessed data from the CPU to the GPU.

[0143] Step 2.2 Customization of reconstruction parameters: Users can adjust key parameters such as the number of iterations, relaxation factors, and geometric settings to optimize the reconstruction effect.

[0144] Step 2.3 CUDA-based parallel computing: Utilize the highly parallel architecture of GPU to efficiently calculate the input data and execute the core reconstruction algorithm according to the user-defined parameters.

[0145] According to the SART algorithm, the reconstruction steps are:

[0146] (1) Forward projection

[0147] In this step, all projections i in each subset s are executed in parallel, using the current volume data x j Calculate the simulated projection:

[0148]

[0149] The purpose of this step is to map the volume data to the projection space through the projection operator to generate simulated projection data.

[0150] (2) Difference calculation

[0151] After the forward projection is completed, the projection error is calculated for each projection i in the same subset (parallel calculation):

[0152]

[0153] This step is done by calculating the difference between the actual projection value and the simulated projection value r i , providing an error signal for subsequent correction.

[0154] (3) Back projection

[0155] In the back-projection stage, each voxel j in the volume space is processed (in parallel), and the error r of all projections in the corresponding subset is used i and weight w {ij}, Accumulate and get the correction value:

[0156]

[0157] This process “back propagates” the projection error from the projection space to the volume space, so that the correction value y of each voxel j It can reflect the combined effect of measurement errors in all relevant projections.

[0158] (4) Update the 3D body

[0159] Finally, the correction value y obtained by back projection j , update the volume data:

[0160]

[0161] This step updates each voxel j synchronously under parallel computing, so that the volume estimate is gradually corrected to be closer to the true value. When the set number of iterations is reached, the reconstructed 3D volume is obtained.

[0162] Step 3: Data return and result combination

[0163] After the calculation of step 3.1 is completed, the reconstructed 3D volume data is transferred from the GPU back to the CPU.

[0164] Step 3.2 uses MPI to merge data to improve the computational efficiency of large-scale reconstruction tasks.

[0165] Step 4: Result output and storage management

[0166] In order to complete the reconstruction more efficiently, a slice-based memory management strategy is adopted, and MPI multi-threaded parallel writing is used.

[0167] Example 6

[0168] ADMM is used for efficient 3D cryo-EM reconstruction. The experimental data sets include: BBb, EMPAIR-10453, EMPAIR-10045, and EMPAIR-10164. The ADMM algorithm regards 3D reconstruction as an optimization problem, and then transforms it into a constrained minimization problem by introducing auxiliary variables. Then, an iterative update strategy is used to reconstruct the volume data. The method includes the following steps:

[0169] Step 1: Data input and preprocessing

[0170] Step 1.1 Input data: including the tilt sequence and its corresponding tilt angle after CTF correction and alignment.

[0171] Step 1.2 MPI acceleration mechanism:

[0172] (1) In this embodiment, MPI technology is used to perform multi-threaded parallel file operations, specifically for reading tilted sequences and writing merged reconstructed sub-volume data, effectively alleviating the I / O bottleneck.

[0173] (2) At the same time, MPI is also used for preprocessing tasks. In the ADMM method, the 3D volume data is aligned according to a customized angle to improve data quality and optimize the subsequent reconstruction process.

[0174] Step 2: Parallel reconstruction on GPU

[0175] Step 2.1 Data transfer: Transfer the preprocessed data from the CPU to the GPU.

[0176] Step 2.2 Customization of reconstruction parameters: Users can adjust key parameters such as the number of iterations, relaxation factors, and geometric settings to optimize the reconstruction effect.

[0177] Step 2.3 CUDA-based parallel computing: Utilize the highly parallel architecture of GPU to efficiently calculate the input data and execute the core reconstruction algorithm according to the user-defined parameters.

[0178] According to the ADMM algorithm, the reconstruction steps are:

[0179] (1) Forward projection

[0180] Add the Lagrange parameter to the forward projection result and perform a soft threshold operation to obtain the ADMM auxiliary variable u:

[0181]

[0182] Where L represents the forward projection operator, x k+1 is the 3D volume data of this iteration, α k is the Lagrangian parameter of the kth iteration, ρ is the regularization penalty term, and λ is the custom penalty parameter. This step uses a soft threshold function to correct the current estimate to achieve a sparse effect and adjust the sparsity of the data.

[0183] (2) Interpolation calculation

[0184] The difference between the forward projection result and the ADMM auxiliary variable u is accumulated to obtain the Lagrangian parameter α k+1 :

[0185] α k+1 =α k +(Lx k+1 -u k+1 )

[0186] This operation updates the Lagrange multiplier α, which is used to track changes in constraints and adjust the optimization process.

[0187] (3) Update the 3D body:

[0188] Update the volume using the conjugate gradient method:

[0189]

[0190] Among them, the first term is the fitting term, which represents the difference between the current volume data x and the target; the second term is the regularization term, which is used to control the sparsity of the solution. After multiple iterations, the optimized volume data is finally obtained.

[0191] In order to compare the efficiency of 3D reconstruction, a variety of data sets and algorithms were used for testing, covering 6 commonly used reconstruction algorithms: WBP, FBP, SIRT, SART, BPT and ADMM. By comparing with IMOD, TOMO3D and RELION, the performance improvement of our invention under different algorithms was evaluated. At the same time, for algorithms such as SART, BPT and ADMM that do not have accelerated versions, the accelerated implementation of the present invention will be compared with the non-accelerated version. By using data sets of different sizes and comparing with commonly used integrated software such as Relion, the acceleration effect and performance of the present invention in various algorithms are comprehensively evaluated. The experimental results are shown in Table 1.

[0192] Table 1 Comparison of running time of different reconstruction methods

[0193]

[0194]

[0195] In Table 1, TiltRec represents the present invention, and TiltRec-cuda is the present invention reconstructed along the y-axis using CUDA;

[0196] TiltRecZ-cuda is the reconstruction using CUDA along the z-axis of the present invention; TiltRec-mpi is the reconstruction using only MPI along the y-axis of the present invention, that is, this method does not use GPU acceleration, and is only used as a comparison scheme; TiltRecZ-mpi is the reconstruction using only MPI along the z-axis of the present invention, that is, this method does not use GPU acceleration, and is only used as a comparison scheme. The tilted data representation in the table is limited by CPU memory and software implementation, and the data set is first downsampled twice before reconstruction.

[0197] According to the runtimes shown in Table 1, TiltRec-cuda consistently maintains the fastest or comparable computation times across almost all datasets and algorithms. This strategy provides significant performance gains, with speedup factors reaching hundreds of times compared to non-accelerated methods.

[0198] The final experimental accuracy is as follows Figure 2 (Comparison of FRC curves of the reconstruction of the present invention (along the y-axis), the reconstruction of the present invention (along the z-axis), the reprojection generated by IMOD and TOMO3D, and the true result is shown. Figure 2It can be seen that the present invention significantly improves the reconstruction efficiency without sacrificing the reconstruction accuracy.

[0199] Finally, although this specification is described according to implementation methods, not every implementation method contains only one independent technical solution. This narrative method of the specification is only for the sake of clarity. Those skilled in the art should regard the specification as a whole. The technical solutions in each embodiment can also be appropriately combined to form other implementation methods that can be understood by those skilled in the art.

Claims

1. A three-dimensional reconstruction method for a cryo-electron microscopy image sequence, characterized in that: The following steps are involved: S1: Acquire cryo-EM image sequences and perform preprocessing; S2: Perform 3D reconstruction on the preprocessed data. First, design the kernel function and use the computing unified device architecture CUDA to handle the core computing-intensive tasks in tomographic reconstruction. Then optimize the data reading efficiency through the GPU memory type and customize the CUDA kernel function based on the voxel scattering strategy. S3: After the calculation is completed, the reconstructed 3D volume data is merged using MPI to improve the computational efficiency of large-scale reconstruction tasks; a slice-based memory management strategy is adopted, and MPI multi-threading is used to write the merged reconstructed sub-volume data in parallel; and the reconstructed data is output.

2. The three-dimensional reconstruction method according to claim 1, characterized in that: In S1: the image sequence includes contrast transfer function correction CTF and aligned tilt sequence and its corresponding tilt angle; MPI technology is used to perform multi-threaded parallel file operations, and MPI is also used for preprocessing tasks, including filtering calculations and geometric adjustments of three-dimensional volume data, to improve data quality and optimize subsequent reconstruction processes.

3. The three-dimensional reconstruction method according to claim 2, characterized in that: The pre-processing comprises: (1) Filter calculation In the weighted back projection method WBP and the filtered back projection method FBP, filtering calculation is involved; WBP performs weighted processing on the projection image before back projection, and the filtering formula is as follows: Where H(·) is a weighting function used to assign weights to each projection; σ Nx is the standard deviation of the Gaussian decay part; weight is the predefined threshold; FBP provides two filters: RamLak filter and Shepp-Logan filter: The filtering formulas of the two methods are: Where H(·) is the filter function, where H RL (·) represents RamLak filter, H SL (·) represents the Shepp-Logan filter, ω represents the frequency, ω c is the cut-off frequency; (2) Geometry adjustment Applying the geometric model, the spatial positioning of the sample in 3D space is precisely adjusted by adjusting the tilt angle θ, the pitch angle φ and the Z-axis offset t. The adjustment formula is as follows:

4. The three-dimensional reconstruction method according to claim 1, characterized in that: In S2, in the custom CUDA kernel function, each voxel is distributed to the detector, and the pixel values ​​of four adjacent points are updated to generate a projection image; Subsequently, these four points are used in the projection data to update the corresponding voxels in the reconstructed volume; both projection and back-projection use transpose operations to enhance the accuracy of the results; assuming that the size of a voxel is Δx×Δy×Δz, and the center of the voxel is located at (x, y, z), the value of the voxel is V(x, y, z), then the contribution of the voxel to a pixel in the image plane is calculated by the interpolation function W(r): Where I(p) is the projection value of pixel p on the detector, V(r i ) is the voxel at position r i The value at W(r i ) is the voxel pixel r on the image plane i The interpolation function uses bilinear difference, that is, W(r i )=(1-α)(1-β)·V(x1,y1)+α(1-β)·V(x2,y1)+(1-α)β·V(x1,y2)+αβ·V(x2,y2) in is the normalized offset of the voxel in the x and y directions, V(x * ,y * ) is the midpoint of the slice image of the three-dimensional body on the z axis (x * ,y * )’s pixel value.

5. The three-dimensional reconstruction method according to claim 1, characterized in that: The reconstruction algorithms in S2 include direct back projection BPT, weighted back projection WBP, filtered back projection FBP, simultaneous iterative reconstruction technology SIRT, synchronous algebraic reconstruction technology SART and ADMM based on alternating multiplier method.

6. The three-dimensional reconstruction method according to claim 1, characterized in that: The kernel function grid configuration of S2 specifically includes: S2-1: Forward projection: A two-dimensional grid is used to generate the projection image; the size of the grid is determined by two parameters: the first dimension is calculated by dividing the number of voxels on the side of the volume by the number of threads per thread block, and the other dimension is the number of projections; to improve computational efficiency, each thread in the grid is responsible for calculating the projection of a specific part of the volume, and these threads iterate along the y direction to calculate the projection value, thereby reducing the time required to calculate the projection coordinates; S2-2: Difference calculation: Use a two-dimensional grid to calculate the difference between the calculated value and the actual projection; the dimension of the grid is determined by two parameters: the number of pixels in one dimension of the projected image and the number of projections in another dimension. Each thread is responsible for calculating the difference between a pixel point and its corresponding actual projection; S2-3: Backprojection: A three-dimensional grid is used to perform the backprojection task; the dimensions of the grid are calculated based on three parameters: the width, thickness, and length of the volume. Each thread in the grid is assigned to a specific voxel and is responsible for calculating the backprojection value of that voxel at all projection angles. S2-4: Update Volume: Use a 1D grid to update the 3D volume, where the number of threads in the grid is equal to the total number of voxels in the volume. Each thread is responsible for updating one voxel.

7. The three-dimensional reconstruction method according to claim 1, characterized in that: CUDA-based parallel computing: Utilizes the highly parallel architecture of the GPU to efficiently compute input data and execute core reconstruction algorithms; Reconstruction parameters can be customized: adjust the number of iterations, relaxation factors, and geometry to set key parameters to optimize the reconstruction effect.

Images