Skip to content

Commit

Permalink
add Fejer1 types and relax sky2pix routines to apply to both
Browse files Browse the repository at this point in the history
  • Loading branch information
@xzackli
xzackli committed Dec 16, 2023
1 parent 8e989e7 commit 691eddc
Show file tree
Hide file tree
Showing 3 changed files with 59 additions and 37 deletions.
4 changes: 2 additions & 2 deletions src/Pixell.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ import Unitful, UnitfulAngles
import Unitful: uconvert, ustrip
using StaticArrays
using DSP: unwrap, unwrap!
import FastTransforms: chebyshevjacobimoments1, clenshawcurtisweights
import FastTransforms: chebyshevjacobimoments1, clenshawcurtisweights, fejerweights1
using Libsharp
import Libsharp: sharp_execute!
import Healpix: Alm, map2alm, alm2map, alm2cl
Expand All @@ -32,7 +32,7 @@ include("plot.jl")
include("transform_distance.jl")
include("utils.jl")

export Enmap, CarClenshawCurtis, Gnomonic, getwcs
export Enmap, CarClenshawCurtis, CarFejer1, Gnomonic, getwcs
export geometry, fullsky_geometry, slice_geometry, pad, SkyBoundingBox
export pix2sky, pix2sky!, sky2pix, sky2pix!, skyarea, pixareamap, pixareamap!, posmap
export read_map, write_map
Expand Down
2 changes: 1 addition & 1 deletion src/enmap.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -216,7 +216,7 @@ function read_map(path::String; hdu::Int=1, sel=(), wcs::Union{WCSTransform,Noth
wcs0 = WCS.from_header(header_str)[1]
@assert wcs0.ctype[1] == "RA---CAR"
@assert wcs0.ctype[2] == "DEC--CAR"
wcs = convert(CarClenshawCurtis{Float64}, wcs0)
wcs = convert(CarClenshawCurtis{Float64}, wcs0) # FIXME select between CC or Fejer1
end
end
Enmap(data, wcs)
Expand Down
90 changes: 56 additions & 34 deletions src/projections/car_proj.jl
View file
Original file line number Diff line number Diff line change
@@ -1,26 +1,43 @@

abstract type AbstractCARWCS{T} <: AbstractWCSTransform end

"""
Fast custom WCS structure.
"""
struct CarClenshawCurtis{T} <: AbstractWCSTransform
struct CarClenshawCurtis{T} <: AbstractCARWCS{T}
cdelt::Tuple{T,T}
crpix::Tuple{T,T}
crval::Tuple{T,T}
unit::T # conversion factor to radians
end

struct CarFejer1{T} <: AbstractCARWCS{T}
cdelt::Tuple{T,T}
crpix::Tuple{T,T}
crval::Tuple{T,T}
unit::T # conversion factor to radians
end
getunit(T::Type{<:Real}, wcs::CarClenshawCurtis) = T(wcs.unit)
getcdelt(wcs::CarClenshawCurtis) = wcs.cdelt
getcrpix(wcs::CarClenshawCurtis) = wcs.crpix
getcrval(wcs::CarClenshawCurtis) = wcs.crval

Base.copy(w::W) where {W<:CarClenshawCurtis} = w # CarClenshawCurtis is fully immutable
getunit(T::Type{<:Real}, wcs::AbstractCARWCS) = T(wcs.unit)
getcdelt(wcs::AbstractCARWCS) = wcs.cdelt
getcrpix(wcs::AbstractCARWCS) = wcs.crpix
getcrval(wcs::AbstractCARWCS) = wcs.crval

Base.copy(w::W) where {W<:AbstractCARWCS} = w # both CAR WCS types are fully immutable

function Base.convert(::Type{CarClenshawCurtis{T}}, w0::WCSTransform) where T
# FIXME add some asserts
return CarClenshawCurtis{T}(
T.(getcdelt(w0)), T.(getcrpix(w0)), T.(getcrval(w0)), getunit(T, w0))
end

function Base.convert(::Type{WCSTransform}, w0::CarClenshawCurtis{T}) where T
function Base.convert(::Type{CarFejer1{T}}, w0::WCSTransform) where T
# FIXME add some asserts
return CarFejer1{T}(
T.(getcdelt(w0)), T.(getcrpix(w0)), T.(getcrval(w0)), getunit(T, w0))
end

function Base.convert(::Type{WCSTransform}, w0::AbstractCARWCS)
return WCSTransform(2;
ctype = ["RA---CAR", "DEC--CAR"],
cdelt = collect(getcdelt(w0)),
Expand All @@ -29,7 +46,7 @@ function Base.convert(::Type{WCSTransform}, w0::CarClenshawCurtis{T}) where T
end

# this kind of WCS only has two spatial dimensions. this check should be constant-propagated
function Base.getproperty(wcs::CarClenshawCurtis, k::Symbol)
function Base.getproperty(wcs::AbstractCARWCS, k::Symbol)
if k == :naxis
return 2
end
Expand All @@ -41,6 +58,11 @@ function Base.show(io::IO, wcs::CarClenshawCurtis{T}) where T
for k in ["naxis","cdelt","crval","crpix"]], ",")
print(io, "CarClenshawCurtis{$(T)}($expr)")
end
function Base.show(io::IO, wcs::CarFejer1{T}) where T
expr = join(["$k=$(getproperty(wcs, Symbol(k)))"
for k in ["naxis","cdelt","crval","crpix"]], ",")
print(io, "CarFejer1{$(T)}($expr)")
end


"""
Expand All @@ -67,7 +89,7 @@ julia> shape, wcs = fullsky_geometry(deg2rad(1))
julia> pix2sky!(shape, wcs, pixcoords, skycoords)
```
"""
function pix2sky!(shape, wcs::CarClenshawCurtis, pixcoords::AbstractArray{TP,2},
function pix2sky!(shape, wcs::AbstractCARWCS, pixcoords::AbstractArray{TP,2},
skycoords::AbstractArray{TS,2}; safe=true) where {TP,TS}

angle_unit = getunit(wcs)
Expand All @@ -93,19 +115,19 @@ function pix2sky!(shape, wcs::CarClenshawCurtis, pixcoords::AbstractArray{TP,2},
end

# the not-in-place version just creates an output array and calls the in-place one above
function pix2sky(shape, wcs::CarClenshawCurtis,
function pix2sky(shape, wcs::AbstractCARWCS,
pixcoords::AbstractArray{TP,2}; safe=true) where TP
skycoords = similar(pixcoords)
return pix2sky!(shape, wcs, pixcoords, skycoords; safe=safe)
end

"""
pix2sky(m::Enmap{T,N,AA,<:CarClenshawCurtis}, ra_pixel, dec_pixel)
pix2sky(m::Enmap{T,N,AA,<:AbstractCARWCS}, ra_pixel, dec_pixel)
Compute the sky position of a single position on the sky.
Only implemented for CAR (Clenshaw-Curtis variant) projections, so
the input map is of type `Enmap{T,N,AA,<:CarClenshawCurtis}`.
Only implemented for CAR projections, so
the input map is of type `Enmap{T,N,AA,<:AbstractCARWCS}`.
This takes pixel indices for RA and DEC, and returns a tuple containing
the corresponding RA and DEC.
Expand All @@ -116,7 +138,7 @@ julia> pix2sky(shape, wcs, 30.0, 80.0)
(151.0, -11.0)
```
"""
function pix2sky(shape, wcs::CarClenshawCurtis, ra_pixel, dec_pixel; safe=true)
function pix2sky(shape, wcs::AbstractCARWCS, ra_pixel, dec_pixel; safe=true)
angle_unit = getunit(wcs)
α₀, δ₀ = getcrval(wcs) .* angle_unit
Δα, Δδ = getcdelt(wcs) .* angle_unit
Expand All @@ -130,7 +152,7 @@ function pix2sky(shape, wcs::CarClenshawCurtis, ra_pixel, dec_pixel; safe=true)
end

# when passing a length-2 vector [ra, dec], return a vector. wraps the pix2sky(shape, wcs, ra_pix, dec_pix)
function pix2sky(shape, wcs::CarClenshawCurtis, pixcoords::AbstractVector; safe=true)
function pix2sky(shape, wcs::AbstractCARWCS, pixcoords::AbstractVector; safe=true)
@assert length(pixcoords) == 2
skycoords = collect(pix2sky(shape, wcs, first(pixcoords), last(pixcoords)))
if safe
Expand All @@ -140,7 +162,7 @@ function pix2sky(shape, wcs::CarClenshawCurtis, pixcoords::AbstractVector; safe=
end


function sky2pix!(shape, wcs::CarClenshawCurtis, skycoords::AbstractArray{TS,2},
function sky2pix!(shape, wcs::AbstractCARWCS, skycoords::AbstractArray{TS,2},
pixcoords::AbstractArray{TP,2}; safe=true) where {TS,TP}

# retrieve WCS info
Expand Down Expand Up @@ -171,15 +193,15 @@ function sky2pix!(shape, wcs::CarClenshawCurtis, skycoords::AbstractArray{TS,2},
end

# the not-in-place version just creates an output array and calls the in-place one above
function sky2pix(shape, wcs::CarClenshawCurtis,
function sky2pix(shape, wcs::AbstractCARWCS,
skycoords::AbstractArray{TS,2}; safe=true) where {TS}
pixcoords = similar(skycoords)
return sky2pix!(shape, wcs, skycoords, pixcoords; safe=safe)
end


"""
sky2pix(m::Enmap{T,N,AA,<:CarClenshawCurtis}, ra, dec)
sky2pix(m::Enmap{T,N,AA,<:AbstractCARWCS}, ra, dec)
Compute 1-indexed pixels into sky coordinates.
Expand All @@ -195,7 +217,7 @@ julia> sky2pix(shape, wcs, deg2rad(30.0), deg2rad(80.0))
(151.0, 171.0)
```
"""
function sky2pix(shape, wcs::CarClenshawCurtis, ra::Number, dec::Number; safe=true)
function sky2pix(shape, wcs::AbstractCARWCS, ra::Number, dec::Number; safe=true)
angle_unit = getunit(wcs)
α₀, δ₀ = getcrval(wcs) .* angle_unit
Δα, Δδ = getcdelt(wcs) .* angle_unit
Expand All @@ -210,7 +232,7 @@ function sky2pix(shape, wcs::CarClenshawCurtis, ra::Number, dec::Number; safe=tr
end
return pix_ra, pix_dec
end
function sky2pix(shape, wcs::CarClenshawCurtis,
function sky2pix(shape, wcs::AbstractCARWCS,
ra::AV, dec::AV; safe=true) where {AV<:AbstractVector}
angle_unit = getunit(wcs)
α₀, δ₀ = getcrval(wcs) .* angle_unit
Expand All @@ -230,28 +252,28 @@ function sky2pix(shape, wcs::CarClenshawCurtis,
end

# when passing a vector [ra, dec], return a vector. wraps the sky2pix(shape, wcs, ra, dec).
function sky2pix(shape, wcs::CarClenshawCurtis, skycoords::AbstractVector; safe=true)
function sky2pix(shape, wcs::AbstractCARWCS, skycoords::AbstractVector; safe=true)
@assert length(skycoords) == 2
return collect(sky2pix(shape, wcs::CarClenshawCurtis,
return collect(sky2pix(shape, wcs::AbstractCARWCS,
first(skycoords), last(skycoords); safe=safe))
end

skyarea(shape, wcs::CarClenshawCurtis) = skyarea_cyl(shape, wcs)
iscyl(wcs::CarClenshawCurtis) = true
skyarea(shape, wcs::AbstractCARWCS) = skyarea_cyl(shape, wcs)
iscyl(wcs::AbstractCARWCS) = true

"""Generate a similar Enmap whose pixel values are the areas of the pixels in steradians."""
function pixareamap(m::Enmap{T,N,AA,W}) where {T,N,AA,W<:CarClenshawCurtis}
function pixareamap(m::Enmap{T,N,AA,W}) where {T,N,AA,W<:AbstractCARWCS}
pixareas = similar(m)
pixareamap!(pixareas)
end

function pixareamap(shape, wcs::CarClenshawCurtis)
function pixareamap(shape, wcs::AbstractCARWCS)
pixareas = Enmap(Array{Float64}(undef, shape), wcs)
pixareamap!(pixareas)
end

function sliced_wcs(wcs::CarClenshawCurtis{T}, cdelt′, crpix′) where T
new_wcs = CarClenshawCurtis{T}(cdelt′, crpix′, wcs.crval, wcs.unit)
function sliced_wcs(wcs::W, cdelt′, crpix′) where {T, W<:AbstractCARWCS{T}}
new_wcs = W(cdelt′, crpix′, wcs.crval, wcs.unit)
return new_wcs
end

Expand All @@ -263,7 +285,7 @@ function pad(m::Enmap, npix_ra::Int, npix_dec::Int; mode=:center)
end
end

function center_pad(m::Enmap{T,N,AA,W}, npix_ra::Int, npix_dec::Int) where {T,N,AA,W<:CarClenshawCurtis}
function center_pad(m::Enmap{T,N,AA,W}, npix_ra::Int, npix_dec::Int) where {T,N,AA,W<:AbstractCARWCS}
new_shape, new_wcs = center_pad(size(m), m.wcs, npix_ra, npix_dec)
arr = AA(undef, new_shape)
fill!(arr, zero(T))
Expand All @@ -274,9 +296,9 @@ function center_pad(m::Enmap{T,N,AA,W}, npix_ra::Int, npix_dec::Int) where {T,N,
return Enmap(arr, new_wcs)
end

function center_pad(shape, wcs::CarClenshawCurtis{T}, npix_ra::Int, npix_dec) where T
function center_pad(shape, wcs::W, npix_ra::Int, npix_dec) where {T,W<:AbstractCARWCS{T}}
new_shape = (shape[1] + 2npix_ra, shape[2] + 2npix_dec, shape[3:end]...)
new_wcs = CarClenshawCurtis{T}(
new_wcs = W(
wcs.cdelt,
wcs.crpix .+ (npix_ra, npix_dec),
wcs.crval,
Expand All @@ -285,7 +307,7 @@ function center_pad(shape, wcs::CarClenshawCurtis{T}, npix_ra::Int, npix_dec) wh
end


function corner_pad(m::Enmap{T,N,AA,W}, npix_ra::Int, npix_dec::Int) where {T,N,AA,W<:CarClenshawCurtis}
function corner_pad(m::Enmap{T,N,AA,W}, npix_ra::Int, npix_dec::Int) where {T,N,AA,W<:AbstractCARWCS}
new_shape, new_wcs = corner_pad(size(m), m.wcs, npix_ra, npix_dec)
arr = AA(undef, new_shape)
fill!(arr, zero(T))
Expand All @@ -296,11 +318,11 @@ function corner_pad(m::Enmap{T,N,AA,W}, npix_ra::Int, npix_dec::Int) where {T,N,
return Enmap(arr, new_wcs)
end

function corner_pad(shape, wcs::CarClenshawCurtis{T}, npix_ra::Int, npix_dec) where T
function corner_pad(shape, wcs::AbstractCARWCS{T}, npix_ra::Int, npix_dec) where T
new_shape = (shape[1] + npix_ra, shape[2] + npix_dec, shape[3:end]...)
new_wcs = deepcopy(wcs)
return new_shape, new_wcs
end

pad(m::Enmap{T,N,AA,W}, npix; kwargs...) where {T,N,AA,W<:CarClenshawCurtis} = pad(m, npix, npix; kwargs...)
pad(m::Enmap{T,N,AA,W}, npix; kwargs...) where {T,N,AA,W<:AbstractCARWCS} = pad(m, npix, npix; kwargs...)

0 comments on commit 691eddc

Please sign in to comment.