Neglecting bad data points; NaN; imresize

Hello Everybody,

I'm working with images where each pixel value represents a number of counts. (X-ray scattering data). Some pixels turn out to malfunction, either with zero counts, or with a saturation. To remove these, we normally create a mask (kind of the inverse of "region of interest", to neglect those bad points/pixels.

Usually, I would make this mask as consisting of "ones" and "NaNs", so that when I multiply my data image/matrix by the mask, I get NaNs for all the bad pixels.

Now, I can sum, or make means using "nansum" or "nanmean", because these operations know that the NaNs don't count for the statistics, which is great. This works for more than half of my work!

However, sometimes, I get data with poor statistics, and to overcome this difficulty, I'm using the "imresize" function, which conveniently "bins" the data automatically to smaller images. The problem, is that "imresize" counts NaN, making a whole box as an invalid square.

To exemplify, say pic1 is my original image, mask is my mask, where I want to hide one part of the data, and pic2 is the image with the bad pixels removed. pic3 is the resized image for improved statistics:

if true
  % %
%
%
pic1=reshape(1:16,4,4)'
%
mask=ones(size(pic1)); mask(2:3,2)=NaN
%
%Now to apply the mask, I get the good data, with the bad pixels removed. 
%
pic2= pic1.*mask
%
%Because pic2 has poor statistics, I want to resize the image to half (or other factors) its %length, to %improve the number of counts, at the cost of loss in spatial resolution:
%
bin=0.5  %0.5, bin factor, it compresses the original matrix to half its size
pic3=imresize(pic2,0.5,'box')
%
%
%
%
%
%
%
end

Unfortunately, I get pic3 =[NaN 5.5000; NaN 13.5000]

instead of pic3=[2.667 5.5000; 12 13.5000]

which would be the result if NaN was neglected.

So, even though I am not counting the bad pixels, which is good, they are propagating in the compression, invalidating some of the good ones.

Do you have any suggestions on how to get around this?

Either using a different strategy in the mask? (replace NaN by something else? zero won't work, as it will affect the statistics)

Maybe there is a symbol I could use to replace for NaN, that the system knows it should not be counted? (again, zero would not work)

Or by using something different than "imresize"?

I could create a loop to manually resize the matrix, and use the "nanmean"… but that would make everything too complicated.

Any suggestion is appreciated!

Thanks!

Answer 1

[Kelly Kearney]

Here's a suggestion to replace the imresize call, using histc and accumarray:

fac = 0.5;
dx = 1./fac;

[r,c] = ndgrid(1:size(pic2,1), 1:size(pic2,2));

[n, ibin] = histc(r(:), 0.5:dx:size(pic2,1)+0.5);
[n, jbin] = histc(c(:), 0.5:dx:size(pic2,2)+0.5);

nr = max(ibin);
nc = max(jbin);

idx = sub2ind([nr nc], ibin, jbin); 
pic3 = accumarray(idx, pic2(:), [nr*nc 1], @nanmean);
pic3 = reshape(pic3, nr, nc)

Cleaner code than looping, though I'm not sure how it compares performance-wise... The initial setup index stuff would only need to be run once, then the imresize could be replaced by a one-liner:

 pic3 = reshape(accumarray(idx, test(:), [nr*nc 1], @nanmean), nr, nc);

[Bruno]

Thanks Kelly, This works really well. It is not as fast has the built in imresize, of course. It takes ca. 5s to convert a 1680 x 1476 image to a quarter of its size, compared to ~1s from imresize, but that is perfectly fine for my needs. Thanks again!

Answer 2

[Bruno]

Hello again,

I ended up coming with a solution that is not entirely satisfactory, but still improves quite a lot. Basically, I use two "reshape" and "transpose" steps, coupled with "nanmean".

So, why is this not a perfect solution? Performance apart, the problem is that "nanmean", just like "mean", has to be performed either along columns or rows. So, it requires two "nanmean" to get the mean of a box. To get the mean of a box, one could use "mean2", but that does not work here because of the NaN entries. That is why we use "nanmean". By using two steps of "nanmean", wherever there is a box with one NaN entry, the relative weights of each number get messed up. Please look at the example below to see what I mean.

In the first box, the mean should be (1+2+9)/3=2.67 (neglecting the NaN). However, because of the two-step mean, what we get is ((1+5)/2+2)/2=2.5… This is not too bad, but still not right…

Any suggestions? maybe a way of using one-step "nanmean"? or some other way of using "reshape", or something completely different?

Thanks! (the code is below)

if true
  % code
pic1=reshape(1:16,4,4)';
%
mask = ones(4,4); mask(2:3,2)=NaN;
%
pic2=pic1.*mask
%
ColSize=size(pic2,1); RowSize=size(pic2,2);
bin=0.5;
  for bbb=1:10
  if ceil(ColSize*bin) ~= floor(ColSize*bin)
      pic2(end,:)=[];
  end
  if ceil(RowSize*bin) ~= floor(RowSize*bin)
      pic2(:,end)=[];
  end
  ColSize=size(pic2,1); RowSize=size(pic2,2);
  end
pic3 = reshape(pic2,1/bin,[]);
pic3 = nanmean(pic3);
pic3 = reshape(pic3,ColSize*bin,[]);
pic3 = pic3';
pic3 = reshape(pic3,1/bin,[]);
pic3 = nanmean(pic3);
pic3 = reshape(pic3,RowSize*bin,[]);
pic3 = pic3'
end

PS: the loop code deletes rows or columns to make sure that the reshape ends up in a integer by integer matrix. (i might change it to add rows and columns instead…). The "bin" parameter is a fraction (e.g. 0.5 instead of 2, which would make more sense), due to "historical reasons".

Answer 3

[Image Analyst]

Unless I misunderstood, I don't see any reason to mess around with nans or reshaping. You have basically salt and pepper noise (either dead pixels at 0 or saturated pixels). A common fix for rthis is the modified median filter. Basically take the median filter. This replaces bad pixels with good values. But only replace the bad pixels, not all of them. So threshold the image to find bad pixels and replace only those like

outputImage(badPixels) = medianFilteredImage(badPixels);

Attached are full demos for both gray scale images and color images.

[Bruno]

Thanks Image Analyst for the great comment! In this particular case we have "salt and pepper noise", but also inactive stripes on the detector (technological limitations) with a width of ~20 pixels. Because of these stripes, we really need to use the mask. But your suggestion and links were really interesting and will definitely be useful for other parts of my work! Thanks again!