Striping noise removal method for ultra-large three-dimensional biological images based on CUDA acceleration

Abstract

The application relates to a CUDA acceleration-based stripe noise removal method for a super-large three-dimensional biological image, which comprises the following steps: constructing a frequency domain Gaussian band wave filter according to two-dimensional sequence slices of the three-dimensional image; dividing the two-dimensional sequence slices of the three-dimensional image into blocks to obtain a plurality of image blocks; realizing fast Fourier transform by using CUDA to obtain a frequency domain image of the image block; filtering; realizing inverse fast Fourier transform by using CUDA to obtain a spatial domain image of the image block; splicing the spatial domain image of each image block into an original image size; and processing the background of the spatial domain image of the original image size. The application utilizes the principle of dividing the original image into blocks, accelerates the calculation process of mutual conversion between the image spatial domain and the frequency domain, removes the stripe noise of the super-large three-dimensional image, effectively improves the operation efficiency, removes the stripe noise of the TB-level super-large three-dimensional image data, and obtains an acceleration ratio of about 9 times, thereby greatly shortening the processing time.

Description

Technical Field

[0001] This invention relates to the field of 3D biomedical image processing and parallel acceleration technology, and in particular to a method for removing stripe noise from ultra-large three-dimensional biological images based on CUDA acceleration. Background Technology

[0002] Biomedical imaging is a crucial intersection of medicine and engineering, providing visualizations of biological structures and functions through imaging techniques and computational methods. However, biomedical images are affected by various factors during the imaging process, including uneven light sources, sensor interference, sampling problems, electromagnetic interference, optical system defects, and environmental vibrations. These factors lead to visible stripe noise in the images, affecting their quality and accuracy. Therefore, stripe noise removal in biomedical images is essential.

[0003] Currently, there are two main methods for removing stripe noise from images: the first is to remove stripe noise directly in the spatial domain by using pixel-level image processing techniques such as median filtering and Gaussian filtering; the second is to convert the image from the spatial domain to the frequency domain, analyze the frequency differences between the signal and noise, and then apply appropriate filtering methods to achieve effective stripe removal. This is currently the most widely used method for stripe noise removal. However, both of these methods only perform well on images with low resolution.

[0004] With the rapid development of microscopic imaging technology, the data size of biological images is gradually increasing, even reaching the terabyte (TB) level. Due to limitations in computing resources and processing speed, stripe noise removal from ultra-large-scale images has become a challenge. Given the limitations of computer memory size, ultra-large-scale images cannot be directly loaded into computer memory, requiring processing methods to address this issue. Furthermore, as the data size increases, the time required for stripe noise removal also increases accordingly, necessitating strategies to accelerate noise removal and reduce time costs. Therefore, how to effectively remove stripe noise from ultra-large 3D images while making the entire process faster and more efficient has become an urgent technical problem to be solved. Summary of the Invention

[0005] To address the problems of ineffective stripe noise removal in ultra-large 3D images and long stripe noise removal time, the present invention aims to provide a CUDA-accelerated method for stripe noise removal in ultra-large 3D biological images that can effectively remove stripe noise from ultra-large 3D images and significantly improve the efficiency of image processing.

[0006] To achieve the above objectives, the present invention adopts the following technical solution: a method for removing stripe noise from ultra-large three-dimensional biological images based on CUDA acceleration, the method comprising the following sequential steps:

[0007] (1) Construct a frequency domain Gaussian strip notch filter based on the two-dimensional sequence slices of the three-dimensional image: For the two-dimensional sequence slices of the three-dimensional image, i.e. the original image, determine the width, rotation angle and cutoff frequency of the stripe noise in the spectrum diagram based on the spectral distribution of the stripe noise, and set the parameters of the frequency domain Gaussian strip notch filter, which include the width, rotation angle and cutoff frequency.

[0008] (2) Divide the two-dimensional sequence slices of the three-dimensional image into blocks according to a fixed size to obtain multiple image blocks: Divide the two-dimensional sequence slices of each three-dimensional image into blocks according to a fixed size, and record the starting position coordinates of all image blocks into which the two-dimensional sequence slices of each three-dimensional image are divided.

[0009] (3) Use CUDA to implement Fast Fourier Transform, and transform the image block from the spatial domain to the frequency domain through Fast Fourier Transform to obtain the frequency domain map of the image block.

[0010] (4) The frequency domain map of each image block is filtered by the constructed frequency domain Gaussian band notch filter to obtain the frequency domain map of the filtered image block.

[0011] (5) Use CUDA to implement inverse fast Fourier transform, and convert the frequency domain map of the filtered image block from the frequency domain back to the spatial domain through inverse fast Fourier transform to obtain the spatial domain map of the image block.

[0012] (6) Stitch the spatial domain maps of each image block to the original image size: Traverse the starting position coordinates of the image blocks recorded during block division, and stitch the spatial domain maps of each image block according to the starting coordinates (P... x ,P y Reassemble the images into a two-dimensional sequence of slices, and stitch all the image blocks together to form a spatial domain map image of the original image size;

[0013] (7) Process the background of the spatial domain map image of the original image size: extract the foreground contour of each spatial domain map image of the original image size to obtain a mask image, then fuse the mask image with the spatial domain map image of the original image size to obtain the foreground image, and refill the background pixel values ​​to remove background noise.

[0014] Step (2) specifically includes the following sequential steps:

[0015] (2a) Read any original image and obtain the dimensions of the X and Y axes of the original image, denoted as X and Y respectively. n Y nThen read the number of original images, denoted as Z. n Z n It is also the Z-axis dimension of the original image;

[0016] (2b) Traverse the Z-axis dimension Z of the original image n Divide the currently traversed two-dimensional sequence slice into blocks, and set the standard image block size [X]. b ,Y b ], so that the size of the image patch is X on the X and Y axes respectively. b Y b And X b =Y b Let N x N y These represent the number of image patches on the X and Y axes, respectively, and are defined as follows:

[0017]

[0018]

[0019] (2c) If the size exceeds the original image size, then the size of the image patch is... The specific formula is as follows:

[0020]

[0021]

[0022] (2d) The final number of image blocks is determined to be N based on the original image size. sum Its definition is as follows:

[0023] N sum =N x *N y *Z n

[0024] (2e) After the original image is divided into blocks, a coordinate system is set according to the origin of the original image. The original image is cropped according to the size of each image block to obtain each image block image, and the starting coordinates (P) of each image block are recorded. x ,P y ) are stored in a list, each image patch having a size of [X i ,Y j ], where i is the block index of the image patch on the X-axis, j is the block index of the image patch on the Y-axis, and X i Y j The definition is as follows:

[0025]

[0026]

[0027] Step (3) specifically includes the following steps:

[0028] (3a) First, let the length of the initial array A for the Fast Fourier Transform be L, and calculate the actual array A for the Fast Fourier Transform. T Length L T Its definition is:

[0029]

[0030] Where ceil represents rounding up, if L T If the value is greater than L, then the initial array A is padded with zeros to obtain the true array A. T The number of zeros is L. T -L;

[0031] (3b) Allocate L for CUDA kernel functions T There are 1024 threads, where the kernel function's thread block size is 1024, and each thread index corresponds to a real array A. T The index number, each thread is responsible for adding the actual array A. T An index number is converted from decimal to binary, and the binary index number is reversed to obtain a new decimal index number. After all threads have completed, the new decimal index is stored in the reversed array A. R middle;

[0032] (3c) Redesign the CUDA kernel function and allocate L to the kernel function. T There are 1024 threads, where the thread block size of the kernel function is 1024, and each thread index m corresponds to the actual array A. T The index number, if the thread index m is less than the index of the reversed array A R If the corresponding index n is used, then the actual array A will be... T The data at indices m and n are swapped. After all threads have finished, a new adjusted array A is obtained. F ;

[0033] (3d) Apply powers of 1 to L with a step size of 2. T Perform the traversal and record the current traversal value M. Calculate the first rotation factor W for each traversal. F Its definition is:

[0034] W F =cos(π÷M)+j*sin(π÷M)

[0035] Where j represents the imaginary part, the CUDA kernel function is redesigned, and L is assigned to the CUDA kernel function. TThere are 1024 threads, where the kernel function's thread block size is 1024. For each thread, index i is used to calculate the index j that needs to exchange data with the current index, and thread index i is updated to j. The second twitch factor W is then calculated. k Its definition is:

[0036]

[0037] Each thread adjusts the array A represented by the current thread index i. F The data in the array is processed using a butterfly operation. Once the traversal is complete, the actual array A is obtained. T The intermediate frequency domain array is obtained by fast Fourier transform. Remove the previously filled zero values ​​from the intermediate frequency domain array. After deleting, we obtain the frequency domain array F of the initial array A. A Complete the fast Fourier transform of the initial array A;

[0038] Step (5) specifically includes the following steps:

[0039] (5a) First, let the length of the initial frequency domain array F for the inverse fast Fourier transform be L. Calculate the actual frequency domain array F for the inverse fast Fourier transform in the same way as in step (3a). T The array length is L T ;

[0040] (5b) The adjusted frequency domain array F is calculated in the same manner as in steps (3b) and (3c). F ;

[0041] (5c) Apply powers of 1 to L with a step size of 2. T Perform the traversal and record the current traversal value M. Calculate the current third rotation factor W in each traversal. R Its definition is:

[0042] W R =cos(π÷M)-j*sin(π÷M)

[0043] Design CUDA kernel functions and allocate L to the kernel functions. T There are 1024 threads, where the kernel function's thread block size is 1024. For each thread, index i is used to calculate the index j that needs to exchange data with the current index, and thread index i is updated to j. The fourth twitch factor W is then calculated. S Its definition is:

[0044]

[0045] Each thread adjusts the frequency domain array F represented by the current thread index i. FThe data in the array is processed using a butterfly operation, and after traversal, the true frequency domain array F is obtained. T The intermediate real space array A T mid ;

[0046] (5d) Redesign the CUDA kernel function and assign L to the kernel function. FT There are 1024 threads, where the kernel function's thread block size is 1024. Each thread is responsible for processing the intermediate real space array A. T mid A data divided by L FT Normalization is performed, and after all threads have finished processing, an intermediate space array is obtained. Remove the previously filled zero value from After deleting, we obtain the spatial domain array A of the initial frequency domain array F. F The inverse fast Fourier transform of the initial frequency domain array F is completed.

[0047] Step (7) specifically includes the following steps:

[0048] (7a) First, calculate the spatial domain map image I of the original image size. o Pixel mean V mean The image is then thresholded, separating pixels with values ​​greater than V. mean The pixel value is set to 255, which is less than V. mean Set to 0, then extract the contour of the thresholded image to obtain contour map I. m ;

[0049] (7b) Calculate the contour map I m Find the area and index of all contours in the dataset, and reorder the contour indices according to their area size. m Keep outlines with an area greater than 50, and delete other outlines;

[0050] (7c) Traversal I m For each column in the image, calculate the number of non-zero pixels in that column. If the number of non-zero pixels in that column exceeds 90% of the number of rows in the image, set all pixel values ​​in that column to 0 to eliminate some vertical lines or noise in the image.

[0051] (7d) Construction and I o Blank image space I of the same size v , to I m Perform pixel traversal; if the current position P (x,y) If the pixel value is 255, then I v At position P (x,y) The pixel value is set to I o At position P (x,y) Pixel values; if P(x,y) If the pixel value is 0, then I v At position P (x,y) The pixel value is set to I o The background pixel values, finally obtained as I v This is the image after stripes have been removed.

[0052] As can be seen from the above technical solution, the beneficial effects of the present invention are as follows: First, the present invention utilizes the principle of dividing the original image into blocks and uses the parallel architecture of the GPU to accelerate the calculation process of mutual conversion between the spatial domain and frequency domain of the image, thereby effectively improving the running efficiency while removing stripe noise from ultra-large three-dimensional images; Second, experiments have shown that the present method can remove stripe noise from ultra-large three-dimensional image data at the TB level on a platform with an i5-10400F CPU, an RTX-3060 GPU, and 16G of memory, and achieves an acceleration ratio of approximately 9 times, greatly shortening the processing time. Attached Figure Description

[0053] Figure 1 This is a flowchart of the method of the present invention;

[0054] Figure 2 Microscopic image before stripe removal;

[0055] Figure 3 This is a microscopic image after stripes have been removed according to the present invention. Detailed Implementation

[0056] like Figure 1 As shown, a method for removing stripe noise from ultra-large 3D biological images based on CUDA acceleration is presented. This method includes the following sequential steps:

[0057] (1) Construct a frequency domain Gaussian strip notch filter based on the two-dimensional sequence slices of the three-dimensional image: For the two-dimensional sequence slices of the three-dimensional image, i.e. the original image, determine the width, rotation angle and cutoff frequency of the stripe noise in the spectrum diagram based on the spectral distribution of the stripe noise, and set the parameters of the frequency domain Gaussian strip notch filter, which include the width, rotation angle and cutoff frequency.

[0058] (2) Divide the two-dimensional sequence slices of the three-dimensional image into blocks according to a fixed size to obtain multiple image blocks: Divide the two-dimensional sequence slices of each three-dimensional image into blocks according to a fixed size, and record the starting position coordinates of all image blocks into which the two-dimensional sequence slices of each three-dimensional image are divided.

[0059] (3) Use CUDA to implement Fast Fourier Transform, and transform the image block from the spatial domain to the frequency domain through Fast Fourier Transform to obtain the frequency domain map of the image block.

[0060] (4) The frequency domain map of each image block is filtered by the constructed frequency domain Gaussian band notch filter to obtain the frequency domain map of the filtered image block.

[0061] (5) Use CUDA to implement inverse fast Fourier transform, and convert the frequency domain map of the filtered image block from the frequency domain back to the spatial domain through inverse fast Fourier transform to obtain the spatial domain map of the image block.

[0062] (6) Stitch the spatial domain maps of each image block to the original image size: Traverse the starting position coordinates of the image blocks recorded during block division, and stitch the spatial domain maps of each image block according to the starting coordinates (P... x ,P y Reassemble the images into a two-dimensional sequence of slices, and stitch all the image blocks together to form a spatial domain map image of the original image size;

[0063] (7) Process the background of the spatial domain map image of the original image size: extract the foreground contour of each spatial domain map image of the original image size to obtain a mask image, then fuse the mask image with the spatial domain map image of the original image size to obtain the foreground image, and refill the background pixel values ​​to remove background noise.

[0064] Step (2) specifically includes the following sequential steps:

[0065] (2a) Read any original image and obtain the dimensions of the X and Y axes of the original image, denoted as X and Y respectively. n Y n Then read the number of original images, denoted as Z. n Z n It is also the Z-axis dimension of the original image;

[0066] (2b) Traverse the Z-axis dimension Z of the original image n Divide the currently traversed two-dimensional sequence slice into blocks, and set the standard image block size [X]. b ,Y b ], so that the size of the image patch is X on the X and Y axes respectively. b Y b And X b =Y b Let N x N y These represent the number of image patches on the X and Y axes, respectively, and are defined as follows:

[0067]

[0068]

[0069] (2c) If the size exceeds the original image size, then the size of the image patch is... The specific formula is as follows:

[0070]

[0071]

[0072] (2d) The final number of image blocks is determined to be N based on the original image size. sum Its definition is as follows:

[0073] N sum =N x *N y *Z n

[0074] (2e) After the original image is divided into blocks, a coordinate system is set according to the origin of the original image. The original image is cropped according to the size of each image block to obtain each image block image, and the starting coordinates (P) of each image block are recorded. x ,P y ) are stored in a list, each image patch having a size of [X i ,Y j ], where i is the block index of the image patch on the X-axis, j is the block index of the image patch on the Y-axis, and X i Y j The definition is as follows:

[0075]

[0076]

[0077] Step (3) specifically includes the following steps:

[0078] (3a) First, let the length of the initial array A for the Fast Fourier Transform be L, and calculate the actual array A for the Fast Fourier Transform. T Length L T Its definition is:

[0079]

[0080] Where ceil represents rounding up, if L T If the value is greater than L, then the initial array A is padded with zeros to obtain the true array A. T The number of zeros is L. T -L;

[0081] (3b) Allocate L for CUDA kernel functions T There are 1024 threads, where the kernel function's thread block size is 1024, and each thread index corresponds to a real array A. T The index number, each thread is responsible for adding the actual array A. TAn index number is converted from decimal to binary, and the binary index number is reversed to obtain a new decimal index number. After all threads have completed, the new decimal index is stored in the reversed array A. R middle;

[0082] (3c) Redesign the CUDA kernel function and allocate L to the kernel function. T There are 1024 threads, where the kernel function's thread block size is 1024, and each thread index m corresponds to the actual array A. T The index number, if the thread index m is less than the index of the reversed array A R If the corresponding index n is used, then the actual array A will be... T The data at indices m and n are swapped. After all threads have finished, a new adjusted array A is obtained. F ;

[0083] (3d) Apply powers of 1 to L with a step size of 2. T Perform the traversal and record the current traversal value M. Calculate the first rotation factor W for each traversal. F Its definition is:

[0084] W F =cos(π÷M)+j*sin(π÷M)

[0085] Where j represents the imaginary part, the CUDA kernel function is redesigned, and L is assigned to the CUDA kernel function. T There are 1024 threads, where the kernel function's thread block size is 1024. For each thread, index i is used to calculate the index j that needs to exchange data with the current index, and thread index i is updated to j. The second twitch factor W is then calculated. k Its definition is:

[0086]

[0087] Each thread adjusts the array A represented by the current thread index i. F The data in the array is processed using a butterfly operation. Once the traversal is complete, the actual array A is obtained. T The intermediate frequency domain array is obtained by fast Fourier transform. Remove the previously filled zero values ​​from the intermediate frequency domain array. After deleting, we obtain the frequency domain array F of the initial array A. A Complete the fast Fourier transform of the initial array A;

[0088] Step (5) specifically includes the following steps:

[0089] (5a) First, let the length of the initial frequency domain array F for the inverse fast Fourier transform be L. Calculate the actual frequency domain array F for the inverse fast Fourier transform in the same way as in step (3a). T The array length is L T ;

[0090] (5b) The adjusted frequency domain array F is calculated in the same manner as in steps (3b) and (3c). F ;

[0091] (5c) Apply powers of 1 to L with a step size of 2. T Perform the traversal and record the current traversal value M. Calculate the current third rotation factor W in each traversal. R Its definition is:

[0092] W R =cos(π÷M)-j*sin(π÷M)

[0093] Design CUDA kernel functions and allocate L to the kernel functions. T There are 1024 threads, where the kernel function's thread block size is 1024. For each thread, index i is used to calculate the index j that needs to exchange data with the current index, and thread index i is updated to j. The fourth twitch factor W is then calculated. S Its definition is:

[0094]

[0095] Each thread adjusts the frequency domain array F represented by the current thread index i. F The data in the array is processed using a butterfly operation, and after traversal, the true frequency domain array F is obtained. T The intermediate real space array A T mid ;

[0096] (5d) Redesign the CUDA kernel function and assign L to the kernel function. FT There are 1024 threads, where the kernel function's thread block size is 1024. Each thread is responsible for processing the intermediate real space array A. T mid A data divided by L FT Normalization is performed, and after all threads have finished processing, an intermediate space array is obtained. Remove the previously filled zero value from After deleting, we obtain the spatial domain array A of the initial frequency domain array F. F The inverse fast Fourier transform of the initial frequency domain array F is completed.

[0097] Step (7) specifically includes the following steps:

[0098] (7a) First, calculate the spatial domain map image I of the original image size. o Pixel mean V mean The image is then thresholded, separating pixels with values ​​greater than V. mean The pixel value is set to 255, which is less than V. mean Set to 0, then extract the contour of the thresholded image to obtain contour map I. m ;

[0099] (7b) Calculate the contour map I m Find the area and index of all contours in the dataset, and reorder the contour indices according to their area size. m Keep outlines with an area greater than 50, and delete other outlines;

[0100] (7c) Traversal I m For each column in the image, calculate the number of non-zero pixels in that column. If the number of non-zero pixels in that column exceeds 90% of the number of rows in the image, set all pixel values ​​in that column to 0 to eliminate some vertical lines or noise in the image.

[0101] (7d) Application and I o Blank image space I of the same size v , to I m Perform pixel traversal; if the current position P (x,y) If the pixel value is 255, then I v At position P (x,y) The pixel value is set to I o At position P (x,y) Pixel values; if P (x,y) If the pixel value is 0, then I v At position P (x,y) The pixel value is set to I o The background pixel values, finally obtained as I v This is the image after stripes have been removed.

[0102] like Figure 2 As shown, the image is a two-dimensional slice of a three-dimensional tree shrew image acquired using a microscope. The image contains vertical stripe noise.

[0103] like Figure 3 As shown, after processing with the present invention, it can be observed that the image after stripe removal is very clear, and no new false stripes appear.

[0104] As shown in Table 1, Table 1 is a schematic table of the running time and processed image size results of the present invention.

[0105] Table 1

[0106]

[0107] In summary, this invention utilizes the principle of segmenting the original image and leverages the parallel architecture of the GPU to accelerate the computational process of converting between the spatial and frequency domains of the image. This effectively improves operational efficiency while removing stripe noise from ultra-large 3D images. Experiments demonstrate that this method, on a platform with an i5-10400F CPU, an RTX-3060 GPU, and 16GB of memory, can remove stripe noise from TB-level ultra-large 3D image data, achieving an approximately 9-fold speedup and significantly reducing processing time.

Claims

1. A method for removing stripe noise from ultra-large three-dimensional biological images based on CUDA acceleration, characterized in that: The method includes the following steps in sequence: (1) Construct a frequency domain Gaussian strip notch filter based on the two-dimensional sequence slices of the three-dimensional image: For the two-dimensional sequence slices of the three-dimensional image, i.e. the original image, determine the width, rotation angle and cutoff frequency of the stripe noise in the spectrum diagram based on the spectral distribution of the stripe noise, and set the parameters of the frequency domain Gaussian strip notch filter, which include the width, rotation angle and cutoff frequency. (2) Divide the two-dimensional sequence slices of the three-dimensional image into blocks according to a fixed size to obtain multiple image blocks: Divide the two-dimensional sequence slices of each three-dimensional image into blocks according to a fixed size, and record the starting position coordinates of all image blocks into which the two-dimensional sequence slices of each three-dimensional image are divided. (3) Use CUDA to implement Fast Fourier Transform, and transform the image block from the spatial domain to the frequency domain through Fast Fourier Transform to obtain the frequency domain map of the image block. (4) The frequency domain map of each image block is filtered by the constructed frequency domain Gaussian band notch filter to obtain the frequency domain map of the filtered image block. (5) Use CUDA to implement inverse fast Fourier transform, and convert the frequency domain map of the filtered image block back to the spatial domain through inverse fast Fourier transform to obtain the spatial domain map of the image block. (6) Stitch the spatial domain maps of each image block to the original image size: Traverse the starting position coordinates of the image blocks recorded during block division, and stitch the spatial domain maps of each image block according to the starting coordinates (P... x ,P y Reassemble the images into a two-dimensional sequence of slices, and stitch all the image blocks together to form a spatial domain map image of the original image size; (7) Process the background of the spatial domain map image of the original image size: extract the foreground contour of each spatial domain map image of the original image size to obtain a mask image, then fuse the mask image with the spatial domain map image of the original image size to obtain the foreground image, and refill the background pixel values ​​to remove background noise.

2. The method for removing stripe noise from ultra-large three-dimensional biological images based on CUDA acceleration according to claim 1, characterized in that: Step (2) specifically includes the following sequential steps: (2a) Read any original image and obtain the dimensions of the X and Y axes of the original image, denoted as X and Y respectively. n Y n Then read the number of original images, denoted as Z. n Z n It is also the Z-axis dimension of the original image; (2b) Traverse the Z-axis dimension Z of the original image n Divide the currently traversed two-dimensional sequence slice into blocks, and set the standard image block size [X]. b ,Y b ], so that the size of the image patch is X on the X and Y axes respectively. b Y b And X b =Y b Let N x N y These represent the number of image patches on the X and Y axes, respectively, and are defined as follows: (2c) If the size exceeds the original image size, then the size of the image patch is [X]. b1 ,Y b1 The specific formula is as follows: (2d) The final number of image blocks is determined to be N based on the original image size. sum Its definition is as follows: N sum =N x *N y *Z n (2e) After the original image is divided into blocks, a coordinate system is set according to the origin of the original image. The original image is cropped according to the size of each image block to obtain each image block image, and the starting coordinates (P) of each image block are recorded. x ,P y ) are stored in a list, each image patch having a size of [X i ,Y j ], where i is the block index of the image patch on the X-axis, j is the block index of the image patch on the Y-axis, and X i Y j The definition is as follows:

3. The method for removing stripe noise from ultra-large three-dimensional biological images based on CUDA acceleration according to claim 1, characterized in that: Step (3) specifically includes the following steps: (3a) First, let the length of the initial array A for the Fast Fourier Transform be L, and calculate the actual array A for the Fast Fourier Transform. T Length L T Its definition is: Where ceil represents rounding up, if L T If the value is greater than L, then the initial array A is padded with zeros to obtain the true array A. T The number of zeros is L. T -L; (3b) Allocate L for CUDA kernel functions T There are 1024 threads, where the thread block size of the kernel function is 1024, and each thread index corresponds to the actual array A. T The index number, each thread is responsible for adding the actual array A. T An index number is converted from decimal to binary, and the binary index number is reversed to obtain a new decimal index number. After all threads have completed, the new decimal index is stored in the reversed array A. R middle; (3c) Redesign the CUDA kernel function and allocate L to the kernel function. T There are 1024 threads, where the kernel function's thread block size is 1024, and each thread index m corresponds to the actual array A. T The index number, if the thread index m is less than the index of the reversed array A R If the corresponding index n is used, then the actual array A will be... T The data at indices m and n are swapped. After all threads have finished, a new adjusted array A is obtained. F ; (3d) Apply powers of 1 to L with a step size of 2. T Perform the traversal and record the current traversal value M. Calculate the first rotation factor W for each traversal. F Its definition is: W F =cos(π÷M)+j*sin(π÷M) Where j represents the imaginary part, the CUDA kernel function is redesigned, and L is assigned to the CUDA kernel function. T There are 1024 threads, where the kernel function's thread block size is 1024. For each thread, index i is used to calculate the index j that needs to exchange data with the current index, and thread index i is updated to j. The second twitch factor W is then calculated. k Its definition is: Each thread adjusts the array A represented by the current thread index i. F The data in the array is processed using a butterfly operation. Once the traversal is complete, the actual array A is obtained. T The intermediate frequency domain array is obtained by fast Fourier transform. Remove the previously filled zero values ​​from the intermediate frequency domain array. After deleting, we obtain the frequency domain array F of the initial array A. A Complete the fast Fourier transform of the initial array A; Step (5) specifically includes the following steps: (5a) First, let the length of the initial frequency domain array F for the inverse fast Fourier transform be L. Calculate the actual frequency domain array F for the inverse fast Fourier transform in the same way as in step (3a). T The array length is L T ; (5b) The adjusted frequency domain array F is calculated in the same manner as in steps (3b) and (3c). F ; (5c) Apply powers of 1 to L with a step size of 2. T Perform the traversal and record the current traversal value M. Calculate the current third rotation factor W in each traversal. R Its definition is: W R =cos(π÷M)-j*sin(π÷M) Design CUDA kernel functions and allocate L to the kernel functions. T There are 1024 threads, where the kernel function's thread block size is 1024. For each thread, index i is used to calculate the index j that needs to exchange data with the current index, and thread index i is updated to j. The fourth twitch factor W is then calculated. S Its definition is: Each thread adjusts the frequency domain array F represented by the current thread index i. F The data in the array undergoes a butterfly operation, and after traversal, the true frequency domain array F is obtained. T The intermediate real space array A T mid ; (5d) Redesign the CUDA kernel function and assign L to the kernel function. FT There are 1024 threads, where the kernel function's thread block size is 1024. Each thread is responsible for processing the intermediate real space array A. T mid A data divided by L FT Normalization is performed, and after all threads have finished processing, an intermediate space array is obtained. Remove the previously filled zero value from After deleting, we obtain the spatial domain array A of the initial frequency domain array F. F The inverse fast Fourier transform of the initial frequency domain array F is completed.

4. The method for removing stripe noise from ultra-large three-dimensional biological images based on CUDA acceleration according to claim 1, characterized in that: Step (7) specifically includes the following steps: (7a) First, calculate the spatial domain map image I of the original image size. o Pixel mean V mean The image is then thresholded, separating pixels with values ​​greater than V. mean The pixel value is set to 255, which is less than V. mean Set to 0, then extract the contour of the thresholded image to obtain contour map I. m ; (7b) Calculate the contour map I m Find the area and index of all contours in the dataset, and reorder the contour indices according to their area size, then assign the I... m Keep outlines with an area greater than 50, and delete other outlines; (7c) Traversal I m For each column in the image, calculate the number of non-zero pixels in that column. If the number of non-zero pixels in that column exceeds 90% of the number of rows in the image, set all pixel values ​​in that column to 0 to eliminate some vertical lines or noise in the image. (7d) Construction and I o Blank image space I of the same size v , to I m Perform pixel traversal; if the current position P (x,y) If the pixel value is 255, then I v At position P (x,y) The pixel value is set to I o At position P (x,y) Pixel values; if P (x,y) If the pixel value is 0, then I v At position P (x,y) The pixel value is set to I o The background pixel values, finally obtained as I v This is the image after stripes have been removed.

Images