failure of sep.extract service for saturated image

[bikramATastro]

Hello,

I have a 2048 X 30000 science image. The image contains few saturated objects (one of them is shown in the link below).

I attempted several times to run the sep.extract method with 100 % failure.
The error message I get is mostly " internal pixel buffer full:"

Then I did few modification to the "sep.extract" arguments have been made such as ,

1. "filter_kernel = None".
2. sep.set_extract_pixstack(600000) [default : 30000]
3. ''mask " parameter is also been changed, in order to avoid the saturated pixel

Different combinations of this three parameters have been implemented.
Still nothing works !!

If you would like, I could share you the link to the science_image as well as dark_frame for your further investigation.
Thanking you in advance for suggesting a solution.
Regards, Bikram

[kbarbary]

Just to check, did you ensure that the array is background-subtracted (using sep.Background), and that the detection threshold is not too low? Those are usually the main culprits.

[bikramATastro]

Yes, I have done that.Actually, after doing photometry for several images, I came across the image containing saturated object. Then the module failed.

[kbarbary]

I guess it's possible that the strongly saturated star has more than 600,000 pixels above threshold. Try setting the pixel stack to something crazy like 50,000,000 with sep.set_extract_pixstack. If that doesn't work, you can share the image with me and I'll investigate further.

[bikramATastro]

I have already increased the "sep.set_extract_pixstack" to very high values and then ran the module. But it ended up with "segmentation fault". I am running the module in HPC with a good memory space.

Here is the link to the science image >>
https://drive.google.com/file/d/0B2V-cmVqYSaEdl91Y2NyLUJzTm8/view?usp=sharing

Here is the link to dark frame >>
https://drive.google.com/file/d/0B2V-cmVqYSaEMWphaHVoQUJlLXc/view?usp=sharing

[kbarbary]

Can you post a minimal script to reproduce one of these errors?

When I read in the science image and run sep.Background I get

>>> f = fitsio.FITS("TDI_Red_12_26_05_2015.fit")
>>> data = f[0].read()
>>> bkg = sep.Background(data)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-3-e95be84d2529> in <module>()
----> 1 bkg = sep.Background(data)

sep.pyx in sep.Background.__cinit__ (sep.c:4241)()

sep.pyx in sep._get_sep_dtype (sep.c:2783)()

ValueError: input array dtype not supported: uint16

which is what I expect because the datatype in the FITS image is indeed uint16, which sep doesn't currently support (without explicitly converting the data array).

[bikramATastro]

I use the "pyfits" module to read the fits files. Following is a sample script.

import sep
import numpy as np
import pyfits
data, header = pyfits.getdata('TDI_Red_12_26_05_2015.fit',0,header=True)

this is a TDI(time delay integration) image so we skip first 2048 columns

...
data = data[:,2048:]
dark, header = pyfits.getdata('super_dark_2_to_7.fit',0,header=True)

dark subtraction

...
dark_subtracted_data = np.subtract(data,dark)
bkg = sep.Background(dark_subtracted_data)

[bikramATastro]

Increasing the "sep.set_extract_pixstack" to very high value gives rise to an error ::
*** Error in `python': double free or corruption (!prev): 0x0000000002805ae0 ***

Is there any solution to this ?

[kbarbary]

I can't reproduce the original problem. I'm using the following script:

from astropy.io.fits import getdata
import sep

data, header = getdata('TDI_Red_12_26_05_2015.fit',0,header=True)

#this is a TDI(time delay integration) image so we skip first 2048 columns
data = data[:,2048:]

dark, header = getdata('super_dark_2_to_7.fit',0,header=True)
dark_subtracted_data = data - dark

bkg = sep.Background(dark_subtracted_data)

bkg_subtracted_data = dark_subtracted_data - bkg

objects = sep.extract(bkg_subtracted_data, 1.5 * bkg.globalrms)

print(len(objects))  # 2679 objects

[kbarbary]

Perhaps you could post a full minimal reproducible example, to make sure we're on the same page?

[bikramATastro]

Now I got the possible reason for my problem.

I do not use the default parameter values in "sep.Background".

bkg = sep.Background(dark_subtracted_data, bw=1000, bh=1000, fw=100, fh=100)

In principle the sky background should be smooth. If we take the default parameter the background it generates is not smooth and local variations can be seen.

So, I conclude that, if I take higher value of the "sep.Background( )" it gives me smoother background but, it can not account for the error given by "sep.extract( )". if I take default value of the "sep.Background( )" it does not give me smoother background but, it solves the problem for the method "sep.extract( )".

[kbarbary]

Note that fw and fh are in units of background boxes. Since your boxes are 1000 x 1000, this is 100000 x 100000, bigger than your image. Did you try much smaller fw, fh values?

[bikramATastro]

Thank you @kbarbary
Yes, I tried with bw=300, bh=300, fw=3, fh=3. I think this is close to optimum combination of these four parameters. And, I got better results. But, still the background is questionable.
I may try fitting 2D polynomial to the sky_background and get a smoother surface.

[kbarbary]

There's definitely some weird stuff going on with sep.Background. It's not giving results I would expect for this image, and I do get a segfault for some parameter combinations.

[bikramATastro]

Yes, I agree. It's the same here. But, I could do 2D polynomial fit to the frame obtained by sep.Background(). The result is encouraging.

[kbarbary]

I was mistaken, sep.Background seems fine to me (though I did find and fix a segfault that occurred when the boxsize is equal to the image size.) For example, using bh=1000, bw=1000, fw=1, fh=1 gives a quite smooth background:

Note that the background uses a cubic spline, so it is already a 2D polynomial. Fitting another one should not be necessary.

The bright spot in the background under the saturated star is somewhat inevitable with the statistical background technique used here. If you want to get rid of that, I suggest a two step procedure where one detects objects with sep.extract, then create a mask array based on the objects and rerun sep.Background() with that mask.

[bikramATastro]

That seems to be a nice idea. Thank you

[kbarbary]

I'm closing this as resolved (as far as I can tell). Feel free to comment or open another more specific issue if you find any other problems.