BM3D is considered state of the art in Image Denoising. I am trying to implement the algorithm in Matlab. Right now, I am done with the 3D wavelet transform in the first step, and to be precise, the equation 3 in the aforementioned paper. What I want to do now is to carry out the equation 5 on my de-noised patches in 3D domain. Can someone please explain the equation 5, the aggregation scheme, a bit more lightly?
The intuition behind BM3D is the same as other non-local algorithms: an image is made by the juxtaposition of many, self-repetitive small image elements called patches. Thus, when looking at small enough pieces, it is possible to find enough similar patches that can be combined and denoised together.
The different non-local algorithms are different in the way used to combine these similar patches. BM3D relies on sparsity. This principle has become widely accepted now, and basically states that for many useful signals (including images) there exists at least one frame where this signal can be decomposed into a small number of meaningful components.
For images, DCT and wavelet basis are obvious candidates for the sparse frame: JPEG / JPEG2k compression rely on the transformation of the input images into 8-by-8 blocks followed by retaining the strongest coefficients, which is indeed a sparsifying transform.
When noise is added to a signal, its frequencies get mixed or added to the frequencies of the noise-free signal. For example, white noise adds non-zero values to the whole spectrum. Thus, finding a sparse approximate of a noisy signal (in a good basis of course) yields usually a good denoising performance.
The 3D transform of BM3D is actually a 2D+1 transform.
Algorithmically, this is what happens (I'll describe the mathematical intuition afterwards):
- make a stack of similar patches. All the patches of the stack are supposed to be small variants of the same ideal (noise free) patch. The variant term is important, you'll see why later;
- take a 1D transform (DCT / wavelet) along one of the horizontal axes;
- take a 1D transform (DCT / wavelet) along the other. What you have now is a stack of transformed patches;
- now, take the transform of choice (DCT / wavelet / Wiener filtering) along the vertical axis. What you are doing here is taking a transform (DCT / wavelet) or directly computing a denoised value (Wiener filtering) for each point in the vertical column;
- (optional) if you used DCT / wavelet transform, apply a sparsifying transform, usually hard thresholding.
- Then, compute the denied patches by inverting the different transforms (in reverse order);
- now, back project the result to the original patch locations. Patches can (and will) overlap when you denies the whole image. This raises an important practical question: how do you fuse the results on overlapping areas? Which denoised patch is the most reliable? BM3D uses Eq. 5 to fuse the results by a weighted averaging: the weight of a back-projected patch is inversely proportional to the variance in the patch stack. When the variance of the stack is high, it means that the patches did not really come from the same ideal one, and thus the result is marked as less reliable.
Mathematically, the intuition is that most patches should be sparse in some good domain (typically, DCT domain for $8 \times 8$ patches). BM3D also supposes that patches in a stack are variants of a same patch, and thus the variations between these transformed patches and the ideal patch (in the transformed domain) should also be sparse. This assumption is enforced by the 3rd transform and the hard thresholding step.