Ideas/Wishlist for AstroImages

[sefffal]

Ideas/Wishlist for AstroImages

As discussed on the JuliaAstro call, I'm interested in developing out AstroImages further. The following is a sort of wishlist of what I would want in such a package. I'm aware several of these features like plot recipes already exist in this package, but I've included them below anyways for completeness.
The purpose of all of this is just to start a discussion and summarize what the community would want.

In terms of previous work in this space, there are to my knowledge:

• AstroImages.jl (of course)
• EasyFits
• DirectImages
• CCDReduction (edit: thanks @mileslucas)

Basics

• Very easy to use: loading and displaying images should be fast and self-explanatory
• It should be possible to display very large images
• Seamless support of multi-dimensional cubes (4+ dimensions) and multiple FITS extensions.
• It should be possible to work with multi-extension cubes as easily as normal cubes

Array types

• We should try and be agnostic to the type of array we wrap in an AstroImage.
• For example, we should support views, static arrays, CuArrays, distributed arrays, etc.

Integration with Images.jl ecosystem

• One of the key features of Images is the use of data types like normed UInt8 to store normalized data in the range [0,1] in 8 bits and color types on the elements to trigger displaying images.
• This approach works well but IMO astronomers will expect their data to stay in the format they load from e.g. FITS files: Ints, Float32s, etc. and isn't often normalized.
• Especially since we will want to handle spectra data as e.g. tagged axes rather than an RGB element type (N.B. this doesn't mean we can't have nice functions to compose RGB images)
• I think it's better to follow the approach this package already uses: create views into the data as e.g. Gray when needed and otherwise preserve the element type.
• Ideally we want our images type to work when passed into any of the Images functions and preserve the wrapper type

FITS Headers

• Headers are a balancing act between convenience and performance.
• We need easy, first class access to header data and comments
• Headers should be preserved through function calls that create new images
• IMO headers should always be copied when new AstroImages are created by applying functions to existing images, unless using a view(). Views should of course still have headers but refer to the parent array.
• We need to support large numbers of headers and changing the headers should not lead to heavy compiler latency, so I think headers should not be stored in a fully type stable way e.g. as named tuples.
• At the same time, we want downstream analysis packages to be able to use values stored in headers performantly
• IMO the best compromise is storing header data in an OrderedDict{Symbol,T} where T is a small union of allowable header types. Perhaps something like Union{Bool,Int,Float32,InlineString128}. Hopefully this will allow the compiler to make reasonably fast code using values stored in headers(?). This might take experimentation.
• One question is if headers should be accessible as properties on the AstroImage. For example: img.totexp=5. This is really convenient but we'll still need another syntax for accessing header comments.
• What I settled on in my experiments was the following:
  img[<:Union{Number,Vector}] for indexing into the image
  img[:head1] for accessing a header value
  img[:head1, /] for accessing a comment (but this may be too punny)
  img.head1 as a shortcut syntax for accessing a header value.
  I implemented this syntax in DirectImages after the discussion we had at Compact syntax for getting, setting header values and comments FITSIO.jl#128 and it has felt very natural. Feedback welcome though.

Axes and WCS

• This section trickier for me as I have less experience here.
• For me personally, I need to be able to use images with relative offset axes so that e.g. [0,0] is the "centre" of the image
• How best should this combine with Unitful and UnitfulAngles?

Displaying Images

• display() of an AstroImages should show an image very quickly without needing Plots or Makie.
• We can do this leveraging the Images support for displaying images by making a view of e.g. Gray(). I have a demo where colormaps can also be used to display an image in e.g. magma very efficiently by using RGB().
• We should have a function like imshow that displays images with more customization, e.g. colorschemes, clipping. I prototyped this and found displaying a large image with full color this way was >10x faster than through plots. E.g. 1600x1600 takes about 150ms.
• IMO scaling like z-scale or log should just be functions users apply to the array before passing to display functions instead of e.g. keyword arguments
• We should also have support for plot()ing images of course.
  • In fact, it's more efficient to plot() an array of Gray for example than making a big heatmap. Hopefully we could leverage such an imshow function to make plotting large images efficient automatically.
• plot should be WCS aware vs imshow that just shows the pixel data.
• Ideally we should also support Makie

External Image Viewing Integration

• For interactive viewing, there is no need to reinvent the wheel. We should have lightweight functions for sending images to programs like DS9
• There is already https://github.com/JuliaAstro/SAOImageDS9.jl and ideas for more features here additional features SAOImageDS9.jl#4
• And ds9show from DirectImages.jl
• It might be possible to use shared memory to show large images.
• Is it better to integrate here, or keep this in SAOImageDS9.jl?

[mileslucas]

To add, there's a struct in CCDReduction.jl (CCDData) that has some mixed image/header operations, including broadcasting rules for headers that copies them automatically, although it isn't perfect.

[mileslucas]

I wrote a zscale implementation in PlotUtils.jl

using Plots
img = # ...
heatmap(img, clims=zscale)

or just

clo, chi = zscale(img)

[sefffal]

Thanks @mileslucas! I'll take a look there as well. How does broadcasting work for headers?

[mileslucas]

How does the broadcasting work for headers?

Awkwardly. https://github.com/JuliaAstro/CCDReduction.jl/blob/dc83e241d54604987377360a799b17a3585e68d2/src/ccddata.jl#L90-L102 is the implementation, based off the docs here: https://docs.julialang.org/en/v1/manual/interfaces/#Selecting-an-appropriate-output-array

[mileslucas]

I'm all for the external viewer incorporation. I've really grown to like using DS9 for diagnostic plotting. I have some open ideas I've yet to pursue here- JuliaAstro/SAOImageDS9.jl#4, JuliaAstro/SAOImageDS9.jl#3

[sefffal]

Nice, I hadn't seen the DS9 package!

I'd like to understand more what you mean by broadcasting headers. Do you mean broadcasting values across multiple headers? Or doing something special with headers when broadcasting images?

[mileslucas]

Sorry for the confusion- what I'm doing with CCDData here is copying the headers when broadcasting the image data. So for example if you had a FITS image loaded and you did img .* 2 the result will have the FITS header copied. If you do this with multiple images (e.g., dark subtraction and flat fielding) it will use the header from the first image semantically in the broadcast statement. Let me know if it's still not clear!

[sefffal]

Ah okay thanks! I understand now.

[sefffal]

Here are the main image types defined by the libraries listed above:

EasyFITS.jl:

struct FitsImage{T,N} <: DenseArray{T,N}
    arr::Array{T,N}
    hdr::FitsHeader
end

AstroImages.jl:

struct AstroImage{T<:Real,C<:Color, N, P}
    data::NTuple{N, Matrix{T}}
    minmax::NTuple{N, Tuple{T,T}}
    wcs::NTuple{N, WCSTransform}
    property::Properties{P}
end

CCDReduction.jl:

abstract type AbstractCCDData{T} <: AbstractMatrix{T} end
struct CCDData{T,M<:AbstractMatrix{T}} <: AbstractCCDData{T}
    data::M
    hdr::FITSHeader
end

DirectImages.jl:

struct DirectImage{T,N,A<:AbstractArray} <: AbstractArray{T,N}
    data::A
    headers::OrderedDict{Symbol,CommentedValue}
end

Observations:

• All packages wrap an array
• EasyFITS, CCDReduction, and DirectImages subtype from AbstractArray (DenseArray, AbstractMatrix, and AbstractArray respectively) but AstroImages does not
• EasyFITS, CCDReduction, and DirectImages also store FITS headers
• EasyFITS and CCDReduction directly store a FITSIO.FITSHeader value, DirectImages stores them in an OrderedDict and converts when reading or writing.
• EasyFITS and DirectImages have a quick syntax for accessing headers: A["HEAD"] and A[:HEAD] respectively.
• To access comments, EasyFITS has get(FitsComment, A, "BITPIX") and DirectImages has A[:BITPIX, /]
• AstroImages stores an NTuple of WCSTransform as well as min/max values for display purposes.

Opinions:

• We should subtype from a DenseArray lilke EasyFITS so that AstroImages can be used transparently as arrays. And we should not assume 2D data as 3+D cubes are common.
• We should minimize the coupling between this package and FITSIO by not directly wrapping FITSIO types. FITS files are the most common but e.g. HDF5 files also exist. It should be possible to construct an AstroImage without importing FITSIO.
• I also expect an OrderedDict is more efficient for constant time access of headers compared to storing a struct of vectors lilke FITSHeader does, but I haven't yet benchmarked it for common cases of a few hundred headers.
• I'm on the fence about storing min/max for brightness and contrast adjustments vs. just modifying the data, or specifying limits in display functions but I'm interested to head other opinions

[mileslucas]

EasyFITS and DirectImages have a quick syntax for accessing headers: A["HEAD"] and A[:HEAD] respectively.

So does CCDData :)
https://github.com/JuliaAstro/CCDReduction.jl/blob/dc83e241d54604987377360a799b17a3585e68d2/src/ccddata.jl#L15-L18

[mileslucas]

Why DenseArray versus AbstractArray?

[mileslucas]

I also expect an OrderedDict is more efficient for constant time access of headers compared to storing a struct of vectors lilke FITSHeader does, but I haven't yet benchmarked it for common cases of a few hundred headers.

Why not leave it parametric? There should be equivalent dict-based interfaces for fits headers and ordered dicts, and if not we need to enable that in FITSIO.jl

[sefffal]

EasyFITS and DirectImages have a quick syntax for accessing headers: A["HEAD"] and A[:HEAD] respectively.

So does CCDData :) https://github.com/JuliaAstro/CCDReduction.jl/blob/dc83e241d54604987377360a799b17a3585e68d2/src/ccddata.jl#L15-L18

Thanks, I missed that!

Why DenseArray versus AbstractArray?

I guess you are right, there's no good reason to prevent people using e.g. sparse arrays.

[sefffal]

Why not leave it parametric?

The three implementations I can think of are: vectors (like FTSHeader), Dicts, and named tuples. Using named tuples might improve type stability in some cases since the type of every header value would be known, but at the cost of a lot of compiler latency.
The list based implementation might be faster for small numbers of headers. I'll do some benchmarking to see where the cross over point is.
My view is that we should try and settle on the best all around implementation for headers which would also keep the type signature simpler. One middle ground could be to support NamedTuple(img::AstroImage) returning the headers as a named tuple and without the comments.

[sefffal]

I realized that FITSHeader from FITSIO already uses a Dict behind the scenes to map header keys to indexes, so I doubt there is a performance difference vs using an OrderedDict directly.

[sefffal]

One question I have for the maintainers of FITSIO is about FITSHeader:

mutable struct FITSHeader
    keys::Vector{String}
    values::Vector{Any}
    comments::Vector{String}
    map::Dict{String, Int}
    #  ...
end

Is values a Vector{Any} on purpose so that users can store arbitrary values, or could it be restricted to a small Union like:

values::Vector{Union{Bool,Int,Float32,Nothing,InlineString127}

where InlineString127 comes from InlineStrings.jl.

If we could change that, then code accessing values from FITSHeaders wouldn't be type unstable. In my testing Julia is able to do union splitting with that element type and doesn't box things.

Edit: I believe inline strings could be used because fits string headers are limited to 68 characters (ASCII, I assume?)
Edit 2: I'm not actually seeing an improvement from InlineStrings vs normal strings. The most important factor seems to be the union splitting, not storing the strings inline in the array.

[sefffal]

I did some testing on restricting the range of header types.

HType = Union{Bool,Int,Float64,Nothing,String}
function test1(v,i)
    v[i]*"abc"
end

dat = HType[1,nothing,false,0.1,"abcd"]
@code_typed test1(dat, 5)

Gives:

CodeInfo(
1 ─ %1 = Base.arrayref(true, v, i)::Union{Nothing, Bool, Float64, Int64, String}
│   %2 = (isa)(%1, String)::Bool
└──      goto #3 if not %2
2 ─ %4 = π (%1, String)
│   %5 = invoke Base.string(%4::String, "2"::String)::String
└──      goto #4
3 ─ %7 = (%1 * "abc")::String
└──      goto #4
4 ┄ %9 = φ (#2 => %5, #3 => %7)::String
└──      return %9
) => String

So the compiler is able to generate fast code consisting of only a isa String and a string append.
If we swap *"abc" for *2 the return type of the function changes to Union{Int64,Float64} which is still a big step up from Any.

On the other hand, if we reuse FITSIO.FITSHeader and the type of values is left as Vector{Any} than we get:

CodeInfo(
1 ─ %1 = Base.arrayref(true, v, i)::Any
│   %2 = (%1 * 2)::Any
└──      return %2
) => Any

which is fully generic and always allocates.

IMO it would be reasonable to restrict the headers to only valid FITS header types to get this improvement.
That would mean that downstream packages could operate directly on AstroImages using data from headers without a performance hit.

Edit: I made that modification to FITSIO and everything seems to work well. Floats and Integers get automatically converted to Float64 and Int and unsupported values display a warning and get stored as a string.
I'll do some benchmarking and open a PR to FITSIO.

[giordano]

Regarding the AstroImage struct as implemented here, it's a subtype of Any because this was leaning more towards Images.jl than anything else, and the conversion to Matrix{Color} was done lazily:

AstroImages.jl/src/AstroImages.jl

Lines 169 to 176 in 1b43abc

	# Lazily reinterpret the image as a Matrix{Color}, upon request.
	function render(img::AstroImage{T,C,N}, header_number = 1) where {T,C,N}
	imgmin, imgmax = extrema(img.minmax[header_number])
	# Add one to maximum to work around this issue:
	# https://github.com/JuliaMath/FixedPointNumbers.jl/issues/102
	f = scaleminmax(_float(imgmin), _float(max(imgmax, imgmax + one(T))))
	return colorview(C, f.(_float.(img.data[header_number])))
	end

And related to this, this render function would always call extrema on the whole image array, which would take a lot every time you want to convert to a Matrix{Color}, thus AstroImage has the minmax field purely for performance reasons, not else: #13. But I'm happy for these details to be revised.

[sefffal]

Thanks @giordano! The lazy conversion / colorview makes a lot of sense.

What happens if the user modifies the array? Maybe there is a way to catch this and recalculate min/max on the next render?

[sefffal]

Thinking aloud, if we use code like https://github.com/JuliaAstro/CCDReduction.jl/blob/dc83e241d54604987377360a799b17a3585e68d2/src/ccddata.jl#L12 from CCDReduction, we could add a flag like minmaxdirty to the struct that gets set on each call to setindex!.

Then, minmax could be recalculated on the next call to render.

I'm not sure what the performance impact of that would be though. I don't think the flag would have to atomic. If another thread modified the data between calculating minmax and using the colorview the scale would just be off.

[sefffal]

Support for Multiple HDUs

• EasyFITS supports this by wrapping an open FITS file handle from which users can access individual FITSImage arrays
• AstroImage stores the data, wcs transforms, etc. in NTuples with an index for each HDU
• DirectImages and CCDReduction only wrap a single HDU at a time.

I personally use muti-extension fits files but each HDU is conceptually a different array. For other people's use cases, are there times grouping together multiple HDUs in a single struct would be better than a vector/tuple of images?

[giordano]

What happens if the user modifies the array? Maybe there is a way to catch this and recalculate min/max on the next render?

That's a good question and I hadn't anticipated such need 😕 This is implicitly assuming the data array won't be modified in-place.

but each HDU is conceptually a different array

I agree with that, each extension is in principle independent. However, the reason why AstroImage allows multiple extensions/images is to nicely handle RGB images, where each extension is a different colour channel. However all colour channels need not to come from the same FITS file, they can also come from different files (situation I had quite frequently).

[sefffal]

Ah of course, RGB 😃
That is one of the really nice features of this package.
I agree this is best handled with a wrapper struct that groups mutiple arrays together like AstroImages does.

Maybe this could be split into a separate wrapper type that groups together N AstroImages and calling render on it could compose the different channels.

For full generality, we could support other color spaces e.g. Lab and mappings like SHO/Hubble palette.

[sefffal]

Of course keeping multiple channels in separate arrays/dimensions gets away from the design choices of the Images ecosystem where they are packed into the array elements. I'm not sure what tradeoff is best.

[giordano]

Maybe this could be split into a separate wrapper type that groups together N AstroImages and calling render on it could compose the different channels.

Yup, I'd be OK with having another struct for handling colours.

For full generality, we could support other color spaces e.g. Lab and mappings like SHO/Hubble palette.

Sure. RGB is the only colour mapping this package is currently specialised for, but it isn't limited to that one. Other options would of course be great!

[sefffal]

For accessing header comments, maybe a clearer syntax than img[:HEAD, /] would be img[:HEAD,Comment].
The former looks really nice but dispatching on / (divide) does feel really sketchy.

[sefffal]

If we use img[:HEAD,Comment] then there is also a nice symmetry for accessing extended COMMENT blocks and HISTORY blocks.

For example, img[Comment] could return a list of all COMMENT lines (not all comments on keywords, just standalone comments) and img[History] could do the same for history lines.

And perhaps push!(img[History], "New history") for appending to the history?

[sefffal]

Looking at render, I was able to adapt it to also support arbitrary color maps from ColorSchemes.jl e.g. viridis and magma, all without allocating a new array.
This takes only ~140ms to display a large 1600x1600 image in VSCode.

Images + MappedArrays are very powerful!

Edit: a preview

[sefffal]

For displaying images, I propose we make an explicit function called something like imshow that takes options for colormaps, limits, etc.
This has the benefit of allowing users to visualize other arbitrary arrays even if not wrapped in an AstroImage.

Wrapping in an AstroImage would opt into automatic display at the REPL (and in Jupyter, etc) which could just fall through to imshow with default values, e.g. greyscale colormap, and limits = extrema.

To configure the default display, we could then have a global configuration (a bit like a plot theme but far simpler).
Something like AstroImages.set_cmap!(:magma) and AstroImages.set_clims!(zscale)

[sefffal]

My current plan is to continue discussions in this thread and build consensus.
At the same time, I'm going to start experimenting in a fork and basically merge in the best implementations from AstroImages, EasyFITS, CCDReductions, and DirectImages which are all MIT licensed.

If/once we're all happy with the basics, I'll open a pull request back here for more detailed review.

[sefffal]

Here is a demonstration with more curved coordinates:

julia> plot(WCSGrid(img), aspectratio=1, framestyle=:box, color=:black)

To make this work, I've created a struct called WCSGrid that wraps a WCSTransform, an extent in pixels e.g. the min and maximum of each pixel coordinate, the axes numbers to plot (typically (1,2)), and the pixel coordinate along any remaining axes. Decoupling the grid from the image should set the stage for overplotting different coordinate systems in future.

Once you have created a grid, you can for example get the ticks for where it intersects the axes via the helper wcsticks:

julia> g = WCSGrid(img)
WCSGrid(WCSTransform(naxis=2), (1, 1600, 1, 1600), (1, 2), (1, 1))
julia> wcsticks(g, 1)
([828.5121645296939], ["-2ʰ"])
julia> wcsticks(g, 2)
([1542.4204273078626, 1012.6018518589875, 496.24868471405944], ["183°", "200°", "216°"])

Or vectors of x and y points defining the gridlines using:

julia> xs, ys = wcsgridlines(g)

I'm still not sure how to handle laying out the grids for images where the coordinates wrap around part way through.

[ptiede]

This looks really great! One question I have the relates to my work is can we make this deal with polarization (Stokes I, Q, U, V and maybe circular polarization as well) in a nice way? This will be important for radio astronomy where the polarization contains a ton of information about the emitting electron populations and can e.g. be used to differentiate between different types of accretion.

You probably also want to change plotting for polarized images. Usually, linear polarization is denoted with ticks overlayed on the stokes I emission.

[sefffal]

Thanks @ptiede! I hope we can. Is that information stored in the WCS headers?
If you have a sample image I could test with, that would help a lot.

[ptiede]

I don't think this will fit in the header, but I am not too familiar with them. I usually think of polarization as being something akin to color, so each image pixel is really four numbers. How I usually dealt with this in the past was to use a StructArray with a custom Stokes vector type that contains I,Q,U,V. So I think to make polarization work you just need to ensure that the images can be stored in some dimension.

I'll try to grab a polarized CLEAN image and send it your way!

[sefffal]

Here are some updates showing different coordinate systems.

Cartesian projected galactic coodinates with rectilinear grid

plot(mws, grid=true, gridcolor=:white, background=:black, framestyle=:box)

Mollweide projected galactic coodinates

plot(ovro, grid=true, gridcolor=:cyan, framestyle=:box, background="#222")

In this case, the x-ticks have gotten all bunched up at the bottom. There's really not much we can do about this other than adding a heuristic to remove plot ticks very close together.

For this reason, I'm keeping the option to directly annotate the grid though it needs a bit more work:

plot(ovro, grid=true, gridcolor=:cyan, framestyle=:box, background="#222", annotategrid=true, ticks=[])

[sefffal]

Summarizing something from the most recent JuliaAstro call, here's the problem I'm currently thinking through.

If a user has some AstroImage, when they slice into it we need to keep track of which index is which from the original parent array, because it is those indices that match up with the WCSTransform (and the FITS headers).

Consider:

X = AstroImage(randn(10,10))
v = X[5:10, :]
plot(v)

We need to know that the first index of v actually corresponds to the 5th element of X when we go to plot the WCS axes.
Besides visualization, this is also important for general image transformations and using e.g. Reproject.jl.

For now what I've done is that indexing into an AstroImage returns a new AstroImage wrapping an OffsetArray. For example:

img = AstroImage("656.fits")
axes(img) # (1:1600,1:1600)
axes( img[800:end, 600:end] ) # (600:1600, 800:1600)

This way Reproject and plotting works smoothly:

plot(img, grid=true, gridcolor=:cyan)
plot(img[800:end, :], grid=true, gridcolor=:cyan)

Note how even though the user took a subset of the first coordinate, the labels along the second coordinate had to change as well?

More dimensions get trickier like when a user selects a slice from a 3D (x, y, wavelength) cube:

A = AstroImage(randn(100,100,10))
b = A[:,:, 10]
c = A[:, 50, :]

We need to know which index b has along the third, "lost" dimension. One reason is because WCS coordinates can be skewed along the third dimension such that for example, moving along spectral index also shifts x and y positions.
We also need to know that c, now a 2D array actually refers to the first and third WCS axis. That way if the user runs plot(c) they get a heatmap of intensity along RA and wavelength.

One way we could do this is to make AstroImage keep track of these indices itself, or we could maybe modify OffsetArrays for our purposes.

astropy completely punts on this and makes the user specify it in their calls to matplotlib, but I think we can do better and make this work really smoothly.

Naturally, getindex with a single element (scalar) shouldn't use any of this machinery and just return a single value. I guess selecting a vector from a cube should keep this info so that you can do e.g.
v = A[50,50,:]; plot(v) and get an automatically labelled spectral axis.

[sefffal]

As Miles has pointed out, some of this machinery already exists when using views. For instance:

julia> A = randn(100,100,10)
julia> v = view(A, :, :, 5)
julia> parentindices(v)
(Base.Slice(Base.OneTo(100)), Base.Slice(Base.OneTo(100)), 5)

So maybe there is a way to re-use this for general slicing (which copies) in addition to creating views.

[sefffal]

Another option is to do this ourselves and store the following fields in the AstroImage struct:

• offset: a tuple of offsets to the indices of the array
• stride: a tuple of strides for when the user does A[1:5:end,:]
• parentaxes: a tuple of indices into the parents axes e.g. (1,3) if they do A[:,5,:]

This is sort of creating our own OffsetArrays, and the OffsetArrays.jl code is quite daunting.

[mileslucas]

@sefffal

Here is an example from imview on the left and plot on the right (with the PyPlot backend)

What are your thoughts on editing the recipe to add

xlim --> extrema(first(axes(image))
ylim --> extrema(second(axes(image))

(or, however you get the indices)

This would eliminate the gridded whitespace in the Plots.jl images. Here's an example of what a plot looks like with these settings: https://juliaastro.github.io/PSFModels.jl/dev/examples/#Fitting-a-PSF

[sefffal]

What are your thoughts on editing the recipe to add...

Great point, yes I think we should have that in the recipe. This is doubly important when world coordinates are in play since they are calculated along the axes of the image. If the plot adds extra white space they can become incorrect.
PSFModels looks really nice BTW.

[sefffal]

There was also some discussion about if images with "flipped" axes e.g. DEC, RA exist in the wild. Here's an image where the rotation matrix almost does this. If you look carefully, the tick labels are where the actual RA and & DEC grid lines intersect the axes.

We're doing the naive thing here and always labelling the x-axis as axis 1 and the y-axis as axis 2, but maybe some future work could do something smarter.

[sefffal]

I've been working steadily on this and am making some good progress.
The biggest features I've implemented are:

• Making AstroImage inherit from DimensionalData's AbstractDimArray instead of wrapping an OffsetArray.
  After comparing a lot of packages, I found this one is pretty much ideal for our use-case. It's well tested and robust and does a lot of heavy lifting for us.
• WCS: pix_to_world(img, coords) and world_to_pix(img, coords) that works correctly even on subarrays or slices of e.g. spectral cubes by automatically handling offsets, hidden dimensions, and axis step sizes.
• Support for many more FITS files with WCS info through better serializing FITS headers in a way that WCSLib can parse.
• implot function for plotting arrays or AstroImages that now tracks WCS coordinates robustly.
• FileIO load and save support for images and tabular data.

I also removed support for just plot(img) in favour of implot(img). In using this package for my work, I realized there are a lot of other plots like histograms we might want to make out of images. Swapping your data to an AstroImage should have no effect on any visualizations you already have.

I'm going to make a few demo notebooks to showcase everything, but here are a few previews.

Loading a cube with axes galactic longitude, galactic latitude, velocity:

fname = download("http://data.astropy.org/tutorials/FITS-cubes/reduced_TAN_C14.fits")
h = load(fname)
implot(h[Z=300], cmap=:turbo)
# Or just implot(h[:,:,300], cmap=:turbo)

Plotting one spectrum extracted from the same cube (this recipe just comes from DimensionalData)

plot(h[X=145,Y=130])

Next up is automatically labeling these axes with WCS axis names and units.

Should have lots to show by the next JuliaAstro call!

[sefffal]

I don't think this will fit in the header, but I am not too familiar with them. I usually think of polarization as being something akin to color, so each image pixel is really four numbers. How I usually dealt with this in the past was to use a StructArray with a custom Stokes vector type that contains I,Q,U,V. So I think to make polarization work you just need to ensure that the images can be stored in some dimension.

I'll try to grab a polarized CLEAN image and send it your way!

Hi @ptiede , would love to test out one of your FITS files!

[mileslucas]

I'm also interested in polarization plotting. My concern is I don't believe there is a standard for arraying polarized data. Perhaps there are standards in the radio community, but I don't know of any in the direct imaging community. I think this is better suited for some different recipes/functions s.t. we aren't tied down to a specific data layout.

For example, given stokes Q and U images, we could plot the polarized intensity (sqrt(Q^2 + U^2)) with the angle of linear polarization (0.5atan(Q, U)) as quivers. There's also the azimuthal stokes parameters, but that's a little bit more complicated (see here, section 4.2), and technically could involve some numerical optimization, which is out of scope for this package. I don't even know if that's a technique used in radio, also I don't know what people do with circular polarization 😅 .

I'd be happy to give you some of my FITS images if you want to play around with them!

[sefffal]

My concern is I don't believe there is a standard for arraying polarized data

I have no experience here, but I did see some allowances for polarization axes in FITS world coordinates (https://www.aanda.org/articles/aa/pdf/2002/45/aah3859.pdf , table 7 for example). Does that have any relevance?

I'd be happy to give you some of my FITS images if you want to play around with them!

Yes please!

[ptiede]

So for linear polarization, the EHT has a few ways to represent images. The one we presented in EHTC VII are of the form:

The tick position angle is the EVPA which is 1/2atan(U/Q), the length of the tick is the absolute linear polarization in Jy (sqrt(Q^2+U^2)) and the color denotes the linear polarization fraction (sqrt(Q^2+U^2)/I). The greyscale image in the background is the total stokes I image. Now this is just for EHT and conventions are a bit of a mess but I think this is a reasonable starting point.

For circular polarization, one thing I have seen is to replace the ticks with little ellipses and then change the ellipse's size depending on the circular polarization. For why circular polarization is interesting, well for the EHT it encodes interesting magnetic field properties around black holes. A recent and quite nice reference for this is https://arxiv.org/pdf/2104.11301.pdf

[sefffal]

I'm working right now on colorbar support. The colormappings we provide through imview and implot are more sophisticated than what Plots.jl can render with it's native colorbars, so I'm implementing them as a second (optional) image series off to the side.

The two options I can think of for displaying a stretched colorbar are to stretch the tick locations (1st image), or stretch the color scheme (2nd image).

implot(img, stretch=logstretch, clims=extrema, cmap=:magma)

What are your preferences?

Edit: The third option would be to do what DS9 does and place tick marks at equal separations along the colorbar, but show the actual pixel level at that location:

My personal preference is option 1.

[sefffal]

Currently implot does not work with the Plotly.jl backend. It has no support for image series, and falls back to treating it as a heatmap.
But I see that in some cases it's possible to use Plotly with images: https://plotly.com/python/images/#zoom-on-static-images
Where is the best place to open a feature request? On the main Plots.jl repo?

[mileslucas]

Where is the best place to open a feature request? On the main Plots.jl repo?

Yes, it seems like this is handled by the logic here.

[mileslucas]

but I did see some allowances for polarization axes in FITS world coordinates (aanda.org/articles/aa/pdf/2002/45/aah3859.pdf , table 7 for example). Does that have any relevance?

I did not know about this, thanks for pointing it out! That being said, I've never used WCS for any high-contrast work, and I haven't personally seen anyone who does. So I wouldn't rely on it.

My personal preference is option 1.

Me too.

[sefffal]

I did not know about this, thanks for pointing it out! That being said, I've never used WCS for any high-contrast work, and I haven't personally seen anyone who does. So I wouldn't rely on it.

Many high contrast imaging instruments I’ve used do have WCS headers attached to their raw data (GPI, NIRC2, HST have them, not sure about VLT-SPHERE or LBT LMIRCam) but I haven’t ever looked at polarization data from any of them in my work. I’ll track down some GPI cubes and see if they have the proper info in the headers.

It’s a good point that we need to make sure all functionality doesn’t assume the data comes from FITS files with well formed headers. Right now the dimensions are labelled as X,Y,Z,Ax{4},Ax{5}, etc automatically.
A clean way solution might be to let users pick from certain predefined dimension labels with special meanings. This way we can mark any axis of a cube as being a polarization axis, a spectral axis, a time axis, etc.
And they could fallback to something like WCS{1}, WCS{2}, … to mean “refer to the header to find out what kind of axis”

[sefffal]

Where is the best place to open a feature request? On the main Plots.jl repo?

Yes, it seems like this is handled by the logic here.

Thanks!

[sefffal]

After your comments @mileslucas and @ptiede, I think I have a solution for files without WCS info where we still want to give certain axes labels and/or special meaning, e.g. polarization, spectral, or time axes.

I've made the AstroImage type accept arbitrary dimensions again from DimensionalData. For instance:

img = AstroImage( randn(101, 101), (-50:50, -50:50)) # Or:
img = AstroImage( randn(101, 101), (X=-50:50, Y=-50:50))

constructs an AstroImage with dimensions centred around zero.

Another example from a real image:

We can use this for arbitrary cubes even if the WCS info isn't set:

cube = AstroImage( randn(15,15,100,45),  (X, Y, Ti, Spec))
# Extracting a spectrum:
p1 = plot(cube[X=1,Y=1,Ti=10,Spec=1:45])
# Plotting a wavelength slice:
p2 = implot(cube[Ti=10,Spec=35])
plot(p1, p2, layout=(2,1), size=(350,700))

Finally, we can use categorical axes to store polarization data in a nice way (edit: or rather, describe how the polarization data is stored in an existing cube)

polcube = AstroImage( randn(15,15,3),  (X=1:15, Y=1:15, Pol=[:I, :Q, :U])) 
# Polarization plotting TBD

To restore the previous behaviour of matching the dimensions to the WCS axes, we can pass wcsdims=true when loading a FITS file. This automatically sets the image dimensions to Wcs{1}, Wcs{2}, ... instead of X, Y, ... which tells downstream code to refer to the appropriate WCS axis.

So now the examples from earlier in this thread look like:

img = load("fits/656nmos.fits", wcsdims=true)
implot(img)

[icweaver]

oop, think my finger slipped, sorry!