ccd2rgb function (#15)

* ccd2rgb function added

* tests for ccd2rgb

* performance fixes and examples

* add JLD to test only deps.

* minor fixes

* tests reduced, travis cache added and optimizations

* stretch function test log replaced by asinh

* reducing tolerance

* Check SHA256 hash of the downloaded files

* tests fix

* atol changed to rtol

* new dispatch

.gitignore

@@ -4,3 +4,6 @@
 
 # generated files
 Manifest.toml
+
+# test files
+/test/data

.travis.yml

@@ -17,6 +17,9 @@ git:
 matrix:
  allow_failures:
  - julia: nightly
+cache:
+  directories:
+  - $TRAVIS_BUILD_DIR/test/data
 
 ## uncomment and modify the following lines to manually install system packages
 #addons:

Project.toml

@@ -10,15 +10,19 @@ ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
 Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
 Interact = "c601a237-2ae4-5e1e-952c-7a85b0c7eef1"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+Reproject = "d1dcc2e6-806e-11e9-2897-3f99785db2ae"
 WCS = "15f3aee2-9e10-537f-b834-a6fb8bdb944d"
 
 [compat]
 julia = "^1.0.0"
+Reproject = "^0.3.0"
 
 [extras]
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Widgets = "cc8bc4a8-27d6-5769-a93b-9d913e69aa62"
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+JLD = "4138dd39-2aa7-5051-a626-17a0bb65d9c8"
+SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 
 [targets]
-test = ["Test", "Random", "Widgets"]
+test = ["Test", "Random", "Widgets", "JLD", "SHA"]

src/AstroImages.jl

@@ -2,9 +2,9 @@ __precompile__()
 
 module AstroImages
 
-using FITSIO, FileIO, Images, Interact, WCS
+using FITSIO, FileIO, Images, Interact, Reproject, WCS
 
-export load, AstroImage
+export load, AstroImage, ccd2rgb
 
 _load(fits::FITS, ext::Int) = read(fits[ext])
 _load(fits::FITS, ext::NTuple{N, Int}) where {N} = ntuple(i-> read(fits[ext[i]]), N)
@@ -149,5 +149,6 @@ Base.length(img::AstroImage{T,C,N}) where {T,C,N} = N
 
 include("showmime.jl")
 include("plot-recipes.jl")
+include("ccd2rgb.jl")
 
 end # module

src/ccd2rgb.jl

@@ -0,0 +1,48 @@
+"""
+    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]))
+
+Converts 3 grayscale ImageHDU into RGB by reprojecting them.
+
+# Arguments
+- `red`: Red channel data.
+- `green`: Green channel data.
+- `blue`: Blue channel data.
+- `stretch`: Stretch function applied.
+- `shape_out`: Shape of output RGB 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)
+```
+"""
+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
+
+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)
+

test/ccd2rgb.jl

@@ -0,0 +1,47 @@
+function download_dep(orig, dest, hash)
+    dest_file = joinpath("data", dest)
+    if isfile(dest_file)
+        dest_hash = open(dest_file, "r") do f
+            bytes2hex(sha256(f))
+        end
+        if dest_hash == hash
+            return nothing
+        end
+    end
+    mkpath("data")
+    download(orig, dest_file)
+    return nothing
+end
+
+@testset "ccd2rgb" begin
+    download_dep("http://chandra.harvard.edu/photo/2009/casa/fits/casa_0.5-1.5keV.fits", "casa_0.5-1.5keV.fits",
+                "5794b9ebced6b991a3e53888d129a38fbf4309250112be530cb6442be812dea6")
+    download_dep("http://chandra.harvard.edu/photo/2009/casa/fits/casa_1.5-3.0keV.fits", "casa_1.5-3.0keV.fits",
+                "a48b2502ceb979dfad0d05fd5ec19bf3e197ff2d1d9c604c9340992d1bf7eec9")
+    download_dep("http://chandra.harvard.edu/photo/2009/casa/fits/casa_4.0-6.0keV.fits", "casa_4.0-6.0keV.fits",
+                "15e90a14515c121c2817e97b255c604ad019c9c2340fda4fb6c5c3da55e1b0c2")
+    download_dep("https://bintray.com/aquatiko/AstroImages.jl/download_file?file_path=ccd2rgb_rounded.jld","ccd2rgb_rounded.jld",
+                "5191e59e527c3667486c680e92c8f77fcdbed1e82d3230317a514d928092107d")
+
+    f(x) = isnan(x) ? RGB.(0.0,0.0,0.0) : x
+    r = FITS(joinpath("data","casa_0.5-1.5keV.fits"))[1]
+    b = FITS(joinpath("data","casa_1.5-3.0keV.fits"))[1]
+    g = FITS(joinpath("data","casa_4.0-6.0keV.fits"))[1]
+    linear_res = ccd2rgb(r, b, g, shape_out = (1000,1000))
+    asinh_res = ccd2rgb(r, b, g, shape_out = (1000,1000), stretch = asinh)
+    linear_res = f.(RGB.(colorview(RGB, round.(red.(linear_res), digits = 8), round.(green.(linear_res), digits = 8), 
+                                                round.(blue.(linear_res), digits = 8))))
+    asinh_res = f.(RGB.(colorview(RGB, round.(red.(asinh_res), digits = 8), round.(green.(asinh_res), digits = 8), 
+                                                round.(blue.(asinh_res), digits = 8))))
+
+    linear_ans = f.(load(joinpath("data","ccd2rgb_rounded.jld"), "linear"))
+    asinh_ans = f.(load(joinpath("data","ccd2rgb_rounded.jld"), "asinh"))
+
+    @test isapprox(red.(linear_res), red.(linear_ans), nans = true, rtol = 3e-5)
+    @test isapprox(blue.(linear_res), blue.(linear_ans), nans = true, rtol = 3e-5)
+    @test isapprox(green.(linear_res), green.(linear_ans), nans = true, rtol = 3e-5)
+
+    @test isapprox(red.(asinh_res), red.(asinh_ans), nans = true, rtol = 3e-5)
+    @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)
+end

test/runtests.jl

@@ -1,6 +1,7 @@
-using AstroImages, FITSIO, Images, Random, Widgets, WCS
+using AstroImages, FITSIO, Images, Random, Widgets, WCS, JLD
 
 using Test, WCS
+using SHA: sha256
 
 import AstroImages: _float, render, _brightness_contrast, brightness_contrast
 
@@ -252,3 +253,4 @@ end
 end
 
 include("plots.jl")
+include("ccd2rgb.jl")