Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

frebin and fshift functions #87

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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"]
1 change: 1 addition & 0 deletions src/AstroLib.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
176 changes: 176 additions & 0 deletions src/frebin.jl
Original file line number Diff line number Diff line change
@@ -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
89 changes: 89 additions & 0 deletions src/fshift.jl
Original file line number Diff line number Diff line change
@@ -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
6 changes: 6 additions & 0 deletions src/utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
Loading