Detecting central blob in image while discarding noisy pixels

TL:DR - I have several sets of images of a detector, and I want to detect the central blob of most intense pixels. However, noisy pixels makes thresholding the image difficult.

Long version: I have a set of 3-d arrays, where the first two dimensions are spatial (corresponding to pixels on a detector) and the 3rd is temporal. The different arrays in the set are replications of the same measurement. After applying fft along the 3rd dimension, I extract a specific frequency and am left with a set of 2d maps of amplitude (actually, it's the amplitude divided by the DC value) and phase. Below is an example of an "easy" set:

Eventually, I want to analyze the phase of the maps. However, I only want to include the phase of pixels who had a significant DC-normalized Fourier amplitude in my analysis. In the maps above, it is fairly simple to detect the largest blob in the dc-normalized maps by thresholding and to take the corresponding pixels in the phase maps:

threshold = 0.75;
peakVal = max(thisAmpMap, [], 'all');
aboveThreshold = thisAmpMap > peakVal*threshold;
largestBlob = double(bwareafilt(aboveThreshold, 1));

However, in some sets of images the maps are noisier; specifically, there are noisy intense pixels along the sides of the images:

Of course, binarizing this image by the maximal value would lead to a detection of an area near the edges, not where the actual signal is.

So far, my attempts at solving this include:

• averaging the dc-normalized maps of a single set. This didn't help in eliminating the effect of the noisy pixels from the thresholding, but it did help with a different case - where the central blob would have side blobs, which are stronger than the central one in some repetitions in the set.
• Smoothing the images before binarizing, by applying medfilt2 several times, which was somewhat effective. Applying medfilt2 with a larger neighborhood helped more, but the detection would still fail in some cases, specifically when the noisy pixels are numerous and not spread-out.
• Cropping the image edges - simply discarding X% of the outer pixels for the calculation of the maximal value for thresholding. This works, but I would like a criteria that's image dependent, since different sets have different margins of noisy pixels.

I have included example data, organized in 3d arrays structured as [repetition X height X width]. They include an "easy" set for reference, and a noisy set where the detection is harder.

Answer 1

[Bjorn Gustavsson]

Perhaps it would be possible to fit a sum of 2-D Gaussians (2-3-4?) to your images using some robust norm (huber-like with an adequate setting to supress the contribution from your salt-noise-pixels), with constraints on how narrow your Gaussian blobs would be allowed to be (2-4 pixels or bigger?). That would require you to get a method to efficiently get sensible start-guesses for the centre-points of the Gaussians, perhaps you can do that manually in some reasonable time? Perhaps you can get those after a significantly wide median-filtering?

HTH

[Yotam Federman]

Thanks for the suggestion! I will look into it and give an update here.

[Bjorn Gustavsson]

It can be a bit fidgety both to make a parameterized multiple-Gaussians function and to get the parameter-optimization to converge. The 2-D Gaussians would become something like this:

function I2DmG = sumof2DGaussians(pars,X,Y)
I2DmG = 0*X;
for i1 = 1:(numel(pars)/6)
    I0 = pars(1+6*(i1-1));
    x0 = pars(2+6*(i1-1));
    y0 = pars(3+6*(i1-1));
    dx = pars(4+6*(i1-1));
    dy = pars(5+6*(i1-1));
    theta = pars(6+6*(i1-1));
    I2DmG = I2DmG + I0*exp(-(X-x0).^2/dx^2-(Y-y0).^2/dy^2);
end

You'll have to fix the rotation of the Gaussians. The key (in my experience) to the fitting is to have good enough start-guesses of the parameters (there is nothing that beats starting close to the goal), but you should find use of either lsqnonlin or fminsearchbnd (on the file exchange), there are also newer optimization-packages there like minimize.