# How one apply correctly FFT in image denoising

I'm writing program (Qt widgets/c++) for removing noise from images. As denoising method, i selected non local means method. This method has incredible quality of restored images (that's why it's the only denoising method in OpenCV), but has huge computation cost, so i made a lot of modified variants of this method (some with multithreading, some algorithmic). But, i'm having problem with the one, involving FFT

I followed all the steps of this article (only one page, 1430) and all works perfectly, except for FFT part, there just 2 lines about it in the paper and i can't understand, HOW should one use fft

This problem has bothered me for months, any help or insight would be greatly appriciated.

Shortened version of question: How can i get summed squared difference of two arrays on the image (the one at top and one in the middle, values are colors) quickly? ( O(n^2) is huge cost, there are lots of this kind of operations, paper above states, that it can be done via FFT with O(n*log n) (says that this 2 arrays forming circular convolution somehow))

• What did you finally end up doing for computing FFT ? Even if FFT is precomputed, the point-wise multiplication and addition of all patch elements takes $O(|P|)$ time where $|P|$ is the size of the patch. How did you overcome this ? Oct 28 '14 at 4:49

The trick inside the paper is the following:

1. What you want to compute is $\sum_{i \in W} |I(x+i)-I(y+i)|^2$, where $I$ is an image, $x$ and $y$ two noisy pixels and $i$ is a 2D offset used to define a patch.
2. Expanding the expression yields: $\sum_i I^2(x+i) + \sum_i I^2(y+i) - 2 \sum_i I(x+i)I(y+i) = A + B - 2C$.
3. $A$ and $B$ are computed using a squared integral image, i.e., an integral image from the squared original image.
4. $C$ is the convolution between the two patches centered on $x$ and $y$. Thus, it can be computed in the Fourier domain, where it becomes a multiplication. You get the value of $C$ by computing the Fourier transform of the patch around $x$, the patch around $y$, pointwise-multiplying these results and taking the inverse Fourier transform of the multplication result.

The Fourier transform is obviously a 2D transform since you are working with 2D data. What you obtain for a given patch is a 2D array of complex values.

In my opinion this article is not the best NL-means speedup strategy. Experiments I did way back in 2007/2008 show that pre-selection of patches are better (both in terms of speed and quality of the results). I have started blogging about these here, but unfortunately I am looking for time to finish the posts.

The original NL-means papers mention blockwise implementations that can be interesting. There are fundamentally 2 ways in implementing NL-means:

1. writing a denoising loop for every pixel in the image
2. writing a denoising loop for each patch, then back-project the patches to form an image.

The first impolementation is the original approach, because in 2005 memory and multicore CPUs were expensive. I chose on the other hand number 2 on recent hardware in the past 2 years. It depends on your typical image size and if you want to be able to compute domain transforms like DFT/DCT (as in the proposed paper and in BM3D).

• Thanks a lot for your answer, that's exactly what i needed, all was ready and working a long time ago, except for 4th item in that list, but now it is a lot clearer. Though one more question, if you don't mind: what will return Fourier transform of patch x or y? Array, vector, or single value? And what is required to use inverse transform? Becouse i am thinking about precomputing fft for every pixel (patches centered around it) and write results to 2d array before denoising and then just use those matrix to get inverse fft, but i don't know if this will be enough for inverse fft
– Shf
May 27 '13 at 20:54
• oh, and should i use 2d fft, or translate patch to 1d array? by the way, i was planning to write after this patchwise implementation anyway, thanks for an advise :) something similar to this also a long time ago- ipol.im/pub/art/2011/bcm_nlm
– Shf
May 27 '13 at 20:58
• I've updated the answer. May 28 '13 at 7:36
• ok, so i can precompute FFT for patches, centered around every pixel before donoising though it will take a lot of memory (mnsize_of_patchsize_of_patchsizeof(double)), but when i will count weights, i would still need to pointwise multiply 2 complex arrays and after that do inverse fft on received 2d array, it's even more then O(n^2) if i'm not mistaken
– Shf
May 28 '13 at 8:38
• Good answer, but how are you deriving that $C$ is a convolution? The way it is written it is a point wise element-by-element product. Where is the convolution? Oct 20 '13 at 3:00