JuliaAstro · abhro · Jun 1, 2025 · Oct 2, 2024 · Oct 2, 2024 · Oct 2, 2024
diff --git a/Project.toml b/Project.toml
@@ -6,6 +6,7 @@ version = "0.4.3"
 [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"
@@ -16,4 +17,5 @@ DelimitedFiles = "1"
 LinearAlgebra = "1"
 Printf = "1"
 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.10"
diff --git 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
@@ -0,0 +1,248 @@
+# 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.
+
+### Examples ###
+
+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> size(image)
+(800, 800)
+
+julia> sum(image)
+512640000
+
+julia> image1 = frebin(image, 850, 900, total=true);
+
+julia> size(image1)
+(850, 900)
+
+julia> sum(image1) ≈ sum(image)
+true
+```
+
+If the new dimensions are integer ratios of the original
+dimensions, some optimizations can be performed.  We can
+also choose to preserve the *average* pixel value rather
+than the *total* pixel value in the image:
+
+```jldoctest
+julia> using AstroLib
+
+julia> image = [x+y for x in 1:800, y in 1:800];
+
+julia> size(image)
+(800, 800)
+
+julia> sum(image)/length(image)
+801.0
+
+julia> image1 = frebin(image, 400, 400, total=false);
+
+julia> size(image1)
+(400, 400)
+
+julia> (sum(image1)/length(image1)) ≈ (sum(image)/length(image))
+true
+
+julia> image2 = frebin(image, 1600, 1600, total=false);
+
+julia> size(image2)
+(1600, 1600)
+
+julia> (sum(image2)/length(image2)) ≈ (sum(image)/length(image))
+true
+```
+
+Of course, this works with 1D arrays as well:
+
+```jldoctest
+julia> using AstroLib
+
+julia> image = [x for x in 1:800];
+
+julia> size(image)
+(800,)
+
+julia> sum(image)
+320400
+
+julia> sum(image)/length(image)
+400.5
+
+julia> image1 = frebin(image, 1600);
+
+julia> size(image1)
+(1600,)
+
+julia> (sum(image1)/length(image1)) ≈ (sum(image)/length(image))
+true
+
+julia> image2 = frebin(image, 410, total=true);
+
+julia> size(image2)
+(410,)
+
+julia> sum(image2) ≈ sum(image)
+true
+```
+
+### 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
@@ -0,0 +1,97 @@
+# 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
+    # note: the idl version has branches checking if fracx and fracy are negative.
+    # these are not necessary here because we're using the floor function rather than
+    # a truncation of a float to an integer, so intx/y will always be <= Δx/y.
+
+    # 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
+
+julia> fshift(image, -3, -1)
+10×10 Matrix{Float64}:
+  6.0   7.0   8.0   9.0  10.0  11.0  12.0  13.0  14.0   5.0
+  7.0   8.0   9.0  10.0  11.0  12.0  13.0  14.0  15.0   6.0
+  8.0   9.0  10.0  11.0  12.0  13.0  14.0  15.0  16.0   7.0
+  9.0  10.0  11.0  12.0  13.0  14.0  15.0  16.0  17.0   8.0
+ 10.0  11.0  12.0  13.0  14.0  15.0  16.0  17.0  18.0   9.0
+ 11.0  12.0  13.0  14.0  15.0  16.0  17.0  18.0  19.0  10.0
+ 12.0  13.0  14.0  15.0  16.0  17.0  18.0  19.0  20.0  11.0
+  3.0   4.0   5.0   6.0   7.0   8.0   9.0  10.0  11.0   2.0
+  4.0   5.0   6.0   7.0   8.0   9.0  10.0  11.0  12.0   3.0
+  5.0   6.0   7.0   8.0   9.0  10.0  11.0  12.0  13.0   4.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
diff --git 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