Question on the implementation of Non Local Means denoising

[TimZaragori]

With the aim of implementing an adapted version of NLM algorithm for dynamic images, I try to understand how the fast version is implemented in 3D in scikit-image (skimage Cython code). From what I see the version implemented corresponds to the algorithm 3 of Froment et al. 2015. In the previous article, the main loop corresponds to a nested loop to apply all the possible shifts, thus in python:

for t_pln in range(-d, d + 1):
    for t_row in range(-d, d + 1):
        for t_col in range(-d, d + 1):
            Do something

However, the scikit-image implementation choose another strategy:

for t_pln in range(-d, d + 1):
    for t_row in range(-d, d + 1):
        for t_col in range(0, d + 1):
            Do something

The shifts done in the third dimension are in the range (0, d) instead of (-d, d). I think that the aim is to reduce the number of shifts done by 2 and thus speeding up the algorithm. From what I understand this missing shifts are compensated during the attribution of the weights with an alpha parameter to take into account "patches on the same column distance is computed twice in this case" (according to the associated comment):

for pln in range(pln_dist_min, pln_dist_max):
    for row in range(row_dist_min, row_dist_max):
        for col in range(col_dist_min, col_dist_max):
            weights[pln, row, col] += alpha * weight
            weights[pln + t_pln, row + t_row, col + t_col] += alpha * weight

However, it's hard for me to visualize how this implementation is similar to the original one. Moreover, I can't find the same results between original implementation and scikit-image one when I try to see how many times each voxel is accessed for the weights assignation. Could someone explain me the thinking behind this implementation?

[jni]

I wonder if @emmanuelle who implemented this can comment? For reference, here are what I think are the most relevant pull requests for that code — there might be some insights in the discussion there:

874 (wow)

2878

2890

[TimZaragori]

Thx for the reading. Just for information find below the code I used to see the number of time each voxel is accessed for the weights (which is for the most part is a copy paste from skimage Cython code) and where I find differences but I may have misunderstood something and do it wrong. Is there something that work differently in Cython and Python and affect that result?

import numpy as np

def _weights_skimage(image, s, d):
    if s % 2 == 0:
        s += 1  # odd value for symmetric patch

    offset = int(s / 2)
    # Image padding: we need to account for patch size, possible shift,
    # + 1 for the boundary effects in finite differences
    pad_size = offset + d + 1
    padded = np.ascontiguousarray(np.pad(image, pad_size, mode='reflect').astype(np.float64))
    weights = np.zeros_like(padded)
    n_pln, n_row, n_col = padded.shape[0], padded.shape[1], padded.shape[2]

    # Outer loops on patch shifts
    # With t2 >= 0, reference patch is always on the left of test patch
    # Iterate over shifts along the plane axis
    for t_pln in range(-d, d + 1):
        pln_dist_min = max(offset, offset - t_pln)
        pln_dist_max = min(n_pln - offset, n_pln - offset - t_pln)
                # alpha is to account for patches on the same column
                # distance is computed twice in this case
        if t_pln == 0:
            alpha = 1.0
        else:
            alpha = 0.5
        # Iterate over shifts along the row axis
        for t_row in range(-d, d + 1):
            row_dist_min = max(offset, offset - t_row)
            row_dist_max = min(n_row - offset, n_row - offset - t_row)
            if t_row == 0:
                alpha = 1.0
            else:
                alpha = 0.5
            # Iterate over shifts along the column axis
            for t_col in range(0, d + 1):
                col_dist_min = offset
                col_dist_max = n_col - offset - t_col
                # Inner loops on pixel coordinates
                # Iterate over planes, taking offset and shift into account
                for pln in range(pln_dist_min, pln_dist_max):
                    # Iterate over rows, taking offset and shift
                    # into account
                    for row in range(row_dist_min, row_dist_max):
                        # Iterate over columns
                        for col in range(col_dist_min, col_dist_max):
                            # Accumulate weights for the different shifts
                            weights[pln, row, col] += alpha * 1
                            weights[pln + t_pln, row + t_row, col + t_col] += alpha * 1
                alpha = 1.0
    return weights

def _weights_froment(image, s, d):
    if s % 2 == 0:
        s += 1  # odd value for symmetric patch

    offset = int(s / 2)
    # Image padding: we need to account for patch size, possible shift,
    # + 1 for the boundary effects in finite differences
    pad_size = offset + d + 1
    padded = np.ascontiguousarray(np.pad(image, pad_size, mode='reflect').astype(np.float64))
    weights = np.zeros_like(padded)
    n_pln, n_row, n_col = padded.shape[0], padded.shape[1], padded.shape[2]

    # Outer loops on patch shifts
    # With t2 >= 0, reference patch is always on the left of test patch
    # Iterate over shifts along the plane axis
    for t_pln in range(-d, d + 1):
        pln_dist_min = max(offset, offset - t_pln)
        pln_dist_max = min(n_pln - offset, n_pln - offset - t_pln)
                # alpha is to account for patches on the same column
                # distance is computed twice in this case
        # Iterate over shifts along the row axis
        for t_row in range(-d, d + 1):
            row_dist_min = max(offset, offset - t_row)
            row_dist_max = min(n_row - offset, n_row - offset - t_row)
            # Iterate over shifts along the column axis
            for t_col in range(-d, d + 1):
                col_dist_min = max(offset, offset - t_col)
                col_dist_max = min(n_col - offset, n_col - offset - t_col)
                # Inner loops on pixel coordinates
                # Iterate over planes, taking offset and shift into account
                for pln in range(pln_dist_min, pln_dist_max):
                    # Iterate over rows, taking offset and shift
                    # into account
                    for row in range(row_dist_min, row_dist_max):
                        # Iterate over columns
                        for col in range(col_dist_min, col_dist_max):
                            # Accumulate weights for the different shifts
                            weights[pln, row, col] += 1
    return weights

image = np.random.rand(50, 50, 25)
weights_skimage = _weights_skimage(image, 3, 5)
weights_froment = _weights_froment(image, 3, 5)
print(np.all(weights_skimage  == weights_froment)) # False

[TimZaragori]

@jni I read with attention the PR you sent me to understand the whole history. It seems that the first implementation was similar to the function _weights_froment above. The acceleration was introduced in #874 with this commit. In this form, I am better able to target what I am struggling to understand in the current implementation : the alpha coefficient. I understand that d(x,y) = d(y,x) and that the computation is done twice if we loop over (-d, d) in the last dimension. So it only loops over (0 ,d) and then compensates in the weights attribution with a coefficient for the combinations that will be done twice. In 3D, the alpha coefficient is defined as alpha = 0.5 if t_col == 0 and (t_pln is not 0 or t_row is not 0) else alpha = 1 (commit). However, voxels are processed twice without alpha=0.5 if (t_pln, t_row, t_col) == (0, 0, 0) right? For this combination, won't the voxels be processed once too often? Next, this part was modified to its actual implementation in #4322 with this commit. I have never done Cython and only a few C but it seems that now the alpha coefficient is defined as alpha = 1 if t_pln == 0 and t_row == 0 else alpha = 0.5 which doesn't seem to be the same as the original implantation. According to all of this I modified the function with the speed trick and alpha = 0.5 if t_col == 0 else 1 and got the same results as weights_froment above assuming this is the ground truth but could you confirm it to me ?

def _weights_skimage_modified(image, s, d):
    if s % 2 == 0:
        s += 1  # odd value for symmetric patch

    offset = int(s / 2)
    # Image padding: we need to account for patch size, possible shift,
    # + 1 for the boundary effects in finite differences
    pad_size = offset + d + 1
    padded = np.ascontiguousarray(np.pad(image, pad_size, mode='reflect').astype(np.float64))
    weights = np.zeros_like(padded)
    n_pln, n_row, n_col = padded.shape[0], padded.shape[1], padded.shape[2]

    # Outer loops on patch shifts
    # With t2 >= 0, reference patch is always on the left of test patch
    # Iterate over shifts along the plane axis
    for t_pln in range(-d, d + 1):
        pln_dist_min = max(offset, offset - t_pln)
        pln_dist_max = min(n_pln - offset, n_pln - offset - t_pln)
                # alpha is to account for patches on the same column
                # distance is computed twice in this case
        # Iterate over shifts along the row axis
        for t_row in range(-d, d + 1):
            row_dist_min = max(offset, offset - t_row)
            row_dist_max = min(n_row - offset, n_row - offset - t_row)
            # Iterate over shifts along the column axis
            for t_col in range(0, d + 1):
                alpha = 0.5 if t_col == 0 else 1
                col_dist_min = offset
                col_dist_max = n_col - offset - t_col
                # Inner loops on pixel coordinates
                # Iterate over planes, taking offset and shift into account
                for pln in range(pln_dist_min, pln_dist_max):
                    # Iterate over rows, taking offset and shift
                    # into account
                    for row in range(row_dist_min, row_dist_max):
                        # Iterate over columns
                        for col in range(col_dist_min, col_dist_max):
                            # Accumulate weights for the different shifts
                            weights[pln, row, col] += alpha * 1
                            weights[pln + t_pln, row + t_row, col + t_col] += alpha * 1
    return weights
weights_skimage_modified = _weights_skimage_modified(image, 3, 5)
print(np.all(weights_skimage_modified  == weights_froment)) # True

[jni]

I'm hoping @emmanuelle or @grlee77 can comment here as I am a bit under the pump and not familiar with this code at all. But thanks for the reading and detailed questions, they make discussion much easier!

[TimZaragori]

@jni Any news on this question ? I saw that the functions had been updated

[grlee77]

comments moved over from Discussion

grlee77 I am taking a look now. I made a minor update to a couple of your prior comments just to fix a link and add code highlighting. I agree that the proposed change here should be made for both the 2d, 3d and 4d cases. I will convert this discussion back to an issue so we can close it. If you view the weights output by the current implementation it is actually very close to the desired ones, but is off slightly. The highest value in the "Froment" weights matrix is 1331 while the maximum difference between the two implementations is only 11, so any difference in output is likely to be very subtle (hopefully not visually detectable). We definitely should still fix it though!

[grlee77]

closed by #5923