LACosmic Differences

[larrybradley]

I just learned that ccdproc has implemented LACosmic and I did a quick comparison with the version in photutils (which I wrote). I wanted to give you some feedback on the differences that I've found.

I had previously compared the photutils version extensively with other implementations and here I ran the same example with ccdproc. I've found that the ccdproc version doesn't produce results that are consistent with any other implementation of the algorithm. ccdproc's version produces only a cosmic-ray mask, so that's what I compare here:

ccdproc.cosmicray_lacosmic cosmic-ray mask:

photutils.lacosmic mask:

IDL version mask

cosmics.py mask:

ccdproc's mask is clearly different than the others. There are actually a few small differences among the other three masks, but only along the top-edge pixels.

I looked into the ccdproc code and the first thing I noticed is that it does not perform any iterations. And, running ccdproc.cosmicray_lacosmic multiple times won't give the same results because each iteration needs the mask of the previous iteration.

I limited photutils.lacosmic to one iteration and it still detects significantly more cosmic rays than ccdproc.cosmicray_lacosmic:

Looking further at the ccdproc code, I think the big difference is the step where neighboring pixels are removed (for one, it should be a two-step process where the neighbors of the neighbors are also checked with a lower threshold). For reference, here's my implementation in photutils: https://github.com/astropy/photutils/blob/master/photutils/detection/lacosmic.py#L150-L156

[mwcraig]

@larrybradley -- thank you very much for the thorough report: this kind of thing (people with more expertise than I have like you and @crawfordsm looking at a package I've had a hand in) is what attracted me to astropy in the first place.

I won't get a chance to look at it in detail until later (am at scipy).

[crawfordsm]

Hi Larry-- Thanks for the comparison. Could you provide a link to the data and the parameters that you used to make the comparison? It will be easier to determine whether or not it is just parameter differences or if there is something inherent in the code.

Previous tests I've run on the algorithm with cosmics.py and lacosmic.cl have not shown such significant differences, so I'm curious whether or not it is more of a documentation issue or a bug in the code.

The cosmic ray tasks were outlined in the initial API set up for ccdprocand they were set up to share some similar infrastructure. So the typically user shouldn't actually be calling cosmicray_lacosmic but should just be calling cosmicray_clean with the cosmic ray identifying algorithm of their choice. It might not be the best way to do it, but it was the way we initially set it up in the API.

Please also see #125 for information about another lacosmic version.

[larrybradley]

@crawfordsm, I've placed the data and resulting mask/cleaned images here: https://www.dropbox.com/sh/yixlek3v4fomj7o/AADDKLTdfZVnH_sMyobHzc5pa The image is one chip of an HST WFC3/UVIS F606W FLT file. I've also placed a script there so that you can exactly replicate my results. Let me know if you need anything else.

I simply looked for "lacosmic" and ran that function. :smiley: If you really don't want anyone to use it, then I'd recommend making it a private function.

[crawfordsm]

Thanks for the data--I set it up to run a test using cosmicray_clean with gbox=2 and for one iteration and the results are very similar to the results from photutils. This is probably the most comparable test as the photutil version. The photutil does the growing automatically, where it is optionally in the ccdproc version.

However at the same time, I can see where the setup currently in ccdproc is pretty confusing. I've set up so that growing and replacing pixels is done in cosmicray_clean whereas there probably be a better way to do this than the current set up. This probably should be something changed to make it more user friendly/obvious to use.

[larrybradley]

Thanks, @crawfordsm. Yes, the results are better with cosmicray_clean, which grows the cosmic rays by gbox pixels. The differences in the images appear to be due to how the neighbors ("growing") are handled and the lack of iterations.

BTW, the background parameter in cosmicray_clean initially confused me before reading the docs because what you really want is a noise/error image. Perhaps a better parameter name is background_rms or error_image or noise_image....?

[crawfordsm]

Hopefully, I've address some of your comments with #147

I've changed background to variance_image, removed cosmicray_clean, adding growing and replacing pixels to cosmicray_lacosmic, and set it so it returns both a data array and a mask.

[larrybradley]

Why do you call this variance_image?

Background variance level.   It should have the same shape as data
          as data. This is the same as the noise array in lacosmic.cl

It's not a variance image but a stddev (or rms or error or uncertainty or sigma; all refer to the same) image. The variance image would be stddev_image**2.

Curiously, I see that this same usage of "variance" is used elsewhere in ccdproc, e.g. create_variance():

var = (gain_value * ccd_data.data + readnoise_value ** 2) ** 0.5

Because you take the sqrt (i.e., ** 0.5), this is the error (sigma), not the variance (sigma**2). Can you please explain why you call this "variance"?

[crawfordsm]

Well, the create_variance frame is historical and that it was originally going to be variance and not error. It got changed and the name never get updated. However, it should actually probably be changed to the correct name.

The variance frame in cosmic ray is because I am currently extremely jet lagged and was rushing a bit to do it and should have known better and used error.

So thanks for pointing this out. I'll update the names in the cosmic ray tasks

[mwcraig]

As long as you are looking, there is a variable maskf defined at https://github.com/astropy/ccdproc/blob/master/ccdproc/core.py#L979 that isn't used in cosmicray_lacosmic. Is it supposed to be?