-
Notifications
You must be signed in to change notification settings - Fork 104
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
Improvement to spatial flexure #1870
base: develop
Are you sure you want to change the base?
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.
Thanks @debora-pe - Overall this looks good, but it will be good to make sure the dev suite passes. Also, it would be good to check what happens when the slit edges are quite diagonal. I also have a few general questions:
(1) Could you mention the reason why this new algorithm is required?
(2) My main concern is that the spectral smashing might smear out the slit edges. Is there are reason that the ross-correlation needs to be done in 1D instead of 2D (along one dimension)?
(3) I quite like that the cross-correlation is done on the sobel images with this algorithm. Would it also work to use the old algorithm, but instead use the sobel edges? This would avoid the smashing.
sci_sobel, _ = trace.detect_slit_edges(_sciimg, bpm=bpm) | ||
slits_sobel, _ = trace.detect_slit_edges(slitmask, bpm=bpm) | ||
# collapse both sobel images along the spectral direction | ||
sci_smash, _, _ = sigma_clipped_stats(sci_sobel, axis=0, mask=bpm) |
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.
Does this work well in cases where the slit edges are significantly tilted wrt detector columns?
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.
Yes, it does seems to work well. I was worried too at first, but if the slits edges are tilted for the trace image they will be tilted also for the science image.
sci_smash, _, _ = sigma_clipped_stats(sci_sobel, axis=0, mask=bpm) | ||
slits_smash, _, _ = sigma_clipped_stats(slits_sobel, axis=0, mask=bpm) | ||
# remove nan values | ||
sci_smash[np.isnan(sci_smash)] = 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.
Is it possible to remove the nan values from the 2D image, instead (before smashing)?
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 nan values are actually created by sigma_clipped_stats
when encounters the masked values. This is why I added it here.
sci_smash[np.isnan(sci_smash)] = 0 | ||
slits_smash[np.isnan(slits_smash)] = 0 | ||
# invert the negative values | ||
sci_smash[sci_smash < 0] *= -1 |
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.
Is this needed? I think keeping the left edges and right edges as +ve/-ve might actually help the cross-correlation.
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.
good point
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 tested this a lot and it seems that inverting the negative value is still better than keeping them negative. In many cases it does not make much difference, but in some cases the cross-correlation is better with the inversion.
|
||
""" | ||
debug = True if outfile is None else False | ||
# TODO: should we use initial or tweaked slits in this plot? |
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 it's the tweaked slits that you're after. As you noted above on line 85, whatever is used here needs to be consistent with that choice. I think once the tweaked slits are available, they get used.
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.
yes, I believe that is correct.
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.
Well in the code above I'm still using the initial slits (not the tweaked ones) since:
slitmask = slits.slit_img(initial=True, exclude_flag=slits.bitmask.exclude_for_flexure)
It seems to me that the initial ones should work better in the cross-correlation, and therefore this plot should be consistent with that. But I do not want users to get confused when seeing this plot and comparing with what is actually use for reduction, i.e., tweaked slits if tweak_slits=True
|
||
|
||
""" |
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.
Add a check to make sure that vrange is a tuple (otherwise warn and set to None?)
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.
Done. Thanks!
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.
a few minor comments from me.
am really happy to see this code advanced.
sci_smash[np.isnan(sci_smash)] = 0 | ||
slits_smash[np.isnan(slits_smash)] = 0 | ||
# invert the negative values | ||
sci_smash[sci_smash < 0] *= -1 |
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.
good point
#tslits_shift = trace_slits.shift_slits(tslits_dict, lag_max) | ||
# Now translate the tilts | ||
# 2D plot | ||
spat_flexure_qa(sciimg, slits, shift, gpm=np.logical_not(bpm), vrange=qa_vrange) |
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.
remove this line? it is duplicated below.
or is the idea to generate 2 PNGs?
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.
No, this one is used only in debug mode and plot in a matplot window the whole image. No PNG is created. I uses the same script, but it does slightly different things.
|
||
""" | ||
debug = True if outfile is None else False | ||
# TODO: should we use initial or tweaked slits in this plot? |
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.
yes, I believe that is correct.
if np.sum(no_nan) > 0: | ||
if np.any(np.absolute(np.diff(spat_flex[no_nan])) > 0.1): | ||
msgs.warn(f'Spatial flexure is not consistent for all images being combined: {spat_flex}.') | ||
comb_spat_flex = np.round(np.mean(spat_flex[no_nan]),3) |
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.
is there a preference to rounding to the nearest pixel instead of staying a 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.
It is only rounding to the 3rd decimal. It's still a float. But the mean calculation can make it to be a very long number.
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.
@rcooke-ast @profxj I have addressed your comments.
I have extensively tested these changes with all the Keck LRIS data in the Dev Suite and all looks good. I will anyway re-run the whole Dev Suite.
@rcooke-ast Here are my answers to your questions:
(1) Could you mention the reason why this new algorithm is required?
The old algorithm was not working well. A peak in the cross-correlation was often not found and the spatial flexure was not corrected.
(2) My main concern is that the spectral smashing might smear out the slit edges. Is there are reason that the ross-correlation needs to be done in 1D instead of 2D (along one dimension)?
As I mentioned in another comment, I was little worried too, but after a lot of testing (also with more tilted slit edges) the correction looks reasonable. I thought about a 2D cross-correlation but I wasn't able to find a way to do it without adding a lot of runtime (it's probably me that I don't know how to deal with it)
(3) I quite like that the cross-correlation is done on the sobel images with this algorithm. Would it also work to use the old algorithm, but instead use the sobel edges? This would avoid the smashing.
I have also tested a lot using the sobel edges, and the correction is not consistently good.
sci_sobel, _ = trace.detect_slit_edges(_sciimg, bpm=bpm) | ||
slits_sobel, _ = trace.detect_slit_edges(slitmask, bpm=bpm) | ||
# collapse both sobel images along the spectral direction | ||
sci_smash, _, _ = sigma_clipped_stats(sci_sobel, axis=0, mask=bpm) |
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.
Yes, it does seems to work well. I was worried too at first, but if the slits edges are tilted for the trace image they will be tilted also for the science image.
sci_smash[np.isnan(sci_smash)] = 0 | ||
slits_smash[np.isnan(slits_smash)] = 0 | ||
# invert the negative values | ||
sci_smash[sci_smash < 0] *= -1 |
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 tested this a lot and it seems that inverting the negative value is still better than keeping them negative. In many cases it does not make much difference, but in some cases the cross-correlation is better with the inversion.
#tslits_shift = trace_slits.shift_slits(tslits_dict, lag_max) | ||
# Now translate the tilts | ||
# 2D plot | ||
spat_flexure_qa(sciimg, slits, shift, gpm=np.logical_not(bpm), vrange=qa_vrange) |
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.
No, this one is used only in debug mode and plot in a matplot window the whole image. No PNG is created. I uses the same script, but it does slightly different things.
|
||
|
||
""" |
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.
Done. Thanks!
|
||
""" | ||
debug = True if outfile is None else False | ||
# TODO: should we use initial or tweaked slits in this plot? |
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.
Well in the code above I'm still using the initial slits (not the tweaked ones) since:
slitmask = slits.slit_img(initial=True, exclude_flag=slits.bitmask.exclude_for_flexure)
It seems to me that the initial ones should work better in the cross-correlation, and therefore this plot should be consistent with that. But I do not want users to get confused when seeing this plot and comparing with what is actually use for reduction, i.e., tweaked slits if tweak_slits=True
if np.sum(no_nan) > 0: | ||
if np.any(np.absolute(np.diff(spat_flex[no_nan])) > 0.1): | ||
msgs.warn(f'Spatial flexure is not consistent for all images being combined: {spat_flex}.') | ||
comb_spat_flex = np.round(np.mean(spat_flex[no_nan]),3) |
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.
It is only rounding to the 3rd decimal. It's still a float. But the mean calculation can make it to be a very long number.
sci_smash, _, _ = sigma_clipped_stats(sci_sobel, axis=0, mask=bpm) | ||
slits_smash, _, _ = sigma_clipped_stats(slits_sobel, axis=0, mask=bpm) | ||
# remove nan values | ||
sci_smash[np.isnan(sci_smash)] = 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.
The nan values are actually created by sigma_clipped_stats
when encounters the masked values. This is why I added it here.
This work makes improvements to the spatial flexure calculation and add a QA plot to check the flexure.
Documentation was also added.
It addresses issue #1843
Dev-Suite PR pypeit/PypeIt-development-suite#339