Side pseudo-reference pixels for subarrays don't do anything in pipeline

JarronL commented 1 year ago

The JWST pipeline doesn't bother to perform any 1/f noise correction with the side pixels for subarray data even if they are manually flagged as reference pixels. We will need to come up with our own 1/f striping correction that doesn't depend on the pipeline. Consider a port of NSClean (https://webb.nasa.gov/content/forScientists/publications.html) for NIRCam purposes.

AarynnCarter commented 12 months ago

Duplicate of #80 and potentially a step towards a solution with #89. Thanks for sharing the NSClean paper!

As an addition, we know from transit observations that ramps are better estimated when correcting 1/f noise at the group level (particularly for shorter ramps). I don't think NSClean does that but it might be possible to adapt it.

JarronL commented 12 months ago

I did some brief tests with NS clean on group-level data. Overall, it seems to work pretty well in cleaning up the spatial 1/f noise, but takes about 1 minute per frame to execute (MacBook Pro w/ M2 Max processor). There are provisions for using cupy built-in, but I don't think that's supported on the Mac, so currently relegated to using numpy. I have not tried on my work desktop (Ubuntu) yet.

The other thing to note is that there is over-subtraction near the edges of the bright signal.

So, instead, I decided to test out a function I wrote a while ago in pynrc called channel_smooth_savgol (pynrc.reduce.channel_smooth_savgol()) that uses a Savitsky-Golay filter to smooth the row data. Turns out that with a window length of 127 pixels, I can pretty effectively reproduce the NSClean model in <10 msec. And it gets rid of the over-subtraction. In the process of implementing into spaceKLIP.

JarronL commented 12 months ago

Upon implementation, one potential problem with these methods is that they fit and subtract background signal offsets in the ND squares. So, I can either try to mask out these regions, or just leave them be with the knowledge that the ND regions are going to be subtracted with PSF references processed in the same way (so it effectively cancels out). So, the background levels across the full array (inside and out of the ND squares) will removed, but any signals from background objects or disks should be preserved.

Left is the original rate.fits image without 1/f noise subtraction, while the middle shows the 1/f noise model subtracted from the data at the group level, then ramp fitting. The final panel shows the difference between these two (left - middle). PSF signal levels are suppressed by about 1%. Not sure how this might affect companion signals after PSF subtraction.

JarronL commented 12 months ago

Just had a thought that this would also be problematic for extended faint signals, like disks. Probably need to think more about how to generate the masks for selection of background pixels. That would likely be something performed at a much later step, then come back to 1/f correction with more informed signal masks to exclude things like disks, planets, ND squares, etc.

AarynnCarter commented 12 months ago

@JarronL Looks really nice, thanks for taking a look at this. It looks to me that both the NSClean method and the SavGol model could do with more aggressive masking of bright pixels. Alternatively we could define regions to compute the 1/f correction over, mask bright pixels within those regions, then mask any pixels outside of those regions.

The right hand plot of your most recent figure is a little concerning to me, and if I squint I can convince myself that there is an oversubtracted halo around the target source in the middle panel.

JarronL commented 12 months ago

Yep, agreed on both points. The Data Mask panel in that first plot was an early attempt, but there is the ability to grow that masked out region and be more aggressive. In the end, I'm not sure how much of an improvement this would even be on this particular dataset (1386). Might be more applicable to full frame data. Anyway, I've added a class that allows for modeling the 1/f noise by adding an arbitrary mask (a default is calculated based on the rate images), and then fitting the background pixels with a variety of models including Savitzky-Golay filter (with interpolation over the masked data), median of each row, or robust mean of each row. Currently adding in exclusion of ND squares based on aperture.

JarronL commented 12 months ago

Okay, updated things to include:

Masking of ND Square regions and COM mask holder that depend on aperture. I already had models of these in pynrc from pre-flight simulations, so just ported them over to webbpsf_ext so that spaceKLIP can access them.
Default masking of circular region at mask aperture location with r=50 pixels (size should maybe depend on wavelength and SW/LW channel?).
Iterative savgol modeling that masks regions seemingly fitting significant signal relative to the median row model. This seems to do a good job of getting around the over-subtraction from faint PSF "glow."

Here are the rate files results relative to no 1/f correction:

Need to try this method on data that exhibits clear 1/f noise in the rate files as this effect is not very apparent in the 1386 dataset.

AarynnCarter commented 11 months ago

Looks great @JarronL! I don't know any datasets with significant 1/f noise, but I think it's possible to run Stage 1 on just the first 2-5 groups of each integration. That should make the 1/f noise less likely to average out in the rate files and the PSF should be much fainter. As a side note, I typically leave things in the _ints.fits format for my reductions.