Repository navigation
Conversation
…o test Nothing that does anything here yet, just bare minimum to test libasdf detection and compile something with it.
Extend the existing (barely used) FitsReader class to implement ImageFileReader Extend FitsImageSource to also allow selecting FITS HDUs by extname (should probably support extver as well, but this is less common)
This will be useful for creating a DetectionImage from some non-FITS file type; mostly as a placeholder until reading a WCS from ASDF is implemented.
...and prepare for the possibility of non-FITS images.
HDUs (even skipping over non-image HDUs) Consequently make FitsReader::get(int) more sane: it returns the N-th image HDU in the file (0-indexed) rather than FITSish 1-indexing. This enables more consistent use of ImageFileReader::get regardless of the underlying file type.
Ndarray objects out of it. There is a bit of a foot-gun here that AsdfFile has to be wrapped in shared_ptr when using this. I may try to get rid of that limitation later.
…e_2d This will be used to implement AsdfImageSource::getImageTile
classes Drop use of enable_shared_from_this on AsdfFile which was misguided. Should clarify that the lifetime of an `AsdfFile::Ndarray` is concurrent with the `AsdfFile` itself (which is kept alive so long as we have a pointer to its `FileAccessor`)
by the ImageFileReader.get() and iterator interfaces Make this work again for FITS--for ASDF we need more work on the libasdf side.
latest updates to libasdf's development API
(I'm not even exactly sure what this is for yet, but it nicely demonstrates how relatively easy it is now to abstract out the file type)
Needed to add setLayer in the ImageSource base class for convenience; by default this can just be a no-op where it's not applicable.
FITS tests--just their mirror images roughly Header metadata isn't copied over yet though.
tests for FitsImageSource_test
ImageFileReader::readImage on the base class--better name and more general Also fixes handling of image_type in FitsReader::get
This caught (and now fixes) some bugs in how the ImageTile::ImageType was being passed through to AsdfImageSource::getImageTile
theory) which will make it easier to avoid having to write different special cases for FITS and ASDF files specifically (even though they're the only applicable cases for now) For the ASDF case specifically support extracting a FitsWCS. Right now just provides a wrapper around the relevant libasdf GWCS extension types--still need to then plug that into WCSLIB and add tests
supported) and initializing the wcsprm, etc. Adds some tests based on the once for FITS but as noted it's a different WCS with different results. I tested the expected against wcslib and GWCS as well.
this will be useful when matching WCS's to images
path for detection images and reference images this just sets up the wiring--still need to implement parsing and interpretation of the WCS path (and figurout where to document the format)
kinda hacky but gets the job done for now; later may change this when the libasdf gwcs module gets moved out to a separate library/plugin
libasdf-gwcs is now a separate package, which was not the case when I last worked on this, so the build system needs to be updated to check for libasdf-gwcs as well, including AST backend support. This also makes a few fixes to existing code that was affected by libasdf API changes so that it's building and the tests passing again before I start work on GWCS suppport.
The getFitsHeaders method is now moved into its own mixin interface class FitsWcsSerializable. For now this is only applicable to the WCS coordinate system, though there are schemes, to be added later, for serializing GWCS to FITS as well. Also adds previously missing test coverage for writing a FITS image with round-tripping of the WCS headers (as used for check images), and the counter example (simulating ASDF+GWCS) where the WCS cannot be output (yet) to the check image.
The ImageSource base class gains ImageSource::getCoordinateSystem for returning any CoordinateSystem (WCS or some future GWCS) as determined appropriate by the ImageSource. This feels like a better level of abstraction, as different ImageSources may have one or more possible ways to represent a coordinate system, and contain the most knowledge about how to select for it. This also fully decouples WCS.h from specific ImageSource implementations (FitsImageSource or the new AsdfImageSource), keeping the WCS implementation more file format agnostic (modulo the parts that are inherent to FITS data structures). Also deleted the generic WCS(ImageSource&) constructor which was actually a bug: Both DetectionImageConfig.cpp and MeasurementImageConfig.cpp had a bug here were a specific ImageSource instance quietly degraded its type to the ImageSource base class, resulting in quietly producing a bogus identity WCS when passed to `WCS(*image_source)`.
This generalizes the notion of "path to the WCS" across supported image formats. This leads to a little bit of awkwardness like supportsWcsPath, as this is not currently meaningful in FITS. However, it's designed with the possibility in mind that it could be meaningful in FITS, e.g. when adding support for ASDF-in-FITS--this is used I think on some JWST data products that embed a full GWCS as ASDF/YAML embedded in a FITS file, but also ships an approximated FITS WCS. Generalizing it also allows deleting several ASDF-specific code paths. Trying to provide a wcs_path to a file that doesn't know what to do with it results in a meaningful warning message in the log, but is otherwise ignored. Note: setWcsPath is implemented on the ImageReader as opposed to ImageSource, as the latter typically represents a single image in a file (like a single FITS HDU) whereas setting the WCS path is more of a file-level operation for now (makes sense when iterating over an ImageFileReader). That said, I'm not convinced myself that that's really the correct reasoning, and might consider revisiting this later.
…gwcs This removes the old code I had for hoisting a WCSLIB WCS from an ASDF file in the narrow case of a lone fitswcs_imaging--might be worth revisiting later so that code is still in the history, but now most of it isn't needed. The old ASDF->WCS tests are rewritten to use the GWCS instead and agrees well with the previous results evaluated through WCSLIB. Still work to be done on figuring out how to handle inverses, and on thread safety.
The only libasdf-gwcs evaluation backend currently supported, AST, holds an AST FrameSet underlying the evaluation context (asdf_gwcs_eval_t). This is pinned to the thread on which it was created, and cannot be shared safely between threads. This uses the newly introduced `asdf_gwcs_eval_copy` routine, which safely copies this context even between threads. Although this may not be necessary for all future hypothetical backends, it's probably best to do in general. libasdf-gwcs hasn't committed yet to whether or not an `asdf_gwcs_eval_t` is always considered immutable or safe to use across threads, so best to copy it before each use which is relatively cheap.
The per-implementation .test() method works on an istream directly so if the stream is not readable the test simply fails; SE++ reports that the file type is not recognized. In fact it should check if the file exists and report more sensible error messages in this case. Adds some tests for ImageFileReader's interface.
This revealed one shortcoming in the libasdf-gwcs implementation in the handling of bad coordinate values in the inverse evaluation particularly (also affects the forward transform in principle but in practice is more likely in the inverse: AST uses its own convention for representing bad values (a special sentinel value AST__BAD). This is a bit of a leak in the abstraction: should be the same bad value regardless of the valuation backend (e.g. NaN). We check for bad values in the results and raise an InvalidCoordinatesException if so, as this is handled in several places I could find like MeasurementFrameRectangleTask when taking vignettes. Tests demonstrate that GWCS::worldToImage round-trips reasonably accurately with GWCS::imageToWorld, and is safe to call from worker threads.
If a specific image within a file is selected via a `[]` suffix on the filename, e.g. `foo.fits[2]` or `foo.asdf[/data]` that means we are selecting a single (for now) specific image within the file. So an `ImageFileReader` created from a filename with such a suffix should return that image and only that image (if it exists) when calling `ImageFileReader::get` (without arguments) or `ImageFileReader::iter`. This also improves the error handling in cases where an explicitly selected image cannot be found (previously it just reported "No measurement or detection image", confusingly)
The default release build configurations set -ffast-math which drops support for subnormals. When astropy is imported, specifically, it does some import time numpy.finfo operations which warns about this. The warning is probably not very informative for users though so better to avoid it.
This can save a noticeable amount of overhead on startup--negligible for a large mosaic run, say, but not for multiple invocations over many files. Astropy is only imported in model_fitting.get_sky_coord, so unless that is used it is extra weight we don't need.
A bit of a pain because libasdf's logging does not yet allow much formatting customization, so it's necessary to capture and parse its log messages to extract its structure manually. libasdf *does* allow setting a custom `FILE *` stream for the log output, so on glibc we can use fopencookie, or funopen on BSD-alikes. This can be removed later when libasdf implements a custom logging callback. This still works quite nicely though, log messages now arrive through the Elements logger and are formatted along with the rest of the application's log messages, making clear that they are libasdf logs.
Two bugs, both the symptom and the cause. The cause was insufficient NaN comparison in ImageMode::getBackGuess--on debug builds a bare assert(!isnan(mean)) caught mean=NaN, but on release builds this fell through. Because of -ffast-math the NaN was simply allowed to propagate through the statistics. In KappaSigmaBinning::computeBins this NaN propagation also meant that a NaN could poison the whole cell, resulting in ``ndata==0``. This was also asserted against on debug builds, but allowed to propagate on release builds. Uses the fastmath_isnan and fastmath_isinf from SEUtils since these were introduced exactly for this reason, that we still need to sweep for non-finite values even of `-ffast-math` builds which are created by default on release.
…lation This turned out to be a more serious issue affecting this code on nearly flat background cells with close adjacent values, resulting in some caes in catastrophic cancellation. I asked Claude to benchmark different alternative algorithms, namely shifted data and Welford's, against similar data that reproduced the problem, both for speed and accuracy. There was also a comment in the code that suggested it was faster that boost::accumulator. It turns out that actually has two algorithms for calculating the variance: `accumlators::lazy_variance` and `accumulators::variance`. That claim was probably against the "lazy" variant, whereas the non-lazy variant is essentially Weldford's. In any case the claim still seems to hold that this is a bit faster. Additionally I realized that the exception message thrown (replacing the previous bare assert()) was maybe not accurate since from the other fixes we should not be getting non-finite values in this code in the first place. They should already be filtered out, and the exception is almost unreachable code otherwise (with this fix). Leave it in anyways for future catches like this: it's still a more helpful error message than just crashing as before.
Otherwise when read on data cubes it could start on arbitrary/junk slices.
As noted in the comments, the previous behavior was to stop enumeration if at least one candidate IMAGE HDU was found, and the next HDU after it is not an HDU. I think this is a bit of a questionable assumption, but it's consistent with the existing behavior and does work well as a heuristic for many real image observation files. Update the tests to match this assumption and document it better for now. Of course users can still explicitly select an HDU that might be missed by this heurstic, either by index or, as now supported by this work, its EXTNAME.
Instead of just enumerating top-level arrays (a weak heuristic, not true for most real files), implement the TODO for breadth-first enumeration of ndarrays in the tree. It might be nice as a later enhancement to list specific tags which should not be descended into during enumeration (e.g. WCS) as a follow-up, as well as add explicit detection image config options for the image path to use (right now it has to be appended to the filename as a `[]` suffix). In the face of ambiguity though, we now refuse to guess which ndarray is the "detection image" and instead report to the user their options and ask them to pick one explicitly.
This is especially useful for ASDF files, in priniciple also for FITS files. For example if one is primarily working on Roman imaging files they can configure detection-image-path=/roman/data once in their config file and leave it. Need to add the same for other image configurations...
--detection-name-name, --detection-image-wcs-name, etc. I realized that calling this "path", although accurate for ASDF files, is at best confusing or unclear for FITS files, and in both cases is too easily confused with "filesystem path", for which there is already an option. Also added support for a --weight-image-name setting. For flags and measurement images this still needs more thinking, and work on the Python config interface...
This cleans up a lot of mess, but now requires libasdf>=0.3. Worth it though--this is much nicer now.
- As discussed in astroramagh-649 this results in buggy, ill-defined behavior - Can still be enabled via -DFINITE_MATH_ONLY=ON if desired, or for comparison - Deletes the IsNan.h workaround which was not a reliable fix anyways, esp. as observed on recent clang versions
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Follow-up to the discussion I opened at #630 to maybe put a little more concreteness behind it.
This primarily does two things:
If it makes it easier I can also split this up into separate PRs, e.g. one for the general refactoring and one with the ASDF-specific stuff.
Writing ASDF files is not supported yet (as it's not yet supported in libasdf, though will be in a later version). So generalization of file output is not part of this.