diff --git a/Project.toml b/Project.toml index 9eff46b..ba4e474 100644 --- a/Project.toml +++ b/Project.toml @@ -6,18 +6,20 @@ version = "0.4.2" [deps] Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] StaticArrays = "0.5, 0.6, 0.7, 0.8, 0.9, 0.10, 0.11, 0.12, 1.0" +Interpolations = "0.13.4" julia = "1" [extras] -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] test = ["Test"] diff --git a/src/AstroLib.jl b/src/AstroLib.jl index 2a09ec2..9c67ae4 100644 --- a/src/AstroLib.jl +++ b/src/AstroLib.jl @@ -6,6 +6,7 @@ __precompile__() module AstroLib using Dates, LinearAlgebra, Printf using StaticArrays +using Interpolations # Note on function definitions in this package. Most functions are defined as # follows: diff --git a/src/frebin.jl b/src/frebin.jl new file mode 100644 index 0000000..6f9a732 --- /dev/null +++ b/src/frebin.jl @@ -0,0 +1,176 @@ +# This file is a part of AstroLib.jl. License is MIT "Expat". +# Copyright (C) 2016 Mosè Giordano. + +function _frebin(image::AbstractArray{T}, nsout::S, nlout::S, total::Bool) where {T<:AbstractFloat,S<:Integer} + + # Determine the size of the input array + ns = size(image, 1) + nl = length(image)/ns + + # Determine if the new sizes are integral factors of the original sizes + sbox = ns/nsout + lbox = nl/nlout + + # Contraction by an integral amount + if (sbox == round(Int, sbox)) && (lbox == round(Int, lbox)) && (ns % nsout == 0) && (nl % nlout == 0) + array_shaped = reshape(image, (Int(sbox), nsout, Int(lbox), nlout)) + return dropdims(sum(array_shaped, dims=(1,3)), dims=(1,3)) .* (total ? 1.0 : 1 ./ (sbox.*lbox)) + end + + # Expansion by an integral amount + if (nsout % ns == 0) && (nlout % nl == 0) + xindex = (1:nsout) / (nsout/ns) + if isone(nl) # 1D case, linear interpolation + interpfunc = extrapolate(interpolate(image, BSpline(Constant())), Flat()) + return interpfunc(xindex) * (total ? sbox : 1.) + end + yindex = (1:nlout) / (nlout/nl) + interpfunc = extrapolate(interpolate(image, BSpline(Constant())), Flat()) + return [interpfunc(x, y) for x in xindex, y in yindex] .* (total ? sbox.*lbox : 1.) + end + + ns1 = ns-1 + nl1 = nl-1 + + # Do 1D case separately + if isone(nl) + result = zeros(eltype(image), nsout) + for i ∈ 0:nsout-1 + rstart = i*sbox # starting position for each box + istart = floor(Int, rstart) + rstop = rstart + sbox # ending position for each box + istop = Int(clamp(floor(rstop), 0, ns1)) + frac1 = rstart-istart + frac2 = 1.0 - (rstop-istop) + + # add pixel values from istart to istop and subtract fractional pixel from istart to start and + # fractional pixel from rstop to istop + + result[i+1] = sum(image[istart+1:istop+1]) - frac1*image[istart+1] - frac2*image[istop+1] + end + return result .* (total ? 1.0 : 1 ./ (sbox.*lbox)) + end + + # Now, do 2D case + # First, bin second dimension + temp = zeros(eltype(image), ns, nlout) + # Loop on output image lines + for i ∈ 0:nlout-1 + rstart = i*lbox # starting position for each box + istart = floor(Int, rstart) + rstop = rstart + lbox + istop = Int(clamp(floor(rstop), 0, nl1)) + frac1 = rstart-istart + frac2 = 1.0 - (rstop-istop) + + # add pixel values from istart to istop and subtract fractional pixel from istart to start and + # fractional pixel from rstop to istop + + if istart == istop + temp[:,i+1] .= (1 .- frac1 .- frac2).*image[:,istart+1] + else + temp[:,i+1] .= dropdims(sum(image[:,istart+1:istop+1], dims=2), dims=2) .- frac1.*image[:,istart+1] .- frac2.*image[:,istop+1] + end + end + temp = temp' + # Bin in first dimension + result = zeros(eltype(image), nlout, nsout) + # Loop on output image samples + for i ∈ 0:nsout-1 + rstart = i*sbox # starting position for each box + istart = floor(Int, rstart) + rstop = rstart + sbox # ending position for each box + istop = Int(clamp(floor(rstop), 0, ns1)) + frac1 = rstart-istart + frac2 = 1.0 - (rstop-istop) + + # add pixel values from istart to istop and subtract fractional pixel from istart to start and + # fractional pixel from rstop to istop + + if istart == istop + result[:,i+1] .= (1 .- frac1 .- frac2).*temp[:,istart+1] + else + result[:,i+1] .= dropdims(sum(temp[:,istart+1:istop+1], dims=2), dims=2) .- frac1.*temp[:,istart+1] .- frac2.*temp[:,istop+1] + end + end + return transpose(result) .* (total ? 1.0 : 1 ./ (sbox.*lbox)) + +end + +""" + frebin(image, nsout, nlout=1, total=false) -> rebinned_image + +### Purpose ### + +Shrink or expand the size of an array an arbitrary amount using interpolation + +### Arguments ### + +* `image`: the array representing the 1D or 2D image to be rebinned. +* `nsout`: number of samples in the output image, numeric scalar. +* `nlout` (optional): number of lines in the output image, numeric scalar (default = 1). +* `total` (optional boolean keyword): if true, the output pixels will be the + sum of pixels within the appropriate box of the input image. + Otherwise they will be the average. Use of the `total` keyword + conserves total counts. + +### Output ### + +The resized image is returned as the function result. + +### Example ### + +Suppose one has an 800 x 800 image array, im, that must be expanded to +a size 850 x 900 while conserving the total counts. The pixel values are +the sum of the x and y coordinates in the original image: + +```jldoctest +julia> using AstroLib + +julia> image = [x+y for x in 1:800, y in 1:800]; + +julia> image1 = frebin(image, 850, 900, total=true); + +julia> sum(image) +512640000 + +julia> sum(image1) +5.126400000000241e8 +``` + +### Notes ### + +If the input image sizes are a multiple of the output image sizes +then `frebin` is equivalent to summing over the grid multiples for +compression, and simple pixel duplication on expansion. + +If the number of output pixels are not integers, the output image +size will be truncated to an integer. The platescale, however, will +reflect the non-integer number of pixels. For example, if you want to +bin a 100 x 100 integer image such that each output pixel is 3.1 +input pixels in each direction use: + n = 100/3.1 # 32.2581 + image_out = frebin(image,n,n) + +The output image will be 32 x 32 and a small portion at the trailing +edges of the input image will be ignored. + +### History ### + +Adapted from May 1998 STIS version, written D. Lindler, ACC +Added /NOZERO, use INTERPOLATE instead of CONGRID, June 98 W. Landsman +Fixed for nsout non-integral but a multiple of image size Aug 98 D.Lindler +DJL, Oct 20, 1998, Modified to work for floating point image sizes when + expanding the image. +Improve speed by addressing arrays in memory order W.Landsman Dec/Jan 2001 + +Code of this function is based on IDL Astronomy User's Library. +""" +function frebin(image::AbstractArray{R}, nsout::Real, nlout::Real=1; total::Bool=false) where {R<:Real} + + # check that the image dimensions are either 1D or 2D + nd = ndims(image) + @assert nd in (1,2) "The input image must be either 1D or 2D!" + + _frebin(float(image), promote(floor(Int, nsout), floor(Int, nlout))..., total) +end diff --git a/src/fshift.jl b/src/fshift.jl new file mode 100644 index 0000000..e33c54c --- /dev/null +++ b/src/fshift.jl @@ -0,0 +1,89 @@ +# This file is a part of AstroLib.jl. License is MIT "Expat". +# Copyright (C) 2016 Mosè Giordano. + +function _fshift(image::AbstractArray{T}, Δx::T, Δy::T) where {T<:AbstractFloat} + + # Separate shift into an integer and fractional shift + intx = floor(Int, Δx) + inty = floor(Int, Δy) + fracx = Δx - intx + fracy = Δy - inty + if fracx < 0 + fracx += 1 + intx -= 1 + end + if fracy < 0 + fracy += 1 + inty -= 1 + end + + # Shift by the integer portion + s = circshift(image, (intx, inty)) + if iszero(fracx) && iszero(fracy) + return s + end + + # Use bilinear interpolation between four pixels + return s .* ((1 .- fracx) .* (1 .- fracy)) .+ + circshift(s, (0,1)) .* ((1 .- fracx) .* fracy) .+ + circshift(s, (1,0)) .* (fracx .* (1 .- fracy)) .+ + circshift(s, (1,1)) .* fracx .* fracy +end + +""" + fshift(image, Δx, Δy) -> shifted_image + +### Purpose ### + +Routine to shift an image by non-integer values + +### Arguments ### + +* `image`: 2D image to be shifted. +* `Δx`: shift in x +* `Δy`: shift in y + +### Output ### + +Shifted image is returned as the function results + +### Example ### + +Suppose we want to shift a 10x10 image by 0.5 pixels in the +x and y directions. The pixel values are the sum of the x and +y coordinates: + +```jldoctest +julia> using AstroLib + +julia> image = [x+y for x in 1:10, y in 1:10]; + +julia> fshift(image, 0.5, 0.5) +10×10 Matrix{Float64}: + 11.0 7.0 8.0 9.0 10.0 11.0 12.0 13.0 14.0 15.0 + 7.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 + 8.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0 + 9.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0 13.0 + 10.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0 13.0 14.0 + 11.0 7.0 8.0 9.0 10.0 11.0 12.0 13.0 14.0 15.0 + 12.0 8.0 9.0 10.0 11.0 12.0 13.0 14.0 15.0 16.0 + 13.0 9.0 10.0 11.0 12.0 13.0 14.0 15.0 16.0 17.0 + 14.0 10.0 11.0 12.0 13.0 14.0 15.0 16.0 17.0 18.0 + 15.0 11.0 12.0 13.0 14.0 15.0 16.0 17.0 18.0 19.0 +``` + +### History ### + +version 2 D. Lindler May, 1992 - rewritten for IDL version 2 +19-may-1992 JKF/ACC - move to GHRS DAF. + +Code of this function is based on IDL Astronomy User's Library. + +""" +function fshift(image::AbstractArray{R}, Δx::Real, Δy::Real) where {R<:Real} + # check that the image dimensions are 2D + nd = ndims(image) + @assert nd == 2 "The input image must be 2D!" + + _fshift(float(image), promote(float(Δx), float(Δy))...) +end \ No newline at end of file diff --git a/src/utils.jl b/src/utils.jl index 01f71a9..ccd645e 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -55,6 +55,12 @@ export euler include("flux2mag.jl") export flux2mag +include("frebin.jl") +export frebin + +include("fshift.jl") +export fshift + include("gal_uvw.jl") export gal_uvw