From 9b426d5531f498a575946b564650a68a824592f4 Mon Sep 17 00:00:00 2001 From: Michael-Reefe Date: Wed, 2 Oct 2024 00:46:40 -0400 Subject: [PATCH 1/5] added frebin and fshift functions --- Project.toml | 3 +- src/AstroLib.jl | 1 + src/frebin.jl | 168 ++++++++++++++++++++++++++++++++++++++++++++++++ src/fshift.jl | 83 ++++++++++++++++++++++++ src/utils.jl | 6 ++ 5 files changed, 260 insertions(+), 1 deletion(-) create mode 100644 src/frebin.jl create mode 100644 src/fshift.jl diff --git a/Project.toml b/Project.toml index 9eff46b..d0fc4c5 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ 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" @@ -15,9 +16,9 @@ StaticArrays = "0.5, 0.6, 0.7, 0.8, 0.9, 0.10, 0.11, 0.12, 1.0" 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..2b7d6e0 --- /dev/null +++ b/src/frebin.jl @@ -0,0 +1,168 @@ +# 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(Linear())), Flat()) + return interpfunc(xindex) * (total ? sbox : 1.) + end + yindex = (1:nlout) / (nlout/nl) + interpfunc = extrapolate(interpolate(image, BSpline(Linear())), 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 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> im = [x+y for x in 1:800, y in 1:800] +julia> im1 = frebin(im, 850, 900, total=true); +julia> sum(im) +512640000 +julia> sum(im1) +5.126400000000241e8 +``` + +### Notes ### + +If the input image sizes are a multiple of the output image sizes +then `frebin` is equivalent to the julia Interpolations.jl functions 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} + _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..23bb9f2 --- /dev/null +++ b/src/fshift.jl @@ -0,0 +1,83 @@ +# 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> im = [x+y for x in 1:10, y in 1:10]; +julia> fshift(im, 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} + _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 From 83f5faa4b898cea62016f055edab81f0fae5c3d4 Mon Sep 17 00:00:00 2001 From: Michael-Reefe Date: Wed, 2 Oct 2024 00:57:20 -0400 Subject: [PATCH 2/5] changed to constant interpolation to agree with the IDL version --- src/frebin.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/frebin.jl b/src/frebin.jl index 2b7d6e0..f408a9c 100644 --- a/src/frebin.jl +++ b/src/frebin.jl @@ -21,11 +21,11 @@ function _frebin(image::AbstractArray{T}, nsout::S, nlout::S, total::Bool) where if (nsout % ns == 0) && (nlout % nl == 0) xindex = (1:nsout) / (nsout/ns) if isone(nl) # 1D case, linear interpolation - interpfunc = extrapolate(interpolate(image, BSpline(Linear())), Flat()) + interpfunc = extrapolate(interpolate(image, BSpline(Constant())), Flat()) return interpfunc(xindex) * (total ? sbox : 1.) end yindex = (1:nlout) / (nlout/nl) - interpfunc = extrapolate(interpolate(image, BSpline(Linear())), Flat()) + interpfunc = extrapolate(interpolate(image, BSpline(Constant())), Flat()) return [interpfunc(x, y) for x in xindex, y in yindex] .* (total ? sbox.*lbox : 1.) end @@ -138,7 +138,7 @@ julia> sum(im1) ### Notes ### If the input image sizes are a multiple of the output image sizes -then `frebin` is equivalent to the julia Interpolations.jl functions for +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 From b72219deb2e1b7778da7caf4d62ceb08a768b761 Mon Sep 17 00:00:00 2001 From: Michael-Reefe Date: Wed, 2 Oct 2024 01:12:10 -0400 Subject: [PATCH 3/5] fixed jldoctests --- src/frebin.jl | 8 ++++---- src/fshift.jl | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/frebin.jl b/src/frebin.jl index f408a9c..df866d3 100644 --- a/src/frebin.jl +++ b/src/frebin.jl @@ -127,11 +127,11 @@ the sum of the x and y coordinates in the original image: ```jldoctest julia> using AstroLib -julia> im = [x+y for x in 1:800, y in 1:800] -julia> im1 = frebin(im, 850, 900, total=true); -julia> sum(im) +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(im1) +julia> sum(image1) 5.126400000000241e8 ``` diff --git a/src/fshift.jl b/src/fshift.jl index 23bb9f2..c7b429c 100644 --- a/src/fshift.jl +++ b/src/fshift.jl @@ -55,8 +55,8 @@ y coordinates: ```jldoctest julia> using AstroLib -julia> im = [x+y for x in 1:10, y in 1:10]; -julia> fshift(im, 0.5, 0.5) +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 From d9cb028fc351bf3436892547061ec6680130f640 Mon Sep 17 00:00:00 2001 From: Michael-Reefe Date: Wed, 2 Oct 2024 01:18:12 -0400 Subject: [PATCH 4/5] fixed jldoctests? --- src/frebin.jl | 5 ++++- src/fshift.jl | 2 ++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/src/frebin.jl b/src/frebin.jl index df866d3..e7561de 100644 --- a/src/frebin.jl +++ b/src/frebin.jl @@ -127,10 +127,13 @@ 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> 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 ``` diff --git a/src/fshift.jl b/src/fshift.jl index c7b429c..3dbecf0 100644 --- a/src/fshift.jl +++ b/src/fshift.jl @@ -55,7 +55,9 @@ 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 From fce96c7c1eeb4407bb32f6fd98cdd7468277402d Mon Sep 17 00:00:00 2001 From: Michael-Reefe Date: Thu, 3 Oct 2024 15:10:50 -0400 Subject: [PATCH 5/5] added compat req and dimensionality checks --- Project.toml | 1 + src/frebin.jl | 7 ++++++- src/fshift.jl | 4 ++++ 3 files changed, 11 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index d0fc4c5..ba4e474 100644 --- a/Project.toml +++ b/Project.toml @@ -13,6 +13,7 @@ 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] diff --git a/src/frebin.jl b/src/frebin.jl index e7561de..6f9a732 100644 --- a/src/frebin.jl +++ b/src/frebin.jl @@ -106,7 +106,7 @@ Shrink or expand the size of an array an arbitrary amount using interpolation ### Arguments ### -* `image`: the array representing the image to be rebinned. +* `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 @@ -167,5 +167,10 @@ 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 index 3dbecf0..e33c54c 100644 --- a/src/fshift.jl +++ b/src/fshift.jl @@ -81,5 +81,9 @@ 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