# Help required with BM3D Image denoising algorithm implementation

This is what I am doing and will be grateful for any help.

1. I selected 16 similar blocks to a reference block from a search block, put them in a 3D array, Wavelet Thresholded them, took inverse transform.

2. These 16 similar blocks belong to one 32 * 32 search block as mentioned above.

3. Now I should multiply each of these blocks with their weight, and then add them all up, creating a 2D array from a 3D array.

4. This new block then replaces all the 16 original blocks that I had selected before for hard thresholding.

Am I correct here? I will be absolutely grateful for a response.

Step 3 is not correct. The first part is ok (each denoised result is weighted by some likelihood of correctness, that was obtained previously from the variance of the data), but you don't have to add them up. Each slice from the 3D array is a 2D patch that is back-projected to the resulting image.

The sum happens when 2 back-projected patches do overlap in the final image domain.

Back-projection is the process of taking a 2D slice from the denoised stack and put it back in image space. This is required because unlike most variants of NL-means BM3D denoises patches and not the central pixel of a patch.

So, you have 2 images when implementing BM3D:

• the original (noisy) one;
• the result (denoised) image.

The spatial location of a patch should not change during the denoising process, otherwise the image content would be destroyed. So, each slice from the denoised patch stack goes back to its original place.

Now, one issue remains: what to do when back-projected patches overlap? What value should we keep for a single pixel? The solution chosen by BM3D is that each pixel contains a weighted mean of the patches that overlap in each pix.

To implement this, you can simply pre-multiply each slice of the denoised stack by its weight, add a third image to store the sum of the weights inside each pixel, sum for each pixel its value from all the denoised patches that contain it, and in a final step (once all the patch stacks haves been processed) divide the result image by the weight image.

Image processing main loop (outline)

There is no "unleft" pixel in the final result. The outline of the main loop is as follows:

1. for every pixel of the image, form a reference patch;
2. for a given reference patch, find n <= N (B=16 in your case) patches whose distance to the reference patch is less than some threshold and make a stack;
3. denoise the stack;
4. back project the patches in the denoised stack;
5. back to 1) for the next pixel).

Thus, you always have one stack for a given patch (degenerate case: only 1 patch in the stack), and a pixel is always involved in several patches.

• Sansuiso, is it possible for you explain it in a bit more detail. A code like way perhaps? Will be terribly grateful. Jan 8, 2014 at 11:41
• Code for BM3D is available for Alessandro Foi's website. If you can be more specific about what is not clear, I'll be happy to reformulate my answer :-) Jan 8, 2014 at 14:58
• thank you. Let me try to make it clearer. When a 2D patch is back-projected to the original patch, what do we do? Or to put in another way, what exactly is final image domain where we add 2 back-projected patches? Jan 8, 2014 at 21:15
• I've added a section on back-projection. Jan 9, 2014 at 7:24
• I am terribly grateful for this help. One final question though for which again the paper I find confusing: I get what we do with pixels that are in 3D array. But what about those that are not in 3D array? See, we take in my case 16 patches to 3D and hard threshold them. Now of course, my search path had many many more patches that were not taken into 3D array and hard thresholded. What do we do about them? From the aggregation statement, I understand that a zero goes to them. Am I correct here? Jan 13, 2014 at 7:06