Merge pull request #35 from mweastwood/query_disc

implement query_disc

deps/build.jl

@@ -29,7 +29,7 @@ provides(AptGet, Dict("libcfitsio3-dev"    => libcfitsio,
 
 if is_apple()
     using Homebrew
-    provides(Homebrew.HB, "homebrew/science/cfitsio", libcfitsio, os=:Darwin)
+    provides(Homebrew.HB, "cfitsio", libcfitsio, os=:Darwin)
     provides(Homebrew.HB, "homebrew/science/healpix", [libchealpix, libhealpix_cxx], os=:Darwin)
     if !has_pkg_config
         Homebrew.add("pkg-config")

deps/src/wrapper/wrapper.cpp

@@ -18,12 +18,14 @@
 using namespace std;
 
 #include <complex>
-#include <xcomplex.h>
-#include <arr.h>
+#include "xcomplex.h"
 
-#include <healpix_map.h>
-#include <alm.h>
-#include <alm_healpix_tools.h>
+#include "arr.h"
+#include "rangeset.h"
+
+#include "healpix_map.h"
+#include "alm.h"
+#include "alm_healpix_tools.h"
 
 // Note:
 //
@@ -43,6 +45,13 @@ Healpix_Map<T> construct_healpix_map(int nside, int order, T* pixels)
     return map;
 }
 
+T_Healpix_Base<int> construct_healpix_base(int nside, int order)
+{
+    Healpix_Ordering_Scheme ordering_scheme = static_cast<Healpix_Ordering_Scheme>(order);
+    auto map = T_Healpix_Base<int>(nside, ordering_scheme, SET_NSIDE);
+    return map;
+}
+
 template <typename T>
 Alm<xcomplex<T> > construct_alm(int lmax, int mmax, complex<T>* coefficients)
 {
@@ -124,5 +133,25 @@ extern "C" {
     {
         return interpolate(nside, order, pixels, theta, phi);
     }
+
+    int* query_disc(int nside, int order,
+                    double theta, double phi, double radius,
+                    bool inclusive, int* output_length)
+    {
+        auto map = construct_healpix_base(nside, order);
+        auto ptg = pointing(theta, phi);
+        auto set = rangeset<int>();
+        if (inclusive)
+            map.query_disc_inclusive(ptg, radius, set);
+        else
+            map.query_disc(ptg, radius, set);
+        // unfortunately we need an extra copy to get this into a standard array that Julia can
+        // take ownership of
+        auto vec = set.toVector();
+        auto arr = new int[vec.size()];
+        copy(vec.begin(), vec.end(), arr);
+        *output_length = vec.size();
+        return arr;
+    }
 }
 

docs/src/library.md

@@ -39,6 +39,8 @@ ang2pix
 pix2ang
 vec2pix
 pix2vec
+LibHealpix.interpolate
+query_disc
 ```
 
 ## Spherical Harmonic Coefficients

src/LibHealpix.jl

@@ -30,6 +30,7 @@ export vec2pix_nest, vec2pix_ring, pix2vec_nest, pix2vec_ring
 export HealpixMap, RingHealpixMap, NestHealpixMap
 export ang2pix, pix2ang, vec2pix, pix2vec
 export isring, isnest
+export query_disc
 
 # alm.jl
 export Alm, @lm, lm

src/io.jl

@@ -18,12 +18,12 @@
 
 Write the `HealpixMap` to disk as a FITS image.
 
-**Arguments**
+**Arguments:**
 
 - `filename` - the name of the output file (eg. `"/path/to/healpix.fits"`)
 - `map` - the Healpix map to write
 
-**Keyword Arguments**
+**Keyword Arguments:**
 
 - `coordsys` - the coordinate system of the map (one of `"G"` galactic, `"E"` ecliptic, or `"C"`
     celestial)
@@ -73,7 +73,7 @@ end
 
 Read a `HealpixMap` (stored as a FITS image) from disk.
 
-**Arguments**
+**Arguments:**
 
 - `filename` - the name of the input file (eg. `"/path/to/healpix.fits"`)
 

src/map.jl

@@ -195,8 +195,8 @@ $(θ, ϕ)$.
 **Arguments:**
 
 - `map` - the input Healpix map
-- `theta` - the inclination angle $θ$
-- `phi` - the azimuthal angle $ϕ$
+- `theta` - the inclination angle $θ$ (in radians)
+- `phi` - the azimuthal angle $ϕ$ (in radians)
 
 **Usage:**
 
@@ -354,15 +354,15 @@ function interpolate(map::HealpixMap{Float64}, θ, ϕ)
 end
 
 doc"""
-    interpolate(map, theta, phi)
+    LibHealpix.interpolate(map, theta, phi)
 
 Linearly interpolate the Healpix map at the given spherical coordinates $(θ, ϕ)$.
 
 **Arguments:**
 
 - `map` - the input Healpix map
-- `theta` - the inclination angle $θ$
-- `phi` - the azimuthal angle $ϕ$
+- `theta` - the inclination angle $θ$ (in radians)
+- `phi` - the azimuthal angle $ϕ$ (in radians)
 
 **Usage:**
 
@@ -379,3 +379,51 @@ julia> healpixmap = RingHealpixMap(Float64, 256)
 """
 interpolate
 
+doc"""
+    query_disc(nside, ordering, theta, phi, radius; inclusive=true)
+    query_disc(map, theta, phi, radius; inclusive=true)
+
+Return a list of all pixels contained within a circular disc of the given radius.
+
+**Arguments:**
+
+- `nside` - the Healpix resolution parameter
+- `ordering` - the ordering of the Healpix map (either `LibHealpix.ring` or `LibHealpix.nest`
+- `theta` - the inclination angle $θ$ (in radians)
+- `phi` - the azimuthal angle $ϕ$ (in radians)
+- `radius` - the radius of the disc (in radians)
+- `map` - the input Healpix map (`nside` and `ordering` will be inferred from the map)
+
+**Keyword Arguments:**
+
+- `inclusive` - if set to `true pixels partially contained within the disc will be included,
+    otherwise they are excluded
+
+**Usage:**
+
+```jldoctest
+julia> query_disc(512, LibHealpix.ring, 0, 0, deg2rad(10/60), inclusive=false)
+4-element Array{Int32,1}:
+ 1
+ 2
+ 3
+ 4
+
+julia> query_disc(512, LibHealpix.ring, 0, 0, deg2rad(10/60), inclusive=true) |> length
+24
+```
+"""
+function query_disc(nside, ordering, θ, ϕ, radius; inclusive=true)
+    len = Ref{Cint}()
+    ptr = ccall(("query_disc", libhealpixwrapper), Ptr{Cint},
+                (Cint, Cint, Cdouble, Cdouble, Cdouble, Bool, Ref{Cint}),
+                nside, ordering, θ, ϕ, radius, inclusive, len)
+    arr = unsafe_wrap(Array, ptr, len[], true)
+    arr .+= Cint(1) # convert to 1-based indexing
+    arr
+end
+
+function query_disc(map::HealpixMap, θ, ϕ, radius; inclusive=true)
+    query_disc(map.nside, ordering(map), θ, ϕ, radius, inclusive=inclusive)
+end
+

src/pixel.jl

@@ -199,61 +199,54 @@ function vec2ang(vec)
     θ, ϕ
 end
 
-types = (Clong == Clonglong)? (:Clong,) : (:Clong, :Clonglong)
-for T in types
-    # Append "64" to the ccall function name when defining for Clonglong
-    # (note we omit this suffix from the Julia function name -- just let dispatch take care of it)
-    funcname(f) = (T == :Clong)? string(f) : string(f)*"64"
-
-    for f in (:nest2ring, :ring2nest)
-        @eval function $f(nside::$T, ipix::$T)
-            ipix -= 1 # Subtract one to convert back to a 0-indexed scheme
-            ipixoutptr = Ref{$T}(0)
-            ccall(($(funcname(f)), libchealpix), Void, ($T, $T, Ref{$T}),
-                  nside, ipix, ipixoutptr)
-            ipixoutptr[] + 1 # Add one to convert to a 1-indexed scheme
-        end
+for f in (:nest2ring, :ring2nest)
+    @eval function $f(nside::Integer, ipix::Integer)
+        ipix -= 1 # Subtract one to convert back to a 0-indexed scheme
+        ipixoutptr = Ref{Clong}(0)
+        ccall(($(string(f)), libchealpix), Void, (Clong, Clong, Ref{Clong}),
+              nside, ipix, ipixoutptr)
+        ipixoutptr[] + 1 # Add one to convert to a 1-indexed scheme
     end
+end
 
-    for f in (:ang2pix_nest, :ang2pix_ring)
-        @eval function $f(nside::$T, θ::Real, ϕ::Real)
-            θ′, ϕ′ = verify_angles(θ, ϕ)
-            ipixptr = Ref{$T}(0)
-            ccall(($(funcname(f)), libchealpix), Void, ($T, Cdouble, Cdouble, Ref{$T}),
-                  nside, θ′, ϕ′, ipixptr)
-            ipixptr[] + 1 # Add one to convert to a 1-indexed scheme
-        end
+for f in (:ang2pix_nest, :ang2pix_ring)
+    @eval function $f(nside::Integer, θ::Real, ϕ::Real)
+        θ′, ϕ′ = verify_angles(θ, ϕ)
+        ipixptr = Ref{Clong}(0)
+        ccall(($(string(f)), libchealpix), Void, (Clong, Cdouble, Cdouble, Ref{Clong}),
+              nside, θ′, ϕ′, ipixptr)
+        ipixptr[] + 1 # Add one to convert to a 1-indexed scheme
     end
+end
 
-    for f in (:pix2ang_nest, :pix2ang_ring)
-        @eval function $f(nside::$T, ipix::$T)
-            ipix -= 1 # Subtract one to convert back to a 0-indexed scheme
-            θptr = Ref{Cdouble}(0.0)
-            ϕptr = Ref{Cdouble}(0.0)
-            ccall(($(funcname(f)), libchealpix), Void, ($T, $T, Ref{Cdouble}, Ref{Cdouble}),
-                  nside, ipix, θptr, ϕptr)
-            θptr[], ϕptr[]
-        end
+for f in (:pix2ang_nest, :pix2ang_ring)
+    @eval function $f(nside::Integer, ipix::Integer)
+        ipix -= 1 # Subtract one to convert back to a 0-indexed scheme
+        θptr = Ref{Cdouble}(0.0)
+        ϕptr = Ref{Cdouble}(0.0)
+        ccall(($(string(f)), libchealpix), Void, (Clong, Clong, Ref{Cdouble}, Ref{Cdouble}),
+              nside, ipix, θptr, ϕptr)
+        θptr[], ϕptr[]
     end
+end
 
-    for f in (:vec2pix_nest, :vec2pix_ring)
-        @eval function $f(nside::$T, vec::AbstractVector)
-            vec′ = UnitVector(vec)
-            ipixptr = Ref{$T}(0)
-            ccall(($(funcname(f)), libchealpix), Void, ($T, Ptr{Cdouble}, Ref{$T}),
-                  nside, vec′, ipixptr)
-            ipixptr[] + 1 # Add one to convert to a 1-indexed scheme
-        end
+for f in (:vec2pix_nest, :vec2pix_ring)
+    @eval function $f(nside::Integer, vec::AbstractVector)
+        vec′ = UnitVector(vec)
+        ipixptr = Ref{Clong}(0)
+        ccall(($(string(f)), libchealpix), Void, (Clong, Ptr{Cdouble}, Ref{Clong}),
+              nside, vec′, ipixptr)
+        ipixptr[] + 1 # Add one to convert to a 1-indexed scheme
     end
+end
 
-    for f in (:pix2vec_nest, :pix2vec_ring)
-        @eval function $f(nside::$T, ipix::$T)
-            ipix -= 1 # Subtract one to convert back to a 0-indexed scheme
-            vec = Ref{UnitVector}(UnitVector(0, 0, 1))
-            ccall(($(funcname(f)), libchealpix), Void, ($T, $T, Ptr{Cdouble}),
-                  nside, ipix, vec)
-            vec[]
-        end
+for f in (:pix2vec_nest, :pix2vec_ring)
+    @eval function $f(nside::Integer, ipix::Integer)
+        ipix -= 1 # Subtract one to convert back to a 0-indexed scheme
+        vec = Ref{UnitVector}(UnitVector(0, 0, 1))
+        ccall(($(string(f)), libchealpix), Void, (Clong, Clong, Ptr{Cdouble}),
+              nside, ipix, vec)
+        vec[]
     end
 end
 

test/map.jl

@@ -181,5 +181,32 @@
             @test LibHealpix.interpolate(map, 0, 0) === expected
         end
     end
+
+    @testset "query_disc" begin
+        function check_query(nside, ordering, θ, ϕ, radius, pixel)
+            vec = ang2vec(θ, ϕ)
+            if ordering == LibHealpix.ring
+                vec′ = pix2vec_ring(nside, pixel)
+            else
+                vec′ = pix2vec_nest(nside, pixel)
+            end
+            distance = acos(clamp(dot(vec, vec′), -1, 1))
+            distance < radius
+        end
+
+        nside = 256
+        for (ordering, map) in ((LibHealpix.ring, RingHealpixMap(Float64, nside)),
+                                (LibHealpix.nest, NestHealpixMap(Float64, nside)))
+            for (θ, ϕ) in ((0, 0), (π, 2π), (π/2, π), (0.1, 0.3), (2.0, 3.0))
+                radius = deg2rad(1)
+                inclusive_pixels = query_disc(nside, ordering, θ, ϕ, radius)
+                exclusive_pixels = query_disc(nside, ordering, θ, ϕ, radius, inclusive=false)
+                @test length(inclusive_pixels) > length(exclusive_pixels)
+                @test inclusive_pixels == query_disc(map, θ, ϕ, radius)
+                @test exclusive_pixels == query_disc(map, θ, ϕ, radius, inclusive=false)
+                @test all(check_query.(nside, ordering, θ, ϕ, radius, exclusive_pixels))
+            end
+        end
+    end
 end