New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add crosstalk to ISR #40
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'd like to see the crosstalk coefficients moved out of the config because I don't understand how this works in the situation where the coefficients depend on sensor. I'd much rather the task be able to get the crosstalk by dataRef.
python/lsst/ip/isr/crosstalk.py
Outdated
minPixelToMask = Field(dtype=float, default=45000, | ||
doc="Set crosstalk mask plane for pixels over this value") | ||
crosstalkMaskPlane = Field(dtype=str, default="CROSSTALK", doc="Name for crosstalk mask plane") | ||
coeffs = ListField( |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I don't think the crosstalk coefficients can be in the config. This means you can only run one task per chip. It needs to be configurable from the dataRef.
python/lsst/ip/isr/crosstalk.py
Outdated
dtype=float, | ||
doc=("Flattened crosstalk coefficient matrix; should have nAmps x nAmps entries. " | ||
"Once 'reshape'-ed, ``coeffs[i][j]`` is the fraction of the j-th amp present on the i-th amp."), | ||
default=[0, 0, 0, 0, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This bakes in an assumption of a 2x2 configuration. I don't think having a default is very useful here.
"""Apply intra-CCD crosstalk correction""" | ||
ConfigClass = CrosstalkConfig | ||
|
||
def run(self, exposure): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think that since the coefficients can't be in the config, this signature needs to take either a dataRef
or an opaque object that contains the details of the crosstalk.
python/lsst/ip/isr/crosstalk.py
Outdated
""" | ||
mi = exposure.getMaskedImage() | ||
mask = mi.getMask() | ||
bg = calculateBackground(mi, badPixels) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Since this is before flatfield, shouldn't the background be calculated per amp?
python/lsst/ip/isr/crosstalk.py
Outdated
|
||
# Set crosstalkStr bit only for those pixels that have been significantly modified | ||
# (not necessarily those that are bright originally). | ||
mask.clearMaskPlane(crosstalkPlane) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm missing something here. Doesn't this clear the mask plane? Where do we set the crosstalk bit?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm curious about @SimonKrughoff 's comment here too.
5f175e8
to
59dc492
Compare
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Some comments on the commits/files you requested I review.
I'm not at all keen on the use of list of list of list...
: is there a better way to do that? I fear that 3-way nested listcomps are going to be measurably slow.
python/lsst/ip/isr/crosstalk.py
Outdated
|
||
# Set crosstalkStr bit only for those pixels that have been significantly modified | ||
# (not necessarily those that are bright originally). | ||
mask.clearMaskPlane(crosstalkPlane) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm curious about @SimonKrughoff 's comment here too.
|
||
Returns | ||
------- | ||
ratios : `list` of `list` of `list` of `float` |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
list of list of list? Isn't there a better way to represent this?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I can't think of a better description. As mentioned in the docstring, it's a matrix of lists of floats.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Modified to "list of list of ndarray", with an explanation that it is a matrix of arrays.
if ii == jj: | ||
continue | ||
jImage = extractAmp(mi.getImage(), jAmp, iAmp.getReadoutCorner()) | ||
pixels = (jImage.getArray()[select] - bg)/iImage.getImage().getArray()[select] |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
How does this select
apply to jImage
, since it was created on iMask
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The amps are assumed to be the same size, and I want to select the same pixels on jImage
as on iImage
.
lo, med, hi = np.percentile(values, [25.0, 50.0, 75.0]) | ||
sigma = 0.741*(hi - lo) | ||
good = np.abs(values - med) < rejSigma*sigma | ||
print(ii, jj, rej, med, sigma, len(values), good.sum()) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
print statement here isn't going to be helpful for debugging: please remove it, or turn it into a log.debug()
, with labels for the various bits.
coeff : `numpy.ndarray` | ||
Crosstalk coefficients. | ||
coeffErr : `numpy.ndarray` | ||
Crosstalk coefficient errors. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Mention that these are calculated as the std-dev of the mean?
values = np.array(ratios[ii][jj]) | ||
values = values[np.abs(values) < 1.0] # Discard unreasonable values | ||
for rej in range(rejIter): | ||
lo, med, hi = np.percentile(values, [25.0, 50.0, 75.0]) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Should this be using afw.stats
or some other standard stats package instead?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
That would be best, but we don't have these robust stats in afw (yet?).
|
||
@classmethod | ||
def parseAndRun(cls, *args, **kwargs): | ||
"""Implement scatter/gather""" |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Should probably document the return value here, since it's a tuple of things.
ebb9bf7
to
433fb85
Compare
Adapted from an implementation long used in obs_subaru for Suprime-Cam and Hyper Suprime-Cam; this was generalised to allow for an arbitrary number of amplifiers in arbitrary orientation, code cleaned up and docs updated. We provide both function and Task interfaces to both measure and remove the crosstalk.
Using the new CrosstalkTask. Crosstalk correction is disabled by default since this is new functionality that requires measured coefficients.
Tests both crosstalk measurement and removal, using both the function and Task interfaces.
433fb85
to
220888b
Compare
No description provided.