0

I have a 2d gaussian whose center has been destroyed by pixel saturation. I need the center to be filled in because a poorly filled in center will confuse a neural network I'm trying to train. See below: enter image description here

The scattered nan values I can handle fairly easily, but the large cluster in the gaussian's center I cannot.

I've tried various methods to correct this, but none seem to work in the sense that the gaussian is filled in correctly.

Here are some other similar answers that I've tried:

Python Image Processing - How to remove certain contour and blend the value with surrounding pixels?

https://docs.astropy.org/en/stable/convolution/index.html

These work well for the small discrete nans floating around the image, but don't adequately address the center cluster.

This is what I get with convolution infilling: enter image description here I've taken slices of the centers as well. enter image description here enter image description here

I do actually have a reference image that does not have nans. However, the scaling of the pixel values are not constant, so I've made a function that takes into account the different scaling of each pixel.

def mult_mean_surround(s_arr, c_arr, coord):
    directions = np.array([[1,0],[-1,0],[0,1],[0,-1],[1,1],[1,-1],[-1,-1],[-1,1]])
    s = np.array([])
    for i in directions:
        try:
            if not np.isnan(s_arr[coord[0]+i[0],coord[1]+i[1]]):
                s=np.append(s,s_arr[coord[0]+i[0],coord[1]+i[1]]/c_arr[coord[0]+i[0],coord[1]+i[1]])
        except IndexError:
            pass
    if len(s)!=0:
        s_arr[coord[0],coord[1]] = c_arr[coord[0],coord[1]] *np.mean(s)

It copies the corresponding pixel values of the reference image and scales it to the correct amount.

Ideally, it would look something like this: enter image description here enter image description here The center is brighter than the rim and it looks more like a gaussian. However, this method is also substantially slower than the rest, so I'm not sure how to get around either of my issues. I've tried boosting speed with cupy to no luck as shown here: Boosting algorithm with cupy

If anyone has any helpful ideas, that would be great.

James Huang
  • 848
  • 1
  • 7
  • 35
  • did you try gaussian mixture model? assuming you have mixture of independent gaussians – qwr Aug 02 '22 at 20:16
  • can you explain your boosting approach more? anyway this is not really a programming question but a statistics question. – qwr Aug 02 '22 at 20:20
  • For each nan you copy the corresponding reference pixel and then multiply it based on the ratio of the surrounding 8 pixels (if they are not nans which are shown in the coords variable of the function). – James Huang Aug 02 '22 at 20:44
  • how is that boosting? anyway is GMM appropriate for your question at hand? – qwr Aug 02 '22 at 20:47
  • I'm talking about using cupy to boost the algorithm I already have or something else. And what is GMM – James Huang Aug 02 '22 at 21:38
  • gaussian mixture model. although from your question it's not clear if you have multiple independent gaussians which is the premise for GMM – qwr Aug 02 '22 at 23:04
  • Normalized convolution can help with this problem https://docs.astropy.org/en/stable/convolution/index.html#using-astropy-s-convolution-to-replace-bad-data but will not recover the central peak – keflavich Aug 26 '22 at 01:19

1 Answers1

1

I am assuming that you are filling the 'hole' with only one gaussian.

First make a mask of all the NaNs, i.e. NaN = 1, not NaN = 0.

You can do a neighbor-count check to remove all mask pixels with no neighbors, then use a clustering algorithm (like DBSCAN) to find the largest cluster of pixels.

Calculate the centroid, width (max x - min x), and height (max y - min y) of the resulting cluster.

You can then use the following code:

import math
def gaussian_fit(query_x,query_y,
centroid_x,centroid_y,
filter_w,filter_h,
sigma_at_edge = 1.0):
    x_coord = (query_x - centroid_x) * 2 / (filter_w * sigma_at_edge)
    y_coord = (query_y - centroid_y) * 2 / (filter_h * sigma_at_edge)
    return math.exp(-1.0*(x_coord**2+y_coord**2))

You may need to rescale the result by some constant.

Michael Sohnen
  • 953
  • 7
  • 15