diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 00000000..69dd5e98 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,37 @@ +name: CI +on: + pull_request: + branches: + - master + push: + branches: + - master + tags: '*' +jobs: + test: + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - '1.6' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'. + - '1' # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia. + - 'nightly' + os: + - ubuntu-latest + arch: + - x64 + steps: + - uses: actions/checkout@v3 + - uses: julia-actions/setup-julia@v1 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v1 + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-runtest@v1 + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v1 + with: + file: lcov.info diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml new file mode 100644 index 00000000..9928b682 --- /dev/null +++ b/.github/workflows/docs.yml @@ -0,0 +1,26 @@ +name: Documentation + +on: + push: + branches: + - master + tags: '*' + pull_request: + +jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + - uses: julia-actions/setup-julia@latest + with: + version: '1.7' + - name: Install dependencies + run: julia --color=yes --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' + - uses: julia-actions/cache@v1 + - name: Build and deploy + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # If authenticating with GitHub Actions token + # DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # If authenticating with SSH deploy key + run: julia --color=yes --project=docs/ docs/make.jl + diff --git a/.gitignore b/.gitignore index 8d931fb1..28387fd8 100644 --- a/.gitignore +++ b/.gitignore @@ -7,3 +7,7 @@ Manifest.toml # test files /test/data + + +# built docs +/docs/build diff --git a/Project.toml b/Project.toml index 03f49a08..ff152f79 100644 --- a/Project.toml +++ b/Project.toml @@ -1,29 +1,51 @@ name = "AstroImages" uuid = "fe3fc30c-9b16-11e9-1c73-17dabf39f4ad" -authors = ["Mosè Giordano", "Rohit Kumar"] -version = "0.2.0" +authors = ["Mosè Giordano", "Rohit Kumar", "William Thompson"] +version = "0.3.0" [deps] +AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" +AstroAngles = "5c4adb95-c1fc-4c53-b4ea-2a94080c53d2" +ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" +DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0" FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" -ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1" -Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0" -Interact = "c601a237-2ae4-5e1e-952c-7a85b0c7eef1" +ImageAxes = "2803e5a7-5153-5ecf-9a86-9b4c37f5f5ac" +ImageBase = "c817782e-172a-44cc-b673-b171935fbb9e" +ImageShow = "4e3cecfd-b093-5904-9786-8bbb286a6a31" +MappedArrays = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900" +PlotUtils = "995b91a9-d308-5afd-9ec6-746e21dbc043" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" -Reproject = "d1dcc2e6-806e-11e9-2897-3f99785db2ae" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" +UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" WCS = "15f3aee2-9e10-537f-b834-a6fb8bdb944d" -MappedArrays = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900" [compat] -julia = "^1.0.0" -Reproject = "^0.3.0" +AbstractFFTs = "1.1" +AstroAngles = "0.1" +ColorSchemes = "3.18" +DimensionalData = "^0.20" +FITSIO = "0.16" +FileIO = "1.14" +ImageAxes = "0.6" +ImageBase = "^0.1.5" +ImageShow = "0.3" +MappedArrays = "0.4" +PlotUtils = "1.2" +RecipesBase = "1.2" +Tables = "1.7" +WCS = "0.6" +julia = "^1.6.0" [extras] +FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb" +ImageBase = "c817782e-172a-44cc-b673-b171935fbb9e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -Widgets = "cc8bc4a8-27d6-5769-a93b-9d913e69aa62" -JLD = "4138dd39-2aa7-5051-a626-17a0bb65d9c8" -SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce" +WCS = "15f3aee2-9e10-537f-b834-a6fb8bdb944d" [targets] -test = ["Test", "Random", "Widgets", "JLD", "SHA"] +test = ["Test", "WCS", "FITSIO", "Random", "Statistics", "ImageBase"] diff --git a/README.md b/README.md index a0ef7e19..2f52584b 100644 --- a/README.md +++ b/README.md @@ -1,24 +1,23 @@ # AstroImages.jl -| **Build Status** | **Code Coverage** | -|:-----------------------------------------:|:-------------------------------:| -| [![Build Status][travis-img]][travis-url] | [![][coveral-img]][coveral-url] | -| [![Build Status][appvey-img]][appvey-url] | [![][codecov-img]][codecov-url] | +| **Documentation** | **Build Status** | **Code Coverage** | +|:------------------|:-----------------------------------------:|:-------------------------------:| +| [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://juliaastro.github.io/AstroImages.jl/dev/) | ![Build status](https://github.com/JuliaAstro/AstroImages/actions/workflows/ci.yml/badge.svg) | [![][codecov-img]][codecov-url] | Introduction ------------ -`AstroImage.jl` allows you to plot images from an -astronomical [`FITS`](https://en.wikipedia.org/wiki/FITS) file using the +`AstroImages.jl` allows you load and visualize images from a +astronomical [`FITS`](https://en.wikipedia.org/wiki/FITS) files using the popular [`Images.jl`](https://github.com/JuliaImages/Images.jl) and [`Plots.jl`](https://github.com/JuliaPlots/Plots.jl) Julia packages. -`AstroImage.jl` uses [`FITSIO.jl`](https://github.com/JuliaAstro/FITSIO.jl) to +`AstroImages.jl` uses [`FITSIO.jl`](https://github.com/JuliaAstro/FITSIO.jl) to read FITS files. Installation ------------ -`AstroImage.jl` is available for Julia 1.0 and later versions, and can be +`AstroImages.jl` is available for Julia 1.6 and later versions, and can be installed with [Julia built-in package manager](https://docs.julialang.org/en/v1/stdlib/Pkg/). @@ -26,6 +25,8 @@ manager](https://docs.julialang.org/en/v1/stdlib/Pkg/). pkg> add AstroImages ``` +You may also need to install `ImageIO.jl` for images to display in certain environments. + Usage ----- @@ -35,6 +36,10 @@ After installing the package, you can start using it with julia> using AstroImages ``` +Images will automatically display in many environments, including VS Code, Jupyter, and Pluto. +If you're using a REPL, you may want to install an external viewer like ImageShow.jl, ElectronDisplay.jl, +or ImageInTerminal.jl. + ## Reading extensions from FITS file You can load and read the the first extension of a FITS file with the `load` @@ -55,6 +60,8 @@ julia> load("file.fits", 3) [...] ``` +If unspecified, the first image HDU will be loaded. + ## AstroImage type The package provides a new type, `AstroImage` to integrate FITS images with @@ -63,31 +70,52 @@ the same syntax as `load`. This command: ```julia julia> img = AstroImage("file.fits") -AstroImages.AstroImage{UInt16,ColorTypes.Gray,1,Float64}[...] ``` -will read the first valid extension from the `file.fits` file and wrap its content in -a `NTuple{N, Matrix{Gray}}`, that can be easily used with `Images.jl` and related packages. +will read the first valid extension from the `file.fits` file. If you are working in a Jupyter notebook, an `AstroImage` object is automatically rendered as a PNG image. `AstroImage` automatically extracts and store `wcs` information of images in a `NTuple{N, WCSTransform}`. -## Forming RGB image -`AstroImage` can automatically construct a RGB image if 3 different colour band data is given. +## Headers +FITS Headers can be accessed directly from an AstroImage: +```julia +julia> img["HEAD1"] = 1.0 +julia> img["HEAD1",Comment] = "A comment describes the meaning of a header keyword" +julia> img["HEAD1"] +1.0 + +julia> push!(img, History, "We can record the history of processes applied to this image in header HISTORY entries.") +``` + +## Visualization + +Any AbstractArray (including an AstroImage) can be displayed using `imview`. This function renders an +arbitrary array into an array of `RGBA` values using a number of parameters. If the input is an AstroImage{<:Number}, +an AstroImage{RGBA} will be returned that retains headers, WCS information, etc. ```julia -julia> img = AstroImage(RGB, ("file1.fits","file2.fits", "file3.fits")) +julia> imview(img; clims=Percent(99.5), cmap=:magma, stretch=identity, contrast=1.0, bias=0.5) ``` -Where 1st index of `file1.fits`, `file2.fits`, `file3.fits` contains band data of red, blue and green channels respectively. -Optionally, `ccd2rgb` method can be used to form a coloured image from 3 bands without creating an `AstroImage`. +Very large Images are automatically downscaled to ensure consistent performance using `restrict` from Images.jl. This function filters the data before downscaling to prevent aliasing, so it may take a moment for truly huge images. In these cases, a faster method that doesn't prevent aliasing would be `imview(img[begin:10:end, begin:10:end])` or similar. + +`imview` is called automatically on `AstroImage{<:Number}` when using a Julia environment with rich graphical IO capabilities (e.g. VSCode, Jupyter, Pluto, etc.). +The defaults for this case can be modified using `AstroImages.set_clims!(...)`, `AstroImages.set_cmap!(...)`, and `AstroImages.set_stretch!(...)`. + +## Forming Color Composite Images -The formed image can be accessed using `img.property.rgb_image`. -`set_brightness!` and `set_contrast!` methods can be used to change brightness and contrast of formed `rgb_image`. -`add_label!` method can be used to add/store Astronomical labels in an `AstroImage`. -`reset!` method resets `brightness`, `contrast` and `label` fields to defaults and construct a fresh `rgb_image` without any brightness, contrast operations. +A color composite image (e.g. RGB) can be constructed using the `composecolors` function. +```julia +julia> rgb = composecolors([img1, img2, img3]) +``` +Where `img1`, `img2`, `img3` are arrays or AstroImages containing data of red, blue and green channels respectively. + +`composecolors` also supports more complex mappings, for example merging two bands according to color schemes from +ColorSchemes.jl. +See the docs for more information. ## Plotting an AstroImage @@ -97,10 +125,26 @@ An `AstroImage` object can be plotted with `Plots.jl` package. Just use ```julia julia> using Plots -julia> plot(img) +julia> implot(img) ``` -and the image will be displayed as a heatmap using your favorite backend. +and the image will be displayed as an image series using your favorite backend. +Plotly, PyPlot, and GR backends have been tested. + +`implot` supports all the same syntax as `imview` in addition to keyword arguments +for controlling axis tick marks, WCS grid lines, and the colorbar. + +## Resolving World Coordinates +If your FITS file contains world coordinate system headers, AstroImages.jl can use WCS.jl to convert between pixel and world coordinates. +This works even if you have sliced or your image to select a region of interest: + +```julia +julia> img_slice = img[100:200,100:200] +julia> coords_world = pix_to_world(img_slice, [5,5]) +[..., ...] +julia> world_to_pix(img_slice, coords_world) +[5.0,5.0] # approximately +``` License ------- @@ -108,14 +152,6 @@ License The `AstroImages.jl` package is licensed under the MIT "Expat" License. The original author is Mosè Giordano. -[travis-img]: https://travis-ci.org/JuliaAstro/AstroImages.jl.svg?branch=master -[travis-url]: https://travis-ci.org/JuliaAstro/AstroImages.jl - -[appvey-img]: https://ci.appveyor.com/api/projects/status/7gaxwe0c8hjx3d1s?svg=true -[appvey-url]: https://ci.appveyor.com/project/giordano/astroimages-jl - -[coveral-img]: https://coveralls.io/repos/JuliaAstro/AstroImages.jl/badge.svg?branch=master&service=github -[coveral-url]: https://coveralls.io/github/JuliaAstro/AstroImages.jl?branch=master [codecov-img]: http://codecov.io/github/JuliaAstro/AstroImages.jl/coverage.svg?branch=master [codecov-url]: http://codecov.io/github/JuliaAstro/AstroImages.jl?branch=master diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 00000000..00f7d4a0 --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,15 @@ +[deps] +AstroImages = "fe3fc30c-9b16-11e9-1c73-17dabf39f4ad" +DemoCards = "311a05b2-6137-4a5a-b473-18580a3d38b5" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5" +ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19" +ImageTransformations = "02fcd773-0e25-5acc-982a-7f6622650795" +Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0" +Photometry = "af68cb61-81ac-52ed-8703-edc140936be4" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Reproject = "d1dcc2e6-806e-11e9-2897-3f99785db2ae" +WCS = "15f3aee2-9e10-537f-b834-a6fb8bdb944d" + +[compat] +Documenter = "0.27" diff --git a/docs/examples/basics/displaying.jl b/docs/examples/basics/displaying.jl new file mode 100644 index 00000000..39496708 --- /dev/null +++ b/docs/examples/basics/displaying.jl @@ -0,0 +1,60 @@ +# --- +# title: Displaying Images +# author: "[William Thompson](https://github.com/sefffal)" +# cover: assets/displaying-images.png +# --- + +# We'll start by downloading a sample image. If you have an image stored locally, +# you would skip this step. +using AstroImages + +AstroImages.set_clims!(Percent(99.5)) #src +AstroImages.set_cmap!(:magma) #src +AstroImages.set_stretch!(identity) #src + + +# Any AbstractArray can be visualized with the `imview` function. +arr = randn(10,10) +imview(arr) + +# Let's load an astronomical image to see how we can tweak its display +fname = download( + "http://www.astro.uvic.ca/~wthompson/astroimages/fits/656nmos.fits", + "eagle-656nmos.fits" +); +img = AstroImage("eagle-656nmos.fits"); +imview(img) + +# We can adjust the color limits manually +imview(img, clims=(0,100)) + +# Or provide a function to calculate them for us +imview(img, clims=extrema) + +# AstroImages includes some handy callables, like Percent and Zscale.flags +# `Percent` sets the limits to include some central percentage of the data range +# For example, 95% sets the color limits to clip the top and bottom 2.5% of pixels. +# Percent(99.5) is the default value of clims. +imview(img, clims=Percent(95)) + + + + +# Arrays wrapped by `AstroImage` are displayed automatically using `imview` +AstroImage(randn(10,10)) + +# The settings for automatic imview are controlled using package defaults that can +# be adjusted to suit your tastes +AstroImages.set_clims!(Zscale()) # Display the full range automatically +AstroImages.set_cmap!(:viridis) +AstroImages.set_stretch!(asinhstretch) +AstroImage(randn(10,10)) + +# --- restore defaults --- #src +AstroImages.set_clims!(Percent(99.5)) #src +AstroImages.set_cmap!(:magma) #src +AstroImages.set_stretch!(identity) #src + +# --- save covers --- #src +mkpath("assets") #src +save("assets/loading-images.png", imview(img)) #src diff --git a/docs/examples/basics/loading.jl b/docs/examples/basics/loading.jl new file mode 100644 index 00000000..fc4c59b9 --- /dev/null +++ b/docs/examples/basics/loading.jl @@ -0,0 +1,23 @@ +# --- +# title: Loading Images +# description: Loading FITS images from files. +# author: "[William Thompson](https://github.com/sefffal)" +# cover: assets/loading-images.png +# --- + +# We'll start by downloading a sample image. If you have an image stored locally, +# you would skip this step. +using AstroImages +fname = download( + "http://www.astro.uvic.ca/~wthompson/astroimages/fits/656nmos.fits", + "eagle-656nmos.fits" +); + +# Load the image by filename. +# If unspecified, the image is loaded from the first image-HDU in the fits file. +img = AstroImage("eagle-656nmos.fits") + + +# --- save covers --- #src +mkpath("assets") #src +save("assets/loading-images.png", imview(img)) #src diff --git a/docs/examples/config.json b/docs/examples/config.json new file mode 100644 index 00000000..02081db5 --- /dev/null +++ b/docs/examples/config.json @@ -0,0 +1,4 @@ +{ + "template": "index.md", + "theme": "list" +} diff --git a/docs/examples/index.md b/docs/examples/index.md new file mode 100644 index 00000000..96ca2426 --- /dev/null +++ b/docs/examples/index.md @@ -0,0 +1,3 @@ +# Examples + +{{{democards}}} \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 00000000..0364239f --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,72 @@ +using Documenter, DemoCards, AstroImages + +# Deps for examples +ENV["GKSwstype"] = "nul" +using Plots, Photometry, ImageTransformations, ImageFiltering, WCS, Reproject + +setup = quote + using AstroImages + using Random + Random.seed!(123456) + + AstroImages.set_clims!(Percent(99.5)) + AstroImages.set_cmap!(:magma) + AstroImages.set_stretch!(identity) +end +DocMeta.setdocmeta!(Photometry, :DocTestSetup, setup; recursive = true) + +# 1. generate demo files +demopage, postprocess_cb, demo_assets = makedemos("examples") # this is the relative path to docs/ + +# if there are generated css assets, pass it to Documenter.HTML +assets = [] +isnothing(demo_assets) || (push!(assets, demo_assets)) + +# 2. normal Documenter usage +format = Documenter.HTML(assets = assets) +makedocs(format = format, + pages = [ + "Home" => "index.md", + demopage, + ], + sitename = "Awesome demos") + + +makedocs( + sitename="AstroImages.jl", + pages = [ + "Home" => "index.md", + "Manual" => [ + "Getting Started" => "manual/getting-started.md", + "Loading & Saving Images" => "manual/loading-images.md", + "Displaying Images" => "manual/displaying-images.md", + "Headers" => "manual/headers.md", + "Dimensions and World Coordinates" => "manual/dimensions-and-world-coordinates.md", + "Polarization" => "manual/polarization.md", + "Spectral Axes" => "manual/spec.md", + "Preserving Wrapper" => "manual/preserving-wrapper.md", + "Conventions" => "manual/conventions.md", + ], + "Guides" => [ + "Blurring & Filtering Images" => "guide/image-filtering.md", + "Transforming Images" => "guide/image-transformations.md", + "Reprojecting Images" => "guide/reproject.md", + "Extracting Photometry" => "guide/photometry.md", + ], + demopage, + "API" => "api.md", + ], + format = Documenter.HTML( + prettyurls = get(ENV, "CI", nothing) == "true" + ), + workdir=".." +) + +# 3. postprocess after makedocs +postprocess_cb() + +deploydocs( + repo = "github.com/sefffal/AstroImages.jl.git", + devbranch = "master", + push_preview = true +) diff --git a/docs/src/api.md b/docs/src/api.md new file mode 100644 index 00000000..f364756c --- /dev/null +++ b/docs/src/api.md @@ -0,0 +1,38 @@ +# API Documentation + + +```@docs +load +save +AstroImage +imview +implot +dims +refdims +Comment +History +pix_to_world +world_to_pix +X +Y +Z +Dim +At +Near +.. +header +wcs +WCSGrid +composecolors +Zscale +Percent +logstretch +powstretch +sqrtstretch +squarestretch +asinhstretch +sinhstretch +powerdiststretch +copyheader +shareheader +``` \ No newline at end of file diff --git a/docs/src/assets/logo.png b/docs/src/assets/logo.png new file mode 100644 index 00000000..46bf5cfc Binary files /dev/null and b/docs/src/assets/logo.png differ diff --git a/docs/src/guide/image-filtering.md b/docs/src/guide/image-filtering.md new file mode 100644 index 00000000..33368c48 --- /dev/null +++ b/docs/src/guide/image-filtering.md @@ -0,0 +1,96 @@ +# Image Filtering + +The package [ImageFiltering.jl](https://juliaimages.org/ImageFiltering.jl/stable/) makes it easy to apply arbitrary filters to images. + +```@setup ex1 +using AstroImages +AstroImages.set_clims!(Percent(99.5)) +AstroImages.set_cmap!(:magma) +AstroImages.set_stretch!(identity) +``` + + +## Gaussian Blurs +Let's start by downloading a radio image of Hercules A: +```@example ex1 +using AstroImages +using ImageFiltering + +fname = download( + "http://www.astro.uvic.ca/~wthompson/astroimages/fits/herca/herca_radio.fits", + "herca-radio.fits" +) + +herca = load("herca-radio.fits") +``` + +Let's now apply a Gaussian blur (aka a low pass filter) using the `imfilter` function: +```@example ex1 +herca_blur_20 = imfilter(herca, Kernel.gaussian(20.0)) +``` +The image has been smoothed out by convolving it with a wide Gaussian. + +Let's now do the opposite and perform a high-pass filter. This will bring out faint variations in structure. +We can do this by subtracting a blurred image from the original: + +```@example ex1 +herca_blur_4 = imfilter(herca, Kernel.gaussian(4.0)) +herca_highpass = herca .- herca_blur_4 +``` +We now see lots of faint structure inside the jets! + + +Finally, let's adjust how the image is displayed and apply a non-linear stretch: +```@example ex1 +imview( + herca_highpass, + cmap=:seaborn_rocket_gradient, + clims=(-50,1500), + stretch=asinhstretch +) +``` + +If you have Plots loaded, we can add a colorbar and coordinate axes by switching to `implot`: +```@example ex1 +using Plots +implot( + herca_highpass, + cmap=:seaborn_rocket_gradient, + clims=(-50,1500), + stretch=asinhstretch +) +``` + + +## Median Filtering +In addition to linear filters using `imfilter`, ImageFiltering.jl also includes a great function called `mapwindow`. This functions allows you to map an arbitrary function over a patch of an image. + +Let's use `mapwindow` to perform a median filter. This is a great way to suppress salt and pepper noise, or remove stars from some images. + +We'll use a Hubble picture of the Eagle nebula: +```@example ex1 +using AstroImages +using ImageFiltering + +fname = download( + "http://www.astro.uvic.ca/~wthompson/astroimages/fits/eagle/673nmos.fits", + "eagle-673nmos.fits" +) + +eagle673 = load("eagle-673nmos.fits") +``` +The data is originally from https://esahubble.org/projects/fits_liberator/eagledata/. + + +We can apply a median filter using `mapwindow`. Make sure the patch size is an odd number in each direction! +```@example ex1 +using Statistics +medfilt = copyheader(eagle673, mapwindow(median, eagle673, (11,11))) +``` + +We use `copyheader` here since `mapwindow` returns a plain array and drops the image meta data. + +We can put this side by side with the original to see how some of the faint stars have been removed from the image: +```@example ex1 +imview([eagle673[1:800,1:800]; medfilt[1:800,1:800]]) +``` \ No newline at end of file diff --git a/docs/src/guide/image-transformations.md b/docs/src/guide/image-transformations.md new file mode 100644 index 00000000..f3461d08 --- /dev/null +++ b/docs/src/guide/image-transformations.md @@ -0,0 +1,57 @@ +# Image Transformations + +The [ImageTransformations.jl](https://juliaimages.org/latest/pkgs/transformations/) package contains many useful functions for manipulating astronomical images. + +Note however that many of these functions drop the AstroImage wrapper and return plain +arrays or OffsetArrays. They can be re-wrapped using `copyheader` or `shareheader` if you'd like to preserve the FITS header, dimension labels, WCS information, etc. + +You can install ImageTransformations by running `] add ImageTransformations` at the REPL. + + +```@setup transforms +using AstroImages +AstroImages.set_clims!(Percent(99.5)) +AstroImages.set_cmap!(:magma) +AstroImages.set_stretch!(identity) +``` + +For these examples, we'll download an image of the Antenae galaxies from Hubble: +```@example transforms +using AstroImages +using ImageTransformations + +fname = download( + "http://www.astro.uvic.ca/~wthompson/astroimages/fits/antenae/blue.fits", + "ant-blue.fits" +) + +antblue = load("ant-blue.fits") + +# We'll change the defaults to avoid setting them each time +AstroImages.set_clims!(Percent(99)) +AstroImages.set_cmap!(:ice) +AstroImages.set_stretch!(asinhstretch) + +imview(antblue) +``` + +## Rotations + +We can rotate images using the `imrotate` function. + +```@example transforms +imrotate(antblue, 3π/4) |> imview +``` +The rotation angle is in radians, but you can use the function `rad2deg` to convert from degrees. + +## Resizing +We can resize images using the `imresize` function: +```@example transforms +imresize(antblue, ratio=0.2) |> imview +``` + +## Arbitrary Transformations +Arbitrary transformations can be performed using ImageTransformation's `warp` function. See the documentation linked above for more details. + +## Mapping from One Coordinate System to Another +For transforming an image from one coordiante system (say, RA & DEC) to another (e.g., galactic lattitude & logitude), see [Reprojecting Images](@ref). \ No newline at end of file diff --git a/docs/src/guide/photometry.md b/docs/src/guide/photometry.md new file mode 100644 index 00000000..58d5a6cc --- /dev/null +++ b/docs/src/guide/photometry.md @@ -0,0 +1,84 @@ +# Photometry + +The following examples are adapted from [Photometry.jl](https://github.com/JuliaAstro/Photometry.jl/) to show the same examples +combined with AstroImages.jl. +To learn how to measure background levels, perform aperture photometry, etc see the [Photometry.jl documentation](https://juliaastro.github.io/Photometry.jl/dev/). + + +## Background Estimation + +From Photometry.jl: +> Estimating backgrounds is an important step in performing photometry. Ideally, we could perfectly describe the background with a scalar value or with some distribution. Unfortunately, it's impossible for us to precisely separate the background and foreground signals. Here, we use mixture of robust statistical estimators and meshing to let us get the spatially varying background from an astronomical photo. +> Let's show an example +> Now let's try and estimate the background using estimate_background. First, we'll si gma-clip to try and remove the signals from the stars. Then, the background is broken down into boxes, in this case of size (50, 50). Within each box, the given statistical estimators get the background value and RMS. By default, we use SourceExtractorBackground and StdRMS. This creates a low-resolution image, which we then need to resize. We can accomplish this using an interpolator, by default a cubic-spline interpolator via ZoomInterpolator. The end result is a smooth estimate of the spatially varying background and background RMS. + +```@setup phot +using AstroImages +AstroImages.set_clims!(Percent(99.5)) +AstroImages.set_cmap!(:magma) +AstroImages.set_stretch!(identity) +``` + +```@example phot +using Photometry +using AstroImages +using Plots # optional, for implot functionality + +# Download our image, courtesy of astropy +image = AstroImage(download("https://rawcdn.githack.com/astropy/photutils-datasets/8c97b4fa3a6c9e6ea072faeed2d49a20585658ba/data/M6707HH.fits")) + +# sigma-clip +clipped = sigma_clip(image, 1, fill=NaN) + +# get background and background rms with box-size (50, 50) +bkg, bkg_rms = estimate_background(clipped, 50) + +imview(image) +imview(clipped) +imview(bkg) +imview(bkg_rms) +``` + +Or, if you have Plots loaded: +```@example phot +using Plots + + AstroImages.set_clims!(Percent(99.5)) + AstroImages.set_cmap!(:magma) + AstroImages.set_stretch!(identity) +plot( + implot(image, title="Original"), + implot(clipped, title="Sigma-Clipped"), + implot(bkg, title="Background"), + implot(bkg_rms, title="Background RMS"), + layout=(2, 2) +) +``` +![](/assets/manual-photometry-2.png) + + +> We could apply a median filter, too, by specifying filter_size +```@example phot +# get background and background rms with box-size (50, 50) and filter_size (5, 5) +bkg_f, bkg_rms_f = estimate_background(clipped, 50, filter_size=5) + +# plot +plot( + implot(bkg, title="Unfiltered", ylabel="Background"), + implot(bkg_f, title="Filtered"), + implot(bkg_rms, ylabel="RMS"), + implot(bkg_rms_f); + layout=(2, 2),) +``` + +> Now we can see our image after subtracting the filtered background and ready for Aperture Photometry! + +```@example phot +subt = image .- bkg_f[axes(image)...] +clims = extrema(vcat(vec(image), vec(subt))) +plot( + implot(image; title="Original", clims), + implot(subt; title="Subtracted", clims), + size=(1600,1000) +) +``` \ No newline at end of file diff --git a/docs/src/guide/reproject.md b/docs/src/guide/reproject.md new file mode 100644 index 00000000..0da3edfc --- /dev/null +++ b/docs/src/guide/reproject.md @@ -0,0 +1 @@ +# Reprojecting Images \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 00000000..5a3b6560 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,11 @@ +# Home + +[GitHub link](https://github.com/JuliaAstro/AstroImages.jl) + + +`AstroImage.jl` allows you to plot images from an +astronomical [`FITS`](https://en.wikipedia.org/wiki/FITS) file using the +popular [`Images.jl`](https://github.com/JuliaImages/Images.jl) +and [`Plots.jl`](https://github.com/JuliaPlots/Plots.jl) Julia packages. +`AstroImage.jl` uses [`FITSIO.jl`](https://github.com/JuliaAstro/FITSIO.jl) to +read FITS files. diff --git a/docs/src/manual/conventions.md b/docs/src/manual/conventions.md new file mode 100644 index 00000000..35881fc9 --- /dev/null +++ b/docs/src/manual/conventions.md @@ -0,0 +1,23 @@ + +# Conventions + +In the Julia Astro ecosystem, images follow the following conventions. + +## Axes +For simple 2D images, the first axis is the horizontal axis and the second axis is the vertical axis. +So images are indexed by `img[xi, yi]`. + +The origin is at the bottom left of the image, so `img[1,1]` refers to the bottom left corner +as does `img[begin,begin]`. +`img[end,end]` is the top right corner, `img[begin,end]` is the top left, etc. + +Note that this is transposed and flipped from how how Julia prints arrays at the REPL, + + +## Pixels +This library considers the exact location of `img[1,1]` to be the center of the pixel in the bottom left corner. +This means that plot limits should have the `1` tick slightly away from the left/bottom spines of the image. +The default plot limits for `implot` are `-0.5` to `end+0.5` along both axes. + +There is a [known bug](https://github.com/JuliaPlots/Plots.jl/issues/4158) with the Plots.jl GR backend that leads ticks to be slightly offset. PyPlot and Plotly backends +show the correct tick locations. \ No newline at end of file diff --git a/docs/src/manual/dimensions-and-world-coordinates.md b/docs/src/manual/dimensions-and-world-coordinates.md new file mode 100644 index 00000000..94c765da --- /dev/null +++ b/docs/src/manual/dimensions-and-world-coordinates.md @@ -0,0 +1,36 @@ +# Dimensions and World Coordinates + +AstroImages are based on [Dimensional Data](https://github.com/rafaqz/DimensionalData.jl). Each axis is assigned a dimension name +and the indices are tracked. +The automatic dimension names are `X`, `Y`, `Z`, `Dim{4}`, `Dim{5}`, and so on; however you can pass in other names or orders to the load function and/or AstroImage contructor: + +```julia +julia> img = load("img.fits",1,(Y=1:1600,Z=1:1600)) +1600×1600 AstroImage{Float32,2} with dimensions: + Y Sampled 1:1600 ForwardOrdered Regular Points, + Z Sampled 1:1600 ForwardOrdered Regular Points +``` +Other useful dimension names are `Spec` for spectral axes, `Pol` for polarization data, and `Ti` for time axes. + +These will be further discussed in Dimensions and World Coordinates. + +For certain applications, it may be useful to use offset axes or axes with different steps. +```julia +julia> img = load("img.fits",1,(X=801:2400,Y=1:2:3200)) +1600×1600 AstroImage{Float32,2} with dimensions: + X Sampled 801:2400 ForwardOrdered Regular Points, + Y Sampled 1:2:3199 ForwardOrdered Regular Points +... +``` + +Unlike OffsetArrays, the usual indexing remains so `img[1,1]` is still the bottom left of the image; +however, data can be looked up according to the offset axes when using specifiers: +```julia +julia> img[X=Near(2000),Y=1..100] +50-element AstroImage{Float32,1} with dimensions: + Y Sampled 1:2:99 ForwardOrdered Regular Points +and reference dimensions: + X Sampled 2000:2000 ForwardOrdered Regular Points + 0.0 +``` + diff --git a/docs/src/manual/displaying-images.md b/docs/src/manual/displaying-images.md new file mode 100644 index 00000000..11a225f7 --- /dev/null +++ b/docs/src/manual/displaying-images.md @@ -0,0 +1,114 @@ +# Displaying Images + +The `imview` and `implot` functions are very similar. +Both allow any abstract array of numbers to be rendered into an image or a Plots.jl image series. +`implot` is largely a superset of `imview` because it also supports colorbars, tick marks, WCS grid lines, overplotting other data & shapes, and automatic +axis and title naming (from the FITS header if available). + +## `imview` + +Any AbstractArray (including an AstroImage) can be displayed using `imview`. This function renders an +arbitrary array into an array of `RGBA` values using a number of parameters. If the input is an AstroImage{<:Number}, +an AstroImage{RGBA} will be returned that retains headers, WCS information, etc. + +```@setup 1 +using AstroImages +using Plots +``` + +The defaults for the `imview` function are: +```@example 1 +img = randn(50,50); +imview(img; clims=Percent(99.5), cmap=:magma, stretch=identity, contrast=1.0, bias=0.5) +``` + +We can adjust the color limits explicitly: +```@example 1 +imview(img; clims=(-1, 1)) +``` + +Or pass a function/callable object to calculate them for us: +```@example 1 +imview(img; clims=Zscale()) +``` + +We turn off the colormap and use it in grayscale mode: +```@example 1 +imview(img; cmap=nothing) +``` + +Pass any color scheme from ColorSchemes.jl: +```@example 1 +imview(img; cmap=:ice) +``` +```@example 1 +imview(img; cmap=:seaborn_rocket_gradient) +``` + +Or an RGB or named color value: +```@example 1 +imview(img; cmap="#F00") +imview(img; cmap="red") +``` + +Let's now switch to an astronomical image: +```@example 1 +fname = download( + "http://www.astro.uvic.ca/~wthompson/astroimages/fits/656nmos.fits", + "eagle-656nmos.fits" +); +eagle = AstroImage("eagle-656nmos.fits") +``` + +We can apply a non-linear stretch like a log-scale, power-scale, or asinh stretch: +```@example 1 +imview(eagle, stretch=asinhstretch) +``` + +Once rendered, we can also tweak the bias and contrast: +```@example 1 +imview(eagle, stretch=asinhstretch, contrast=1.5) +``` +```@example 1 +imview(eagle, stretch=asinhstretch, contrast=1.5, bias=0.6) +``` +These are the parameters that change when you click and drag in some applications like DS9. + +Once rendered via `imview`, the resulting image can be saved in traditional image formats like PNG, JPG, GIF, etc: +```julia +save("out.png", imview(eagle, cmap=:viridis)) +``` + +Very large Images are automatically downscaled to ensure consistent performance using `restrict` from Images.jl. This function filters the data before downscaling to prevent aliasing, so it may take a moment for truly huge images. In these cases, a faster method that doesn't prevent aliasing would be `imview(img[begin:10:end, begin:10:end])` or similar. + +`imview` is called automatically on `AstroImage{<:Number}` when using a Julia environment with rich graphical IO capabilities (e.g. VSCode, Jupyter, Pluto, etc.). +The defaults for this case can be modified using `AstroImages.set_clims!(...)`, `AstroImages.set_cmap!(...)`, and `AstroImages.set_stretch!(...)`. + +## Note on Views +The function `imview` has its name because it produces a "view" into the image. The result from calling `imview` is an object that lazily maps data values into RGBA colors on the fly. +This means that if you change the underlying data array, the view will update (the next time it is shown). +If you have many data files to render, you may find it faster to create a single `imview` and then mutate the data in the underlying array. This is faster since `imview` only has to resolve colormaps and compute limits once. + +For example: +```julia +data = randn(100,100) +iv = imview(data) +display(iv) +data[1:50,1:50] .= 0 +display(iv) +``` +`iv` will reflect the changes to `data` when it is displayed the second time. + +## `implot` + +`implot` is a Plots.jl recipe, which means before you can use it you first have to load `Plots.jl`: + +```@example 1 +using Plots +``` + +`implot` accepts all the arguments `imview` does for controlling how data is rendered to the screen. + +```@example +implot(img; clims=Percent(99.5), cmap=:magma, stretch=identity, contrast=1.0, bias=0.5) +``` \ No newline at end of file diff --git a/docs/src/manual/getting-started.md b/docs/src/manual/getting-started.md new file mode 100644 index 00000000..b4c406c9 --- /dev/null +++ b/docs/src/manual/getting-started.md @@ -0,0 +1,24 @@ +# Getting Started + +To get started, you will first need to install AstroImages. +After starting Julia, enter package-mode by typing `]` and then +```julia-repl +pkg> add AstroImages +``` + +To display images and save them in traditional graphics formats like PNG, JPG, GIF, etc., you +will also need to add the `ImageIO` package. Once installed, this package doesn't need to be +loaded explicitly. + + +For some of the more advanced visualizations you may also want `Plots`: +```julia-repl +pkg> add Plots +``` + +To load the package, run: +```julia +using AstroImages +# And if desired: +using Plots +``` \ No newline at end of file diff --git a/docs/src/manual/headers.md b/docs/src/manual/headers.md new file mode 100644 index 00000000..81805f43 --- /dev/null +++ b/docs/src/manual/headers.md @@ -0,0 +1,47 @@ +# Headers + +FITS files consist of one or more HDUs (header data units), and each HDU can contain an N-dimensional image or table. +Before the data is a *header*. Headers contain (key, value, comment) groups as well as dedicated long-form COMMENT and HISTORY sections used to document, for example, the series of post-processing steps applied to an image. + +## Accessing Headers + +Here are some examples of how to set and read keys, comments, and history. + +Well start by making a blank image. +```julia +img = AstroImage(zeros(10,10)) +# Set keys to values with different data types +img["KEY1"] = 2 # Integer +img["KEY2"] = 2.0 # Float +img["KEY3"] = "STRING" +img["KEY4"] = true +img["KEY5"] = false +img["KEY6"] = nothing + +# Set comments +img["KEY1", Comment] = "A key with an integer value" + +# Read keys +a = img["KEY3"] + +# Read comment +com = img["KEY1", Comment] + +# Add long-form COMMENT +push!(img, Comment, """ +We now describe how to add a long form comment to the end of a header. +""") + +# Add HISTORY entry +push!(img, History, """ +We now describe how to add a long form history to the end of a header. +""") + +# Retrieve long form comments/ history +comment_strings = img[Comment] +history_strings = img[History] +``` + +Note that floating point values are formatted as ASCII strings when written to the FITS files, so the precision may be limited. + +`AstroImage` objects wrap a FITSIO.jl `FITSHeader`. If necessary, you can recover it using `header(img)`; however, in most cases you can access header keywords directly from the image. \ No newline at end of file diff --git a/docs/src/manual/loading-images.md b/docs/src/manual/loading-images.md new file mode 100644 index 00000000..edb3d9a4 --- /dev/null +++ b/docs/src/manual/loading-images.md @@ -0,0 +1,133 @@ +# Loading Images + +FITS (Flexible Image Transport System) files can be loaded and saved using AstroImages thanks to the FITSIO package. + +AstroImages is registered with [FileIO](https://juliaio.github.io/FileIO.jl/stable/), so if you have FileIO and AstroImages +installed you can get started with the `load` function. When you pass a file name with the appropriate file extension (".fits", ".fit", etc.) +FileIO will import AstroImages automatically. + +Alternatively, you can use the `AstroImage` contructor instead of load. This will work on fits files with any file extension, including compressed +files (e.g. ".fits.gz"). + +```julia +julia> img = load("myfitsimg.fits") +1600×1600 AstroImage{Float32,2} with dimensions: + X Sampled Base.OneTo(1600) ForwardOrdered Regular Points, + Y Sampled Base.OneTo(1600) ForwardOrdered Regular Points + 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 + ⋮ ⋱ + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +``` + +Note: if you are in an interactive environment like VSCode, Jupyter, or Pluto, instead of a REPL, AstroImages are automatically +rendered to images and displayed. You can see this plain text output by explicitly calling: +`show(stdout, MIME("text/plain"), img)`. + +Or: +```julia + julia> img = AstroImage("myfitsimg.fits.gz") +1600×1600 AstroImage{Float32,2} with dimensions: + X Sampled Base.OneTo(1600) ForwardOrdered Regular Points, + Y Sampled Base.OneTo(1600) ForwardOrdered Regular Points + 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 + ⋮ ⋱ + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +``` + +A FITS file can contain multiple N-dimensional images and tables. +When you call load or AstroImage with a file name and no other arguments, the package will search through the file +and return the first image HDU. That is, it will skip any FITS tables or empty HDUs with only headers. + +You can also specify an HDU number explicitly: +```julia +julia> img = load("myfitsimg.fits",1) +1600×1600 AstroImage{Float32,2} with dimensions: + X Sampled Base.OneTo(1600) ForwardOrdered Regular Points, + Y Sampled Base.OneTo(1600) ForwardOrdered Regular Points +... +``` +This way, you can load specific images from multi-extension files. + +You can load all HDUs simultaneously by passing `:`: + +```julia +julia> hdus = load("multiext.fits", :); +julia> hdus[2] # Second HDU as an AstroImage +10×10 AstroImage{Float64,2} with dimensions: + X Sampled Base.OneTo(10) ForwardOrdered Regular Points, + Y Sampled Base.OneTo(10) ForwardOrdered Regular Points + -0.777315 -1.36683 -0.580179 1.39629 … -2.14298 0.450059 0.432065 + -1.09619 0.789249 0.938415 0.959903 -0.88995 -1.29406 -0.4291 + 0.47427 -1.41855 0.814823 -1.15975 0.0427149 -1.20116 -0.0920709 + -0.179858 -1.60228 1.09648 -0.497927 -1.31824 -0.156529 -0.0223846 + 2.64162 0.131437 0.320476 0.331197 -0.914713 -1.55162 -0.18862 + 0.209669 -1.17923 -0.656512 0.000775311 … 0.377461 -0.24278 0.967202 + 1.01442 -0.762895 -2.13238 -0.456932 -0.415733 -1.21416 -1.6108 + 0.385626 0.389335 -0.00726015 0.309936 -0.533175 0.157878 0.100876 + -1.24799 0.461216 -0.868826 -0.255654 -0.37151 0.49479 -1.87129 + 1.39356 2.29254 0.0548325 1.50674 -0.0880865 0.580978 -1.81629 +julia> # Or: +julia> hdu1, hdu2, hdu3 = load("multiext.fits", :); +``` + +There is also limited support for table HDUs. In this case, a bare-bones Tables.jl compatible +object is returned. + +## Dimension Names +You may have noticed the entries above the image array: +``` +10×10 AstroImage{Float64,2} with dimensions: + X Sampled Base.OneTo(10) ForwardOrdered Regular Points, + Y Sampled Base.OneTo(10) ForwardOrdered Regular Points +``` + +AstroImages are based on [Dimensional Data](https://github.com/rafaqz/DimensionalData.jl). Each axis is assigned a dimension name +and the indices are tracked. +The automatic dimension names are `X`, `Y`, `Z`, `Dim{4}`, `Dim{5}`, and so on; however you can pass in other names or orders to the load function and/or AstroImage contructor: + +```julia +julia> img = load("img.fits",1,(Y=1:1600,Z=1:1600)) +1600×1600 AstroImage{Float32,2} with dimensions: + Y Sampled 1:1600 ForwardOrdered Regular Points, + Z Sampled 1:1600 ForwardOrdered Regular Points +``` +Other useful dimension names are `Spec` for spectral axes, `Pol` for polarization data, and `Ti` for time axes. + +These will be further discussed in Dimensions and World Coordinates. + + +## Saving Images +You can save one or more AstroImages and tables to a FITS file using the `save` function: + +```julia +julia> save("abc.fits", astroimage1, astroimage2, table1) +``` + +You can also save individual images to traditional graphics formats by first rendering them with `imview` (for more on imview, see Displaying Images). +```julia +julia> save("abc.png", imview(astroimage1)) +``` + +You can save animated GIFs by saving a 3D datacube that has been rendered with imview: +```julia +julia> cube = imview(AstroImage(randn(100,100,10))); +julia> save("abc.gif", cube, fps=10) + +julia> # Or a more complex example (changing color schemes each frame) +julia> img = randn(10,10) +julia> cube2 = [imview(img1, cmap=:magma) ;;; imview(img2, cmap=:plasma) ;;; imview(img3, cmap=:viridis)] +julia> # Alternative syntax: +julia> cube2 = cat(imview(img1, cmap=:magma), imview(img2, cmap=:plasma), imview(img3, cmap=:viridis), dims=3) +julia> save("abc.gif", cube, fps=10) +``` + diff --git a/docs/src/manual/polarization.md b/docs/src/manual/polarization.md new file mode 100644 index 00000000..e69de29b diff --git a/docs/src/manual/preserving-wrapper.md b/docs/src/manual/preserving-wrapper.md new file mode 100644 index 00000000..59df3fd8 --- /dev/null +++ b/docs/src/manual/preserving-wrapper.md @@ -0,0 +1,32 @@ +# Preserving the AstroImage Wrapper + +Wherever possible, overloads have been added to DimensionalData and AstroImages so that common operations retain the `AstroImage` wrapper with associated dimensions, FITS header, and WCS information. +Most of the time this works automatically if libraries follow good patterns like allocating outputs using `Base.similar`. +However, some other library functions may follow patterns like allocating a plain `Array` of the correct size and then filling it. + +To make it easier to work with these libraries, AstroImages exports two functions `copyheader` and `shareheader`. +These functions wrap an AbstractArray in an AstroImage while copying over the header, dimensions, and WCS info. + +Consider the function: +```julia +function badfunc(arr) + out = zeros(size(arr)) # instead of similar(arr) + out .= arr.^2 + return out +end +``` + +Calling `badfunc(astroimg)` will return a plain `Array` . + +We can use `copyheader` to retain the `AstroImage` wrapper: +```julia +copyheader(astroimg, badfunc(astroimg)) +``` + +For particularly incompatible functions that require an Array (not subtype of AbstractArray) we can go one step further: +```julia +copyheader(astroimg, worsefunc(parent(astroimg))) +# Or: +copyheader(astroimg, worsefunc(collect(astroimg))) +``` + diff --git a/docs/src/manual/spec.md b/docs/src/manual/spec.md new file mode 100644 index 00000000..e69de29b diff --git a/examples/imview-pluto.jl b/examples/imview-pluto.jl new file mode 100644 index 00000000..9c429eb4 --- /dev/null +++ b/examples/imview-pluto.jl @@ -0,0 +1,72 @@ +### A Pluto.jl notebook ### +# v0.19.6 + +using Markdown +using InteractiveUtils + +# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error). +macro bind(def, element) + quote + local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end + local el = $(esc(element)) + global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el) + el + end +end + +# ╔═╡ 685479e8-1ad5-48d8-b9fe-f2cf8a672700 +using AstroImages, PlutoUI + +# ╔═╡ 59e1675f-9426-4bc4-88cc-e686ed90b6b5 +md""" +Download a FITS image and open it. +Apply `restrict` to downscale 2x for faster rendering. +""" + +# ╔═╡ d1e5947b-2c1a-46fc-ab8f-feeba03453e7 +img = AstroImages.restrict( + AstroImage(download("http://www.astro.uvic.ca/~wthompson/astroimages/fits/656nmos.fits")) +); + +# ╔═╡ c9ebe984-4630-47c1-a941-795293f5b3c1 +md""" +Display options +""" + +# ╔═╡ a3e81f3f-203b-47b7-ac60-b4267eddfad4 +md""" + +| parameter | value | +|-----------|-------| +|`cmap` | $( @bind cmap Select([:magma, :turbo, :ice, :viridis, :seaborn_icefire_gradient, "red"]) ) | +|`clims`| $( @bind clims Select([Percent(99.5), Percent(95), Percent(80), Zscale(), (0, 400)]) ) | +| `stretch` | $( @bind stretch Select([identity, asinhstretch, logstretch, sqrtstretch, powstretch, powerdiststretch, squarestretch])) | +| `contrast` | $(@bind contrast Slider(0:0.1:2.0, default=1.0)) | +| `bias` | $(@bind bias Slider(0:0.1:1.0, default=0.5)) | +""" + +# ╔═╡ 2315ffec-dc49-413a-b0d6-1bcce2addd76 +imview(img; cmap, clims, stretch, contrast, bias) + +# ╔═╡ d2bd2f13-ed23-42c5-9317-5b48ec3a8bb7 +md""" +## `implot` +Uncomment the following cells to use `Plots` instead. +""" + +# ╔═╡ fe6b5b76-8b77-4bfc-a2e8-bcc0b78ad764 +# using Plots + +# ╔═╡ f557784e-828c-415e-abb0-964b3a9fe8ef +# implot(img; cmap, clims, stretch, contrast, bias) + +# ╔═╡ Cell order: +# ╠═685479e8-1ad5-48d8-b9fe-f2cf8a672700 +# ╟─59e1675f-9426-4bc4-88cc-e686ed90b6b5 +# ╠═d1e5947b-2c1a-46fc-ab8f-feeba03453e7 +# ╟─c9ebe984-4630-47c1-a941-795293f5b3c1 +# ╟─a3e81f3f-203b-47b7-ac60-b4267eddfad4 +# ╠═2315ffec-dc49-413a-b0d6-1bcce2addd76 +# ╟─d2bd2f13-ed23-42c5-9317-5b48ec3a8bb7 +# ╠═fe6b5b76-8b77-4bfc-a2e8-bcc0b78ad764 +# ╠═f557784e-828c-415e-abb0-964b3a9fe8ef diff --git a/src/AstroImages.jl b/src/AstroImages.jl index 6699e330..a30900b3 100644 --- a/src/AstroImages.jl +++ b/src/AstroImages.jl @@ -1,71 +1,65 @@ -__precompile__() - module AstroImages -using FITSIO, FileIO, Images, Interact, Reproject, WCS, MappedArrays - -export load, AstroImage, ccd2rgb, set_brightness!, set_contrast!, add_label!, reset! +using FITSIO +using FileIO +# Rather than pulling in all of Images.jl, just grab the packages +# we need to extend to our basic functionality. +# We also need ImageShow so that user's images appear automatically. +using ImageBase, ImageShow#, ImageAxes + +using WCS +using Statistics +using MappedArrays +using ColorSchemes +using DimensionalData +using Tables +using RecipesBase +using AstroAngles +using Printf +using PlotUtils: PlotUtils +using PlotUtils: optimize_ticks, AbstractColorList +using UUIDs # can remove once reigstered with FileIO + + +export load, + save, + AstroImage, + AstroImageVec, + AstroImageMat, + WCSGrid, + composecolors, + Zscale, + Percent, + logstretch, + powstretch, + sqrtstretch, + squarestretch, + asinhstretch, + sinhstretch, + powerdiststretch, + imview, + render, # deprecated + header, + copyheader, + shareheader, + wcs, + Comment, + History, + # Dimensions + Centered, + Spec, + Pol, + Ti, + X, Y, Z, Dim, + At, Near, .., + dims, refdims, + recenter, + pix_to_world, + pix_to_world!, + world_to_pix, + world_to_pix! -_load(fits::FITS, ext::Int) = read(fits[ext]) -_load(fits::FITS, ext::NTuple{N, Int}) where {N} = ntuple(i-> read(fits[ext[i]]), N) -_load(fits::NTuple{N, FITS}, ext::NTuple{N, Int}) where {N} = ntuple(i -> _load(fits[i], ext[i]), N) - -_header(fits::FITS, ext::Int) = WCS.from_header(read_header(fits[ext], String))[1] -_header(fits::FITS, ext::NTuple{N, Int}) where {N} = - ntuple(i -> WCS.from_header(read_header(fits[ext[i]], String))[1], N) -_header(fits::NTuple{N, FITS}, ext::NTuple{N, Int}) where {N} = - ntuple(i -> _header(fits[i], ext[i]), N) -""" - load(fitsfile::String, n=1) -Read and return the data from `n`-th extension of the FITS file. - -Second argument can also be a tuple of integers, in which case a -tuple with the data of each corresponding extension is returned. -""" -function FileIO.load(f::File{format"FITS"}, ext::Int=1) - fits = FITS(f.filename) - out = _load(fits, ext) - header = _header(fits,ext) - close(fits) - return out, header -end - -function FileIO.load(f::File{format"FITS"}, ext::NTuple{N,Int}) where {N} - fits = FITS(f.filename) - out = _load(fits, ext) - header = _header(fits, ext) - close(fits) - return out, header -end - -function FileIO.load(f::NTuple{N, String}) where {N} - fits = ntuple(i-> FITS(f[i]), N) - ext = indexer(fits) - out = _load(fits, ext) - header = _header(fits, ext) - for i in 1:N - close(fits[i]) - end - return out, header -end - -function indexer(fits::FITS) - ext = 0 - for (i, hdu) in enumerate(fits) - if hdu isa ImageHDU && length(size(hdu)) >= 2 # check if Image is atleast 2D - ext = i - break - end - end - if ext > 1 - @info "Image was loaded from HDU $ext" - elseif ext == 0 - error("There are no ImageHDU extensions in '$(fits.filename)'") - end - return ext -end -indexer(fits::NTuple{N, FITS}) where {N} = ntuple(i -> indexer(fits[i]), N) # Images.jl expects data to be either a float or a fixed-point number. Here we define some # utilities to convert all data types supported by FITS format to float or fixed-point: @@ -85,163 +79,478 @@ for n in (8, 16, 32, 64) end end -mutable struct Properties{P <: Union{AbstractFloat, FixedPoint}} - rgb_image::MappedArrays.MultiMappedArray{RGB{P},2,Tuple{Array{P,2},Array{P,2},Array{P,2}},Type{RGB{P}},typeof(ImageCore.extractchannels)} - contrast::Float64 - brightness::Float64 - label::Array{Tuple{Tuple{Float64,Float64},String},1} - function Properties{P}(;kvs...) where P - obj = new{P}() - obj.contrast = 1.0 - obj.brightness = 0.0 - obj.label = Array{Tuple{Tuple{Float64,Float64},String}}(undef,0) - for (k,v) in kvs - setproperty!(obj, k, v) - end - return obj - end + +""" +Provides access to a FITS image along with its accompanying +header and WCS information, if applicable. +""" +struct AstroImage{T,N,D<:Tuple,R<:Tuple,A<:AbstractArray{T,N},W<:Tuple} <: AbstractDimArray{T,N,D,A} + # Parent array we are wrapping + data::A + # Fields for DimensionalData + dims::D + refdims::R + # FITS Heads beloning to this image, if any + header::FITSHeader + # cached WCSTransform objects for this data. + wcs::Vector{WCSTransform} + # A flag that is set when a user modifies a WCS header. + # The next access to the wcs object will regenerate from + # the new header on demand. + wcs_stale::Base.RefValue{Bool} + # Correspondance between dims & refdims -> WCS Axis numbers + wcsdims::W end +# Provide type aliases for 1D and 2D versions of our data structure. +const AstroImageVec{T,D} = AstroImage{T,1} where {T} +const AstroImageMat{T,D} = AstroImage{T,2} where {T} + + +""" + Centered() -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} +Pass centered as a dimesion range to automatically center a dimension +along that axis. + +Example: +```julia +cube = load("abc.fits", (X=Centered(), Y=Centered(), Pol=[:I, :Q, :U])) +``` + +In that case, cube will have dimsions with the centre of the image at 0 +in both the X and Y axes. +""" +struct Centered end + +# Default dimension names if none are provided +const dimnames = ( + X, Y, Z, + (Dim{i} for i in 4:10)... +) + +const Spec = Dim{:Spec} +const Pol = Dim{:Pol} + +""" + wcsax(img, dim) + +Return the WCS axis number associated with a dimension, even if the image +has been slices or otherwise transformed. + +Internally, the order is stored in the field `wcsdims`. +""" +function wcsax(img::AstroImage, dim) + return findfirst(di->name(di)==name(dim), img.wcsdims) end +# Accessors +""" + header(img::AstroImage) + +Return the underlying FITSIO.FITSHeader object wrapped by an AstroImage. +Note that this object has less flexible getindex and setindex methods. +Indexing by symbol, Comment, History, etc are not supported. +""" +header(img::AstroImage) = getfield(img, :header) """ - AstroImage([color=Gray,] data::Matrix{Real}) - AstroImage(color::Type{<:Color}, data::NTuple{N, Matrix{T}}) where {T<:Real, N} + header(array::AbstractArray) -Construct an `AstroImage` object of `data`, using `color` as color map, `Gray` by default. +Returns an empty FITSIO.FITSHeader object when called with a non-AstroImage +abstract array. """ -AstroImage(color::Type{<:Color}, data::Matrix{T}, wcs::WCSTransform) where {T<:Real} = - AstroImage{T,color, 1, Float64}((data,), (extrema(data),), (wcs,), Properties{Float64}()) -function AstroImage(color::Type{<:AbstractRGB}, data::NTuple{N, Matrix{T}}, wcs::NTuple{N, WCSTransform}) where {T <: Union{AbstractFloat, FixedPoint}, N} - if N == 3 - img = ccd2rgb((data[1], wcs[1]),(data[2], wcs[2]),(data[3], wcs[3])) - return AstroImage{T,color,N, widen(T)}(data, ntuple(i -> extrema(data[i]), N), wcs, Properties{widen(T)}(rgb_image = img)) +header(::AbstractArray) = emptyheader() + +""" + wcs(img) + +Computes and returns a list of World Coordinate System WCSTransform objects from WCS.jl. +The resultss are cached after the first call, so subsequent calls are fast. +Modifying a WCS header invalidates this cache automatically, so users should call `wcs(...)` +each time rather than keeping the WCSTransform object around. +""" +function wcs(img::AstroImage) + if getfield(img, :wcs_stale)[] + empty!(getfield(img, :wcs)) + append!(getfield(img, :wcs), wcsfromheader(img)) + getfield(img, :wcs_stale)[] = false end + return getfield(img, :wcs) end -function AstroImage(color::Type{<:AbstractRGB}, data::NTuple{N, Matrix{T}}, wcs::NTuple{N, WCSTransform}) where {T<:Real, N} - if N == 3 - img = ccd2rgb((data[1], wcs[1]),(data[2], wcs[2]),(data[3], wcs[3])) - return AstroImage{T,color,N, Float64}(data, ntuple(i -> extrema(data[i]), N), wcs, Properties{Float64}(rgb_image = img)) - end +""" + wcs(img, index) + +Computes and returns a World Coordinate System WCSTransform objects from WCS.jl by index. +This is to support files with multiple WCS transforms specified. +`wcs(img,1)` is useful for selecting selecting the first WCSTranform object. +The resultss are cached after the first call, so subsequent calls are fast. +Modifying a WCS header invalidates this cache automatically, so users should call `wcs(...)` +each time rather than keeping the WCSTransform object around. +""" +wcs(img, ind) = wcs(img)[ind] +""" +wcs(array) + +Returns a list with a single basic WCSTransform object when called with a non-AstroImage +abstract array. +""" +wcs(arr::AbstractArray) = [emptywcs(arr)] + +# Implement DimensionalData interface +Base.parent(img::AstroImage) = getfield(img, :data) +DimensionalData.dims(A::AstroImage) = getfield(A, :dims) +DimensionalData.refdims(A::AstroImage) = getfield(A, :refdims) +DimensionalData.data(A::AstroImage) = getfield(A, :data) +DimensionalData.name(::AstroImage) = DimensionalData.NoName() +DimensionalData.metadata(::AstroImage) = DimensionalData.Dimensions.LookupArrays.NoMetadata() + +@inline function DimensionalData.rebuild( + img::AstroImage, + data, + # Fields for DimensionalData + dims::Tuple=DimensionalData.dims(img), + refdims::Tuple=DimensionalData.refdims(img), + name::Union{Symbol,DimensionalData.AbstractName,Nothing}=nothing, + metadata::Union{DimensionalData.LookupArrays.AbstractMetadata,Nothing}=nothing, + # FITS Header beloning to this image, if any + header::FITSHeader=deepcopy(header(img)), + # A cached WCSTransform object for this data + wcs::Vector{WCSTransform}=getfield(img, :wcs), + wcs_stale::Bool=getfield(img, :wcs_stale)[], + wcsdims::Tuple=(dims...,refdims...), +) + return AstroImage(data, dims, refdims, header, wcs, Ref(wcs_stale), wcsdims) +end +# Keyword argument version. +# We have to define this since our struct contains additional field names. +@inline function DimensionalData.rebuild( + img::AstroImage; + data, + # Fields for DimensionalData + dims::Tuple=DimensionalData.dims(img), + refdims::Tuple=DimensionalData.refdims(img), + name::Union{Symbol,DimensionalData.AbstractName,Nothing}=nothing, + metadata::Union{DimensionalData.LookupArrays.AbstractMetadata,Nothing}=nothing, + # FITS Header beloning to this image, if any + header::FITSHeader=deepcopy(header(img)), + # A cached WCSTransform object for this data + wcs::Vector{WCSTransform}=getfield(img, :wcs), + wcs_stale::Bool=getfield(img, :wcs_stale)[], + wcsdims::Tuple=(dims...,refdims...), +) + return AstroImage(data, dims, refdims, header, wcs, Ref(wcs_stale), wcsdims) end -function AstroImage(color::Type{<:Color}, data::NTuple{N, Matrix{T}}, wcs::NTuple{N, WCSTransform}) where {T<:Real, N} - return AstroImage{T,color, N, Float64}(data, ntuple(i -> extrema(data[i]), N), wcs, Properties{Float64}()) +@inline DimensionalData.rebuildsliced( + f::Function, + img::AstroImage, + data, + I, + header=deepcopy(header(img)), + wcs=getfield(img, :wcs), + wcs_stale=getfield(img, :wcs_stale)[], + wcsdims=getfield(img, :wcsdims), +) = rebuild(img, data, DimensionalData.slicedims(f, img, I)..., nothing, nothing, header, wcs, wcs_stale, wcsdims) + +# Return result wrapped in AstroImage +# For these functions that return lazy wrappers, we want to share header +for f in [ + :(Base.adjoint), + :(Base.transpose), + :(Base.view) +] + # TODO: these functions are copying headers + @eval ($f)(img::AstroImage) = shareheader(img, $f(parent(img))) end -AstroImage(data::Matrix{T}) where {T<:Real} = AstroImage{T,Gray,1, Float64}((data,), (extrema(data),), (WCSTransform(2),), Properties{Float64}()) -AstroImage(data::NTuple{N, Matrix{T}}) where {T<:Real, N} = AstroImage{T,Gray,N, Float64}(data, ntuple(i -> extrema(data[i]), N), ntuple(i-> WCSTransform(2), N), Properties{Float64}()) -AstroImage(data::Matrix{T}, wcs::WCSTransform) where {T<:Real} = AstroImage{T,Gray,1, Float64}((data,), (extrema(data),), (wcs,), Properties{Float64}()) -AstroImage(data::NTuple{N, Matrix{T}}, wcs::NTuple{N, WCSTransform}) where {T<:Real, N} = AstroImage{T,Gray,N, Float64}(data, ntuple(i -> extrema(data[i]), N), wcs, Properties{Float64}()) -""" - AstroImage([color=Gray,] filename::String, n::Int=1) - AstroImage(color::Type{<:Color}, file::String, n::NTuple{N, Int}) where {N} -Create an `AstroImage` object by reading the `n`-th extension from FITS file `filename`. +""" + AstroImage(img::AstroImage) -Use `color` as color map, this is `Gray` by default. +Returns its argument. Useful to ensure an argument is converted to an +AstroImage before proceeding. """ -AstroImage(color::Type{<:Color}, file::String, ext::Int) = - AstroImage(color, file, (ext,)) -AstroImage(color::Type{<:Color}, file::String, ext::NTuple{N, Int}) where {N} = - AstroImage(color, load(file, ext)...) +AstroImage(img::AstroImage) = img -AstroImage(file::String, ext::Int) = AstroImage(Gray, file, ext) -AstroImage(file::String, ext::NTuple{N, Int}) where {N} = AstroImage(Gray, file, ext) -AstroImage(color::Type{<:Color}, fits::FITS, ext::Int) = - AstroImage(color, _load(fits, ext), WCS.from_header(read_header(fits[ext],String))[1]) -AstroImage(color::Type{<:Color}, fits::FITS, ext::NTuple{N, Int}) where {N} = - AstroImage(color, _load(fits, ext), ntuple(i -> WCS.from_header(read_header(fits[ext[i]], String))[1], N)) -AstroImage(color::Type{<:Color}, fits::NTuple{N, FITS}, ext::NTuple{N, Int}) where {N} = - AstroImage(color, ntuple(i -> _load(fits[i], ext[i]), N), ntuple(i -> WCS.from_header(read_header(fits[i][ext[i]], String))[1], N)) +""" + AstroImage(data::AbstractArray, [header::FITSHeader,] [wcs::WCSTransform,]) -AstroImage(files::NTuple{N,String}) where {N} = - AstroImage(Gray, load(files)...) -AstroImage(color::Type{<:Color}, files::NTuple{N,String}) where {N} = - AstroImage(color, load(files)...) -AstroImage(file::String) = AstroImage((file,)) +Create an AstroImage from an array, and optionally header or header and a +WCSTransform. +""" +function AstroImage( + data::AbstractArray{T,N}, + dims::Union{Tuple,NamedTuple}=(), + refdims::Union{Tuple,NamedTuple}=(), + header::FITSHeader=emptyheader(), + wcs::Union{WCSTransform,Nothing}=nothing; + wcsdims=nothing +) where {T, N} + wcs_stale = isnothing(wcs) + if isnothing(wcs) + wcs = [emptywcs(data)] + end + # If the user passes in a WCSTransform of their own, we use it and mark + # wcs_stale=false. It will be kept unless they manually change a WCS header. + # If they don't pass anything, we start with empty WCS information regardless + # of what's in the header but we mark it as stale. + # If/when the WCS info is accessed via `wcs(img)` it will be computed and cached. + # This avoids those computations if the WCS transform is not needed. + # It also allows us to create images with invalid WCS header, + # only erroring when/if they are used. + + # Fields for DimensionalData. + # TODO: cleanup logic + if dims == () + # if wcsdims + # ourdims = Tuple(Wcs{i} for i in 1:ndims(data)) + # else + ourdims = dimnames[1:ndims(data)] + # end + dims = map(ourdims, axes(data)) do dim, ax + dim(ax) + end + end + # Replace any occurences of Centered() with an automatic range + # from the data. + dimvals = map(dims, axes(data)) do dim, ax + if dim isa Centered + ax .- mean(ax) + else + dim + end + end + if dims isa NamedTuple + dims = NamedTuple{keys(dims)}(dimvals) + elseif !(dims isa NTuple{N,Dimensions.Dimension} where N) && + !(all(d-> d isa Union{UnionAll,DataType} && d <: Dimensions.Dimension, dims)) + k = name.(dimnames[1:ndims(data)]) + dims = NamedTuple{k}(dimvals) + end + dims = DimensionalData.format(dims, data) + if length(dims) != ndims(data) + error("Number of dims does not match the shape of the data") + end -# 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]))) + if isnothing(wcsdims) + wcsdims = (dims...,refdims...) + end + + return AstroImage(data, dims, refdims, header, wcs, Ref(wcs_stale), wcsdims) end +function AstroImage( + darr::AbstractDimArray, + header::FITSHeader=emptyheader(), + wcs::Union{Vector{WCSTransform},Nothing}=nothing; +) + wcs_stale = isnothing(wcs) + if isnothing(wcs) + wcs = [emptywcs(darr)] + end + wcsdims = (dims(darr)..., refdims(darr)...) + return AstroImage(parent(darr), dims(darr), refdims(darr), header, wcs, Ref(wcs_stale), wcsdims) +end +AstroImage( + data::AbstractArray, + dims::Union{Tuple,NamedTuple}, + header::FITSHeader, + wcs::Union{Vector{WCSTransform},Nothing}=nothing; +) = AstroImage(data, dims, (), header, wcs) +AstroImage( + data::AbstractArray, + header::FITSHeader, + wcs::Union{Vector{WCSTransform},Nothing}=nothing; +) = AstroImage(data, (), (), header, wcs) + + +# TODO: ensure this gets WCS dims. +AstroImage(data::AbstractArray, wcs::Vector{WCSTransform}) = AstroImage(data, emptyheader(), wcs) -""" - set_brightness!(img::AstroImage, value::AbstractFloat) -Sets brightness of `rgb_image` to value. """ -function set_brightness!(img::AstroImage, value::AbstractFloat) - if isdefined(img.property, :rgb_image) - diff = value - img.property.brightness - img.property.brightness = value - img.property.rgb_image .+= RGB{typeof(value)}(diff, diff, diff) - else - throw(DomainError(value, "Can't apply operation. AstroImage dosen't contain :rgb_image")) - end -end +Index for accessing a comment associated with a header keyword +or COMMENT entry. + +Example: +``` +img = AstroImage(randn(10,10)) +img["ABC"] = 1 +img["ABC", Comment] = "A comment describing this key" + +push!(img, Comment, "The purpose of this file is to demonstrate comments") +img[Comment] # ["The purpose of this file is to demonstrate comments")] +``` +""" +struct Comment end """ - set_contrast!(img::AstroImage, value::AbstractFloat) +Allows accessing and setting HISTORY header entries -Sets contrast of rgb_image to value. +img = AstroImage(randn(10,10)) +push!(img, History, "2023-04-19: Added history entry.") +img[History] # ["2023-04-19: Added history entry."] """ -function set_contrast!(img::AstroImage, value::AbstractFloat) - if isdefined(img.property, :rgb_image) - diff = (value / img.property.contrast) - img.property.contrast = value - img.property.rgb_image = colorview(RGB, red.(img.property.rgb_image) .* diff, green.(img.property.rgb_image) .* diff, - blue.(img.property.rgb_image) .* diff) - else - throw(DomainError(value, "Can't apply operation. AstroImage dosen't contain :rgb_image")) +struct History end + + +# We might want getproperty for header access in future. +# function Base.getproperty(img::AstroImage, ::Symbol) +# error("getproperty reserved for future use.") +# end + +# Getting and setting comments +Base.getindex(img::AstroImage, inds::AbstractString...) = getindex(header(img), inds...) # accesing header using strings +function Base.setindex!(img::AstroImage, v, ind::AbstractString) # modifying header using a string + setindex!(header(img), v, ind) + # Mark the WCS object as being out of date if this was a WCS header keyword + if ind ∈ WCS_HEADERS + getfield(img, :wcs_stale)[] = true end end +Base.getindex(img::AstroImage, inds::Symbol...) = getindex(img, string.(inds)...) # accessing header using symbol +Base.setindex!(img::AstroImage, v, ind::Symbol) = setindex!(img, v, string(ind)) +Base.getindex(img::AstroImage, ind::AbstractString, ::Type{Comment}) = get_comment(header(img), ind) # accesing header comment using strings +Base.setindex!(img::AstroImage, v, ind::AbstractString, ::Type{Comment}) = set_comment!(header(img), ind, v) # modifying header comment using strings +Base.getindex(img::AstroImage, ind::Symbol, ::Type{Comment}) = get_comment(header(img), string(ind)) # accessing header comment using symbol +Base.setindex!(img::AstroImage, v, ind::Symbol, ::Type{Comment}) = set_comment!(header(img), string(ind), v) # modifying header comment using Symbol + +# Support for special HISTORY and COMMENT entries +function Base.getindex(img::AstroImage, ::Type{History}) + hdr = header(img) + ii = findall(==("HISTORY"), hdr.keys) + return view(hdr.comments, ii) +end +function Base.getindex(img::AstroImage, ::Type{Comment}) + hdr = header(img) + ii = findall(==("COMMENT"), hdr.keys) + return view(hdr.comments, ii) +end +# Adding new comment and history entries +function Base.push!(img::AstroImage, ::Type{Comment}, history::AbstractString) + hdr = header(img) + push!(hdr.keys, "COMMENT") + push!(hdr.values, nothing) + push!(hdr.comments, history) +end +function Base.push!(img::AstroImage, ::Type{History}, history::AbstractString) + hdr = header(img) + push!(hdr.keys, "HISTORY") + push!(hdr.values, nothing) + push!(hdr.comments, history) +end """ - add_label!(img::AstroImage, x::Real, y::Real, label::String) + copyheader(img::AstroImage, data) -> imgnew +Create a new image copying the header of `img` but +using the data of the AbstractArray `data`. Note that changing the +header of `imgnew` does not affect the header of `img`. +See also: [`shareheader`](@ref). +""" +copyheader(img::AstroImage, data::AbstractArray) = + AstroImage(data, dims(img), refdims(img), deepcopy(header(img)), copy(getfield(img, :wcs)), Ref(getfield(img, :wcs_stale)[]), getfield(img,:wcsdims)) -Stores label to coordinates (x,y) in AstroImage's property label. """ -function add_label!(img::AstroImage, x::Real, y::Real, label::String) - push!(img.property.label, ((x,y), label)) -end + shareheader(img::AstroImage, data) -> imgnew +Create a new image reusing the header dictionary of `img` but +using the data of the AbstractArray `data`. The two images have +synchronized header; modifying one also affects the other. +See also: [`copyheader`](@ref). +""" +shareheader(img::AstroImage, data::AbstractArray) = AstroImage(data, dims(img), refdims(img), header(img), getfield(img, :wcs), Ref(getfield(img, :wcs_stale)[]), getfield(img,:wcsdims)) +# Share header if an AstroImage, do nothing if AbstractArray +maybe_shareheader(img::AstroImage, data) = shareheader(img, data) +maybe_shareheader(::AbstractArray, data) = data +maybe_copyheader(img::AstroImage, data) = copyheader(img, data) +maybe_copyheader(::AbstractArray, data) = data + + +Base.promote_rule(::Type{AstroImage{T}}, ::Type{AstroImage{V}}) where {T,V} = AstroImage{promote_type{T,V}} + + + +Base.copy(img::AstroImage) = rebuild(img, copy(parent(img))) +Base.convert(::Type{AstroImage}, A::AstroImage) = A +Base.convert(::Type{AstroImage}, A::AbstractArray) = AstroImage(A) +Base.convert(::Type{AstroImage{T}}, A::AstroImage{T}) where {T} = A +Base.convert(::Type{AstroImage{T}}, A::AstroImage) where {T} = shareheader(A, convert(AbstractArray{T}, parent(A))) +Base.convert(::Type{AstroImage{T}}, A::AbstractArray{T}) where {T} = AstroImage(A) +Base.convert(::Type{AstroImage{T}}, A::AbstractArray) where {T} = AstroImage(convert(AbstractArray{T}, A)) +Base.convert(::Type{AstroImage{T,N,D,R,AT}}, A::AbstractArray{T,N}) where {T,N,D,R,AT} = AstroImage(convert(AbstractArray{T}, A)) + +# TODO: share headers in View. Needs support from DimensionalData. """ - reset!(img::AstroImage) + emptyheader() + +Convenience function to create a FITSHeader with no keywords set. +""" +emptyheader() = FITSHeader(String[],[],String[]) -Resets AstroImage property fields. -Sets brightness to 0.0, contrast to 1.0, empties label -and form a fresh rgb_image without any brightness, contrast operations on it. """ -function reset!(img::AstroImage{T,C,N}) where {T,C,N} - img.property.contrast = 1.0 - img.property.brightness = 0.0 - img.property.label = [] - if N == 3 && C == RGB - shape_out = size(img.property.rgb_image) - img.property.rgb_image = ccd2rgb((img.data[1], img.wcs[1]),(img.data[2], img.wcs[2]),(img.data[3], img.wcs[3]), - shape_out = shape_out) - end -end + recenter(img::AstroImage, newcentx, newcenty, ...) + +Adjust the dimensions of an AstroImage so that they are centered on the pixel locations given by `newcentx`, .. etc. -Images.colorview(img::AstroImage) = render(img) +This does not affect the underlying array, it just updates the dimensions associated with it by the AstroImage. -Base.size(img::AstroImage) = Base.size.(img.data) +Example: +```julia +a = AstroImage(randn(11,11)) +a[1,1] # Bottom left +a[At(1),At(1)] # Bottom left +r = recenter(a, 6, 6) +r[1,1] # Still bottom left +r[At(1),At(1)] # Center pixel +``` +""" +function recenter(img::AstroImage, centers::Number...) + newdims = map(dims(img), axes(img), centers) do d, a, c + return AstroImages.name(d) => a .- c + end + newdimsformatted = AstroImages.DimensionalData.format(NamedTuple(newdims), parent(img)) + l = length(newdimsformatted) + if l < ndims(img) + newdimsformatted = (newdimsformatted..., dims(img)[l+1:end]...) + end + AstroImages.rebuild(img, parent(img), newdimsformatted) +end -Base.length(img::AstroImage{T,C,N}) where {T,C,N} = N +include("wcs.jl") +include("io.jl") +include("imview.jl") include("showmime.jl") include("plot-recipes.jl") + +include("contrib/images.jl") +include("contrib/abstract-ffts.jl") +# include("contrib/reproject.jl") + include("ccd2rgb.jl") +include("precompile.jl") + +function __init__() + + # You can only `imview` 2D slices. Add an error hint if the user + # tries to display a cube. + if isdefined(Base.Experimental, :register_error_hint) + Base.Experimental.register_error_hint(MethodError) do io, exc, argtypes, kwargs + if exc.f == imview && first(argtypes) <: AbstractArray && ndims(first(argtypes)) != 2 + print(io, "\nThe `imview` function only supports 2D images. If you have a cube, try viewing one slice at a time: imview(cube[:,:,1])\n") + end + end + end + + # 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")] + ) + +end end # module diff --git a/src/ccd2rgb.jl b/src/ccd2rgb.jl index 8ecf15a7..835836b3 100644 --- a/src/ccd2rgb.jl +++ b/src/ccd2rgb.jl @@ -1,48 +1,77 @@ """ - ccd2rgb(red::ImageHDU, green::ImageHDU, blue::ImageHDU; stretch = identity, shape_out = size(red)) - ccd2rgb(red::Tuple{AbstractMatrix, WCSTransform}, green::Tuple{AbstractMatrix, WCSTransform}, - blue::Tuple{AbstractMatrix, WCSTransform}; stretch = identity, shape_out = size(red[1])) + composecolors( + images, + cmap=["#F00", "#0F0", "#00F"]; + clims, + stretch, + contrast, + bias, + multiplier + ) -Converts 3 grayscale ImageHDU into RGB by reprojecting them. +Create a color composite of multiple images by applying `imview` and blending +the results. This function can be used to create RGB composites using any number of channels +(e.g. red, green, blue, and hydrogen alpha) as well as more exotic images like blending +radio and optical data using two different colormaps. -# Arguments -- `red`: Red channel data. -- `green`: Green channel data. -- `blue`: Blue channel data. -- `stretch`: Stretch function applied. -- `shape_out`: Shape of output RGB image. +`cmap` should be a list of colorants, named colors (see Colors.jl), or colorschemes (see ColorSchemes.jl). +`clims`, `stretch`, `contrast`, and `bias` are passed on to `imview`. They can be a single value or +a list of different values for each image. -# Examples -```julia-repl -julia> ccd2rgb(r, b, g, shape_out = (1000,1000)) - -julia> ccd2rgb(r, b, g, shape_out = (1000,1000), stretch = log) - -julia> ccd2rgb(r, b, g, shape_out = (1000,1000), stretch = sqrt) - -julia> ccd2rgb(r, b, g, shape_out = (1000,1000), stretch = asinh) +Examples: +```julia +# Basic RGB +composecolors([redimage, greenimage, blueimage]) +# Non-linear stretch before blending +composecolors([redimage, greenimage, blueimage], stretch=asinhstretch) +# More than three channels are allowed (H alpha in pink) +composecolors( + [antred, antgreen, antblue, anthalp], + ["red", "green", "blue", "maroon1"], + multiplier=[1,2,1,1] +) +# Can mix +composecolors([radioimage, xrayimage], [:ice, :magma], clims=extrema) +composecolors([radioimage, xrayimage], [:magma, :viridis], clims=[Percent(99), Zscale()]) ``` """ -function ccd2rgb(red::Tuple{AbstractMatrix, WCSTransform}, green::Tuple{AbstractMatrix, WCSTransform}, - blue::Tuple{AbstractMatrix, WCSTransform}; stretch = identity, shape_out = size(red[1])) - red_rp = reproject(red, red[2], shape_out = shape_out)[1] - green_rp = reproject(green, red[2], shape_out = shape_out)[1] - blue_rp = reproject(blue, red[2], shape_out = shape_out)[1] - - I = (red_rp .+ green_rp .+ blue_rp) ./ 3 - I .= (x -> stretch(x)/x).(I) - - red_rp .*= I - green_rp .*= I - blue_rp .*= I - - m1 = maximum(x->isnan(x) ? -Inf : x, red_rp) - m2 = maximum(x->isnan(x) ? -Inf : x, green_rp) - m3 = maximum(x->isnan(x) ? -Inf : x, blue_rp) - return colorview(RGB, red_rp./m1 , green_rp./m2, blue_rp./m3) -end +function composecolors( + images, + cmap=nothing; + clims=_default_clims[], + stretch=_default_stretch[], + contrast=1.0, + bias=0.5, + multiplier=1.0 +) + if isempty(images) + error("At least one image is required.") + end + if !allequal(size.(images)) + error("Images must have the same dimensions to compose them.") + end + if length(images) == 3 && isnothing(cmap) + cmap = ["red", "green", "blue"] + end + if length(cmap) < length(images) + error("Please provide a color channel for each image") + end -ccd2rgb(red::ImageHDU, green::ImageHDU, blue::ImageHDU; stretch = identity, shape_out = size(red)) = - ccd2rgb((read(red), WCS.from_header(read_header(red, String))[1]), (read(green), WCS.from_header(read_header(green, String))[1]), - (read(blue), WCS.from_header(read_header(blue, String))[1]), stretch = stretch, shape_out = shape_out) + # Use imview to render each channel to RGBA + images_rendered = broadcast(images, cmap, clims, stretch, contrast, bias) do image, cmap, clims, stretch, contrast, bias + imview(image; cmap, clims, stretch, contrast, bias) + end + # Now blend, ensuring each color channel never exceeds [0,1] + combined = mappedarray(images_rendered...) do channels... + pxblended = sum(channels .* multiplier) + return typeof(pxblended)( + clamp(pxblended.r,0,1), + clamp(pxblended.g,0,1), + clamp(pxblended.b,0,1), + clamp(pxblended.alpha,0,1) + ) + end + return combined +end +export composechannels diff --git a/src/contrib/abstract-ffts.jl b/src/contrib/abstract-ffts.jl new file mode 100644 index 00000000..51118018 --- /dev/null +++ b/src/contrib/abstract-ffts.jl @@ -0,0 +1,27 @@ +using AbstractFFTs + +for f in [ + :(AbstractFFTs.fft), + :(AbstractFFTs.bfft), + :(AbstractFFTs.ifft), + :(AbstractFFTs.fft!), + :(AbstractFFTs.bfft!), + :(AbstractFFTs.ifft!), + :(AbstractFFTs.rfft), +] + # TODO: should we try to alter the image headers to change the units? + @eval ($f)(img::AstroImage, args...; kwargs...) = copyheader(img, $f(parent(img), args...; kwargs...)) +end + +for f in [ + :(AbstractFFTs.fftshift), +] + # TODO: should we try to alter the image headers to change the units? + @eval ($f)(img::AstroImage, args...; kwargs...) = shareheader(img, $f(parent(img), args...; kwargs...)) +end + + + + +# AbstractFFTs.complexfloat(img::AstroImage{T,N,D,R,StridedArray{T}}) where {T<:Complex{<:BlasReal}} = img +# AbstractFFTs.realfloat(img::AstroImage{T,N,D,R,StridedArray{T}}) where {T<:BlasReal} = img diff --git a/src/contrib/images.jl b/src/contrib/images.jl new file mode 100644 index 00000000..e66e9c76 --- /dev/null +++ b/src/contrib/images.jl @@ -0,0 +1,83 @@ +# function warp(img::AstroImageMat, args...; kwargs...) +# out = warp(parentt(img), args...; kwargs...) +# return copyheaders(img, out) +# end + + + +# Instead of using a datatype like N0f32 to interpret integers as fixed point values in [0,1], +# we use a mappedarray to map the native data range (regardless of type) to [0,1] +ImageBase.normedview(img::AstroImageMat{<:FixedPoint}) = img +function ImageBase.normedview(img::AstroImageMat{T}) where T + imgmin, imgmax = extrema(skipmissingnan(img)) + Δ = abs(imgmax - imgmin) + # Do not introduce NaNs if limits are identical + if Δ == 0 + Δ = oneunit(imgmin) + end + normeddata = mappedarray( + pix -> (pix - imgmin)/Δ, + pix_norm -> convert(T, pix_norm*Δ + imgmin), + img + ) + return shareheader(img, normeddata) +end + +""" + clampednormedview(arr, (min, max)) + +Given an AbstractArray and limits `min,max` return a view of the array +where data between [min, max] are scaled to [0, 1] and datat outside that +range are clamped to [0, 1]. + +See also: normedview +""" +function clampednormedview(img::AbstractArray{T}, lims) where T + imgmin, imgmax = lims + Δ = imgmax - imgmin + # Do not introduce NaNs if colorlimits are identical + if Δ == 0 + Δ = oneunit(imgmin) + end + normeddata = mappedarray( + pix -> clamp((pix - imgmin)/Δ, zero(pix), oneunit(pix)), + img + ) + return maybe_shareheader(img, normeddata) +end + +# Restrict downsizes images by roughly a factor of two. +# We want to keep the wrapper but downsize the underlying array +# TODO: correct dimensions after restrict. +ImageBase.restrict(img::AstroImage, ::Tuple{}) = img +function ImageBase.restrict(img::AstroImage, region::Dims) + restricted = restrict(parent(img), region) + steps = cld.(size(img), size(restricted)) + newdims = Tuple(d[begin:s:end] for (d,s) in zip(dims(img),steps)) + return rebuild(img, restricted, newdims) +end + + +ImageBase.pixelspacing(img::AstroImage) = step.(dims(img)) + + +# ImageContrastAdjustment +# function ImageContrastAdjustment.adjust_histogram(::Type{T}, +# img::AstroImageMat, +# f::Images.ImageContrastAdjustment.AbstractHistogramAdjustmentAlgorithm, +# args...; kwargs...) where T +# out = similar(img, axes(img)) +# adjust_histogram!(out, img, f, args...; kwargs...) +# return out +# end + + + + +# """ +# ImageMetadata.parent(img::AstroImage) + +# Returns the underlying wrapped array of `img`. +# """ +# ImageMetadata.parent(img::AstroImage) = getfield(img, :data) + diff --git a/src/contrib/reproject.jl b/src/contrib/reproject.jl new file mode 100644 index 00000000..a4179cbd --- /dev/null +++ b/src/contrib/reproject.jl @@ -0,0 +1,37 @@ + + + +#= +Additional methods to allow Reproject to work. +=# + +using Reproject + +""" + img_proj, mask = reproject(img_in::AstroImageMat, img_out::AstroImageMat) + +Reprojects the AstroImageMat `img_in` to the coordinates of `img_out` +according to the WCS information/headers using interpolation. +""" +function Reproject.reproject(img_in::AstroImageMat, img_out::AstroImageMat) + data_out, mask = reproject(img_in, img_out) + # TODO: should copy the WCS headers from img_out and the remaining + # headers from img_in. + return copyheaders(img_in, data_out) +end + +function Reproject.parse_input_data(input_data::AstroImageMat, hdu) + input_data, input_data.wcs +end +function Reproject.parse_output_projection(output_data::AstroImageMat, hdu) + output_data.wcs, size(output_data) +end +function Reproject.pad_edges(array_in::AstroImageMat{T}) where {T} + image = Matrix{T}(undef, size(array_in)[1] + 2, size(array_in)[2] + 2) + image[2:end-1,2:end-1] = array_in + image[2:end-1,1] = array_in[:,1] + image[2:end-1,end] = array_in[:,end] + image[1,:] = image[2,:] + image[end,:] = image[end-1,:] + return AstroImageMat(image, headers(array_in), wcs(array_in)) +end diff --git a/src/imview.jl b/src/imview.jl new file mode 100644 index 00000000..77153b17 --- /dev/null +++ b/src/imview.jl @@ -0,0 +1,407 @@ +# These reproduce the behaviour of DS9 according to http://ds9.si.edu/doc/ref/how.html + +""" + logstretch(num,a=1000) + +A log-stretch as defined by the SAO DS9 application: http://ds9.si.edu/doc/ref/how.html +""" +logstretch(x,a=1000) = log(a*x+1)/log(a) +""" + powstretch(num, a=1000) + +A power-stretch as defined by the SAO DS9 application: http://ds9.si.edu/doc/ref/how.html +""" +powstretch(x,a=1000) = (a^x - 1)/a +""" + sqrtstretch(num) + +A square root stretch (simply defined as Base.sqrt) +""" +sqrtstretch(x) = sqrt(x) +""" + squarestretch(num) + +A squarestretch-stretch as defined by the SAO DS9 application: http://ds9.si.edu/doc/ref/how.html +""" +squarestretch(x) = x^2 +""" + asinhstretch(num) + +A hyperbolic arcsin stretch as defined by the SAO DS9 application: http://ds9.si.edu/doc/ref/how.html. +""" +asinhstretch(x) = asinh(10x)/3 +""" + sinhstretch(num) + +A hyperbolic sin stretch as defined by the SAO DS9 application: http://ds9.si.edu/doc/ref/how.html +""" +sinhstretch(x) = sinh(3x)/10 +# These additional stretches reproduce behaviour from astropy +""" + logstretch(num,a=1000) + +A power distance stretch as defined by astropy. +""" +powerdiststretch(x, a=1000) = (a^x - 1) / (a - 1) + +""" + Percent(99.5) + +Returns a callable that calculates display limits that include the given +percent of the image data. +Reproduces the behaviour of the SAO DS9 scale menu. + +Example: +```julia +julia> imview(img, clims=Percent(90)) +``` +This will set the limits to be the 5th percentile to the 95th percentile. +""" +struct Percent + perc::Float64 + trim::Float64 + Percent(percentage::Number) = new(Float64(percentage), (1 - percentage/100)/2) +end +(p::Percent)(data::AbstractArray) = quantile(vec(data), (p.trim, 1-p.trim)) +(p::Percent)(data) = p(collect(data)) +Base.broadcastable(p::Percent) = Ref(p) +Base.show(io::IO, p::Percent; kwargs...) = print(io, "Percent($(p.perc))", kwargs...) + +""" + Zscale(options)(data) + +Wraps PlotUtils.zscale in a callable with default parameters. +This is a common algorithm for agressively stretching astronomical data +to see faint structure that originated in IRAF: https://iraf.net/forum/viewtopic.php?showtopic=134139 +but is now seen in many other applications/libraries (DS9, Astropy, etc.) + +Usage: +``` +imview(img, clims=Zscale()) +implot(img, clims=Zscale(contrast=0.1)) +``` + +Default parameters: +``` +nsamples::Int=1000 +contrast::Float64=0.25 +max_reject::Float64=0.5 +min_npixels::Float64=5 +k_rej::Float64=2.5 +max_iterations::Int=5 +``` +""" +Base.@kwdef struct Zscale + nsamples::Int=1000 + contrast::Float64=0.25 + max_reject::Float64=0.5 + min_npixels::Float64=5 + k_rej::Float64=2.5 + max_iterations::Int=5 +end +(z::Zscale)(data::AbstractArray) = PlotUtils.zscale(vec(data), z.nsamples; z.contrast, z.max_reject, z.min_npixels, z.k_rej, z.max_iterations) +(z::Zscale)(data) = z(collect(data)) +Base.show(io::IO, z::Zscale; kwargs...) = print(io, "Zscale()", kwargs...) +Base.broadcastable(z::Zscale) = Ref(z) + +const _default_cmap = Base.RefValue{Union{Symbol,Nothing}}(:magma) +const _default_clims = Base.RefValue{Any}(Percent(99.5)) +const _default_stretch = Base.RefValue{Any}(identity) + +""" + set_cmap!(cmap::Symbol) + set_cmap!(cmap::Nothing) + +Alter the default color map used to display images when using +`imview` or displaying an AstroImageMat. +""" +function set_cmap!(cmap) + # Ensure it's valid + _lookup_cmap(cmap) + _default_cmap[] = cmap +end +""" + set_clims!(clims::Tuple) + set_clims!(clims::Callable) + +Alter the default limits used to display images when using +`imview` or displaying an AstroImageMat. +""" +function set_clims!(clims) + _default_clims[] = clims +end +""" + set_stretch!(stretch::Function) + +Alter the default value stretch functio used to display images when using +`imview` or displaying an AstroImageMat. +""" +function set_stretch!(stretch) + _default_stretch[] = stretch +end + + + +# Helper to iterate over data skipping missing and non-finite values. +skipmissingnan(itr) = Iterators.filter(el->!ismissing(el) && isfinite(el), itr) + +# Convert argument into a colorscheme or AbstractColorList which allow converting +# from numerical data into colors. +function _lookup_cmap(cmap::Symbol) + if cmap ∉ keys(ColorSchemes.colorschemes) + error("$cmap not found in ColorSchemes.colorschemes. See: https://juliagraphics.github.io/ColorSchemes.jl/stable/catalogue/") + end + return ColorSchemes.colorschemes[cmap] +end +_lookup_cmap(::Nothing) = nothing +_lookup_cmap(acl::AbstractColorList) = acl +_lookup_cmap(colorant::Colorant) = PlotUtils.cgrad([:black, colorant]) +_lookup_cmap(colorant::String) = PlotUtils.cgrad([:black, colorant]) + +function _resolve_clims(img::AbstractArray, clims) + # Tuple or abstract array + if typeof(clims) <: AbstractArray || typeof(clims) <: Tuple + if length(clims) != 2 + error("clims must have exactly two values if provided.") + end + imgmin = first(clims) + imgmax = last(clims) + # Or as a callable that computes them given an iterator + else + imgmin, imgmax = clims(skipmissingnan(img)) + end + + return imgmin, imgmax +end + + +""" + imview(img; clims=Percent(99.5), stretch=identity, cmap=:magma, contrast=1.0, bias=0.5) + +Create a read only view of an array or AstroImageMat mapping its data values +to Colors according to `clims`, `stretch`, and `cmap`. + +The data is first clamped to `clims`, which can either be a tuple of (min, max) +values or a function accepting an iterator of pixel values that returns (min, max). +By default, `clims=extrema` i.e. the minimum and maximum of `img`. +Convenient functions to use for `clims` are: +`extrema`, `zscale`, and `percent(p)` + +Next, the data is rescaled to [0,1] and remapped according to the function `stretch`. +Stretch can be any monotonic fuction mapping values in the range [0,1] to some range [a,b]. +Note that `log(0)` is not defined so is not directly supported. +For a list of convenient stretch functions, see: +`logstretch`, `powstretch`, `squarestretch`, `asinhstretch`, `sinhstretch`, `powerdiststretch` + +Finally the data is mapped to RGB values according to `cmap`. If cmap is `nothing`, +grayscale is used. ColorSchemes.jl defines hundreds of colormaps. A few nice ones for +images include: `:viridis`, `:magma`, `:plasma`, `:thermal`, and `:turbo`. + +Crucially, this function returns a view over the underlying data. If `img` is updated +then those changes will be reflected by this view with the exception of `clims` which +is not recalculated. + +Note: if clims or stretch is a function, the pixel values passed in are first filtered +to remove non-finite or missing values. + +### Defaults +The default values of `clims`, `stretch`, and `cmap` are `extrema`, `identity`, and `nothing` +respectively. +You may alter these defaults using `AstroImages.set_clims!`, `AstroImages.set_stretch!`, and +`AstroImages.set_cmap!`. + +### Automatic Display +Arrays wrapped by `AstroImageMat()` get displayed as images automatically by calling +`imview` on them with the default settings when using displays that support showing PNG images. + +### Missing data +Pixels that are `NaN` or `missing` will be displayed as transparent when `cmap` is set +or black if. ++/- Inf will be displayed as black or white respectively. + +### Exporting Images +The view returned by `imview` can be saved using general `FileIO.save` methods. +Example: +```julia +v = imview(data, cmap=:magma, stretch=asinhstretch, clims=percent(95)) +save("output.png", v) +``` +""" +function imview( + img::AbstractArray{T}; + clims=_default_clims[], + stretch=_default_stretch[], + cmap=_default_cmap[], + contrast=1.0, + bias=0.5 +) where {T} + + # Create flipped view of to match conventions of other programs. + # Origin is centre of pixel (1,1) at bottom left. + if ndims(img) == 2 + imgT = view( + permutedims(img,(2,1)), + reverse(axes(img,2)), + :, + ) + elseif ndims(img) >= 3 + newdims = (2,1, 3:ndims(img)...) + ds = Tuple(((:) for _ in 2:ndims(img))) + imgT = view( + permutedims(img,newdims), + reverse(axes(img,2)), + ds..., + ) + else + imgT = img + end + isempt = isempty(imgT) + if isempt + @warn "imview called with empty argument" + return fill(RGBA{N0f8}(0,0,0,0), 1,1) + end + # Users will occaisionally pass in data that is 0D, filled with NaN, or filled with missing. + # We still need to do something reasonable in those caes. + nonempty = any(x-> !ismissing(x) && isfinite(x), imgT) + if !nonempty + @warn "imview called with all missing or non-finite values" + return map(px->RGBA{N0f8}(0,0,0,0), imgT) + end + + imgmin, imgmax = _resolve_clims(imgT, clims) + normed = clampednormedview(imgT, (imgmin, imgmax)) + return _imview(imgT, normed, stretch, _lookup_cmap(cmap), contrast, bias) +end + +# Unwrap AstroImages before view, then rebuild. +# We have to permute the dimensions of the image to get the origin at the bottom left. +# But we don't want this to affect the dimensions of the array. +# Also, this reduces the number of methods we need to compile for imview by standardizing types +# earlier on. The compiled code for showing an array is the same as an array wrapped by an +# AstroImage, except for one unwrapping step. +function imview( + img::AstroImage; + kwargs... +) + return shareheader(img, imview(parent(img); kwargs...)) +end + +# Special handling for complex images +""" + imview(img::AbstractArray{<:Complex}; ...) + +When applied to an image with complex values, display the magnitude +of the pixels using `imview` and display the phase angle as a panel below +using a cyclical color map. +For more customatization, you can create a view like this yourself: +```julia +vcat( + imview(abs.(img)), + imview(angle.(img)), +) +``` +""" +function imview(img::AbstractArray{T}; kwargs...) where {T<:Complex} + mag_view = imview(abs.(img); kwargs...) + angle_view = imview(angle.(img), clims=(-pi, pi), stretch=identity, cmap=:cyclic_mygbm_30_95_c78_n256_s25) + vcat(mag_view,angle_view) +end + +function _imview(img, normed::AbstractArray{T}, stretch, cmap, contrast, bias) where T + + function colormap(pixr, pixn)::RGBA{N0f8} + if ismissing(pixr) || !isfinite(pixr) || ismissing(pixn) || !isfinite(pixn) + # We check pixr in addition to pixn because we want to preserve if the pixels + # are +-Inf + stretched = pixr + else + stretched = (stretch(pixn) - bias)*contrast+0.5 + end + + # We treat NaN/missing values as transparent + pix= if ismissing(stretched) || isnan(stretched) + RGBA{N0f8}(0,0,0,0) + # We treat Inf values as white / -Inf as black + elseif isinf(stretched) + if stretched > 0 + RGBA{N0f8}(1,1,1,1) + else + RGBA{N0f8}(0,0,0,1) + end + else + if isnothing(cmap) + # true/false used as numerical values to prevent unucessary promotion + s = clamp(stretched, false, true) + RGBA{N0f8}(s,s,s,1) + else + # Look up colormap + RGBA{N0f8}(get(cmap, stretched, (false, true))) + end + end + return pix + end + mapper = mappedarray(colormap, img, normed) + + return maybe_copyheader(img, mapper) +end + + +""" + imview_colorbar(img; orientation=:vertical) +Create a colorbar for a given image matching how it is displayed by +`imview`. Returns an image. + +`orientation` can be `:vertical` or `:horizontal`. +""" +function imview_colorbar( + img::AbstractArray; + orientation=:vertical, + clims=_default_clims[], + stretch=_default_stretch[], + cmap=_default_cmap[], + contrast=1, + bias=0.5 +) + imgmin, imgmax = _resolve_clims(img, clims) + cbpixlen = 100 + data = repeat(range(imgmin, imgmax, length=cbpixlen), 1,10) + if orientation == :vertical + data = data' + elseif orientation == :horizontal + data = data + else + error("Unsupported orientation for colorbar \"$orientation\"") + end + + # # Stretch the colors: + # # Construct the image to use as a colorbar + # cbimg = imview(data; clims=(imgmin,imgmax), stretch, cmap, k_min=3) + # # And the colorbar tick locations & labels + # ticks, cbmin, cbmax = optimize_ticks(imgmin, imgmax) + # # Now map these to pixel locations through streching and colorlimits: + # stretchmin = stretch(zero(eltype(data))) + # stretchmax = stretch(one(eltype(data))) + # normedticks = clampednormedview(ticks, (imgmin, imgmax)) + # ticksloc = map(ticks,normedticks) do tick, tickn + # return cbpixlen * tickn + # end + + # Strech the ticks + # Construct the image to use as a colorbar + cbimg = imview(data; clims=(imgmin,imgmax), stretch=identity, cmap, contrast, bias) + # And the colorbar tick locations & labels + ticks, _, _ = optimize_ticks(Float64(imgmin), Float64(imgmax), k_min=3) + # Now map these to pixel locations through streching and colorlimits: + stretchmin = stretch(zero(eltype(data))) + stretchmax = stretch(oneunit(eltype(data))) + normedticks = clampednormedview(ticks, (imgmin, imgmax)) + ticksloc = map(normedticks) do tickn + stretched = stretch(tickn) + stretchednormed = (stretched - stretchmin) * (stretchmax - stretchmin) + return cbpixlen * stretchednormed + end + ticklabels = map(ticks) do t + @sprintf("%4g", t) + end + return cbimg, (ticksloc, ticklabels) +end diff --git a/src/io.jl b/src/io.jl new file mode 100644 index 00000000..651a16c6 --- /dev/null +++ b/src/io.jl @@ -0,0 +1,180 @@ + +""" +AstroImage(fits::FITS, ext::Int=1) + +Given an open FITS file from the FITSIO library, +load the HDU number `ext` as an AstroImage. +""" +AstroImage(fits::FITS, ext::Int=1, args...; kwargs...) = AstroImage(fits[ext], args...; kwargs...) + +""" +AstroImage(hdu::HDU) + +Given an open FITS HDU, load it as an AstroImage. +""" +AstroImage(hdu::HDU, args...; kwargs...) = AstroImage(read(hdu), args..., read_header(hdu); kwargs...) + +""" +img = AstroImage(filename::AbstractString, ext::Integer=1) + +Load an image HDU `ext` from the FITS file at `filename` as an AstroImage. +""" +function AstroImage(filename::AbstractString, ext::Integer, args...; kwargs...) + return FITS(filename, "r") do fits + return AstroImage(fits[ext], args...; kwargs...) + end +end +""" +img1, img2 = AstroImage(filename::AbstractString, exts) + +Load multiple image HDUs `exts` from an FITS file at `filename` as an AstroImage. +`exts` must be a tuple, range, :, or array of Integers. +All listed HDUs in `exts` must be image HDUs or an error will occur. + +Example: +```julia +img1, img2 = AstroImage("abc.fits", (1,3)) # loads the first and third HDU as images. +imgs = AstroImage("abc.fits", 1:3) # loads the first three HDUs as images. +imgs = AstroImage("abc.fits", :) # loads all HDUs as images. +``` +""" +function AstroImage(filename::AbstractString, exts::Union{NTuple{N,<:Integer},AbstractArray{<:Integer}}, args...; kwargs...) where {N} + return FITS(filename, "r") do fits + return map(exts) do ext + return AstroImage(fits[ext], args...; kwargs...) + end + end +end +function AstroImage(filename::AbstractString; kwargs...) where {N} + return FITS(filename, "r") do fits + ext = indexer(fits) + return AstroImage(fits[ext]; kwargs...) + end +end +function AstroImage(filename::AbstractString, ::Colon, args...; kwargs...) where {N} + return FITS(filename, "r") do fits + return map(fits) do hdu + return AstroImage(hdu, args...; kwargs...) + end + end +end + + +""" +load(fitsfile::String) + +Read and return the data from the first ImageHDU in a FITS file +as an AstroImage. If no ImageHDUs are present, an error is returned. + +load(fitsfile::String, ext::Int) + +Read and return the data from the HDU `ext`. If it is an ImageHDU, +as AstroImage is returned. If it is a TableHDU, a plain Julia +column table is returned. + +load(fitsfile::String, :) + +Read and return the data from each HDU in an FITS file. ImageHDUs are +returned as AstroImage, and TableHDUs are returned as column tables. + +load(fitsfile::String, exts::Union{NTuple, AbstractArray}) + +Read and return the data from the HDUs given by `exts`. ImageHDUs are +returned as AstroImage, and TableHDUs are returned as column tables. + +!!! Currently any header on TableHDUs are not supported and are ignored. +""" +function fileio_load(f::File{format"FITS"}, ext::Union{Int,Nothing}=nothing, args...; kwargs...) where {N} + return FITS(f.filename, "r") do fits + if isnothing(ext) + ext = indexer(fits) + end + _loadhdu(fits[ext], args...; kwargs...) + end +end +function fileio_load(f::File{format"FITS"}, exts::Union{NTuple{N,<:Integer},AbstractArray{<:Integer}}, args...; kwargs...) where {N} + return FITS(f.filename, "r") do fits + map(exts) do ext + _loadhdu(fits[ext], args...; kwargs...) + end + end +end +function fileio_load(f::File{format"FITS"}, ::Colon, args...; kwargs...) where {N} + return FITS(f.filename, "r") do fits + exts_resolved = 1:length(fits) + map(exts_resolved) do ext + _loadhdu(fits[ext], args...; kwargs...) + end + end +end + +_loadhdu(hdu::FITSIO.TableHDU) = Tables.columntable(hdu) +function _loadhdu(hdu::FITSIO.ImageHDU, args...; kwargs...) + if size(hdu) != () + return AstroImage(hdu, args...; kwargs...) + else + # Sometimes files have an empty data HDU that shows up as an image HDU but has headers. + # Fallback to creating an empty AstroImage with those headers. + emptydata = fill(missing) + return AstroImage(emptydata, (), (), read_header(hdu), [emptywcs(emptydata)], Ref(false), ()) + end +end +function indexer(fits::FITS) + ext = 0 + for (i, hdu) in enumerate(fits) + if hdu isa ImageHDU && length(size(hdu)) >= 1 # check if Image is atleast 1D + ext = i + break + end + end + if ext > 1 + @info "Image was loaded from HDU $ext" + elseif ext == 0 + error("There are no ImageHDU extensions in '$(fits.filename)'") + end + return ext +end +indexer(fits::NTuple{N,FITS}) where {N} = ntuple(i -> indexer(fits[i]), N) + + +# Fallback for saving arbitrary arrays +function fileio_save(f::File{format"FITS"}, args...) + return writefits(f.filename, args...) +end +""" + writefits("abc.fits", img1, img2, table1,...) + +Write arguments to a FITS file. + +See also `Fileio.save` +""" +function writefits(fname, args...) + FITS(fname, "w") do fits + for arg in args + writearg(fits, arg) + end + end +end +writearg(fits, img::AstroImage) = write(fits, parent(img), header=header(img)) +# Fallback for writing plain arrays +writearg(fits, arr::AbstractArray) = write(fits, arr) +# For table compatible data. +# This allows users to round trip: dat = load("abc.fits", :); write("abc", dat) +# when it contains FITS tables. +function writearg(fits, table) + if !Tables.istable(table) + error("Cannot save argument to FITS file. Value is not an AbstractArray or table.") + end + # FITSIO has fairly restrictive input types for writing tables (assertions for documentation only) + colname_strings = string.(collect(Tables.columnnames(table)))::Vector{String} + columns = collect(Tables.columns(table))::Vector + write( + fits, + colname_strings, + columns; + hdutype=TableHDU + # TODO: In future, we want to be able to access and round-trip coments + # on table HDUs + # header=nothing + ) +end diff --git a/src/plot-recipes.jl b/src/plot-recipes.jl index a1132e69..bc5bd477 100644 --- a/src/plot-recipes.jl +++ b/src/plot-recipes.jl @@ -1,76 +1,1128 @@ -using RecipesBase - -@recipe function f(img::AstroImage, header_number::Int) - seriestype := :heatmap - aspect_ratio := :equal - # Right now we only support single frame images, - # gray scale is a good choice. - color := :grays - img.data[header_number] +@userplot ImPlot +@recipe function f(h::ImPlot) + if length(h.args) != 1 || !(typeof(h.args[1]) <: AbstractArray) + error("Image plots require an arugment that is a subtype of AbstractArray. Got: $(typeof(h.args))") + end + data = only(h.args) + if !(typeof(data) <: AstroImage) + data = AstroImage(only(h.args)) + end + T = eltype(data) + if ndims(data) != 2 + error("Image passed to `implot` must be two-dimensional. Got ndims(img)=$(ndims(data))") + end + + wcsn = get(plotattributes, :wcsn, 1) + # Show WCS coordinates if wcsticks is true or unspecified, and has at least one WCS axis present. + showwcsticks = get(plotattributes, :wcsticks, true) && !all(==(""), wcs(data, wcsn).ctype) + showwcstitle = get(plotattributes, :wcstitle, true) && length(refdims(data)) > 0 && !all(==(""), wcs(data, wcsn).ctype) + + + + minx = first(parent(dims(data,1))) + maxx = last(parent(dims(data,1))) + miny = first(parent(dims(data,2))) + maxy = last(parent(dims(data,2))) + extent = (minx-0.5, maxx+0.5, miny-0.5, maxy+0.5) + if haskey(plotattributes, :xlims) + extent = (plotattributes[:xlims]..., extent[3:4]...) + end + if haskey(plotattributes, :ylims) + extent = (extent[1:2]..., plotattributes[:ylims]...) + end + if showwcsticks + wcsg = WCSGrid(data, Float64.(extent), wcsn) + gridspec = wcsgridspec(wcsg) + end + + # Use package defaults if not user provided. + clims --> _default_clims[] + stretch --> _default_stretch[] + cmap --> _default_cmap[] + + bias = get(plotattributes, :bias, 0.5) + contrast = get(plotattributes, :contrast, 1) + platescale = get(plotattributes, :platescale, 1) + + grid := false + # In most cases, a grid framestyle is a nicer looking default for images + # but the user can override. + framestyle --> :box + + + if T <: Colorant + imgv = data + else + clims = plotattributes[:clims] + stretch = plotattributes[:stretch] + cmap = plotattributes[:cmap] + if T <: Complex + img = abs.(data) + img["UNIT"] = "magnitude" + else + img = data + end + imgv = imview(img; clims, stretch, cmap, contrast, bias) + end + + # Reduce large images using the same heuristic as Images.jl + maxpixels = get(plotattributes, :maxpixels, 10^6) + _length1(A::AbstractArray) = length(eachindex(A)) + _length1(A) = length(A) + while _length1(imgv) > maxpixels + imgv = restrict(imgv) + end + + # We have to do a lot of flipping to keep the orientation corect + yflip := false + xflip := false + + + # Disable equal aspect ratios if the scales are totally different + displayed_data_ratio = (extent[2]-extent[1])/(extent[4]-extent[3]) + if displayed_data_ratio >= 7 + aspect_ratio --> :none + end + + + # we have a wcs flag (from the image by default) so that users can skip over + # plotting in physical coordinates. This is especially important + # if the WCS headers are mallformed in some way. + showgrid = get(plotattributes, :xgrid, true) && get(plotattributes, :ygrid, true) + # Display a title giving our position along unplotted dimensions + if length(refdims(imgv)) > 0 + if showwcstitle + refdimslabel = join(map(refdims(imgv)) do d + # match dimension with the wcs axis number + i = wcsax(imgv, d) + ct = wcs(imgv, wcsn).ctype[i] + label = ctype_label(ct, wcs(imgv, wcsn).radesys) + if label == "NONE" + label = name(d) + end + value = pix_to_world(imgv, [1,1]; wcsn, all=true, parent=true)[i] + unit = wcs(imgv, wcsn).cunit[i] + if ct == "STOKES" + return _stokes_name(_stokes_symbol(value)) + else + return @sprintf("%s = %.5g %s", label, value, unit) + end + end, ", ") + else + refdimslabel = join(map(d->"$(name(d))= $(d[1])", refdims(imgv)), ", ") + end + title --> refdimslabel + end + + # To ensure the physical axis tick labels are correct the axes must be + # tight to the image + xl = (first(dims(imgv,1))-0.5)*platescale, (last(dims(imgv,1))+0.5)*platescale + yl = (first(dims(imgv,2))-0.5)*platescale, (last(dims(imgv,2))+0.5)*platescale + ylims --> yl + xlims --> xl + + subplot_i = 0 + # Actual image series (RGB pixels by this point) + @series begin + subplot_i += 1 + subplot := subplot_i + colorbar := false + aspect_ratio --> 1 + + # Note: if the axes are on unusual sides (e.g. y-axis at right, x-axis at top) + # then these coordinates are not correct. They are only correct exactly + # along the axis. + # In astropy, the ticks are actually tilted to reflect this, though in general + # the transformation from pixel to coordinates can be non-linear and curved. + + if showwcsticks + xticks --> (gridspec.tickpos1x, wcslabels(wcs(imgv, wcsn), wcsax(imgv, dims(imgv,1)), gridspec.tickpos1w)) + xguide --> ctype_label(wcs(imgv, wcsn).ctype[wcsax(imgv, dims(imgv,1))], wcs(imgv, wcsn).radesys) + + yticks --> (gridspec.tickpos2x, wcslabels(wcs(imgv, wcsn), wcsax(imgv, dims(imgv,2)), gridspec.tickpos2w)) + yguide --> ctype_label(wcs(imgv, wcsn).ctype[wcsax(imgv, dims(imgv,2))], wcs(imgv, wcsn).radesys) + end + + + ax1 = collect(parent(dims(imgv,1))) .* platescale + ax2 = collect(parent(dims(imgv,2))) .* platescale + # Views of images are not currently supported by plotly() so we have to collect them. + # ax1, ax2, view(parent(imgv), reverse(axes(imgv,1)),:) + ax1, ax2, parent(imgv)[reverse(axes(imgv,1)),:] + end + + # If wcs=true (default) and grid=true (not default), overplot a WCS + # grid. + if showgrid && showwcsticks + + # Plot the WCSGrid as a second series (actually just lines) + @series begin + subplot := 1 + # Use a default grid color that shows up across more + # color maps + if !haskey(plotattributes, :xforeground_color_grid) && !haskey(plotattributes, :yforeground_color_grid) + gridcolor --> :lightgray + end + + wcsg, gridspec + end + end + + + # Disable the colorbar. + # Plots.jl does not give us sufficient control to make sure the range and ticks + # are correct after applying a non-linear stretch. + # We attempt to make our own colorbar using a second plot. + showcolorbar = !(T <: Colorant) && get(plotattributes, :colorbar, true) != :none + if T <: Complex + layout := @layout [ + imgmag{0.5h} + imgangle{0.5h} + ] + end + if showcolorbar + if T <: Complex + layout := @layout [ + imgmag{0.95w, 0.5h} colorbar{0.5h} + imgangle{0.95w, 0.5h} colorbarangle{0.5h} + ] + else + layout := @layout [ + img{0.95w} colorbar + ] + end + colorbar_title = get(plotattributes, :colorbar_title, "") + if !haskey(plotattributes, :colorbar_title) + if haskey(header(img), "UNIT") + colorbar_title = string(img[:UNIT]) + elseif haskey(header(img), "BUNIT") + colorbar_title = string(img[:BUNIT]) + end + end + + subplot_i += 1 + @series begin + subplot := subplot_i + aspect_ratio := :none + colorbar := false + cbimg, cbticks = imview_colorbar(img; clims, stretch, cmap, contrast, bias) + xticks := [] + ymirror := true + yticks := cbticks + yguide := colorbar_title + xguide := "" + xlims := Tuple(axes(cbimg, 2)) + ylims := Tuple(axes(cbimg, 2)) + title := "" + # Views of images are not currently supported by plotly so we have to collect them + # view(cbimg, reverse(axes(cbimg,1)),:) + cbimg[reverse(axes(cbimg,1)),:] + end + end + + + # TODO: refactor to reduce duplication + if T <: Complex + img = angle.(data) + img["UNIT"] = "angle (rad)" + imgv = imview(img, clims=(-1pi, 1pi),stretch=identity, cmap=:cyclic_mygbm_30_95_c78_n256_s25) + @series begin + subplot_i += 1 + subplot := subplot_i + colorbar := false + title := "" + aspect_ratio --> 1 + + + # Note: if the axes are on unusual sides (e.g. y-axis at right, x-axis at top) + # then these coordinates are not correct. They are only correct exactly + # along the axis. + # In astropy, the ticks are actually tilted to reflect this, though in general + # the transformation from pixel to coordinates can be non-linear and curved. + + if showwcsticks + xticks --> (gridspec.tickpos1x, wcslabels(wcs(imgv, wcsn), wcsax(imgv, dims(imgv,1)), gridspec.tickpos1w)) + xguide --> ctype_label(wcs(imgv, wcsn).ctype[wcsax(imgv, dims(imgv,1))], wcs(imgv, wcsn).radesys) + + yticks --> (gridspec.tickpos2x, wcslabels(wcs(imgv, wcsn), wcsax(imgv, dims(imgv,2)), gridspec.tickpos2w)) + yguide --> ctype_label(wcs(imgv, wcsn).ctype[wcsax(imgv, dims(imgv,2))], wcs(imgv, wcsn).radesys) + end + + ax1 = collect(parent(dims(imgv,1))) .* platescale + ax2 = collect(parent(dims(imgv,2))) .* platescale + # Views of images are not currently supported by plotly() so we have to collect them. + # ax1, ax2, view(parent(imgv), reverse(axes(imgv,1)),:) + ax1, ax2, parent(imgv)[reverse(axes(imgv,1)),:] + end + + if showcolorbar + colorbar_title = get(plotattributes, :colorbar_title, "") + if !haskey(plotattributes, :colorbar_title) && haskey(header(img), "UNIT") + colorbar_title = string(img[:UNIT]) + end + + + @series begin + subplot_i += 1 + subplot := subplot_i + aspect_ratio := :none + colorbar := false + cbimg, _ = imview_colorbar(img; stretch=identity, clims=(-pi, pi), cmap=:cyclic_mygbm_30_95_c78_n256_s25) + xticks := [] + ymirror := true + ax = axes(cbimg,1) + yticks := ([first(ax), mean(ax), last(ax)], ["-π", "0", "π"]) + yguide := colorbar_title + xguide := "" + xlims := Tuple(axes(cbimg, 2)) + ylims := Tuple(axes(cbimg, 2)) + title := "" + view(cbimg, reverse(axes(cbimg,1)),:) + end + end + + end + + + return end -@recipe function f(img::AstroImage) - seriestype := :heatmap - aspect_ratio := :equal - # Right now we only support single frame images, - # gray scale is a good choice. - color := :grays - img.data[1] + +""" + implot( + img::AbstractArray; + clims=Percent(99.5), + stretch=identity, + cmap=:magma, + bias=0.5, + contrast=1, + wcsticks=true, + grid=true, + platescale=1 + ) + +Create a read only view of an array or AstroImageMat mapping its data values +to an array of Colors. Equivalent to: + + implot( + imview( + img::AbstractArray; + clims=Percent(99.5), + stretch=identity, + cmap=:magma, + bias=0.5, + contrast=1, + ), + wcsn=1, + wcsticks=true, + wcstitle=true, + grid=true, + platescale=1 + ) + +### Image Rendering +See `imview` for how data is mapped to RGBA pixel values. + +### WCS & Image Coordinates +If provided with an AstroImage that has WCS headers set, the tick marks and plot grid +are calculated using WCS.jl. By default, use the first WCS coordinate system. +The underlying pixel coordinates are those returned by `dims(img)` multiplied by `platescale`. +This allows you to overplot lines, regions, etc. using pixel coordinates. +If you wish to compute the pixel coordinate of a point in world coordinates, see `world_to_pix`. + +* `wcsn` (default `1`) select which WCS transform in the headers to use for ticks & grid +* `wcsticks` (default `true` if WCS headers present) display ticks and labels, and title using world coordinates +* `wcstitle` (default `true` if WCS headers present and `length(refdims(img))>0`). When slicing a cube, display the location along unseen axes in world coordinates instead of pixel coordinates. +* `grid` (default `true`) show a grid over the plot. Uses WCS coordinates if `wcsticks` is true, otherwise pixel coordinates multiplied by `platescale`. +* `platescale` (default `1`). Scales the underlying pixel coordinates to ease overplotting, etc. If `wcsticks` is false, the displayed pixel coordinates are also scaled. + + +### Defaults +The default values of `clims`, `stretch`, and `cmap` are `extrema`, `identity`, and `nothing` +respectively. +You may alter these defaults using `AstroImages.set_clims!`, `AstroImages.set_stretch!`, and +`AstroImages.set_cmap!`. +""" +implot + +struct WCSGrid + img::AstroImage + extent::NTuple{4,Float64} + wcsn::Int end -@recipe function f(img::AstroImage, wcs::WCSTransform) - seriestype := :heatmap - aspect_ratio := :equal - color := :grays - xformatter := x -> pix2world_xformatter(x, wcs) - yformatter := y -> pix2world_yformatter(y, wcs) - xlabel := labler_x(wcs) - ylabel := labler_y(wcs) - img.data + +""" + wcsticks(img, axnum) + +Generate nice tick labels for an AstroImageMat along axis `axnum` +Returns a vector of pixel positions and a vector of strings. + +Example: +plot(img, xticks=wcsticks(WCSGrid(img), 1), yticks=wcsticks(WCSGrid(img), 2)) +""" +function wcsticks(wcsg::WCSGrid, axnum, gs = wcsgridspec(wcsg)) + tickposx = axnum == 1 ? gs.tickpos1x : gs.tickpos2x + tickposw = axnum == 1 ? gs.tickpos1w : gs.tickpos2w + return tickposx, wcslabels( + wcs(wcsg.img, wcsg.wcsn), + axnum, + tickposw + ) end -function pix2world_xformatter(x, wcs::WCSTransform) - res = round(pix_to_world(wcs, [float(x), float(x)])[1], digits = 2) - if wcs.cunit[1] == "deg" # TODO: add symbols for more units - return string(res)*"°" +# Function to generate nice string coordinate labels given a WCSTransform, axis number, +# and a vector of tick positions in world coordinates. +# This is used for labelling ticks and for annotating grid lines. +function wcslabels(w::WCSTransform, axnum, tickposw) + + if length(tickposw) == 0 + return String[] + end + + # Select a unit converter (e.g. 12.12 -> (a,b,c,d)) and list of units + if w.cunit[axnum] == "deg" + if startswith(uppercase(w.ctype[axnum]), "RA") + converter = deg2hms + units = hms_units + else + converter = deg2dmsmμ + units = dmsmμ_units + end else - return res[1] + converter = x->(x,) + units = ("",) + end + + # Format inital ticklabel + ticklabels = fill("", length(tickposw)) + # We only include the part of the label that has changed since the last time. + # Split up coordinates into e.g. sexagesimal + parts = map(tickposw) do w + vals = converter(w) + return vals + end + + # Start with something impossible of the same size: + last_coord = Inf .* converter(first(tickposw)) + zero_coords_i = maximum(map(parts) do vals + changing_coord_i = findfirst(vals .!= last_coord) + if isnothing(changing_coord_i) + changing_coord_i = 1 + end + last_coord = vals + return changing_coord_i + end) + + + # Loop through using only the relevant part of the label + # Start with something impossible of the same size: + last_coord = Inf .* converter(first(tickposw)) + for (i,vals) in enumerate(parts) + changing_coord_i = findfirst(vals .!= last_coord) + if isnothing(changing_coord_i) + changing_coord_i = 1 + end + # Don't display just e.g. 00" when we should display 50'00" + if changing_coord_i > 1 && vals[changing_coord_i] == 0 + changing_coord_i = changing_coord_i -1 + end + val_unit_zip = zip(vals[changing_coord_i:zero_coords_i],units[changing_coord_i:zero_coords_i]) + ticklabels[i] = mapreduce(*, enumerate(val_unit_zip)) do (coord_i,(val,unit)) + # If the last coordinate we print if the last coordinate we have available, + # display it with decimal places + if coord_i + changing_coord_i - 1== length(vals) + str = @sprintf("%.2f", val) + else + str = @sprintf("%02d", val) + end + if length(str) > 0 + return str * unit + else + return str + end + end + last_coord = vals end + + return ticklabels +end + +# Extended form of deg2dms that further returns mas, microas. +function deg2dmsmμ(deg) + d,m,s = deg2dms(deg) + s_f = floor(s) + mas = (s - s_f)*1e3 + mas_f = floor(mas) + μas = (mas - mas_f)*1e3 + return (d,m,s_f,mas_f,μas) end +const dmsmμ_units = [ + "°", + "'", + "\"", + "mas", + "μas", +] +const hms_units = [ + "ʰ", + "ᵐ", + "ˢ", +] -function pix2world_yformatter(x, wcs::WCSTransform) - res = round(pix_to_world(wcs, [float(x), float(x)])[2], digits = 2) - if wcs.cunit[2] == "deg" # TODO: add symbols for more units - return string(res)*"°" +function ctype_label(ctype,radesys) + if length(ctype) == 0 + return radesys + elseif startswith(ctype, "RA") + return "Right Ascension ($(radesys))" + elseif startswith(ctype, "GLON") + return "Galactic Longitude" + elseif startswith(ctype, "TLON") + return "ITRS" + elseif startswith(ctype, "DEC") + return "Declination ($(radesys))" + elseif startswith(ctype, "GLAT") + return "Galactic Latitude" + # elseif startswith(ctype, "TLAT") + elseif ctype == "STOKES" + return "Polarization" else - return res[1] + return ctype end end -function labler_x(wcs::WCSTransform) - if length(wcs.ctype[1]) == 0 - return wcs.radesys - elseif wcs.ctype[1][1:2] == "RA" - return "Right Ascension (ICRS)" - elseif wcs.ctype[1][1:4] == "GLON" - return "Galactic Coordinate" - elseif wcs.ctype[1][1:4] == "TLON" - return "ITRS" + + +""" + WCSGrid(img::AstroImageMat, ax=(1,2), coords=(first(axes(img,ax[1])),first(axes(img,ax[2])))) + +Given an AstroImageMat, return information necessary to plot WCS gridlines in physical +coordinates against the image's pixel coordinates. +This function has to work on both plotted axes at once to handle rotation and general +curvature of the WCS grid projected on the image coordinates. + +""" +function WCSGrid(img::AstroImageMat, wcsn=1) + minx = first(dims(img,2)) + maxx = last(dims(img,2)) + miny = first(dims(img,1)) + maxy = last(dims(img,1)) + extent = (minx-0.5, maxx+0.5, miny-0.5, maxy+0.5) + @show extent + return WCSGrid(img, extent, wcsn) +end + + + +# Recipe for a WCSGrid with lines, optional ticks (on by default), +# and optional grid labels (off by defaut). +# The AstroImageMat plotrecipe uses this recipe for grid lines if `grid=true`. +@recipe function f(wcsg::WCSGrid, gridspec=wcsgridspec(wcsg)) + label --> "" + xs, ys = wcsgridlines(gridspec) + + if haskey(plotattributes, :foreground_color_grid) + color --> plotattributes[:foreground_color_grid] + elseif haskey(plotattributes, :foreground_color) + color --> plotattributes[:foreground_color] + else + color --> :black + end + if haskey(plotattributes, :foreground_color_text) + textcolor = plotattributes[:foreground_color_text] else - return wcs.ctype[1] + textcolor = plotattributes[:color] + end + annotate = haskey(plotattributes, :gridlabels) && plotattributes[:gridlabels] + + xguide --> ctype_label(wcs(wcsg.img, wcsg.wcsn).ctype[wcsax(wcsg.img, dims(wcsg.img,1))], wcs(wcsg.img, wcsg.wcsn).radesys) + yguide --> ctype_label(wcs(wcsg.img, wcsg.wcsn).ctype[wcsax(wcsg.img, dims(wcsg.img,2))], wcs(wcsg.img, wcsg.wcsn).radesys) + + xlims --> wcsg.extent[1], wcsg.extent[2] + ylims --> wcsg.extent[3], wcsg.extent[4] + + grid := false + tickdirection := :none + + xticks --> wcsticks(wcsg, 1, gridspec) + yticks --> wcsticks(wcsg, 2, gridspec) + + @series xs, ys + + # We can optionally annotate the grid with their coordinates. + # These come after the grid lines so they appear overtop. + if annotate + @series begin + # TODO: why is this reverse necessary? + rotations = reverse(rad2deg.(gridspec.annotations1θ)) + ticklabels = wcslabels(wcs(wcsg.img), 1, gridspec.annotations1w) + seriestype := :line + linewidth := 0 + # TODO: we need to use requires to load in Plots for the necessary text control. Future versions of RecipesBase might fix this. + series_annotations := [ + Main.Plots.text(" $l", :right, :bottom, textcolor, 8, rotation=(-95 <= r <= 95) ? r : r+180) + for (l, r) in zip(ticklabels, rotations) + ] + gridspec.annotations1x, gridspec.annotations1y + end + @series begin + rotations = rad2deg.(gridspec.annotations2θ) + ticklabels = wcslabels(wcs(wcsg.img), 2, gridspec.annotations2w) + seriestype := :line + linewidth := 0 + series_annotations := [ + Main.Plots.text(" $l", :right, :bottom, textcolor, 8, rotation=(-95 <= r <= 95) ? r : r+180) + for (l, r) in zip(ticklabels, rotations) + ] + gridspec.annotations2x, gridspec.annotations2y + end + end + + return end -function labler_y(wcs::WCSTransform) - if length(wcs.ctype[2]) == 0 - return wcs.radesys - elseif wcs.ctype[2][1:3] == "DEC" - return "Declination (ICRS)" - elseif wcs.ctype[2][1:4] == "GLAT" - return "Galactic Coordinate" - elseif wcs.ctype[2][1:4] == "TLAT" - return "ITRS" +# Helper: true if all elements in vector are equal to each other. +allequal(itr) = all(==(first(itr)), itr) + +# This function is responsible for actually laying out grid lines for a WCSGrid, +# ensuring they don't exceed the plot bounds, finding where they intersect the axes, +# and picking tick locations at the appropriate intersections with the left and +# bottom axes. +function wcsgridspec(wsg::WCSGrid) + # Most of the complexity of this function is making sure everything + # generalizes to N different, possiby skewed axes, where a change in + # the opposite coordinate or even an unplotted coordinate affects + # the grid. + + # x and y denote pixel coordinates (along `ax`), u and v are world coordinates roughly along same. + minx, maxx, miny, maxy = wsg.extent + + # Find the extent of this slice in world coordinates + posxy = [ + minx minx maxx maxx + miny maxy miny maxy + ] + posuv = pix_to_world(wsg.img, posxy; wsg.wcsn, parent=true) + (minu, maxu), (minv, maxv) = extrema(posuv, dims=2) + + # In general, grid can be curved when plotted back against the image, + # so we will need to sample multiple points along the grid. + # TODO: find a good heuristic for this based on the curvature. + N_points = 50 + urange = range(minu, maxu, length=N_points) + vrange = range(minv, maxv, length=N_points) + + # Find nice grid spacings using PlotUtils.optimize_ticks + # These heuristics can probably be improved + # TODO: this does not handle coordinates that wrap around + Q=[(1.0,1.0), (3.0, 0.8), (2.0, 0.7), (5.0, 0.5)] + k_min = 3 + k_ideal = 5 + k_max = 10 + + tickpos2x = Float64[] + tickpos2w = Float64[] + gridlinesxy2 = NTuple{2,Vector{Float64}}[] + # Not all grid lines will intersect the x & y axes nicely. + # If we don't get enough valid tick marks (at least 2) loop again + # requesting more locations up to three times. + local tickposv + j = 5 + while length(tickpos2x) < 2 && j > 0 + k_min += 2 + k_ideal += 2 + k_max += 2 + j -= 1 + + tickposv = optimize_ticks(6minv, 6maxv; Q, k_min, k_ideal, k_max)[1]./6 + + empty!(tickpos2x) + empty!(tickpos2w) + empty!(gridlinesxy2) + for tickv in tickposv + # Make sure we handle unplotted slices correctly. + griduv = repeat(posuv[:,1], 1, N_points) + griduv[1,:] .= urange + griduv[2,:] .= tickv + posxy = world_to_pix(wsg.img, griduv; wsg.wcsn, parent=true) + + # Now that we have the grid in pixel coordinates, + # if we find out where the grid intersects the axes we can put + # the labels in the correct spot + + # We can use these masks to determine where, and in what direction + # the gridlines leave the plot extent + in_horz_ax = minx .<= posxy[1,:] .<= maxx + in_vert_ax = miny .<= posxy[2,:] .<= maxy + in_axes = in_horz_ax .& in_vert_ax + if count(in_axes) < 2 + continue + elseif all(in_axes) + point_entered = [ + posxy[1,begin] + posxy[2,begin] + ] + point_exitted = [ + posxy[1,end] + posxy[2,end] + ] + elseif allequal(posxy[1,findfirst(in_axes):findlast(in_axes)]) + point_entered = [ + posxy[1,max(begin,findfirst(in_axes)-1)] + # posxy[2,max(begin,findfirst(in_axes)-1)] + miny + ] + point_exitted = [ + posxy[1,min(end,findlast(in_axes)+1)] + # posxy[2,min(end,findlast(in_axes)+1)] + maxy + ] + # Vertical grid lines + elseif allequal(posxy[2,findfirst(in_axes):findlast(in_axes)]) + point_entered = [ + minx #posxy[1,max(begin,findfirst(in_axes)-1)] + posxy[2,max(begin,findfirst(in_axes)-1)] + ] + point_exitted = [ + maxx #posxy[1,min(end,findlast(in_axes)+1)] + posxy[2,min(end,findlast(in_axes)+1)] + ] + else + # Use the masks to pick an x,y point inside the axes and an + # x,y point outside the axes. + i = findfirst(in_axes) + x1 = posxy[1,i] + y1 = posxy[2,i] + x2 = posxy[1,i+1] + y2 = posxy[2,i+1] + if x2-x1 ≈ 0 + @warn "undef slope" + end + + # Fit a line where we cross the axis + m1 = (y2-y1)/(x2-x1) + b1 = y1-m1*x1 + # If the line enters via the vertical axes... + if findfirst(in_vert_ax) <= findfirst(in_horz_ax) + # Then we simply evaluate it at that axis + x = abs(x1-maxx) < abs(x1-minx) ? maxx : minx + x = clamp(x,minx,maxx) + y = m1*x+b1 + else + # We must find where it enters the plot from + # bottom or top + x = abs(y1-maxy) < abs(y1-miny) ? (maxy-b1)/m1 : (miny-b1)/m1 + x = clamp(x,minx,maxx) + y = m1*x+b1 + end + + # From here, do a linear fit to find the intersection with the axis. + point_entered = [ + x + y + ] + + + # Use the masks to pick an x,y point inside the axes and an + # x,y point outside the axes. + i = findlast(in_axes) + x1 = posxy[1,i-1] + y1 = posxy[2,i-1] + x2 = posxy[1,i] + y2 = posxy[2,i] + if x2-x1 ≈ 0 + @warn "undef slope" + end + + # Fit a line where we cross the axis + m2 = (y2-y1)/(x2-x1) + b2 = y2-m2*x2 + if findlast(in_vert_ax) > findlast(in_horz_ax) + # Then we simply evaluate it at that axis + x = abs(x1-maxx) < abs(x1-minx) ? maxx : minx + x = clamp(x,minx,maxx) + y = m2*x+b2 + else + # We must find where it enters the plot from + # bottom or top + x = abs(y1-maxy) < abs(y1-miny) ? (maxy-b2)/m2 : (miny-b2)/m2 + x = clamp(x,minx,maxx) + y = m2*x+b2 + end + + # From here, do a linear fit to find the intersection with the axis. + point_exitted = [ + x + y + ] + end + + + if point_entered[1] == minx + push!(tickpos2x, point_entered[2]) + push!(tickpos2w, tickv) + end + if point_exitted[1] == minx + push!(tickpos2x, point_exitted[2]) + push!(tickpos2w, tickv) + end + + + posxy_neat = [point_entered posxy[[1,2],in_axes] point_exitted] + # posxy_neat = posxy + # TODO: do unplotted other axes also need a fit? + + gridlinexy = ( + posxy_neat[1,:], + posxy_neat[2,:] + ) + push!(gridlinesxy2, gridlinexy) + end + end + + # Then do the opposite coordinate + k_min = 3 + k_ideal = 5 + k_max = 10 + tickpos1x = Float64[] + tickpos1w = Float64[] + gridlinesxy1 = NTuple{2,Vector{Float64}}[] + # Not all grid lines will intersect the x & y axes nicely. + # If we don't get enough valid tick marks (at least 2) loop again + # requesting more locations up to three times. + local tickposu + j = 5 + while length(tickpos1x) < 2 && j > 0 + k_min += 2 + k_ideal += 2 + k_max += 2 + j -= 1 + + tickposu = optimize_ticks(6minu, 6maxu; Q, k_min, k_ideal, k_max)[1]./6 + + empty!(tickpos1x) + empty!(tickpos1w) + empty!(gridlinesxy1) + for ticku in tickposu + # Make sure we handle unplotted slices correctly. + griduv = repeat(posuv[:,1], 1, N_points) + griduv[1,:] .= ticku + griduv[2,:] .= vrange + posxy = world_to_pix(wsg.img, griduv; wsg.wcsn, parent=true) + + # Now that we have the grid in pixel coordinates, + # if we find out where the grid intersects the axes we can put + # the labels in the correct spot + + # We can use these masks to determine where, and in what direction + # the gridlines leave the plot extent + in_horz_ax = minx .<= posxy[1,:] .<= maxx + in_vert_ax = miny .<= posxy[2,:] .<= maxy + in_axes = in_horz_ax .& in_vert_ax + + + if count(in_axes) < 2 + continue + elseif all(in_axes) + point_entered = [ + posxy[1,begin] + posxy[2,begin] + ] + point_exitted = [ + posxy[1,end] + posxy[2,end] + ] + # Horizontal grid lines + elseif allequal(posxy[1,findfirst(in_axes):findlast(in_axes)]) + point_entered = [ + posxy[1,findfirst(in_axes)] + miny + ] + point_exitted = [ + posxy[1,findlast(in_axes)] + maxy + ] + # push!(tickpos1x, posxy[1,findfirst(in_axes)]) + # push!(tickpos1w, ticku) + # Vertical grid lines + elseif allequal(posxy[2,findfirst(in_axes):findlast(in_axes)]) + point_entered = [ + minx + posxy[2,findfirst(in_axes)] + ] + point_exitted = [ + maxx + posxy[2,findfirst(in_axes)] + ] + else + # Use the masks to pick an x,y point inside the axes and an + # x,y point outside the axes. + i = findfirst(in_axes) + x1 = posxy[1,i] + y1 = posxy[2,i] + x2 = posxy[1,i+1] + y2 = posxy[2,i+1] + if x2-x1 ≈ 0 + @warn "undef slope" + end + + # Fit a line where we cross the axis + m1 = (y2-y1)/(x2-x1) + b1 = y1-m1*x1 + # If the line enters via the vertical axes... + if findfirst(in_vert_ax) < findfirst(in_horz_ax) + # Then we simply evaluate it at that axis + x = abs(x1-maxx) < abs(x1-minx) ? maxx : minx + x = clamp(x,minx,maxx) + y = m1*x+b1 + else + # We must find where it enters the plot from + # bottom or top + x = abs(y1-maxy) < abs(y1-miny) ? (maxy-b1)/m1 : (miny-b1)/m1 + x = clamp(x,minx,maxx) + y = m1*x+b1 + end + + # From here, do a linear fit to find the intersection with the axis. + point_entered = [ + x + y + ] + + # Use the masks to pick an x,y point inside the axes and an + # x,y point outside the axes. + i = findlast(in_axes) + x1 = posxy[1,i-1] + y1 = posxy[2,i-1] + x2 = posxy[1,i] + y2 = posxy[2,i] + if x2-x1 ≈ 0 + @warn "undef slope" + end + + # Fit a line where we cross the axis + m2 = (y2-y1)/(x2-x1) + b2 = y2-m2*x2 + if findlast(in_vert_ax) > findlast(in_horz_ax) + # Then we simply evaluate it at that axis + x = abs(x1-maxx) < abs(x1-minx) ? maxx : minx + x = clamp(x,minx,maxx) + y = m2*x+b2 + else + # We must find where it enters the plot from + # bottom or top + x = abs(y1-maxy) < abs(y1-miny) ? (maxy-b2)/m2 : (miny-b2)/m2 + x = clamp(x,minx,maxx) + y = m2*x+b2 + end + + # From here, do a linear fit to find the intersection with the axis. + point_exitted = [ + x + y + ] + end + + posxy_neat = [point_entered posxy[[1,2],in_axes] point_exitted] + # TODO: do unplotted other axes also need a fit? + + if point_entered[2] == miny + push!(tickpos1x, point_entered[1]) + push!(tickpos1w, ticku) + end + if point_exitted[2] == miny + push!(tickpos1x, point_exitted[1]) + push!(tickpos1w, ticku) + end + + gridlinexy = ( + posxy_neat[1,:], + posxy_neat[2,:] + ) + push!(gridlinesxy1, gridlinexy) + end + end + + # Grid annotations are simpler: + annotations1w = Float64[] + annotations1x = Float64[] + annotations1y = Float64[] + annotations1θ = Float64[] + for ticku in tickposu + # Make sure we handle unplotted slices correctly. + griduv = posuv[:,1] + griduv[1] = ticku + griduv[2] = mean(vrange) + posxy = world_to_pix(wsg.img, griduv; wsg.wcsn, parent=true) + if !(minx < posxy[1] < maxx) || + !(miny < posxy[2] < maxy) + continue + end + push!(annotations1w, ticku) + push!(annotations1x, posxy[1]) + push!(annotations1y, posxy[2]) + + # Now find slope (TODO: stepsize) + # griduv[ax[2]] -= 1 + griduv[2] += 0.1step(vrange) + posxy2 = world_to_pix(wsg.img, griduv; wsg.wcsn, parent=true) + θ = atan( + posxy2[2] - posxy[2], + posxy2[1] - posxy[1], + ) + push!(annotations1θ, θ) + end + annotations2w = Float64[] + annotations2x = Float64[] + annotations2y = Float64[] + annotations2θ = Float64[] + for tickv in tickposv + # Make sure we handle unplotted slices correctly. + griduv = posuv[:,1] + griduv[1] = mean(urange) + griduv[2] = tickv + posxy = world_to_pix(wsg.img, griduv; wsg.wcsn, parent=true) + if !(minx < posxy[1] < maxx) || + !(miny < posxy[2] < maxy) + continue + end + push!(annotations2w, tickv) + push!(annotations2x, posxy[1]) + push!(annotations2y, posxy[2]) + + griduv[1] += 0.1step(urange) + posxy2 = world_to_pix(wsg.img, griduv; wsg.wcsn, parent=true) + θ = atan( + posxy2[2] - posxy[2], + posxy2[1] - posxy[1], + ) + push!(annotations2θ, θ) + end + + return (; + gridlinesxy1, + gridlinesxy2, + tickpos1x, + tickpos1w, + tickpos2x, + tickpos2w, + + annotations1w, + annotations1x, + annotations1y, + annotations1θ, + + annotations2w, + annotations2x, + annotations2y, + annotations2θ, + ) +end + +# From a WCSGrid, return just the grid lines as a single pair of x & y coordinates +# suitable for plotting. +function wcsgridlines(wcsg::WCSGrid) + return wcsgridlines(wcsgridspec(wcsg)) +end +function wcsgridlines(gridspec::NamedTuple) + # Unroll grid lines into a single series separated by NaNs + xs1 = mapreduce(vcat, gridspec.gridlinesxy1, init=Float64[]) do gridline + return vcat(gridline[1], NaN) + end + ys1 = mapreduce(vcat, gridspec.gridlinesxy1, init=Float64[]) do gridline + return vcat(gridline[2], NaN) + end + xs2 = mapreduce(vcat, gridspec.gridlinesxy2, init=Float64[]) do gridline + return vcat(gridline[1], NaN) + end + ys2 = mapreduce(vcat, gridspec.gridlinesxy2, init=Float64[]) do gridline + return vcat(gridline[2], NaN) + end + + xs = vcat(xs1, NaN, xs2) + ys = vcat(ys1, NaN, ys2) + return xs, ys +end + + + + +@userplot PolQuiver +@recipe function f(h::PolQuiver) + cube = only(h.args) + bins = get(plotattributes, :bins, 4) + ticklen = get(plotattributes, :ticklen, nothing) + minpol = get(plotattributes, :minpol, 0.1) + + i = cube[Pol=At(:I)] + q = cube[Pol=At(:Q)] + u = cube[Pol=At(:U)] + polinten = @. sqrt(q^2 + u^2) + linpolfrac = polinten ./ i + + binratio=1/bins + xs = imresize([x for x in dims(cube,1), y in dims(cube,2)], ratio=binratio) + ys = imresize([y for x in dims(cube,1), y in dims(cube,2)], ratio=binratio) + qx = imresize(q, ratio=binratio) + qy = imresize(u, ratio=binratio) + qlinpolfrac = imresize(linpolfrac, ratio=binratio) + qpolintenr = imresize(polinten, ratio=binratio) + + + # We want the longest ticks to be around 1 bin long by default. + qmaxlen = quantile(filter(isfinite,qpolintenr), 0.98) + if isnothing(ticklen) + a = bins / qmaxlen else - return wcs.ctype[2] + a = ticklen / qmaxlen + end + # Only show arrows where the data is finite, and more than a couple pixels + # long. + mask = (isfinite.(qpolintenr)) .& (qpolintenr .>= minpol.*qmaxlen) + pointstmp = map(xs[mask],ys[mask],qx[mask],qy[mask]) do x,y,qxi,qyi + return ([x, x+a*qxi, NaN], [y, y+a*qyi, NaN]) + end + xs = reduce(vcat, getindex.(pointstmp, 1)) + ys = reduce(vcat, getindex.(pointstmp, 2)) + + colors = qlinpolfrac[mask] + if !isnothing(colors) + line_z := repeat(colors, inner=3) + end + + label --> "" + color --> :turbo + framestyle --> :box + aspect_ratio --> 1 + linewidth --> 1.5 + colorbar --> true + colorbar_title --> "Linear polarization fraction" + + xl = first(dims(i,2)), last(dims(i,2)) + yl = first(dims(i,1)), last(dims(i,1)) + ylims --> yl + xlims --> xl + + @series begin + xs, ys end end + +""" + polquiver(polqube::AstroImage) + +Given a data cube (of at least 2 spatial dimensions, plus a polarization axis), +plot a vector field of polarization data. +The tick length represents the polarization intensity, sqrt(q^2 + u^2), +and the color represents the linear polarization fraction, sqrt(q^2 + u^2) / i. + +There are several ways you can adjust the appearance of the plot using keyword arguments: +* `bins` (default = 1) By how much should we bin down the polarization data before drawing the ticks? This reduced clutter from higher resolution datasets. Can be fractional. +* `ticklen` (default = bins) How long the 98th percentile arrow should be. By default, 1 bin long. Make this larger to draw longer arrows. +* `color` (default = :turbo) What colorscheme should be used for linear polarization fraction. +* `minpol` (default = 0.2) Hides arrows that are shorter than `minpol` times the 98th percentile arrow to make a cleaner image. Set to 0 to display all data. + +Use `implot` and `polquiver!` to overplot polarization data over an image. +""" +polquiver diff --git a/src/precompile.jl b/src/precompile.jl new file mode 100644 index 00000000..8edfd32b --- /dev/null +++ b/src/precompile.jl @@ -0,0 +1,44 @@ + + + +for T in [Float32, Float64, Int, Int8, UInt8, N0f8] + + a = rand(T, 5, 5) + i = AstroImage(a) + for stretch in [ + logstretch, + powstretch, + sqrtstretch, + squarestretch, + asinhstretch, + sinhstretch, + powerdiststretch + ] + # Easiest way to precompile everything we need is just to call these functions. + # They have no side-effects. + imview(a; stretch) + + # And precompile on an astroimage + imview(i; stretch) + end + TI = typeof(i) + precompile(parent, (TI,)) + precompile(header, (TI,)) + precompile(wcs, (TI,)) + precompile(getindex, (TI, Symbol)) + precompile(getindex, (TI, String)) + precompile(getindex, (TI, Int)) + precompile(getindex, (TI, Int, Int)) + precompile(getindex, (TI, Vector{Int})) + precompile(getindex, (TI, Vector{Bool})) + precompile(getindex, (TI, Matrix{Bool})) + precompile(setindex!, (TI, Matrix{Bool})) + precompile(world_to_pix, (TI, Vector{Float64})) + precompile(pix_to_world, (TI, Vector{Float64})) +end + + +# From trace-compile: +precompile(Tuple{typeof(AstroImages.imview), Array{Float64, 2}}) +precompile(Tuple{AstroImages.var"#imview##kw", NamedTuple{(:clims, :stretch, :cmap, :contrast, :bias), Tuple{AstroImages.Percent, typeof(Base.identity), Symbol, Int64, Float64}}, typeof(AstroImages.imview), AstroImages.AstroImage{Float64, 2, Tuple{DimensionalData.Dimensions.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}}, DimensionalData.Dimensions.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{}, Array{Float64, 2}, Tuple{DimensionalData.Dimensions.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}}, DimensionalData.Dimensions.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}}}}}) +precompile(Tuple{AstroImages.var"#imview_colorbar##kw", NamedTuple{(:clims, :stretch, :cmap, :contrast, :bias), Tuple{AstroImages.Percent, typeof(Base.identity), Symbol, Int64, Float64}}, typeof(AstroImages.imview_colorbar), AstroImages.AstroImage{Float64, 2, Tuple{DimensionalData.Dimensions.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}}, DimensionalData.Dimensions.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{}, Array{Float64, 2}, Tuple{DimensionalData.Dimensions.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}}, DimensionalData.Dimensions.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}}}}}) diff --git a/src/showmime.jl b/src/showmime.jl index c6dd36a8..f77aecee 100644 --- a/src/showmime.jl +++ b/src/showmime.jl @@ -1,20 +1,19 @@ -_brightness_contrast(color, matrix::AbstractMatrix{T}, brightness, contrast) where {T} = - @. color(matrix / 255 * T(contrast) + T(brightness) / 255) -""" - brightness_contrast(image::AstroImage; brightness_range = 0:255, contrast_range = 1:1000, header_number = 1) +# This is used in VSCode and others -Visualize the fits image by changing the brightness and contrast of image. +# If the user displays a AstroImageMat of colors (e.g. one created with imview) +# fal through and display the data as an image +Base.show(io::IO, mime::MIME"image/png", img::AstroImageMat{T}; kwargs...) where {T<:Colorant} = + show(io, mime, parent(img), kwargs...) -Users can also provide their own range as keyword arguments. -""" -function brightness_contrast(img::AstroImage{T,C,N}; brightness_range = 0:255, - contrast_range = 1:1000, header_number = 1) where {T,C,N} - @manipulate for brightness in brightness_range, contrast in contrast_range - _brightness_contrast(C, img.data[header_number], brightness, contrast) - end -end +# Otherwise, call imview with the default settings. +Base.show(io::IO, mime::MIME"image/png", img::AstroImageMat{T}; kwargs...) where {T<:Union{Number,Missing}} = + show(io, mime, imview(img), kwargs...) -# This is used in Jupyter notebooks -Base.show(io::IO, mime::MIME"text/html", img::AstroImage; kwargs...) = - show(io, mime, brightness_contrast(img), kwargs...) + +# Deprecated +# Lazily reinterpret the AstroImageMat as a Matrix{Color}, upon request. +# By itself, Images.colorview works fine on AstroImages. But +# AstroImages are not normalized to be between [0,1]. So we override +# colorview to first normalize the data using scaleminmax +@deprecate render(img::AstroImageMat) imview(img, clims=extrema, cmap=nothing) diff --git a/src/wcs.jl b/src/wcs.jl new file mode 100644 index 00000000..80541a74 --- /dev/null +++ b/src/wcs.jl @@ -0,0 +1,606 @@ +const WCS_HEADERS_TEMPLATES = [ + "DATE", + "MJD", + + + "WCSAXESa", + "WCAXna", + "WCSTna", + "WCSXna", + "CRPIXja", + "jCRPna", + "jCRPXn", + "TCRPna", + "TCRPXn", + "PCi_ja", + "ijPCna", + "TPn_ka", + "TPCn_ka", + "CDi_ja", + "ijCDna", + "TCn_ka", + "TCDn_ka", + "CDELTia", + "iCDEna", + "iCDLTn", + "TCDEna", + "TCDLTn", + "CROTAi", + "iCROTn", + "TCROTn", + "CUNITia", + "iCUNna", + "iCUNIn", + "TCUNna", + "TCUNIn", + "CTYPEia", + "iCTYna", + "iCTYPn", + "TCTYna", + "TCTYPn", + "CRVALia", + "iCRVna", + "iCRVLn", + "TCRVna", + "TCRVLn", + "LONPOLEa", + "LONPna", + "LATPOLEa", + "LATPna", + "RESTFREQ", + "RESTFRQa", + "RFRQna", + "RESTWAVa", + "RWAVna", + "PVi_ma", + "iVn_ma", + "iPVn_ma", + "TVn_ma", + "TPVn_ma", + "PROJPm", + "PSi_ma", + "iSn_ma", + "iPSn_ma", + "TSn_ma", + "TPSn_ma", + "VELREF", + "CNAMEia", + "iCNAna", + "iCNAMn", + "TCNAna", + "TCNAMn", + "CRDERia", + "iCRDna", + "iCRDEn", + "TCRDna", + "TCRDEn", + "CSYERia", + "iCSYna", + "iCSYEn", + "TCSYna", + "TCSYEn", + "CZPHSia", + "iCZPna", + "iCZPHn", + "TCZPna", + "TCZPHn", + "CPERIia", + "iCPRna", + "iCPERn", + "TCPRna", + "TCPERn", + "WCSNAMEa", + "WCSNna", + "TWCSna", + "TIMESYS", + "TREFPOS", + "TRPOSn", + "TREFDIR", + "TRDIRn", + "PLEPHEM", + "TIMEUNIT", + "DATEREF", + "MJDREF", + "MJDREFI", + "MJDREFF", + "JDREF", + "JDREFI", + "JDREFF", + "TIMEOFFS", + "DATE-OBS", + "DOBSn", + "DATE-BEG", + "DATE-AVG", + "DAVGn", + "DATE-END", + "MJD-OBS", + "MJDOBn", + "MJD-BEG", + "MJD-AVG", + "MJDAn", + "MJD-END", + "JEPOCH", + "BEPOCH", + "TSTART", + "TSTOP", + "XPOSURE", + "TELAPSE", + "TIMSYER", + "TIMRDER", + "TIMEDEL", + "TIMEPIXR", + "OBSGEO-X", + "OBSGXn", + "OBSGEO-Y", + "OBSGYn", + "OBSGEO-Z", + "OBSGZn", + "OBSGEO-L", + "OBSGLn", + "OBSGEO-B", + "OBSGBn", + "OBSGEO-H", + "OBSGHn", + "OBSORBIT", + "RADESYSa", + "RADEna", + "RADECSYS", + "EPOCH", + "EQUINOXa", + "EQUIna", + "SPECSYSa", + "SPECna", + "SSYSOBSa", + "SOBSna", + "VELOSYSa", + "VSYSna", + "VSOURCEa", + "VSOUna", + "ZSOURCEa", + "ZSOUna", + "SSYSSRCa", + "SSRCna", + "VELANGLa", + "VANGna", + "RSUN_REF", + "DSUN_OBS", + "CRLN_OBS", + "HGLN_OBS", + "HGLT_OBS", + "NAXISn", + "CROTAn", + "PROJPn", + "CPDISja", + "CQDISia", + "DPja", + "DQia", + "CPERRja", + "CQERRia", + "DVERRa", + "A_ORDER", + "B_ORDER", + "AP_ORDER", + "BP_ORDER", + "A_DMAX", + "B_DMAX", + "A_p_q", + "B_p_q", + "AP_p_q", + "BP_p_q", + "CNPIX1", + "PPO3", + "PPO6", + "XPIXELSZ", + "YPIXELSZ", + "PLTRAH", + "PLTRAM", + "PLTRAS", + "PLTDECSN", + "PLTDECD", + "PLTDECM", + "PLTDECS", + "PLATEID", + "AMDXm", + "AMDYm", + "WATi_m" +] + +# Expand the headers containing lower case specifers into N copies +# Find all lower case templates +const WCS_HEADERS = Set(mapreduce(vcat, WCS_HEADERS_TEMPLATES) do template + if any(islowercase, template) + template_chars = Vector{Char}(template) + chars = template_chars[islowercase.(template_chars)] + out = String[template] + for replace_target in chars + newout = String[] + for template in out + for i in [""; string.(1:4); string.('a':'d')] + push!(newout, replace(template, replace_target=>i)) + end + end + append!(out, newout) + end + out + else + template + end +end) + + + +""" + emptywcs() + +Given an AbstractArray, return a blank WCSTransform of the appropriate +dimensionality. +""" +emptywcs(data::AbstractArray) = WCSTransform(ndims(data)) +emptywcs(img::AstroImage) = WCSTransform(length(dims(img))+length(refdims(img))) + + + +# """ +# filterwcsheader(hdrs::FITSHeader) + +# Return a new FITSHeader containing WCS header from `hdrs`. +# This is useful for creating a new image with the same coordinates +# as another. +# """ +# function filterwcsheader(hdrs::FITSHeader) +# include_keys = intersect(keys(hdrs), WCS_HEADERS) +# return FITSHeader( +# include_keys, +# map(key -> hdrs[key], include_keys), +# map(key -> get_comment(hdrs, key), include_keys), +# ) +# end + +""" + wcsfromheader(img::AstroImage; relax=WCS.HDR_ALL, ignore_rejected=true) + +Helper function to create a WCSTransform from an array and +FITSHeaders. +""" +function wcsfromheader(img::AstroImage; relax=WCS.HDR_ALL) + # We only need to stringify WCS header. This might just be 4-10 header keywords + # out of thousands. + local wcsout + + io = IOBuffer() + serializeheader(io, header(img)) + hdrstr = String(take!(io)) + + + # Load the headers without ignoring rejected to get error messages + try + wcsout = WCS.from_header( + hdrstr; + ignore_rejected=false, + relax + ) + catch err + # Load them again ignoring error messages + wcsout = WCS.from_header( + hdrstr; + ignore_rejected=true, + relax + ) + # If that still fails, the use gets the stack trace here + # If not, print a warning about rejected header + @warn "WCSTransform was generated by ignoring rejected header. It may not be valid." exception=err + end + + if length(wcsout) == 0 + return [emptywcs(img)] + else + return wcsout + end +end + +# Map FITS stokes numbers to a symbol +function _stokes_symbol(i) + return if i == 1 + :I + elseif i == 2 + :Q + elseif i == 3 + :U + elseif i == 4 + :V + elseif i == -1 + :RR + elseif i == -2 + :LL + elseif i == -3 + :RL + elseif i == -4 + :LR + elseif i == -5 + :XX + elseif i == -6 + :YY + elseif i == -7 + :XY + elseif i == -8 + :YX + else + @warn "unknown FITS stokes number $i. See \"Representations of world coordinates in FITS\", Table 7." + nothing + end +end +function _stokes_name(symb) + return if symb == :I + "Stokes Unpolarized" + elseif symb == :Q + "Stokes Linear Q" + elseif symb == :U + "Stokes Linear U" + elseif symb == :V + "Stokes Circular" + elseif symb == :RR + "Right-right cicular" + elseif symb == :LL + "Left-left cicular" + elseif symb == :RL + "Right-left cross-cicular" + elseif symb == :LR + "Left-right cross-cicular" + elseif symb == :XX + "X parallel linear" + elseif symb == :YY + "Y parallel linear" + elseif symb == :XY + "XY cross linear" + elseif symb == :YX + "YX cross linear" + else + @warn "unknown FITS stokes key $symb. See \"Representations of world coordinates in FITS\", Table 7." + "" + end +end + + + + +# Smart versions of pix_to_world and world_to_pix +""" + pix_to_world(img::AstroImage, pixcoords; all=false) + +Given an astro image, look up the world coordinates of the pixels given +by `pixcoords`. World coordinates are resolved using WCS.jl and a +WCSTransform calculated from any FITS header present in `img`. If +no WCS information is in the header, or the axes are all linear, this will +just return pixel coordinates. + +`pixcoords` should be the coordinates in your current selection +of the image. For example, if you select a slice like this: +```julia +julia> cube = load("some-3d-cube.fits") +julia> slice = cube[10:20, 30:40, 5] +``` + +Then to look up the coordinates of the pixel in the bottom left corner of +`slice`, run: +```julia +julia> world_coords = pix_to_world(img, [1, 1]) +[10, 30] +``` + +If WCS information was present in the header of `cube`, then those coordinates +would be resolved using axis 1, 2, and 3 respectively. + +To include world coordinates in all axes, pass `all=true` +```julia +julia> world_coords = pix_to_world(img, [1, 1], all=true) +[10, 30, 5] +``` + +!! Coordinates must be provided in the order of `dims(img)`. If you transpose +an image, the order you pass the coordinates should not change. +""" +function WCS.pix_to_world(img::AstroImage, pixcoords; wcsn=1, all=false, parent=false) + if pixcoords isa Array{Float64} + pixcoords_prepared = pixcoords + else + pixcoords_prepared = [Float64(c) for c in pixcoords] + end + D_out = wcs(img,wcsn).naxis + if ndims(pixcoords_prepared) > 1 + worldcoords_out = similar(pixcoords_prepared, Float64, D_out, size(pixcoords_prepared,2)) + else + worldcoords_out = similar(pixcoords_prepared, Float64, D_out) + end + + # Find the coordinates in the parent array. + # Dimensional data + # pixcoords_floored = floor.(Int, pixcoords) + # pixcoords_frac = (pixcoords .- pixcoords_floored) .* step.(dims(img)) + # parentcoords = getindex.(dims(img), pixcoords_floored) .+ pixcoords_frac + if parent + parentcoords = pixcoords + else + parentcoords = pixcoords .* step.(dims(img)) .+ first.(dims(img)) + end + # WCS.jl is very restrictive. We need to supply a Vector{Float64} + # as input, not any other kind of collection. + # TODO: avoid allocation in case where refdims=() and pixcoords isa Array{Float64} + if ndims(pixcoords_prepared) > 1 + parentcoords_prepared = zeros(wcs(img,wcsn).naxis, size(pixcoords_prepared,2)) + else + parentcoords_prepared = zeros(wcs(img,wcsn).naxis) + end + # out = zeros(Float64, wcs(img,wcsn).naxis, size(pixcoords,2)) + for (i, dim) in enumerate(dims(img)) + j = wcsax(img, dim) + parentcoords_prepared[j,:] .= parentcoords[i,:] .- 1 + end + for dim in refdims(img) + j = wcsax(img, dim) + # Non numeric reference dims can be used, e.g. a polarization axis of symbols I, Q, U, etc. + if eltype(dim) <: Number + z = dim[1] - 1 + else + # Find the index of the symbol into the parent cube + parentrefdim = img.wcsdims[findfirst(d->name(d)==name(dim), img.wcsdims)] + z = findfirst(==(first(dim)), collect(parentrefdim)) - 1 + end + parentcoords_prepared[j,:] .= z + end + + # Get world coordinates along all slices + WCS.pix_to_world!(wcs(img, wcsn), parentcoords_prepared, worldcoords_out) + + # If user requested world coordinates in all dims, not just selected + # dims of img + if all + return worldcoords_out + end + + # Otherwise filter to only return coordinates along selected dims. + if ndims(pixcoords_prepared) > 1 + world_coords_of_these_axes = zeros(length(dims(img)), size(pixcoords_prepared,2)) + else + world_coords_of_these_axes = zeros(length(dims(img))) + end + for (i, dim) in enumerate(dims(img)) + j = wcsax(img, dim) + world_coords_of_these_axes[i,:] .= worldcoords_out[j,:] + end + + return world_coords_of_these_axes +end + + +## +function WCS.world_to_pix(img::AstroImage, worldcoords; parent=false, wcsn=1) + if worldcoords isa Array{Float64} + worldcoords_prepared = worldcoords + else + worldcoords_prepared = [Float64(c) for c in worldcoords] + end + D_out = wcs(img,wcsn).naxis + if ndims(worldcoords_prepared) > 1 + out = similar(worldcoords_prepared, Float64, D_out, size(worldcoords_prepared,2)) + else + out = similar(worldcoords_prepared, Float64, D_out) + end + return WCS.world_to_pix!(out, img, worldcoords_prepared; parent) +end +function WCS.world_to_pix!(pixcoords_out, img::AstroImage, worldcoords; wcsn=1, parent=false) + # # Find the coordinates in the parent array. + # # Dimensional data + # worldcoords_floored = floor.(Int, worldcoords) + # worldcoords_frac = (worldcoords .- worldcoords_floored) .* step.(dims(img)) + # parentcoords = getindex.(dims(img), worldcoords_floored) .+ worldcoords_frac + # WCS.jl is very restrictive. We need to supply a Vector{Float64} + # as input, not any other kind of collection. + # TODO: avoid allocation in case where refdims=() and worldcoords isa Array{Float64} + if ndims(worldcoords) > 1 + worldcoords_prepared = zeros(wcs(img,wcsn).naxis,size(worldcoords,2)) + else + worldcoords_prepared = zeros(wcs(img,wcsn).naxis) + end + # TODO: we need to pass in ref dims locations as well, and then filter the + # output to only include the dims of the current slice? + # out = zeros(Float64, wcs(img,wcsn).naxis, size(worldcoords,2)) + for (i, dim) in enumerate(dims(img)) + j = wcsax(img, dim) + worldcoords_prepared[j,:] = worldcoords[i,:] + end + for dim in refdims(img) + j = wcsax(img, dim) + # Non numeric reference dims can be used, e.g. a polarization axis of symbols I, Q, U, etc. + if eltype(dim) <: Number + z = dim[1] + else + # Find the index of the symbol into the parent cube + parentrefdim = img.wcsdims[findfirst(d->name(d)==name(dim), img.wcsdims)] + z = findfirst(==(first(dim)),collect(parentrefdim)) -1 + end + worldcoords_prepared[j,:] .= z + end + + # This returns the parent pixel coordinates. + # TODO: switch to non-allocating version. + pixcoords_out .= WCS.world_to_pix(wcs(img, wcsn), worldcoords_prepared) + + if !parent + coordoffsets = zeros(wcs(img,wcsn).naxis) + coordsteps = zeros(wcs(img,wcsn).naxis) + for (i, dim) in enumerate(dims(img)) + j = wcsax(img, dim) + coordoffsets[j] = first(dims(img)[i]) + coordsteps[j] = step(dims(img)[i]) + end + for dim in refdims(img) + j = wcsax(img, dim) + # Non numeric reference dims can be used, e.g. a polarization axis of symbols I, Q, U, etc. + if eltype(dim) <: Number + coordoffsets[j] = first(dim) + coordsteps[j] = step(dim) + else + # Find the index of the symbol into the parent cube + parentrefdim = img.wcsdims[findfirst(d->name(d)==name(dim), img.wcsdims)] + z = findfirst(==(first(dim)),collect(parentrefdim)) + coordoffsets[j] = z - 1 + coordsteps[j] = 1 + end + end + + pixcoords_out .-= coordoffsets + pixcoords_out .= (pixcoords_out .+ 1) ./ coordsteps + end + return pixcoords_out +end + + + + + + +## For now, we use a copied version of FITSIO's show method for FITSHeader. +# We have to be careful to format things in a way WCSLib will like. +# In particular, we can't put newlines after each 80 characters. +# FITSIO has to do this so users can see the header. + +# functions for displaying header values in show(io, header) +hdrval_repr(v::Bool) = v ? "T" : "F" +hdrval_repr(v::String) = @sprintf "'%-8s'" v +hdrval_repr(v::Union{AbstractFloat, Integer}) = string(v) + +function serializeheader(io, hdr::FITSHeader) + n = length(hdr) + for i=1:n + if hdr.keys[i] == "COMMENT" || hdr.keys[i] == "HISTORY" + lastc = min(72, lastindex(hdr.comments[i])) + @printf io "%s %s" hdr.keys[i] hdr.comments[i][1:lastc] + print(io, " "^(72-lastc)) + else + @printf io "%-8s" hdr.keys[i] + if hdr.values[i] === nothing + print(io, " ") + rc = 50 # remaining characters on line + elseif hdr.values[i] isa String + val = hdrval_repr(hdr.values[i]) + @printf io "= %-20s" val + rc = length(val) <= 20 ? 50 : 70 - length(val) + else + val = hdrval_repr(hdr.values[i]) + @printf io "= %20s" val + rc = length(val) <= 20 ? 50 : 70 - length(val) + end + if length(hdr.comments[i]) > 0 + lastc = min(rc-3, lastindex(hdr.comments[i])) + @printf io " / %s" hdr.comments[i][1:lastc] + rc -= lastc + 3 + end + print(io, " "^rc) + end + if i == n + print(io, "\nEND"*(" "^77)) + else + print(io) + end + end +end diff --git a/test/ccd2rgb.jl b/test/ccd2rgb.jl index 70ee6739..f65fccd4 100644 --- a/test/ccd2rgb.jl +++ b/test/ccd2rgb.jl @@ -45,8 +45,8 @@ end @test isapprox(blue.(asinh_res), blue.(asinh_ans), nans = true, rtol = 3e-5) @test isapprox(green.(asinh_res), green.(asinh_ans), nans = true, rtol = 3e-5) - @testset "AstroImage using ccd2rgb and properties" begin - img = AstroImage(RGB, (joinpath("data","casa_0.5-1.5keV.fits"), joinpath("data","casa_1.5-3.0keV.fits"), + @testset "AstroImageMat using ccd2rgb and properties" begin + img = AstroImageMat(RGB, (joinpath("data","casa_0.5-1.5keV.fits"), joinpath("data","casa_1.5-3.0keV.fits"), joinpath("data","casa_4.0-6.0keV.fits"))) @test RGB.(img.property.rgb_image) isa Array{RGB{Float64},2} @@ -85,7 +85,7 @@ end @test isapprox(blue.(img.property.rgb_image), blue.(ans), nans = true) - img = AstroImage(joinpath("data","casa_0.5-1.5keV.fits")) + img = AstroImageMat(joinpath("data","casa_0.5-1.5keV.fits")) @test_throws DomainError set_brightness!(img, 1.2) @test_throws DomainError set_contrast!(img, 1.2) end diff --git a/test/plots.jl b/test/plots.jl index 1ea2bcba..8678272c 100644 --- a/test/plots.jl +++ b/test/plots.jl @@ -3,7 +3,7 @@ using AstroImages: pix2world_xformatter, pix2world_yformatter @testset "Plot recipes" begin data = randn(10, 10) - img = AstroImage(data) + img = AstroImageMat(data) wcs1 = WCSTransform(2; ctype = ["RA---AIR", "DEC--AIR"]) wcs2 = WCSTransform(2; ctype = ["GLON--", "GLAT--"]) wcs3 = WCSTransform(2; ctype = ["TLON--", "TLAT--"]) diff --git a/test/runtests.jl b/test/runtests.jl index 98545eaf..43c9bd31 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,9 +1,6 @@ -using AstroImages, FITSIO, Images, Random, Widgets, WCS, JLD +using AstroImages, FITSIO, Random, WCS, ImageBase, Test, WCS, Statistics -using Test, WCS -using SHA: sha256 - -import AstroImages: _float, render, _brightness_contrast, brightness_contrast +import AstroImages: _float, render @testset "Conversion to float and fixed-point" begin @testset "Float" begin @@ -37,29 +34,15 @@ end FITS(fname, "w") do f write(f, data) end - @test load(fname, 1)[1] == data - @test load(fname, (1, 1))[1] == (data, data) + @test load(fname, 1) == data + @test load(fname, (1, 1)) == (data, data) img = AstroImage(fname) - rendered_img = colorview(img) - @test iszero(minimum(rendered_img)) + rendered_img = imview(img) + @test eltype(rendered_img) <: RGBA end rm(fname, force = true) end -@testset "Control contrast" begin - @test @inferred(_brightness_contrast(Gray, ones(Float32, 2, 2), 100, 100)) isa - Array{Gray{Float32},2} - Random.seed!(1) - M = rand(2, 2) - @test @inferred(_brightness_contrast(Gray, M, 100, 100)) ≈ - Gray.([0.48471895908315565 0.514787046406301; - 0.5280458879184159 0.39525854250610226]) rtol = 1e-12 - @test @inferred(_brightness_contrast(Gray, M, 0, 255)) == Gray.(M) - @test @inferred(_brightness_contrast(Gray, M, 255, 0)) == Gray.(ones(size(M))) - @test @inferred(_brightness_contrast(Gray, M, 0, 0)) == Gray.(zeros(size(M))) - @test brightness_contrast(AstroImage(M)) isa Widgets.Widget{:manipulate,Any} -end - @testset "default handler" begin fname = tempname() * ".fits" @testset "less dimensions than 2" begin @@ -67,7 +50,7 @@ end FITS(fname, "w") do f write(f, data) end - @test_throws ErrorException AstroImage(fname) + @test ndims(AstroImage(fname)) == 1 end @testset "no ImageHDU" begin @@ -86,23 +69,20 @@ end FITS(fname, "w") do f write(f, indata; varcols=["vcol", "VCOL"]) - @test_throws MethodError AstroImage(f) + @test_throws Exception AstroImage(f) end end @testset "Opening AstroImage in different ways" begin data = rand(2,2) - wcs = WCSTransform(2;) FITS(fname, "w") do f write(f, data) end f = FITS(fname) + header = read_header(f[1]) @test AstroImage(fname, 1) isa AstroImage - @test AstroImage(Gray ,fname, 1) isa AstroImage - @test AstroImage(Gray, f, 1) isa AstroImage - @test AstroImage(data, wcs) isa AstroImage - @test AstroImage((data,data), (wcs,wcs)) isa AstroImage - @test AstroImage(Gray, data, wcs) isa AstroImage + @test AstroImage(f, 1) isa AstroImage + @test AstroImage(data, header) isa AstroImage close(f) end @@ -131,126 +111,246 @@ end end @testset "Utility functions" begin - @test size(AstroImage((rand(10,10), rand(10,10)))) == ((10,10), (10,10)) - @test length(AstroImage((rand(10,10), rand(10,10)))) == 2 -end - -@testset "multi image AstroImage" begin - data1 = rand(10,10) - data2 = rand(10,10) - fname = tempname() * ".fits" - FITS(fname, "w") do f - write(f, data1) - write(f, data2) - end - - img = AstroImage(fname, (1,2)) - @test length(img.data) == 2 - @test img.data[1] == data1 - @test img.data[2] == data2 - - f = FITS(fname) - img = AstroImage(Gray, f, (1,2)) - @test length(img.data) == 2 - @test img.data[1] == data1 - @test img.data[2] == data2 - close(f) + @test size(AstroImage(rand(10,10))) == (10,10) + @test length(AstroImage(rand(10,10))) == 100 end @testset "multi wcs AstroImage" begin fname = tempname() * ".fits" f = FITS(fname, "w") - inhdr = FITSHeader(["CTYPE1", "CTYPE2", "RADESYS", "FLTKEY", "INTKEY", "BOOLKEY", "STRKEY", "COMMENT", - "HISTORY"], - ["RA---TAN", "DEC--TAN", "UNK", 1.0, 1, true, "string value", nothing, nothing], - ["", - "", - "", - "floating point keyword", - "", - "boolean keyword", - "string value", - "this is a comment", - "this is a history"]) + inhdr = FITSHeader([ + "FLTKEY", "INTKEY", "BOOLKEY", "STRKEY", "COMMENT", "HISTORY", + "CRVAL1a", + "CRVAL2a", + "CRPIX1a", + "CRPIX2a", + "CDELT1a", + "CDELT2a", + "CTYPE1a", + "CTYPE2a", + "CUNIT1a", + "CUNIT2a", + + "CRVAL1b", + "CRVAL2b", + "CRPIX1b", + "CRPIX2b", + "CDELT1b", + "CDELT2b", + "CTYPE1b", + "CTYPE2b", + "CUNIT1b", + "CUNIT2b", + ], + [ + 1.0, 1, true, "string value", nothing, nothing, + 0.5, + 89.5, + 1, + 1, + 1, + -1, + "RA---TAN", + "DEC--TAN", + "deg ", + "deg ", + + 0.5, + 89.5, + 1, + 1, + 1, + -1, + "RA---TAN", + "DEC--TAN", + "deg ", + "deg ", + ], + [ + "floating point keyword", "", "boolean keyword", "string value", "this is a comment", "this is a history", + "", + "", + "", + "", + "", + "", + "Terrestrial East Longitude", + "Terrestrial North Latitude", + "", + "", + + "", + "", + "", + "", + "", + "", + "Terrestrial East Longitude", + "Terrestrial North Latitude", + "", + "", + ]) indata = reshape(Float32[1:100;], 5, 20) write(f, indata; header=inhdr) - write(f, indata; header=inhdr) close(f) - img = AstroImage(fname, (1,2)) + img = AstroImage(fname) f = FITS(fname) - @test length(img.wcs) == 2 - @test WCS.to_header(img.wcs[1]) === WCS.to_header(WCS.from_header(read_header(f[1], String))[1]) - @test WCS.to_header(img.wcs[2]) === WCS.to_header(WCS.from_header(read_header(f[2], String))[1]) - - img = AstroImage(Gray, f, (1,2)) - @test length(img.wcs) == 2 - @test WCS.to_header(img.wcs[1]) === WCS.to_header(WCS.from_header(read_header(f[1], String))[1]) - @test WCS.to_header(img.wcs[2]) === WCS.to_header(WCS.from_header(read_header(f[2], String))[1]) + @test length(wcs(img)) == 2 + @test WCS.to_header(wcs(img,1)) === WCS.to_header(WCS.from_header(read_header(f[1], String))[1]) + @test WCS.to_header(wcs(img,2)) === WCS.to_header(WCS.from_header(read_header(f[1], String))[2]) + + img = AstroImage(f) + @test length(wcs(img)) == 2 + @test WCS.to_header(wcs(img,1)) === WCS.to_header(WCS.from_header(read_header(f[1], String))[1]) + @test WCS.to_header(wcs(img,2)) === WCS.to_header(WCS.from_header(read_header(f[1], String))[2]) close(f) end -@testset "multi file AstroImage" begin - fname1 = tempname() * ".fits" - f = FITS(fname1, "w") - inhdr = FITSHeader(["CTYPE1", "CTYPE2", "RADESYS", "FLTKEY", "INTKEY", "BOOLKEY", "STRKEY", "COMMENT", - "HISTORY"], - ["RA---TAN", "DEC--TAN", "UNK", 1.0, 1, true, "string value", nothing, nothing], - ["", - "", - "", - "floating point keyword", - "", - "boolean keyword", - "string value", - "this is a comment", - "this is a history"]) - - indata1 = reshape(Int[1:100;], 5, 20) - write(f, indata1; header=inhdr) - close(f) +## +@testset "imview" begin - fname2 = tempname() * ".fits" - f = FITS(fname2, "w") - indata2 = reshape(Int[1:100;], 5, 20) - write(f, indata2; header=inhdr) - close(f) + arr1 = permutedims(reshape(1:9,3,3)) + img = AstroImage(arr1) - fname3 = tempname() * ".fits" - f = FITS(fname3, "w") - indata3 = reshape(Int[1:100;], 5, 20) - write(f, indata3; header=inhdr) - close(f) + @test imview(arr1) == imview(img) + @test imview(img) isa AstroImage + @test !(imview(arr1) isa AstroImage) + + img_rendered_1 = imview(img, clims=(1,9), stretch=identity, contrast=1, bias=0.5, cmap=nothing) - img = AstroImage((fname1, fname2, fname3)) - f1 = FITS(fname1) - f2 = FITS(fname2) - f3 = FITS(fname3) - - @test length(img.data) == length(img.wcs) == 3 - @test img.data[1] == indata1 - @test img.data[2] == indata2 - @test img.data[3] == indata3 - @test WCS.to_header(img.wcs[1]) == WCS.to_header(img.wcs[2]) == - WCS.to_header(img.wcs[3]) == WCS.to_header(WCS.from_header(read_header(f1[1], String))[1]) - @test eltype(eltype(img.data)) == Int - - img = AstroImage(Gray, (f1, f2, f3), (1,1,1)) - @test length(img.data) == length(img.wcs) == 3 - @test img.data[1] == indata1 - @test img.data[2] == indata2 - @test img.data[3] == indata3 - @test WCS.to_header(img.wcs[1]) == WCS.to_header(img.wcs[2]) == - WCS.to_header(img.wcs[3]) == WCS.to_header(WCS.from_header(read_header(f1[1], String))[1]) - @test eltype(eltype(img.data)) == Int - close(f1) - close(f2) - close(f3) - rm(fname1, force = true) - rm(fname2, force = true) - rm(fname3, force = true) + # Image Orientation + @test CartesianIndex(3,1) == argmin(Gray.(img_rendered_1)) + @test CartesianIndex(1,3) == argmax(Gray.(img_rendered_1)) + + # Rendering Basics + @test allunique(img_rendered_1) + # It is intended that the rendered image is flipped vs it's data + @test img_rendered_1[3,1] == RGBA(0,0,0,1) + @test img_rendered_1[1,3] == RGBA(1,1,1,1) + @test all(p -> p.r==p.g==p.b && p.alpha==1, img_rendered_1) + + # Limits + img_rendered_2 = imview(img, clims=(3,7), stretch=identity, contrast=1, bias=0.5, cmap=nothing) + @test length(unique(img_rendered_2)) == 5 + @test count(==(RGBA(0,0,0,1)), img_rendered_2) == 3 + @test count(==(RGBA(1,1,1,1)), img_rendered_2) == 3 + + # Calculated limits + @test img_rendered_1 == imview(img, clims=extrema, stretch=identity, contrast=1, bias=0.5, cmap=nothing) + img_rendered_3 = imview(img, clims=Zscale(), stretch=identity, contrast=1, bias=0.5, cmap=nothing) + img_rendered_4 = imview(img, clims=Percent(100), stretch=identity, contrast=1, bias=0.5, cmap=nothing) + @test img_rendered_1 == img_rendered_3 + @test img_rendered_1 == img_rendered_4 + + # Stretching + for stretchfunc in (sqrtstretch, asinhstretch, powerdiststretch, logstretch, powstretch, squarestretch, sinhstretch) + img_rendered_5 = imview(arr1, clims=(1,9), stretch=stretchfunc, contrast=1, bias=0.5, cmap=nothing) + @test extrema(Gray.(img_rendered_5)) == (0,1) + manual_stretch = stretchfunc.(AstroImages.clampednormedview(arr1,(1,9))) + @test Gray.(img_rendered_5) ≈ + N0f8.((manual_stretch.-minimum(manual_stretch)) ./ + (maximum(manual_stretch)-minimum(manual_stretch)))'[end:-1:begin,:] + end + + # Contrast/Bias + @test Gray.(imview(img, clims=extrema, stretch=identity, contrast=1, bias=0.6, cmap=nothing)) == + N0f8.(clamp.(N0f8.(Gray.(img_rendered_1)) .- 0.1,false,true)) + + img_rendered_5 = imview(arr1, clims=(1,9), stretch=sqrtstretch, contrast=0.5, bias=0.5, cmap=nothing) + + # Missing/NaN + for m in (NaN, missing) + arr2 = [ + 1 2 3 + 4 m 6 + 7 8 9 + ] + @test imview(arr2) == imview(AstroImage(arr2)) + @test imview(arr2)[2,2].alpha == 0 + @test 8 == count(img_rendered_1 .== imview(arr2, clims=(1,9), stretch=identity, contrast=1, bias=0.5, cmap=nothing)) + end + + img_rendered_6 = imview([1, 2, NaN, missing, -Inf, Inf], clims=extrema) + img_rendered_6b = imview([1, 2], clims=extrema) + + @test img_rendered_6[1] == img_rendered_6b[1] + @test img_rendered_6[2] == img_rendered_6b[2] + @test img_rendered_6[1].alpha == 1 + @test img_rendered_6[2].alpha == 1 + @test img_rendered_6[3].alpha == 0 + @test img_rendered_6[4].alpha == 0 + @test img_rendered_6[5] == RGBA(0,0,0,1) + @test img_rendered_6[6] == RGBA(1,1,1,1) end -include("plots.jl") -include("ccd2rgb.jl") +## + +1 + +## + +# @testset "multi file AstroImage" begin +# fname1 = tempname() * ".fits" +# f = FITS(fname1, "w") +# inhdr = FITSHeader(["CTYPE1", "CTYPE2", "RADESYS", "FLTKEY", "INTKEY", "BOOLKEY", "STRKEY", "COMMENT", +# "HISTORY"], +# ["RA---TAN", "DEC--TAN", "UNK", 1.0, 1, true, "string value", nothing, nothing], +# ["", +# "", +# "", +# "floating point keyword", +# "", +# "boolean keyword", +# "string value", +# "this is a comment", +# "this is a history"]) + +# indata1 = reshape(Int[1:100;], 5, 20) +# write(f, indata1; header=inhdr) +# close(f) + +# fname2 = tempname() * ".fits" +# f = FITS(fname2, "w") +# indata2 = reshape(Int[1:100;], 5, 20) +# write(f, indata2; header=inhdr) +# close(f) + +# fname3 = tempname() * ".fits" +# f = FITS(fname3, "w") +# indata3 = reshape(Int[1:100;], 5, 20) +# write(f, indata3; header=inhdr) +# close(f) + +# img = AstroImage((fname1, fname2, fname3)) +# f1 = FITS(fname1) +# f2 = FITS(fname2) +# f3 = FITS(fname3) + +# @test length(img.data) == length(img.wcs) == 3 +# @test img.data[1] == indata1 +# @test img.data[2] == indata2 +# @test img.data[3] == indata3 +# @test WCS.to_header(img.wcs[1]) == WCS.to_header(img.wcs[2]) == +# WCS.to_header(img.wcs[3]) == WCS.to_header(WCS.from_header(read_header(f1[1], String))[1]) +# @test eltype(eltype(img.data)) == Int + +# img = AstroImage(Gray, (f1, f2, f3), (1,1,1)) +# @test length(img.data) == length(img.wcs) == 3 +# @test img.data[1] == indata1 +# @test img.data[2] == indata2 +# @test img.data[3] == indata3 +# @test WCS.to_header(img.wcs[1]) == WCS.to_header(img.wcs[2]) == +# WCS.to_header(img.wcs[3]) == WCS.to_header(WCS.from_header(read_header(f1[1], String))[1]) +# @test eltype(eltype(img.data)) == Int +# close(f1) +# close(f2) +# close(f3) +# rm(fname1, force = true) +# rm(fname2, force = true) +# rm(fname3, force = true) +# end + +# include("plots.jl") +# include("ccd2rgb.jl")