Repository navigation
LACosmic Differences #130
Description
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
Activity
mwcraig commentedon Jul 9, 2014
@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 commentedon Jul 10, 2014
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 commentedon Jul 10, 2014
@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. 😃 If you really don't want anyone to use it, then I'd recommend making it a private function.
crawfordsm commentedon Jul 10, 2014
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 thephotutilversion. Thephotutildoes the growing automatically, where it is optionally in theccdprocversion.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_cleanwhereas 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 commentedon Jul 18, 2014
Thanks, @crawfordsm. Yes, the results are better with
cosmicray_clean, which grows the cosmic rays bygboxpixels. The differences in the images appear to be due to how the neighbors ("growing") are handled and the lack of iterations.BTW, the
backgroundparameter incosmicray_cleaninitially confused me before reading the docs because what you really want is a noise/error image. Perhaps a better parameter name isbackground_rmsorerror_imageornoise_image....?crawfordsm commentedon Jul 22, 2014
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 commentedon Jul 23, 2014
Why do you call this
variance_image?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.5Because 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 commentedon Jul 23, 2014
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 commentedon Jul 23, 2014
As long as you are looking, there is a variable
maskfdefined at https://github.com/astropy/ccdproc/blob/master/ccdproc/core.py#L979 that isn't used incosmicray_lacosmic. Is it supposed to be?