JuliaAstro / AstroImages.jl

Visualization of astronomical images
Other
29 stars 9 forks source link

New Features and Improvements to AstroImages #30

Closed sefffal closed 2 years ago

sefffal commented 2 years ago

History: https://github.com/JuliaAstro/AstroImages.jl/issues/29

Status: T-2 weeks until ready to review & merge

This PR is an overhaul of AstroImages to address the following issues:

It implements a new AstroImage struct based on DimensionalData.jl that:

It implements a new imview function for rendering arrays to static images that supports:

It also implements plot recipes that:

Still TODO:

Known issues:

Gallery

imview(img; clims=percent(99), cmap=:seaborn_rocket_gradient)
implot(img, contrast=2, stretch=asinhstretch)

More:

sefffal commented 2 years ago

I realized this morning that this implementation assumes there is only one WCS system. Multiple systems are common so I think I'll need to change the wcs field into a collection the way it was in the original implementation.

timholy commented 2 years ago

Also, the large-image issue with ImageView is fixed in v0.11 (https://github.com/JuliaImages/ImageView.jl/pull/265).

sefffal commented 2 years ago

Hi @timholy, thanks for your comments.

Main thing might be to consider whether some of this should either move to JuliaImages (parts are not very astro-specific) or whether JuliaImages needs to enhance support for some features you might need.

It occurs to me that the imview function is not astronomy specific and could potentially be moved into one of the more general images packages, if you’re interested. It has fairly light dependencies: mostly just mappedarrays, colors, colorschemes, and plotutils.

It seems to mix well with the rest of the Images ecosystem. I’ve had a lot of fun using Images functions to manipulate arrays once they’ve been rendered to RGBA via imview, eg. creating mosaics with MosaicViews, saving animated GIFs, etc.

On the other side of things, I’ve also hit a few snags. One is preserving the AstroImage wrapper, for example through warp or filter operations. For now I’ve just been adding methods in contrib/images.jl as I go. Do you have any suggestions?

timholy commented 2 years ago

If I understand correctly, imview is partly like https://github.com/JuliaImages/ImageCore.jl/blob/master/src/map.jl on steroids. I would be very happy to have further elementwise transformations in ImageCore. I guess the question is whether it's possible to factor just the contrast-generation out?

sefffal commented 2 years ago

If I understand correctly, imview is partly like https://github.com/JuliaImages/ImageCore.jl/blob/master/src/map.jl on steroids. I would be very happy to have further elementwise transformations in ImageCore. I guess the question is whether it's possible to factor just the contrast-generation out?

It absolutely would be. Currently if you pass cmap=nothing you get a grayscale image in return. The pixel values are still RGBA though since I have missing and NaN mapped to transparent.

sefffal commented 2 years ago

Update on creating RGB composites

I've created the function composecolors which aims to be a more general way to create color composites. It supports creating RGB images of course, but it can also work with any number of channels. Like implot, it leverages imshow so that as much as possible goes through a single code path.

Basic usage:

composecolors([red, green, blue])

A more complicated example of mapping red, green, blue, and hydrogen alpha channels into a single image:

composecolors([antred, antgreen, antblue, anthalp], ["red", "green", "blue", "maroon1"], multiplier=[1,2,1,1])

image

It can also use colorshemes for each channel which IMO works well for combining instruments. Here's a radio and xray composite of Hercules A:

composecolors([chandra, vla], [:ice, :magma])

image

Unlike ccd2rgb it accepts any abstractarray (or AstroImage of course) as input, but not a FITSHdu or filepath. It also doesn't apply reproject automatically since this is quite expensive. In my experience, one wants to tweak the color mappings a few times to get the best result, so IMO better to reproject/align the images first as a separate step.

Edit: stretches are supported with the stretch keyword argument which can be a single function, or a list of functions to pair with each image. All the arguments of imview are supported in this way (e.g. contrast, bias).

Edit 2: this view is created lazily so you can update the underlying data and the image will recompute. This is also good for very large images since viewing it doesn’t create a copy (Unless the display machinery does, but we can’t control that).

rafaqz commented 2 years ago

Great to see this use of DimensionalData.jl, pretty much how I hoped it could be plugged in. One comment: dont export Between. It clashes with DataFrames.jl and I'm hoping to remove it in a future release, in favour of using .. intervals.

sefffal commented 2 years ago

Some comments about CI

Thanks @giordano I've made your suggested CI changes.

On the topic of docs, I would like the docs to have lots of examples showing images. Do you know a way to have images created in code blocks get displayed in the docs pages automatically?

Right now I have to run the examples locally, save the images, and then embed them manually. It would be nice if they could be kept in sync automatically.

sefffal commented 2 years ago

Note that I had suggested using the julia-actions/cache action for caching 🙂

Thanks, fixed!

giordano commented 2 years ago

On the topic of docs, I would like the docs to have lots of examples showing images. Do you know a way to have images created in code blocks get displayed in the docs pages automatically?

I don't know if Documenter.jl has a way to automatically output an image, but in docs/make.jl you can run all the code you want. Including automatically generating images.

sefffal commented 2 years ago

Looks like DemoCards.jl is used by a lot of packages for this sort of thing (AlgebraOfGraphics, Images, Plots, etc).

sefffal commented 2 years ago

This sort of works with ImageInTerminal.jl which is neat!

I used collect since the type signatures can get pretty gnarly:

julia> typeof(imview(img))
AstroImageMat{RGBA{N0f8}, Tuple{X{Sampled{Int64, OneTo{Int64}, ForwardOrdered, Regular{Int64}, Points, NoMetadata}}, Y{Sampled{Int64, OneTo{Int64}, ForwardOrdered, Regular{Int64}, Points, NoMetadata}}}, Tuple{}, ReadonlyMultiMappedArray{RGBA{N0f8}, 2, Tuple{SubArray{Float32, 2, Matrix{Float32}, Tuple{StepRange{Int64, Int64}, Slice{OneTo{Int64}}}, false}, ReadonlyMappedArray{Any, 2, SubArray{Float32, 2, Matrix{Float32}, Tuple{StepRange{Int64, Int64}, Slice{OneTo{Int64}}}, false}, var"#129#130"{Float64}}}, var"#colormap#80"{typeof(identity), ColorScheme{Vector{RGB{Float64}}, String, String}, Float64, Float64}}, Tuple{X{Sampled{Int64, OneTo{Int64}, ForwardOrdered, Regular{Int64}, Points, NoMetadata}}, Y{Sampled{Int64, OneTo{Int64}, ForwardOrdered, Regular{Int64}, Points, NoMetadata}}}} (alias for AstroImage{ColorTypes.RGBA{FixedPointNumbers.Normed{UInt8, 8}}, 2, Tuple{X{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}, Y{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, Tuple{}, MappedArrays.ReadonlyMultiMappedArray{ColorTypes.RGBA{FixedPointNumbers.Normed{UInt8, 8}}, 2, Tuple{SubArray{Float32, 2, Array{Float32, 2}, Tuple{StepRange{Int64, Int64}, Base.Slice{Base.OneTo{Int64}}}, false}, MappedArrays.ReadonlyMappedArray{Any, 2, SubArray{Float32, 2, Array{Float32, 2}, Tuple{StepRange{Int64, Int64}, Base.Slice{Base.OneTo{Int64}}}, false}, AstroImages.var"#129#130"{Float64}}}, AstroImages.var"#colormap#80"{typeof(identity), ColorSchemes.ColorScheme{Array{ColorTypes.RGB{Float64}, 1}, String, String}, Float64, Float64}}, Tuple{X{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}, Y{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}})

Any suggestions on how to cut this down? I doubt that level of detail is interesting to anyone.

sefffal commented 2 years ago

@timholy and @rafaqz I'm hitting a number of errors when I try to use functions like imfilter or imrotate on AstroImages or any DimArray.

For example with imfilter:

julia> using DimensionalData, ImageTransformations, ImageFiltering
julia> da = DimArray(randn(10,10), (X,Y))
...
julia> imfilter(da, Kernel.gaussian(1.0))
ERROR: ArgumentError: Invalid index `CartesianIndex(1,)`. Did you mean `At(CartesianIndex(1,))`? Use stardard indices, `Selector`s, or `Val` for compile-time `At`.
Stacktrace:
  [1] selectindices
    @ C:\Users\William\.julia\packages\DimensionalData\jl2xf\src\LookupArrays\selector.jl:637 [inlined]
  [2] (::DimensionalData.Dimensions.LookupArrays.var"#37#38")(l::DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}, s::CartesianIndex{1})
    @ DimensionalData.Dimensions.LookupArrays C:\Users\William\.julia\packages\DimensionalData\jl2xf\src\LookupArrays\selector.jl:630
  [3] map
    @ .\tuple.jl:247 [inlined]
  [4] selectindices
    @ C:\Users\William\.julia\packages\DimensionalData\jl2xf\src\LookupArrays\selector.jl:629 [inlined]
  [5] dims2indices
    @ C:\Users\William\.julia\packages\DimensionalData\jl2xf\src\Dimensions\indexing.jl:17 [inlined]
  [6] dims2indices
    @ C:\Users\William\.julia\packages\DimensionalData\jl2xf\src\Dimensions\indexing.jl:11 [inlined]
  [7] getindex
    @ C:\Users\William\.julia\packages\DimensionalData\jl2xf\src\array\indexing.jl:16 [inlined]
  [8] copydata!
    @ C:\Users\William\.julia\packages\ImageFiltering\cs86y\src\border.jl:679 [inlined]
  [9] copydata!
    @ C:\Users\William\.julia\packages\ImageFiltering\cs86y\src\border.jl:686 [inlined]
 [10] _copy!(dst::OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}, src::BorderArray{Float64, 2, DimArray{Float64, 2, Tuple{X{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}, Y{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}}, Tuple{}, Matrix{Float64}, DimensionalData.NoName, DimensionalData.Dimensions.LookupArrays.NoMetadata}, Pad{2}}, #unused#::Pad{2})
    @ ImageFiltering C:\Users\William\.julia\packages\ImageFiltering\cs86y\src\borderarray.jl:177
 [11] copy!
    @ C:\Users\William\.julia\packages\ImageFiltering\cs86y\src\borderarray.jl:163 [inlined]
 [12] padarray(#unused#::Type{Float64}, img::DimArray{Float64, 2, Tuple{X{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}, Y{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}}, Tuple{}, Matrix{Float64}, DimensionalData.NoName, DimensionalData.Dimensions.LookupArrays.NoMetadata}, border::Pad{2})
    @ ImageFiltering C:\Users\William\.julia\packages\ImageFiltering\cs86y\src\border.jl:664
 [13] imfilter!
    @ C:\Users\William\.julia\packages\ImageFiltering\cs86y\src\imfilter.jl:712 [inlined]
 [14] imfilter!(r::ComputationalResources.CPUThreads{ImageFiltering.Algorithm.FIRTiled{2}}, out::DimArray{Float64, 2, Tuple{X{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}, Y{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}}, Tuple{}, Matrix{Float64}, DimensionalData.NoName, DimensionalData.Dimensions.LookupArrays.NoMetadata}, img::DimArray{Float64, 2, Tuple{X{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}, Y{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}}, Tuple{}, Matrix{Float64}, DimensionalData.NoName, DimensionalData.Dimensions.LookupArrays.NoMetadata}, kernel::Tuple{OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}, OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}}, border::Pad{0})
    @ ImageFiltering C:\Users\William\.julia\packages\ImageFiltering\cs86y\src\imfilter.jl:703
 [15] imfilter!(out::DimArray{Float64, 2, Tuple{X{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}, Y{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}}, Tuple{}, Matrix{Float64}, DimensionalData.NoName, DimensionalData.Dimensions.LookupArrays.NoMetadata}, img::DimArray{Float64, 2, Tuple{X{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}, Y{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}}, Tuple{}, Matrix{Float64}, DimensionalData.NoName, DimensionalData.Dimensions.LookupArrays.NoMetadata}, kernel::Tuple{OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}, OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}}, border::Pad{0}, alg::ImageFiltering.Algorithm.FIRTiled{2})
    @ ImageFiltering C:\Users\William\.julia\packages\ImageFiltering\cs86y\src\imfilter.jl:612
 [16] imfilter!
    @ C:\Users\William\.julia\packages\ImageFiltering\cs86y\src\imfilter.jl:606 [inlined]
 [17] imfilter
    @ C:\Users\William\.julia\packages\ImageFiltering\cs86y\src\imfilter.jl:27 [inlined]
 [18] imfilter(::Type{Float64}, ::DimArray{Float64, 2, Tuple{X{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}, Y{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}}, Tuple{}, Matrix{Float64}, DimensionalData.NoName, DimensionalData.Dimensions.LookupArrays.NoMetadata}, ::Tuple{OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}, OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}}, ::String)
    @ ImageFiltering C:\Users\William\.julia\packages\ImageFiltering\cs86y\src\imfilter.jl:22
 [19] imfilter
    @ C:\Users\William\.julia\packages\ImageFiltering\cs86y\src\imfilter.jl:18 [inlined]
 [20] imfilter
    @ C:\Users\William\.julia\packages\ImageFiltering\cs86y\src\imfilter.jl:10 [inlined]
 [21] imfilter(::DimArray{Float64, 2, Tuple{X{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}, Y{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}}, Tuple{}, Matrix{Float64}, DimensionalData.NoName, DimensionalData.Dimensions.LookupArrays.NoMetadata}, ::OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}})
    @ ImageFiltering C:\Users\William\.julia\packages\ImageFiltering\cs86y\src\imfilter.jl:5
 [22] top-level scope
    @ REPL[34]:1

And imrotate:

julia> imrotate(da, 0.4)
ERROR: ArgumentError: Invalid index `Interpolations.WeightedAdjIndex{2, Float64}(9, (0.9766076262980938, 0.023392373701906166))`. Did you mean `At(Interpolations.WeightedAdjIndex{2, Float64}(9, (0.9766076262980938, 0.023392373701906166)))`? Use stardard indices, `Selector`s, or `Val` for compile-time `At`.
Stacktrace:
  [1] selectindices
    @ C:\Users\William\.julia\packages\DimensionalData\jl2xf\src\LookupArrays\selector.jl:637 [inlined]
  [2] (::DimensionalData.Dimensions.LookupArrays.var"#37#38")(l::DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}, s::Interpolations.WeightedAdjIndex{2, Float64})
    @ DimensionalData.Dimensions.LookupArrays C:\Users\William\.julia\packages\DimensionalData\jl2xf\src\LookupArrays\selector.jl:630
  [3] map
    @ .\tuple.jl:247 [inlined]
  [4] selectindices
    @ C:\Users\William\.julia\packages\DimensionalData\jl2xf\src\LookupArrays\selector.jl:629 [inlined]
  [5] dims2indices
    @ C:\Users\William\.julia\packages\DimensionalData\jl2xf\src\Dimensions\indexing.jl:17 [inlined]
  [6] dims2indices
    @ C:\Users\William\.julia\packages\DimensionalData\jl2xf\src\Dimensions\indexing.jl:11 [inlined]
  [7] getindex
    @ C:\Users\William\.julia\packages\DimensionalData\jl2xf\src\array\indexing.jl:16 [inlined]
  [8] BSplineInterpolation
    @ C:\Users\William\.julia\packages\Interpolations\y4lLj\src\b-splines\indexing.jl:8 [inlined]
  [9] FilledExtrapolation
    @ C:\Users\William\.julia\packages\Interpolations\y4lLj\src\extrapolation\filled.jl:35 [inlined]
 [10] _getindex
    @ C:\Users\William\.julia\packages\ImageTransformations\kQuat\src\ImageTransformations.jl:77 [inlined]
 [11] warp!(out::OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}, img::Interpolations.FilledExtrapolation{Float64, 2, Interpolations.BSplineInterpolation{Float64, 2, DimArray{Float64, 2, Tuple{X{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}, Y{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}}, Tuple{}, Matrix{Float64}, DimensionalData.NoName, DimensionalData.Dimensions.LookupArrays.NoMetadata}, Interpolations.BSpline{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, Interpolations.BSpline{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Float64}, tform::CoordinateTransformations.AffineMap{StaticArrays.SMatrix{2, 2, Float64, 4}, StaticArrays.SVector{2, Float64}})  
    @ ImageTransformations C:\Users\William\.julia\packages\ImageTransformations\kQuat\src\warp.jl:174
 [12] warp
    @ C:\Users\William\.julia\packages\ImageTransformations\kQuat\src\warp.jl:162 [inlined]
 [13] #warp#8
    @ C:\Users\William\.julia\packages\ImageTransformations\kQuat\src\warp.jl:181 [inlined]
 [14] warp
    @ C:\Users\William\.julia\packages\ImageTransformations\kQuat\src\warp.jl:180 [inlined]
 [15] imrotate(img::DimArray{Float64, 2, Tuple{X{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}, Y{DimensionalData.Dimensions.LookupArrays.NoLookup{Base.OneTo{Int64}}}}, Tuple{}, Matrix{Float64}, DimensionalData.NoName, DimensionalData.Dimensions.LookupArrays.NoMetadata}, θ::Float64, inds::Nothing; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ImageTransformations C:\Users\William\.julia\packages\ImageTransformations\kQuat\src\warp.jl:250
 [16] imrotate (repeats 2 times)
    @ C:\Users\William\.julia\packages\ImageTransformations\kQuat\src\warp.jl:247 [inlined]
 [17] top-level scope
    @ REPL[33]:1

Do you have any suggestions for how I could resolve these? If needed I could add overrides everywhere that unwrap the astroimage, but I'm wondering if there is an easier way.

Thanks!

rafaqz commented 2 years ago

This is a method I've missed on the lookup index in DimensionalData.jl. I will add and test in the next few days, although apologies I'm too busy today.

rafaqz commented 2 years ago

@sefffal this should work with DimensionalData.jl v0.20.6 when it's registered. The problem was not properly handling the mixed Int/CartesianIndex that ImageFiltering.jl uses.

sefffal commented 2 years ago

Thanks so much for your help @rafaqz !

I didn’t actually know mixed int and CatesianIndexes were allowed. Corner cases like that are another reason why an AbstractArrayTests.jl would be useful (https://github.com/rafaqz/DimensionalData.jl/issues/374#issuecomment-1141668819). I would actually argue that there should be tests for all the documented interfaces in Base.

rafaqz commented 2 years ago

Absolutely. We're discussing ways to implement a package like that over at: https://github.com/lorenzoh/Invariants.jl/issues/2

sefffal commented 2 years ago

The docs aren't ready yet & there are still some rough edges but I think this is about ready to merge.

The last major outstanding thing is FileIO.jl registration. I currently have

    # TODO: This should be registered correctly with FileIO
    del_format(format"FITS")
    add_format(format"FITS",
        # See https://www.loc.gov/preservation/digital/formats/fdd/fdd000317.shtml#sign
        [0x53,0x49,0x4d,0x50,0x4c,0x45,0x20,0x20,0x3d,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x54],
        [".fit", ".fits", ".fts", ".FIT", ".FITS", ".FTS", ".fit",],
        [:FITSIO => UUID("525bcba6-941b-5504-bd06-fd0dc1a4d2eb")],
        [:AstroImages => UUID("fe3fc30c-9b16-11e9-1c73-17dabf39f4ad")]
    )

in AstroImages.__init__() which is obviously a hack.

When running load("img.fits") FileIO goes down the list, so AstroImages will be used if installed and otherwise fallback to FITSIO.

Does anyone know the procedure to request that change in the FileIO registry? Do I need to run this the (other) FITSIO maintainers?

Edit: it would be nice for this to work with .fits.gz but I wasn't able to make that work. I'm not sure FileIO can differentiate .fits.gz from .gz.

sefffal commented 2 years ago

Even with the FileIO registration question outstanding, I think this is ready for a review and merge. After merge, I can keep adding tests and continue progress on the docs.

sefffal commented 2 years ago

Thanks @giordano and everyone else for your help and guidance along the way! Will merge now and prep to tag a new release.