mweastwood / LibHealpix.jl

A Julia wrapper of the Healpix library.
http://mweastwood.info/LibHealpix.jl/stable/
Other
11 stars 7 forks source link

gaussian blur #32

Open mweastwood opened 6 years ago

mweastwood commented 6 years ago

Could be a useful function to have. It's common to want to degrade maps to a lower resolution. Below is a prototype implementation (would definitely want to tidy this up first).

function smooth(map, width, output_nside=nside(map))
    σ = width/(2sqrt(2log(2)))
    kernel = HealpixMap(Float64, output_nside)
    for pix = 1:length(kernel)
        θ, ϕ = LibHealpix.pix2ang_ring(nside(kernel), pix)
        θ = rad2deg(θ)
        kernel[pix] = exp(-θ^2/(2σ^2))
    end
    dΩ = 4π/length(kernel)
    kernel = HealpixMap(kernel.pixels / (sum(kernel.pixels)*dΩ))

    lmax = mmax = 1000
    map_alm = map2alm(map, lmax, mmax, iterations=10)
    kernel_alm = map2alm(kernel, lmax, mmax, iterations=10)
    output_alm = Alm(Complex128, lmax, mmax)
    for m = 0:mmax, l = m:lmax
        output_alm[l, m] = sqrt((4π)/(2l+1))*map_alm[l, m]*kernel_alm[l, 0]
    end

    alm2map(output_alm, output_nside)
end