A real-time image denoising method based on GPU

Abstract

A real-time image denoising method based on GPU belongs to the field of image processing. The input picture is read by using opencv, and the picture pixel value is normalized to [0, 255]. The image data is transmitted from the host end to the device end, and the image data is stored in the shared memory. Each thread reads the corresponding 9 points in the shared memory in turn and stores them in the 3x3 window window. The number of salt and pepper noise points in the 3x3 window is judged, and the corresponding denoising method is used for processing. The data processed by the device end is transmitted to the host end, and the opencv is used to read and display the image after denoising. The method solves the problems of long filtering time, poor removal effect of mixed salt and pepper noise and Gaussian noise in the current nonlinear filtering field, and the processed image has the characteristics of less detail information loss, high image clarity and obvious texture characteristics. The denoised image can be processed and displayed in real time.

Description

Technical Field

[0001] This invention relates to the field of image processing, and more particularly to a GPU-based real-time image denoising method. Background Technology

[0002] During image acquisition and transmission, noise is inevitably generated due to interference from random signals. Common noises include dark current noise and non-uniform photoresponse noise, which can be mainly classified into four types: Gaussian noise, salt-and-pepper noise, multiplicative noise, and Poisson noise. Among them, Gaussian noise and salt-and-pepper noise are the most common types and have a significant impact on the quality of the final acquired image.

[0003] Common noise reduction methods include linear filtering and nonlinear filtering. Linear filtering, such as mean filtering and block filtering, can effectively remove noise, but it severely damages the texture features of the image, resulting in a blurry image and poor image processing performance. Nonlinear filtering, including median filtering and bilateral filtering, can effectively remove noise without damaging image details, but its slow computation speed makes it difficult to meet real-time requirements. Therefore, we urgently need a noise reduction method that satisfies the requirements of good noise reduction effect, no significant loss of image details, and real-time performance.

[0004] Fortunately, the advent of GPUs has made real-time image processing possible through software programming. Based on NVIDIA's CUDA platform, parallel programming using CUDA is easy; it can be combined with the OpenCV library to directly read images in formats such as JPG, and can also directly display the output images. Because nonlinear filtering algorithms operate essentially the same on each pixel, they are very suitable for parallel implementation.

[0005] Therefore, we adopted an improved filtering algorithm based on the CUDA platform to process Gaussian noise and salt-and-pepper noise to meet real-time requirements. Summary of the Invention

[0006] The purpose of this invention is to provide a GPU-based image denoising method. This method is based on the CUDA platform and can achieve real-time image denoising. It solves the problems of long filtering time and poor removal effect when salt-and-pepper noise and Gaussian noise are mixed in the current field of nonlinear filtering. The processed image has the characteristics of less loss of detail information, high image clarity and obvious texture features. The denoised image can be processed and displayed in real time.

[0007] This invention provides a GPU-based real-time image denoising method, which mainly includes the following steps:

[0008] S1. Use OpenCV to read the input image and normalize the image pixel values ​​to [0, 255].

[0009] S2. Transmit the image data from the host to the device and store the image data in shared memory. Each thread reads the corresponding 9 points in the shared memory in turn and stores them in a 3x3 window. Determine the number of salt-and-pepper noise points in the 3x3 window and use the corresponding noise reduction method for processing.

[0010] S3. Transmit the data processed by the device to the host, and use OpenCV to read and display the noise-reduced image.

[0011] Furthermore, step S2 mainly includes the following steps:

[0012] S21. Allocate device storage space memory. Since global memory read speed is slow, requiring dozens of cycles to read a single piece of data, while shared memory only requires one to two cycles to read data, each thread sequentially reads the corresponding 9 pieces of data stored in shared memory and stores the data in a 3x3 window.

[0013] S22. Determine the number of salt-and-pepper noise points in the window. If the pixel value of a point is equal to 0 or 255, then this point is determined to be a salt-and-pepper noise point, and the device function related to median filtering is directly called to perform salt-and-pepper noise denoising. Otherwise, continue to determine until it is determined that there are no salt-and-pepper noise points in the 3x3 window, then the device function related to bilateral filtering is called to perform filtering operations.

[0014] S23, Improved bilateral filtering algorithm, such as Figure 2 As shown, research revealed that the structural similarity coefficient (ssim) between the ideal image and the image after bilateral filtering is close to 97%. This indicates that the noise-free image experiences almost no loss of detail after bilateral filtering, and bilateral filtering itself also has a good effect on processing Gaussian noise. Therefore, bilateral filtering can be directly used to process noise-free or Gaussian noise-containing regions.

[0015] S24. Improved median filtering algorithm, such as Figure 3 As shown, median filtering is performed by repeatedly calling the device function `sort` (a function that can only be called on the device side) to sort the three values ​​using a relevant algorithm.

[0016] Further improved bilateral filtering algorithms based on S23 include the following steps:

[0017] S231. Calculate the average value of pixels within a 3x3 window.

[0018] S232. Based on the mean calculated in the previous step, select the corresponding spatial domain standard deviation (space_coeff) and value domain standard deviation (value_coeff).

[0019] S233. Based on the value of space_coeff, call the 3x3 space field template space stored in constant memory; based on the value of value_coeff, determine the array value_lut containing 256 data items stored in constant memory.

[0020] S234. Calculate the absolute value of the difference between each pixel in the window and the center pixel of the 3x3 window, window[1][1].

[0021] S235. The absolute value of the difference calculated in the previous step is used as the address of the value_lut array (lookup table), and the value corresponding to the value_lut array is output and stored in the 3x3 value range template value.

[0022] S236. Multiply the corresponding positions of the 3x3 space and value and perform normalization to obtain the final 3x3 convolution template for bilateral filtering.

[0023] S237. Multiply the 3x3 convolution template template with the corresponding element in the template window containing the data pixel values, and add the products together to obtain the pixel values ​​after bilateral filtering.

[0024] Further improved median filtering algorithms based on S24 mainly include the following steps:

[0025] S241. Find the sorting results of a11, a12, a13 from largest to smallest: max1, mid1, min1; the sorting results of a21, a22, a23 from largest to smallest: max2, mid2, min2; and the sorting results of a31, a32, a33 from largest to smallest: max3, mid3, min3.

[0026] S242. Find the minimum value min_max among max1, max2, and max3; find the middle value mid_mid among mid1, mid2, and mid3; and find the maximum value max_min among min1, min2, and min3.

[0027] S243. Find the median value mid_data among min_max, mid_mid, and max_min. mid_data is the median value of the nine numbers from a11 to a33, which is the output of the median filter.

[0028] The beneficial effects of this invention are:

[0029] 1. Good image processing effect. The above-mentioned noise reduction method can achieve good noise reduction effect for images containing salt-and-pepper noise and Gaussian noise. It can remove the relevant noise relatively well, and the texture information of the processed image is more complete, and the image blur is greatly reduced.

[0030] 2. It can directly read or display images. It can directly read images in many formats such as JPG, PNG, and BMP using the OpenCV library, and can directly display the output images. This allows for a clearer observation of the image processing effects, eliminating the hassle of saving the output data to a TXT file and then reading the TXT file through software like MATLAB to observe the images.

[0031] 3. Images can be directly compressed or saved. OpenCV can be used to directly compress output images and specify the compression ratio or save them as images in many formats such as JPG and PNG.

[0032] 4. Real-time processing is possible. Processing speed is significantly improved, achieving more than 50 times the acceleration compared to CPU processing alone. For 512x512 images, processing time using NVIDIA's Jetson Nano development board is less than 10ms. Attached Figure Description

[0033] Figure 1 This is a schematic diagram of the overall process of a GPU-based real-time image denoising method provided in an embodiment of the present invention.

[0034] Figure 2 This is a flowchart illustrating an improved fast bilateral filtering algorithm based on shared memory provided in an embodiment of the present invention.

[0035] Figure 3 This is a flowchart illustrating an improved fast median filtering algorithm based on shared memory provided in an embodiment of the present invention.

[0036] Figure 4 This is a diagram showing the spatial domain and value domain relationship of each point within a 3x3 window provided in an embodiment of the present invention. Detailed Implementation

[0037] Reference Figure 1 This invention provides a GPU-based real-time image denoising method, which includes the following steps:

[0038] S1. Use OpenCV to read the input image and normalize the pixel values ​​to [0, 255].

[0039] Step S1 can be detailed as follows:

[0040] S11. Configure the OpenCV library on the CUDA platform. After configuration, you can use the `imread` function in OpenCV to read images in formats such as PNG, JPG, and BMP. If you read the image using the `IMREAD_DEPTH` method, you will directly read the actual pixel values ​​in the image without any processing. If you read the image using the default method, if the image is larger than 8 bits, the `imread` function will directly normalize the pixel values ​​to [0, 255]. Therefore, we will directly read the image using the default method of the `imread` function.

[0041] S12. Since the 8-bit images read by OpenCV are generally of type uchar, type conversion is required. The data can be converted to int or short type. The convertTo function in the OpenCV library can be used directly for type conversion. For example, image1.convertTo(image2, CV_32F) can directly convert the data type of the read image image1 to a 32-bit single-precision floating-point number and store the converted data in image2.

[0042] S2. Transmit the image data from the host to the device, determine the number of salt-and-pepper noise points in the 3x3 window, and process it using the appropriate noise reduction method.

[0043] Further step S2 mainly includes the following:

[0044] S21. Allocate host and device storage space memory. Since the image resolution may be very high, static memory allocation might cause the program to fail. Dynamic memory allocation avoids memory-related issues; therefore, host-side memory is dynamically allocated. Because data access speed in shared memory is tens of times faster than in global memory, image data is stored in shared memory, and each thread reads the corresponding 9 points from shared memory and stores them in a 3x3 window. Assuming the image height is `height` and width is `width`, the size of the allocated memory needs to be considered when allocating device-side memory. If the image data of type `int` is stored in a one-dimensional array `data`, the expression for dynamic host-side memory allocation is as follows:

[0045] int*data=(int*)malloc(height*width*sizeof(int));

[0046] Among them, malloc is the memory allocation function;

[0047] S22. Determine the number of salt-and-pepper noise points in the 3x3 window. Points with a pixel value of 0 or 255 are considered salt-and-pepper noise points. Record the count of salt-and-pepper noise points in the window. If count = 1, stop counting the salt-and-pepper noise points and immediately call the device function fast_median_filter. If it is determined that none of the 9 points are salt-and-pepper noise points, i.e., count = 0, then immediately call the device function fast_bilateral_filter.

[0048] S23. Describe the algorithm of the device function fast_bilateral_filter as follows: Figure 2 As shown, fast bilateral filtering is performed according to the algorithm described.

[0049] S24. Describe the algorithm of the device function fast_meidan_filter as follows: Figure 3 As shown, according to the algorithm, fast median filtering is performed by calling the device function sort, which sorts the three values ​​by size.

[0050] Further step S21 can be detailed as follows:

[0051] S211. Data is transferred from the host to the device using the cudaMemcpy function.

[0052] S212. Allocate the number of threads per block and the number of threads per grid. Since image data is generally stored in a two-dimensional matrix, and hardware resources can be fully utilized when the number of threads per block is a multiple of 16, we set the block and grid as follows:

[0053] dim3ThreadsPerBlock(16,16); (Each block is 16*16 in size)

[0054] dim3BlocksPerGrid((width+ThreadsPerBlock.x-1) / ThreadsPerBlock.x,(height+

[0055] ThreadsPerBlock.y-1) / ThreadsPerBlock.y);

[0056] Where height is the image height and width is the image width;

[0057] ThreadsPerBlock represents the number of threads contained in a predefined block. These are defined in two dimensions, including both the x and y directions.

[0058] ThreadsPerBlock.x represents the number of threads contained in each block (thread block) in the x-dimensional dimension;

[0059] ThreadsPerBlock.y represents the number of threads contained in each block (thread block) along the y-axis;

[0060] (width+ThreadsPerBlock.x-1) / ThreadsPerBlock.x means calculating the number of blocks (thread blocks) contained in the grid (thread grid) corresponding to the x-dimensional dimension of the image (the width of the image) based on the preset value of ThreadsPerBlock.x.

[0061] (height+ThreadsPerBlock.y-1) / ThreadsPerBlock.y means calculating the number of blocks (thread blocks) contained in the grid (thread grid) corresponding to the y-axis of the image (the height of the image) based on the preset value of ThreadsPerBlock.y.

[0062] S213. The pixel coordinates corresponding to the thread index in the actual image can be expressed by the following formula (the coordinates of the first pixel are (0,0):

[0063] row=threadIdx.y+blockIdx.y*blockDim.y;

[0064] col=threadIdx.x+blockIdx.x*blockDim.x;

[0065] Wherein, `blockDim.y` represents the number of threads contained in each block (thread block) along the y-axis. `blockIdx.y` represents the current index position of the block (thread block) along the y-axis. `threadIdx.y` represents the current index position of the thread (thread) along the y-axis within the block. `blockDim.x` represents the number of threads contained in each block (thread block) along the x-axis. `blockIdx.x` represents the current index position of the block (thread block) along the x-axis. `threadIdx.x` represents the current index position of the thread (thread) along the x-axis within the block.

[0066] S214. Quickly traverse the pixels in the image from left to right and top to bottom, and store the pixels at the corresponding image positions (row-1,col-1), (row-1,col), (row-1,col+1), (row,col-1), (row,col), (row,col+1), (row+1,col-1), (row+1,col+1), (row+1,col+1) in shared memory into a 3x3 window. Since pixels at the image boundaries have negligible impact on the overall image, they are directly assigned values ​​without further processing. Because shared memory is an on-chip cache shared by all threads within the same thread block, it is much faster than global memory (accessible to all threads in the thread grid) and local memory (unique to each thread). Therefore, we perform operations in shared memory.

[0067] Further step S22 can be detailed as follows:

[0068] Because filtering methods such as Gaussian filtering and mean filtering significantly damage the texture information of images (regions with drastic pixel value changes), they can easily cause negative effects such as image blurring. Bilateral filtering, however, considers both the value domain similarity and spatial domain similarity between pixel values ​​within a window, effectively preserving the texture information of the image during the filtering process. The spatial domain expression for a typical bilateral filtering 3x3 window is shown below:

[0069] d(i,j)=exp(-((i-0)^2+(j-0)^2) / (2*σ d ^2))(-1<=i<=1,-1<=j<=1,σ d (Smoothing coefficient in the spatial domain)

[0070] The value range expression for its 3x3 window is shown in the following formula:

[0071] r(i,j)=exp(-((I(i,j)–I(1,1))^2) / (2*σ r ^2))(0<=i<=2,0<=j<=2,σ r I(i,j) is the range smoothing coefficient, representing the pixel value at (i,j).

[0072] The positions of each point in d(i,j) are as follows Figure 4 As shown on the left, observing the expression d(i,j) reveals that for a specific position in the 3x3 template, once σ d Once the size of σ is determined, d(i,j) is also determined. Therefore, we select several typical σ... d Given the value, we can directly calculate the 3x3 spatial domain template corresponding to σd. To speed up the calculation, we can enlarge the spatial domain template by a certain factor to make the processed result closer to an integer or simply round the calculated result to an integer.

[0073] The positions of each point in I(i,j) are as follows: Figure 4 As shown on the right, observing the expression for r(i,j), we can see that the range of values ​​for abs(I(i,j)–I(1,1)) is [0,255]. Once σ r Once the size is determined, we can completely predetermine σ. r The corresponding r(i,j) value is calculated, and the calculated value is multiplied by 2^10 and stored in an array of length 256. This array is used as a lookup table. The value of abs(I(i,j)–I(1,1)) is input into this array as the address, and the output value is the corresponding r(i,j) value at this time.

[0074] Since constant memory is on-chip memory and read-only, its access speed is comparable to shared memory, making it ideal for storing lookup table data. Therefore, we can transfer the data of the storage space field template and value field template to the relevant array defined on constant memory by calling cudaMemcpyToSymbol.

[0075] Let window be a 3x3 window that reads data in shared memory. If window[i][j] == 0 or

[0076] If window[i][j] == 255 (0 <= i <= 2, 0 <= j <= 2), then record count = 1, and immediately stop judging the salt-and-pepper noise points, directly enter the next step, and call the relevant device functions. Otherwise, keep counting the salt-and-pepper noise points in window until it is judged that there are no salt-and-pepper noise points in the window, record count = 0, enter the next step, and call the relevant device functions.

[0077] The improved bilateral filtering algorithm in step S23 includes the following specific steps:

[0078] S231. Calculate the average value of the pixels in the 3x3 window window in the shared memory.

[0079] S232. According to the mean value calculated in the previous step, select the corresponding spatial domain coefficient space_coeff and range coefficient value_coeff.

[0080] S233. According to the value of space_coeff, call the corresponding 3x3-sized spatial domain template space stored in the constant memory; according to the value of value_coeff, call the corresponding array value_lut storing 256 data stored in the constant memory.

[0081] S234. Calculate the absolute value abs_value of the difference between each pixel point window[i][j] (0 <= i <= 2, 0 <= j <= 2) in window and the central pixel point window[1][1] of window.

[0082] S235. Use the abs_value calculated in the previous step as the address of the value_lut array, and output the corresponding value and store it in the 3x3-sized range template value.

[0083] S236. Multiply the corresponding positions of the 3x3-sized space and value and perform a normalization operation to obtain the final 3x3-sized convolution template template of the bilateral filtering. The calculation formula is as follows.

[0084] template_value = sum(space[i][j] * value[i][j]) (0 <= i <= 2, 0 <= j <= 2, sum represents summation)

[0085] template[i][j] = space[i][j] * value[i][j] / template_value (0 <= i <= 2, 0 <= j <= 2)

[0086] S237. Multiply the corresponding elements of the 3x3 convolution template "template" and the template "window" containing element pixel values, and add the products to obtain the pixel value "bilateral_value" after bilateral filtering.

[0087] bilateral_data = sum(template[i][j] * window[i][j]) (0 <= i <= 2, 0 <= j <= 2, sum represents summation, row and col respectively represent the image position indices corresponding to the central element of window under the current thread)

[0088] The improved median filtering algorithm based on S24 mainly includes the following content and steps:

[0089] For the common algorithm of taking the median value of 9 elements, it is necessary to compare the sizes of 9 elements, which requires a relatively large number of comparisons and takes a long time. Commonly used ones include quicksort, bubble sort, etc. In the worst case, it needs to compare 36 times. However, when observing median filtering, we can find that it is not necessary to sort the entire sequence, but only to find the median value of 9 elements. If we adopt the Figure 3 shown fast median filtering method, only 21 comparisons are needed to find the median value, which greatly reduces the sorting time.

[0090] First, we write a device function "sort" to sort the three numbers a1, a2, and a3, and output the sorting result to an array. Next, perform operations according to the following steps:

[0091] S241. By calling the device function "sort", obtain the sorting results max1, mid1, min1 of a11, a12, a13 from largest to smallest, the sorting results max2, mid2, min2 of a21, a22, a23 from largest to smallest, and the sorting results max3, mid3, min3 of a31, a32, a33 from largest to smallest.

[0092] S242. By calling the device function "sort", obtain the minimum value min_max among max1, max2, max3, obtain the median value mid_mid among mid1, mid2, mid3, and obtain the maximum value max_min among min1, min2, min3.

[0093] S243. Obtain the median value mid_data among min_max, mid_mid, max_min. mid_data is the median value of the 9 numbers from a11 to a33, that is, the output result of median filtering.

[0094] After obtaining the filtering result, we store it in the corresponding position of the result array. Once the noise reduction process for the entire image is complete, we output the result array.

[0095] S3. Transfer the data processed on the device to the host, and use OpenCV to read and display the denoised image. Use cudaMemcpy to transfer the data in the array obtained after denoising on the device to the host. Note that since OpenCV typically displays 8-bit images, the data type usually needs to be uchar, so we need to perform type conversion on the output data.

Claims

1. A GPU-based real-time image denoising method, characterized in that, The method mainly includes the following steps: S1. Use OpenCV to read the input image and normalize the image pixel values ​​to [0, 255]. S2. Transmit the image data from the host to the device, store the image data in shared memory, each thread reads the corresponding 9 pixels in the shared memory in turn and stores them in a 3x3 window, determine the number of salt and pepper noise points in the 3x3 window, and process them using the corresponding noise reduction method. S3. Transmit the data processed by the device to the host, and use OpenCV to read and display the noise-reduced image; In step S2: the number of salt-and-pepper noise points in the window is determined. If the pixel value of a point is equal to 0 or 255, then the point is determined to be a salt-and-pepper noise point, and the device function related to median filtering is directly called to denoise the salt-and-pepper noise; otherwise, the determination continues until it is determined that there are no salt-and-pepper noise points in the 3x3 window, and then the device function related to bilateral filtering is called to perform the filtering operation. Further improvements to the bilateral filtering algorithm include the following steps: S231. Calculate the average value of pixels within a 3x3 window. S232. Based on the average value calculated in the previous step, select the corresponding spatial domain standard deviation (space_coeff) and value domain standard deviation (value_coeff). S233. Based on the value of space_coeff, call the 3x3 space field template space stored in constant memory; based on the value of value_coeff, determine the array value_lut containing 256 data items stored in constant memory; S234. Calculate the absolute value of the difference between each pixel in the window and the center pixel of the 3x3 window, window[1][1]. S235. The absolute value of the difference calculated in the previous step is used as the address of the value_lut array. The value corresponding to the value_lut array is output and stored in the 3x3 value range template value. S236. Multiply the corresponding positions of the space and value in the 3x3 size and perform normalization to obtain the final 3x3 size convolution template for bilateral filtering. S237. Multiply the 3x3 convolution template template with the corresponding element in the template window containing the data pixel values, and add the products together to obtain the pixel values ​​after bilateral filtering.

2. The GPU-based real-time image denoising method according to claim 1, characterized in that, The improved median filtering algorithm includes the following steps: S241. Find the sorting results of a11, a12, a13 from largest to smallest: max1, mid1, min1; the sorting results of a21, a22, a23 from largest to smallest: max2, mid2, min2; and the sorting results of a31, a32, a33 from largest to smallest: max3, mid3, min3. S242. Find the minimum value min_max among max1, max2, and max3; find the middle value mid_mid among mid1, mid2, and mid3; and find the maximum value max_min among min1, min2, and min3. S243. Find the median value mid_data among min_max, mid_mid, and max_min. mid_data is the median value of the nine numbers from a11 to a33, which is the output of the median filter.

Images