CODEWITHSUNDEEP
image-processingscipy
pootle
pootle

Reputation: 547

Using scipy.ndimage.uniform_filter to find stars in astro photo, but puzzled by results

I am searching for stars in an umage of night sky, after making a mask I use scipy.ndimage.uniform_filter at different sizes to find the stars. It looks to work reasonably well, but I expected once I used a small enough size, I would just get more hits as I reduced the size further, but it doesn't do this consistently, I am just a bit baffled by this.

There is an extract from around one of the hit areas at the bottom

The code below gives me:

size: 3, len: 621
size: 4, len: 340
size: 5, len: 200
size: 6, len: 0
size: 7, len: 0
size: 8, len: 24
size: 9, len: 8
size: 10, len: 0
size: 11, len: 0
size: 12, len: 0

Why do size 6 & 7 give zero hits? This seems totally bizarre to me!

 def __init__(self, filename):
        self.good=False
        self.img = scipy.ndimage.imread(filename, flatten=True)

    def checkcandidates(self, meanfact=3.0, maxwindow=25):
        mask = self.img > self.img.mean()*meanfact
        for wsize in range(3,maxwindow):
            m2 = scipy.ndimage.uniform_filter(mask, size=wsize)
            xc,yc = m2.nonzero()
            print("size: %d, len: %d" %(wsize, len(xc)))

Here's part of the mask centred on one of the stars:

>>> sc1.showCoords(1360,493,10,usemask=True)    
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0
0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0
0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0
0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0
0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0
0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0
0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0
0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Upvotes: 0

Views: 872

Answers (1)

Warren Weckesser
Warren Weckesser

Reputation: 114811

This looks like a bug, or at least a nasty implementation detail that will result in bugs in users' code.

First, read the note in the uniform_filter docstring:

The multi-dimensional filter is implemented as a sequence of
one-dimensional uniform filters. The intermediate arrays are stored
in the same data type as the output. Therefore, for output types
with a limited precision, the results may be imprecise because
intermediate results may be stored with insufficient precision.

So let's look at how one row of your input array is processed by uniform_filter1d for different size filters.

Here's a small one-dimensional input:

In [416]: x
Out[416]: array([0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0])

Apply uniform_filter1d with increasing sizes:

In [417]: from scipy.ndimage.filters import uniform_filter1d

In [418]: uniform_filter1d(x, 3)
Out[418]: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0])

In [419]: uniform_filter1d(x, 4)
Out[419]: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0])

In [420]: uniform_filter1d(x, 5)
Out[420]: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0])

In [421]: uniform_filter1d(x, 6)
Out[421]: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

In [422]: uniform_filter1d(x, 7)
Out[422]: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

In [423]: uniform_filter1d(x, 8)
Out[423]: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0])

In [424]: uniform_filter1d(x, 9)
Out[424]: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

Like your example, the output is all zeros when the size is 6 or 7.

I suspect this is a floating point precision problem. Note what happens when we make the input an array of floating point values:

In [439]: f = uniform_filter1d(x.astype(float), 6)

In [440]: f
Out[440]: 
array([  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         1.66666667e-01,   3.33333333e-01,   5.00000000e-01,
         6.66666667e-01,   8.33333333e-01,   1.00000000e+00,
         1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
         8.33333333e-01,   6.66666667e-01,   5.00000000e-01,
         3.33333333e-01,   1.66666667e-01,   5.55111512e-17,
         5.55111512e-17,   5.55111512e-17])

In [441]: f.max()
Out[441]: 0.99999999999999989

So the intermediate values computed using floating point do not give the expected value of 1 in the "middle" of that output. When this array is converted back to the input data type (int), the result is all zeros:

In [442]: f.astype(int)
Out[442]: array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

Given that behavior, I recommend converting your input array to floating point before calling uniform_filter, and adding a final step that converts the result back to integers in a way that you control, and that matches how you want to classify a "hit". Or even use a different function altogether.

Upvotes: 1

Related Questions