A CUDA-based deconvolution image restoration method

Abstract

The application discloses a CUDA-based deconvolution image restoration method, and has the characteristics that fixed shared memory is allocated in a host, Device memory is allocated in a graphics processor, and collected original image data is stored in the fixed shared memory and the Device memory; the cut-off frequency radius of the original image is calculated; the required Block thread block number is designed; according to the cut-off frequency radius, an imaging system optical transfer function template is generated, and after Fourier transform, the optical transfer function template is stored in the Device memory; Fourier transform is performed on original image data to be deconvoluted, according to the Block thread block number obtained above and the optical transfer function template after Fourier transform, deconvolution processing calculation is performed on the original image data after Fourier transform, and finally, a restored image is obtained, and the advantage is that compared with traditional deconvolution processing, the processing time is greatly shortened under the same number of convolutions.

Description

TECHNICAL FIELD

[0001] The present application relates to a deconvolution image restoration method, in particular to a CUDA-based deconvolution image restoration method. BACKGROUND

[0002] With the rapid development of dynamic microscopic imaging technology and modern life science research, it is essential to explore the dynamic process of biology from the perspective of cell dynamics. Live cell imaging technology can capture the dynamic changes of single cells, subcellular structures and organisms in vivo, and provide effective time and spatial information, which provides an important technical means for the research of cell biology, medicine, developmental biology and neuroscience. In the image acquisition process of live cell imaging, the image sensor is affected by various complex factors, such as the continuous dynamic changes of cells, so that the image acquired by the image sensor will inevitably have blurred area information, which greatly affects the image imaging effect and subsequent result analysis. Therefore, it is essential to restore the blurred image to its original appearance.

[0003] The current common image restoration methods mainly include Wiener filtering, inverse filtering restoration, constrained least squares method restoration, Lucy-Richardson deconvolution algorithm restoration and blind deconvolution restoration. The Lucy-Richardson deconvolution algorithm is referred to as (R-L) algorithm, which is one of the classical algorithms of image restoration method. Its outstanding advantages are low requirement for prior knowledge of image degradation process, superior restoration effect and the like. The R-L algorithm has the best restoration effect on low-noise motion blur images, and is widely used in the image post-processing part of live cell imaging technology. However, the R-L algorithm requires complex matrix calculation and large equation solving, which requires a large amount of operation time and memory for processing, and the traditional CPU architecture is difficult to meet the demand of real-time image acquisition-image processing-image display of live cells.

[0004] CUDA (Compute Unified Device Architecture) is a platform and programming model launched by NVIDIA, a graphics card manufacturer. By utilizing the powerful graphics processing capability of GPU, the computing performance and speed are greatly improved, which is especially suitable for image computing field. SUMMARY

[0005] The technical problem to be solved by the present application is to provide a CUDA-based deconvolution image restoration method with short time consumption and obvious image restoration effect.

[0006] The technical scheme adopted by the present application to solve the above technical problem is: a CUDA-based deconvolution image restoration method, comprising the following steps:

[0007] Step (1): the host assigns a fixed shared memory as a temporary area for raw image data and final result image;

[0008] Step (2): a Device memory is assigned in the graphics processor for saving optical transfer function template, image data, complex plane required by Fourier transform and result image data;

[0009] Step (3): raw image is collected by using laser as light source and image collection device matched with the optical microscope, and the collected raw image data is stored in the fixed shared memory and the Device memory of the graphics processor;

[0010] Step (4): the cut-off frequency radius of the raw image is calculated according to the numerical aperture of the objective lens, the physical size of a single pixel point of the raw image, the side length of the raw image and the wavelength of the laser;

[0011] Step (5): the number of Thread threads required in the operation process of the graphics processor is calculated according to the data size, and the number of Block thread blocks required is designed;

[0012] Step (6): the imaging system optical transfer function template is generated according to the cut-off frequency radius calculated in step (4), and the template image is stored in the Device memory of the graphics processor after Fourier transform;

[0013] Step (7): the raw image data to be deconvoluted is obtained from the fixed shared memory and Fourier transformed, and the Fourier transformed raw image data is deconvoluted and calculated in the graphics processor, wherein the Fourier transformed optical transfer function template obtained in step (6) is used in the deconvolution calculation, the number of Block thread blocks used in the CUDA operation is calculated in step (5), until the result image data is obtained, the result image data is converted from the frequency domain to the spatial domain, and finally the result image in the real plane is obtained and stored in the Device memory of the graphics processor;

[0014] Step (8): the result image in the Device memory of the graphics processor is copied to the fixed shared memory to obtain the final restored image.

[0015] Compared with the prior art, the CUDA parallel acceleration R-L deconvolution method used in the application greatly improves the operation speed, and can guarantee a processing speed of about 71 frames / s in a large field of view (2048x2048). As can be seen from Table 1, the deconvolution processing method used in the embodiment of the application greatly shortens the processing time compared with the traditional deconvolution processing under the condition of the same number of convolutions.

[0016] Table 1

[0017]

[0018]

[0019] The table takes the iteration number as two times for comparison.

[0020] Since the operation speed of the application is extremely fast, sufficient time is reserved for further optimization processing of the image or other processing operations such as region recognition, thereby ensuring the continuity of the video stream mode, and the traditional method cannot achieve the above effects.

[0021] The calculation method of the cut-off frequency radius of the original image in step (4) is as follows:

[0022] Let the cut-off frequency radius of the original image be ρ c , and the calculation formula is as follows:

[0023] ρ c = ρ max × P size × N

[0024] Wherein, ρ max is the theoretical cut-off frequency of the microscope optical system, NA is the numerical aperture, λ is the wavelength, P size is the physical size of a single pixel point of the image, Wherein, U Pixel is the physical size of a single pixel point of the camera unit, M obj is the magnification of the objective lens, M eye is the magnification of the microscope eyepiece, and N is the side length of the original image.

[0025] The calculation method of the Block thread block number designed in step (5) is as follows:

[0026] Step (5-1): Calculate the number of parallel threads:

[0027] Thread number = blockIdx.x*blockDim.x + threadIdx.x

[0028] Step (5-2): Obtain the Thread number and Block thread block number used by querying the performance of the current graphics card, let the maximum number of threads contained in each Block thread block be S, and set the original image to be square, the total number of image pixel points N×N as the required total thread number, then the Block thread block number is

[0029] The specific steps of step (6) are as follows:

[0030] Step (6-1): according to the length of the original image data to be processed, a plurality of thread blocks are divided in the Device memory grid of the graphic processor, and a kernel function is started to calculate the image cut-off frequency radius ρ c is brought into the real plane allocated by the graphic processor, and after processing and calculation, a binary optical transfer function OTF plane template is obtained, denoted as Mask otf ; the plane Mask otf is traversed, and if the distance between any pixel point and the center is less than the cut-off frequency radius ρ c , the pixel value of the pixel point is assigned as 1; if the distance between any pixel point and the center is greater than the cut-off frequency radius ρ c , the pixel value of the pixel point is assigned as 0; the basic method and formula are as follows:

[0031]

[0032] In the above formula, represents the gray value of the pixel point with the horizontal coordinate value x otf and the vertical coordinate value y j in the plane Mask j , j = 0, 1, 2, …, Len-1; Len is the length of the original image data, so as to obtain the optical transfer function OTF template;

[0033] Step (6-2): Fourier transform the optical transfer template Mask otf in the graphic processor to obtain its frequency domain template, denoted as Mask otf_cpl , and store it in the Device memory of the graphic processor.

[0034] The specific method of step (7) is as follows:

[0035] Step (7-1): use the cudaMemcpyAsync function to copy the original image data in the fixed shared memory to the Device memory of the graphic processor, denoted as Buf device ;

[0036] Step (7-2): use the fast Fourier library in CUDA to perform Fourier transform on the original image data Buf device to obtain the complex plane expression of the original image data in the frequency domain, denoted as Buf device_cpl ;

[0037] Step (7-3): use the optical transfer function frequency domain template Mask otf_cpl as the optical transfer function template for deconvolution operation in the graphic processor to perform Fourier transform on the original image data Buf device_cplThe deconvolution operation is performed:

[0038] First, the image degradation model is described according to the deconvolution principle as follows:

[0039]

[0040] where f(x, y) represents the original image, h(x, y) represents the point spread function of the optical system, n(x, y) represents the additive noise, and g(x, y) represents the degraded blurred image;

[0041] Then, the original image f(x, y) is iteratively operated using the following formula,

[0042]

[0043] to obtain the result image data, where is a convolution operator, T is a matrix transpose symbol, n represents the number of iterations, and the convolution iteration number n is preferably 2-3 times to obtain the best image quality.

[0044] Step (7-4): The result image data is converted from the frequency domain to the spatial domain to finally obtain the result image in the real plane, denoted as DM real and stored in the Device memory of the graphic processor. BRIEF DESCRIPTION OF DRAWINGS

[0045] Figure 1 is a flowchart of the CUDA-based deconvolution image restoration method of the present application;

[0046] Figure 2 is a flowchart of the GPU partial correlation calculation method of the present application;

[0047] Figure 3 is the original image of a biological sample according to an embodiment of the present application;

[0048] Figure 4 is the restored image after the original image of a biological sample is processed using the restoration method of the present application. Figure 3 is the restored image after the original image of a biological sample is processed using the restoration method of the present application. DETAILED DESCRIPTION

[0049] The present application is further described in detail below with reference to the embodiments and the accompanying drawings.

[0050] Embodiment:

[0051] A CUDA-based deconvolution image restoration method, comprising the following steps:

[0052] Step (1): The host Host allocates fixed shared memory as a staging area for raw image data and final result images; the image data type is ushort, and the corresponding type length is sizeof(ushort); the total length of the fixed shared memory is sizeof(ushort) x N x N x m, denoted as Buf host [m], where m = 3; mainly includes the following contents:

[0053] Step (1-1): The host Host memory is used for fixed shared memory in the real-time process of image acquisition-image processing-image display, forming a circulating shared memory module, i.e., a circular buffer area;

[0054] Step (1-2): Use the cudaMallocHost() function in CUDA C / C++ to allocate fixed memory, and the image data type is consistent with the above image data type, which is ushort type, and the corresponding type length is sizeof(ushort);

[0055] Step (1-3): If the image side length N is 1024, then the length corresponding to a single image is:

[0056] L img = sizeof(ushort) x 1024 x 1024; three fixed memory spaces of the same size are applied, Buf host [0], Buf host [1], and Buf host [2], which is used for GPU to prepare input data, and its advantage is not to participate in the page exchange in memory, improving system efficiency and further improving image processing speed;

[0057] Step (2): Allocate Device memory in the graphics processor for saving optical transfer function templates used in the deconvolution operation process, image data, complex plane required for Fourier transform, and result image data; the detailed steps mainly include the following contents:

[0058] Step (2-1): Use the cudaMalloc() function in CUDA C / C++ to allocate GPU memory space;

[0059] Step (2-2): The calculation of the optical transfer function (OTF) template used in image processing requires real plane and complex plane, and corresponding size memory is applied for subsequent OTF template calculation, denoted as Mask otf and Mask otf_cpl ;

[0060] Step (2-3): The original image data needs Fourier transform in the deconvolution operation process, therefore, the Fourier plane handle in CUDA C / C++ is applied in advance, denoted as IM fft ; the memory length is sizeof(ushort) x 1024 x 1024;

[0061] Step (3): The inverted fluorescence microscope is used as the collection hardware device of the application, the laser is used as the excitation light source, the CMOS camera is used for image collection, and the collected original image data is stored in the fixed shared memory and the Device memory of the graphics processor; the collected image original data is 16-bit single-channel image data, the image data type is ushort type,

[0062] Step (4): According to the numerical aperture of the objective lens, the physical size of the original image single pixel point, the original image side length, and the wavelength of the laser, the cutoff frequency radius of the original image is calculated; the calculation method of the cutoff frequency radius of the original image is as follows:

[0063] Let the cutoff frequency radius of the original image be p c , and the calculation formula is as follows:

[0064] p c = p max x P size x N

[0065] Where p max is the theoretical cutoff frequency of the microscope optical system, NA is the numerical aperture of the objective lens during shooting, l is the wavelength of the laser, P size is the physical size of the original image single pixel point, Where U Pixel is the physical size of the camera unit pixel point, M obj is the objective lens magnification, M eye is the magnification of the microscope eyepiece, and N is the original image side length.

[0066] A specific example is: the numerical aperture is NA = 1.49; the original image size is 1024 x 1024, that is, the image side length N = 1024; the physical size of the camera unit pixel point U Pixel = 6.5 pm, the objective lens magnification M obj = 100, the magnification of the microscope eyepiece M eye = 1.5 times, and the image pixel point physical size P size = 43.33 nm is calculated from the above formula; the cutoff frequency radius p c is about 224, which is calculated from this.

[0067] Step (5): According to the data size, the number of threads required in the operation process of the graphics processor is calculated, and the number of block threads required is designed;

[0068] Step (5-1): The number of threads in parallel is calculated:

[0069] Thread number = blockIdx.x*blockDim.x + threadIdx.x

[0070] Step (5-2): The number of threads and block threads used is obtained by querying the performance of the current graphics card. The maximum number of threads contained in each block thread is denoted as S. The original image is square, and the total number of image pixels N*N is used as the total number of threads required. Then the number of block threads is

[0071] According to the query, the number of GPUs used by the computer in this embodiment is 1, the number of blocks is (2147483647, 65535, 65535), and the number of threads in each block is 1024. The size of the BlockNum and ThreadNum of the kernel function is designed. In the embodiment, the original image data is a one-dimensional matrix of 1024*1024. Since the number of blocks on the machine is sufficient, one-dimensional can be used. The maximum number of threads of the machine is 1024, so the calculation method of BlockNum is:

[0072] int blockNum = (Len-0.5) / threadNum+1

[0073] Len is the length of the image data, which is equal to 1024*1024. The minus 0.5 and the plus 1 in the formula are to prevent integer division and consider the next bit.

[0074] Step (6): According to the cut-off frequency radius calculated in step (4), generate an imaging system optical transfer function template, and store the Fourier transform of the template image in the Device memory of the graphics processor;

[0075] The specific steps are as follows:

[0076] Step (6-1): According to the length of the original image data to be processed, divide a plurality of thread blocks in the Device memory grid of the graphics processor, start the kernel function Kernel function 1 to calculate, and input the image cut-off frequency radius p c into the real plane allocated by the graphics processor, and after processing and calculation, a binary optical transfer function plane template is obtained, denoted as Mask otf ; traverse the plane Maskotf , the distance between any pixel point and the center is less than the cutoff frequency radius p c , the pixel value of the pixel point is assigned as 1; the distance between any pixel point and the center is greater than the cutoff frequency radius p c , the pixel value of the pixel point is assigned as 0; the basic method and formula are as follows:

[0077]

[0078] In the above formula, denotes the optical transfer function plane template Mask otf , the horizontal coordinate value is x j , the vertical coordinate value is y j , j = 0, 1, 2, …, Len-1; Len is the length of the original image data;

[0079] Step (6-2): Then continue in the GPU, call Kernel function 2, the optical transfer function plane template Mask otf is subjected to Fourier transform in the graphics processor to obtain the optical transfer function frequency domain template, denoted as Mssk otf_cpl , and stored in the Device memory of the graphics processor.

[0080] Step (7): Obtain the original image data to be deconvolution processed from the fixed shared memory for Fourier transform, and perform deconvolution processing calculation on the original image data subjected to Fourier transform in the graphics processor, wherein the optical transfer function template subjected to Fourier transform obtained in step (6) is used in the deconvolution processing calculation, the number of Block thread blocks used in the CUDA operation is calculated in step (5), until the result image data is obtained, the result image data is converted from the frequency domain to the spatial domain, and finally the result image data is obtained. The result image of the real plane is stored in the Device memory of the graphics processor; the specific method is as follows:

[0081] Step (7-1): Use the cudaMemcpyAsync(..., cudaMemcpyHostToDevice) function to copy the original image data in the fixed shared memory to the Device memory of the graphics processor, denoted as Buf device ;

[0082] Step (7-2): Use the fast Fourier library in CUDA to perform Fourier transform on the original image data Buf device , to obtain the complex plane expression of the original image data in the frequency domain, denoted as Buf deuice_cpl ;

[0083] Step (7-3): The optical transfer function frequency domain template Maskotf_cpl As an optical transfer function template for deconvolution operations, the original image data Buf after Fourier transform is processed in the graphics processor. device_cpl Perform deconvolution operation:

[0084] First, based on the deconvolution principle, the image degradation model can be described as follows:

[0085]

[0086] Where f(x,y) represents the original image, h(x,y) represents the point spread function of the optical system, n(x,y) represents additive noise, and g(x,y) represents the degraded blurred image;

[0087] Then, the original image f(x, y) is iteratively calculated using the following formula.

[0088]

[0089] The resulting image data is obtained, where... This is the convolution operator, where T is the matrix transpose symbol and n represents the number of iterations.

[0090] Step (7-4): Transform the resulting image data from the frequency domain to the spatial domain to obtain the final real-plane result image, denoted as DM. real And store it in the graphics processor's Device memory.

[0091] Step (8): Transfer the resulting image DM from the graphics processor's device memory. real The image DM is copied from the Device memory to the fixed shared memory to obtain the final restored image DM. Specifically, the cudaMemcpy() function in CUDAC / C++ is used, expressed as follows: cudaMemcpy(DM, DM...) real ->data,DM real ->size, cudaMemcpyDeviceToHost); where cudaMemcpy is the data copy function, and cudaMemcpyDeviceToHost indicates the transfer type of copying memory from device to host.

[0092] The comparison chart of biological sample treatment effects used in this embodiment is shown below. Figure 3 and Figure 4 As shown. Figure 3 For the original image, Figure 4The post-image is processed. In this embodiment, the number of iterations of the deconvolution is 2, and the processing time is about 4 milliseconds. The results show that the method of the application has extremely fast calculation and processing speed, the texture of the processed image is restored, the deblurring effect is superior, the calculation efficiency of the deconvolution is effectively improved, sufficient processing time is reserved for subsequent image processing, and the video stream mode of real-time image acquisition-image deconvolution calculation-image display can be realized.

Claims

1. A CUDA-based deconvolution image restoration method, characterized by It comprises the following steps: Step (1): the host assigns a fixed shared memory as a temporary area for the original image data and the final result image; Step (2): a device memory is allocated in the graphics processor for saving the optical transfer function template, image data, complex plane required for Fourier transform and result image data; Step (3): the original image is collected by the image collection device of the optical microscope with laser as the light source, and the collected original image data is stored in the fixed shared memory and the device memory of the graphics processor; Step (4): the cutoff frequency radius of the original image is calculated according to the numerical aperture of the objective lens, the physical size of the single pixel point of the original image, the side length of the original image and the wavelength of the laser; Step (5): the number of threads required in the operation process of the graphics processor is calculated according to the data size, and the number of block thread blocks required is designed; Step (6): the imaging system optical transfer function template is generated according to the cutoff frequency radius calculated in step (4), and the Fourier transform of the imaging system optical transfer function template is stored in the device memory of the graphics processor; Step (7): the original image data to be deconvoluted is obtained from the fixed shared memory for Fourier transform, and the deconvoluted original image data is processed and calculated in the graphics processor, wherein the Fourier transformed optical transfer function template obtained in step (6) is used in the deconvolution calculation, the number of block thread blocks used in the CUDA operation is calculated in step (5), until the result image data is obtained, the result image data is converted from the frequency domain to the spatial domain, and the result image of the real plane is finally obtained and stored in the device memory of the graphics processor; Step (8): the result image in the device memory of the graphics processor is copied to the fixed shared memory to obtain the final restored image.

2. The CUDA-based deconvolution image restoration method of claim 1, wherein: The calculation method of the cutoff frequency radius of the original image in step (4) is as follows: Let the cutoff frequency radius of the original image be denoted by p c The formula for calculating it is as follows: p c = p max × P size × N wherein ρ max is the theoretical cut-off frequency of the microscope optical system, NA is the numerical aperture of the objective, λ is the wavelength, P size is the physical size of a single pixel of the original image, wherein U Pixel is the physical size of a single pixel of the camera unit, M obj is the magnification of the objective, M eye is the magnification of the microscope eyepiece, and N is the side length of the original image.

3. The CUDA-based deconvolution image restoration method of claim 2, wherein The calculation method of the number of block thread blocks designed in step (5) is as follows: Step (5-1): calculate the number of parallel threads: Thread number = blockIdx.x*blockDim.x + threadIdx.x Step (5-2): Obtain the number of Thread and Block by querying the performance of the current graphics card in use, record the maximum number of threads contained in each Block as S, set the original image as a square, and set the total number of image pixels N x N as the total number of threads required, then the number of Block is 4. The CUDA-based deconvolution image restoration method of claim 3, wherein The specific steps of step (6) are as follows: Step (6-1): According to the length of the original image data to be processed, a plurality of thread blocks are divided in the Device memory grid of the graphic processor, and the kernel function is started for calculation. The image cutoff frequency radius ρ obtained in step (4) is taken as the input of the kernel function, and the image cutoff frequency radius ρ is taken as the input of the kernel function. c The image cutoff frequency radius ρ is taken as the input of the kernel function, and the image cutoff frequency radius ρ is taken as the input of the kernel function. otf ; the plane MasK otf is traversed, and any pixel point is assigned a pixel value of 1 if the distance between the position of the pixel point and the center is less than the cutoff frequency radius ρ c , and a pixel value of 0 if the distance between the position of the pixel point and the center is greater than the cutoff frequency radius ρ c . The basic method and formula are as follows: In the above formula, represents the optical transfer function plane template Mask otf The horizontal coordinate value is x j , and the vertical coordinate value is y j The gray value of the pixel point, j=0, 1, 2, …, Len-1; Len is the length of the original image data; Step (6-2): The optical transfer function plane template Mask otf In the described graphic processor, Fourier transform is carried out to obtain an optical transfer function frequency domain template, denoted as Mask otf_cpl , which is stored in the Device memory of the described graphic processor.

5. The CUDA-based deconvolution image restoration method of claim 4, wherein The specific method of step (7) is as follows: Step (7-1): using cudaMemcpyAsync function to copy the raw image data in the fixed shared memory to the Device memory of the graphic processor, denoted as Buf device ; Step (7-2): Perform Fourier transform on the original image data Buf device using the fast Fourier library in CUDA to obtain the complex plane expression of the original image data in the frequency domain, denoted as Buf device_cpl ; Step (7-3): The optical transfer function frequency domain template Mask otf_cpl As an optical transfer function template for deconvolution operation, the original image data Buf device_cpl deconvolution operation processing: Firstly, according to the deconvolution principle, the image degradation model is described as: Wherein, f(x, y) represents the original image, h(x, y) represents the point spread function of the optical system, n(x, y) represents the additive noise, and g(x, y) represents the degraded blurred image; Then, the original image f(x, y) is iteratively operated using the following formula, , resultant image data is obtained, where is a convolution operator, T is a matrix transpose symbol, and n represents the number of iterations; Step (7-4): The resulting image data is converted from frequency domain to spatial domain, and the final real plane resulting image is obtained, denoted as DM real and stored to the Device memory of the graphic processor.

Images