How to Denoise images using Curvelet Transform in Matlab

Curvelet transform is a powerful tool that can capture details along the curvature in images. It is useful when it comes to feature extraction and pattern recognition. <!--more--> Curvelet transform is also efficient in image denoising. This is the removal of noise signals in an image.

This tutorial will look into the Curvelet transform analysis and denoising of images using Matlab. We will also discuss the application of the Curvelet toolbox.

Table of contents #

• Prerequisites
• Theoretical background of Curvelet transform
• Curvelet Toolbox
• Curvelet transform of an image
• Curvelet based denoising of an image
• Matlab code for denoising
• Conclusion

Prerequisite #

To follow along, you need to have:

• Matlab installed in your computer.
• A clear understanding of Matlab basics.

Theoretical background of Curvelet transform #

Curvelet transform is a geometric scale transform, used to represent edges and curves more efficiently than traditional wavelets. A huge disadvantage of wavelets is that they are not directional.

Wavelets are effective in capturing coefficients along the horizontal, vertical, and diagonal axes but poor in determining coefficients along the curvature.

The complex transform is one way to improve directional selectivity. However, it is difficult to design complex wavelets with perfect reconstruction properties and filter characteristics.

Ridgelet transform was proposed by Candes and Donoho in 1999. It was regarded as an anisotropic geometric wavelet transform.

In 2000, they introduced 1st and 2nd generation Curvelet transform, called Ridgelet transform-based Curvelet transform.

The 2nd generation Curvelet transform is an efficient tool in numerous applications including image processing, seismic data exploration, fluid mechanics, and solving partial differential equations. It is also efficient in representing curve-like edges.

However, this generation suffered two main drawbacks:

• It is not optimal for sparse approximation of curve features beyond C-square singularity.
• It is highly redundant, thus giving a poor output.

Later, a fast and less redundant Curvelet transform version was introduced. This newly implemented version, also known as Fast discrete curvelet transform (FDCT), is implemented in two ways:

1. Unequally spaced Fast Fourier transform(USFFT).

2. Wrapping function.

Wrapping-based FDCT is faster than USFFT, which makes it widely used.

In general, Discrete Curvelet transform can be expressed as follows:

Where: C: is the Curvelet coefficients which are the function of j, $\theta, k_1$, and $k_2$. j: is for scale. These are the layers of decomposition. For example, if j is 3, we have three levels of decomposition of the image using a Curvelet. $\theta$: This is the Curvelet orientation. It is because the Curvelet is arranged in the spatial domain. $k_1$ and $k_2$: They are the spatial location of the Curvelet. $\varphi[x,y]$: Is the Curvelet function. f[x,y]: This is the input image with a size of MxN.

This curvelet transform is usually expressed in the frequency domain, as shown below:

curvelet transform = IFFT{FFT(curvelet) x FFT(image)}

Curvelet toolbox #

Mathworks Matlab has no function to implement Curvelet transform. However, they can be computed using a third-party toolbox created by Curvelab.

Fortunately, the software is free. All you need is to create an account, then download it.

The Curvelab has two functions:

c=fdct_wrapping(x, is_real, finest, nbscales, nbangles_coarse)  % for FDCT
x=ifdct_wrapping(c, is_real, M, N)

In the function above; x: is the input image which is an MxN matrix.

is_real: Is the type of transform. Here, 0 is the default and is used for complex-valued curvelet. 1 is for a real-valued curvelet.

Finest: This can be Curvelet or wavelet. The value 1 is for curvelet, and 2 is for the wavelet.

nbscales: This shows the number of decomposition levels. If this value is not defined, the default value is given by:

ceil[log2(min(M, N))-3]

nbangles_coarse: This is the orientation of the angles. This value must be a multiple of 4. The minimum value is 8, and the default is 16.

C: This is the output cell array of the Curvelet coefficients. It is in the form of C{j}{k1, k2}.

Curvelet transform of an image #

We first need to select our image from the stored folder and then read the image. In this tutorial, we use the images in the CurveLab (Lena).

[filename, pathname] = uigetfile('*.*', 'select your input grayscale image');
filewithpath = strcat(pathname, filename);
img = imread(filewithpath);

The function uigetfile allows one to select an image from different folders. It takes in the filename and path as arguments.

Strcat converts this path to a string, and then the image is read using the imread function. All pixel values are converted to double, as shown below:

imgd = double(img);

Now, we are taking the FDCT of the input image, and extracting the approximation coefficients.

Since these coefficients are so large, we normalize them for display using the code below:

C = fdct_wrapping(imgd, 1, 2, 3); %Taking FDCT of input images.
cA = C{1, 1}{1, 1}; %Extracting approximation coefficients.
cAd = 225.*(cA./max(cA));  %Normalizing coefficients for display only.
figure(1)
imshow(uint8(cAd), []);   % Showing approximation coefficients.

To get the FDCT of the input image, we use the curve lab function fdct_wrapping. This function takes the image pixel values that are stored in the imgd variable, real-valued curvelets (1), wavelet at finest level (2), and the decomposition level (3).

We then extract approximate coefficients using C{1,1}{1,1}. To normalize the coefficients, we use the formula:

225.*(cA./max(cA))
Where;
cA: Are the extracted coefficients.

Then the display uses the imshow function. For instance, let's show all the detailed coefficients at scale 2 using a for loop:

%Showing all the detailed coefficients at scale 2.

figure(2)
for k=1:16
    x=C{1, 2}{1, k};
    xr=imresize(x, [512, 512]);
    imshow(uint8(xr), []);
    pause(1);
end

k is the number of edges, ranging from 1-16. c Is the FDCT of the input image. The detailed coefficient at the second decomposition level is extracted using the varying k. It is done by C{1, 2}{1, k}.

These extracted coefficients are resized to a square dimension of 512 using the imresize function.

The square coefficients are then displayed using the imshow function one by one:

The white colors are the Curvelet coefficients of scale 2 at different orientations. They kept changing since the number of edges was 16. This is how Curvelets are efficient in capturing coefficients along the curvature.

The image below shows the approximated coefficients:

Curvelet based denoising of an image #

Below is the curvelet-based denoising scheme:

This scheme takes the noisy images as the input, performs FDCT to obtain the approximation coefficients cA and detailed coefficients cD.

It then analyzes the threshold for the detailed coefficients. The denoised image is obtained when IFDCT is performed on these thresholded coefficients.

Matlab code for denoising #

We first retrieve the input image from our folder using the uigetfile function and read it using the imread function.

We then determine the size of this image using the size function, as shown below:

[filename, pathname] = uigetfile('*.*', 'select your input grayscale image');
filewithpath = strcat(pathname, filename);
img = imread(filewithpath);   %reading the image.
imgd = double(img); %converting the image matrix to double datatype.

n = size(img, 1);   %getting the size of the image.

Define the noise variance (sigma) and then add noise to the image using the code below:

sigma = 20;

noisy_img = imgd + sigma*randn(n);   %Adding noise

sigma defines the noise variance. You get noise when you multiply this variance by a random matrix of equal distribution using the function randn.

Let'calculate the norm of the curvelet. The norm gives a measure of the element's magnitude:

%Compute norm of curvelets(exact)
F = ones(n);
X = fftshift(ifft2(F)) * sqrt(numel(F));
Cn = fdct_wrapping(X, 0, 2);
E = cell(size(Cn));

for s = 1:length(Cn)
    E{s} = cell(size(Cn{s}));

    for w = 1:length(Cn{s})
        A = Cn{s}{w};
        E{s}{w} = sqrt(sum(sum(A.*conj(A))) / numel(A));
    end
end

The code above is for obtaining the norm of the Curvelet. The code remains unchanged in all programs. It is the general way of obtaining the Curvelet's norm.

Let's now find the Curvelet transform of the noisy image using the fdct_wrapping function:

%Taking curvelet transform
C = fdct_wrapping(noisy_img, 1, 2);

In this case, the decomposition level is not given. Thus, the default one is used. To apply the threshold, we use the code below:

%Applying thresholding
Ct = C;

for s = 2:length(C)
    thresh = 3*sigma + sigma*(s == length(C));

    for w = 1: length(C{s})
        Ct{s}{w} = C{s}{w}.*(abs(C{s}{w}) > thresh*E{s}{w}); %hard thresholding
    end
end

All the curvelet coefficients are stored in the variable Ct. We use a for loop to read all the thresholded coefficients.

The loop starts from two (2: length(c)). This is because the approximated coefficients are not thresholded. We obtain the threshold using the following formula:

3*sigma + sigma*(s == length(C))   %hard thresholding

We also perform IFDCT to reconstruct this image using the threshold coefficients:

% Taking inverse curvelet transform
restored_img = real(ifdct_wrapping(Ct, 1));

To get the IFDCT of the threshold coefficients, we use the ifdct_wrapping function. After this, we visualize the output using the imshow function:

%plotting results
subplot(131); imshow(uint8(img), []); title('Original image');
subplot(132); imshow(uint8(noisy_img), []); title('Noisy image');
subplot(133); imshow(uint8(restored_img), []); title('Denoised image');

When performing a denoising operation, it is good to get the Signal to noise ratio(SNR) to see how well your program is performing. This SNR can be obtained using the code below:

%Finding SNR
orig_vs_Noisy = 20*log10(norm(imgd(:))/norm(imgd(:)-noisy_img(:)));
orig_vs_Denoised = 20*log10(norm(imgd(:))/norm(imgd(:)-restored_img(:)));

Orig_vs_Noisy stores the SNR of the original image to the noisy image while orig_vs_Denoised stores the SNR of the original to the denoised image.

imgd is the matrix forming the image but in the double form.

When we run our program, we get:

We can see the SNR if the analysis is well performed. To get the SNR of the original image to the noise, we execute orig_vs_Noisy in the command window.

To get that of the original image to the denoised, we execute orig_vs_Denoised:

orig_vs_Noisy

orig_vs_Noisy =

   16.2380
orig_vs_Denoised =

   19.5566

As we can see, there is an improvement in the SNR from 16 to 19. Our program can now be described as effective.

Conclusion #

Curvelet transform is efficient for denoising and transforming images.

In the Curvelab, we can access Matlab functions and C++ code. This makes it an efficient tool for performing several operations.

Happy coding!

---

Peer Review Contributions by: Monica Masae