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; } } }