The code generates blocks of -nan for some input images

[ajeb78]

nantest.zip
I've found some images (typically mono images with mostly low pixel values and high dynamic range, such as unstretched astronomical images, when processed with BM3D, generate blocks of -nan in the output. Looking at the basic output, the -nan values appear to start off in the first stage of the algorithm usually as single pixels and then presumably get spread to a whole block in the second stage where processing of a block encounters the pixel with -nan value. An example file that creates -nan results is attached: the sigma for this file is 0.000105 (using bgnoise measurement from Siril).

[ajeb78]

I haven't got to the bottom of why this happens, but if it's causing you problems a quick and dirty fix is given below. This catches any pixels that are -nan at the end of the 1st and 2nd steps and eliminates the regions of -nan in the output image and replaces them based on nearest neighbours.

Change the end of bm3d_1st_step() to:
//! Final reconstruction
for (unsigned k = 0; k < width * height * chnls; k++) {
img_basic[k] = numerator[k] / denominator[k];
if (isnan(img_basic[k])) {
unsigned count = 0;
img_basic[k] = 0.f;
if (k > chnls-1) {
img_basic[k] += (isnan(img_basic[k-chnls]) ? 0 : img_basic[k-chnls]);
count = count + (isnan(img_basic[k-chnls]) ? 0 : 1);
}
if (k < count-chnls) {
img_basic[k] += (isnan(img_basic[k+chnls]) ? 0 : img_basic[k+chnls]);
count = count + (isnan(img_basic[k+chnls]) ? 0 : 1);
}
if (count != 0)
img_basic[k] /= count;
}
}
}

Change the end of bm3d_2nd_step() to:
//! Final reconstruction
for (unsigned k = 0; k < width * height * chnls; k++) {
img_denoised[k] = numerator[k] / denominator[k];
if (isnan(img_denoised[k])) {
unsigned count = 0;
img_denoised[k] = 0.f;
if (k > chnls-1) {
img_denoised[k] += (isnan(img_denoised[k-chnls]) ? 0 : img_denoised[k-chnls]);
count = count + (isnan(img_denoised[k-chnls]) ? 0 : 1);
}
if (k < count-chnls) {
img_denoised[k] += (isnan(img_denoised[k+chnls]) ? 0 : img_denoised[k+chnls]);
count = count + (isnan(img_denoised[k+chnls]) ? 0 : 1);
}
if (count != 0)
img_denoised[k] /= count;
}
}
}