+ cat2mask

python/spherelib/__init__.py

@@ -23,6 +23,6 @@
 """
 
 import map
-from spherelib_map import regress
+from spherelib_map import regress, fastfill, cat2mask
 from fitsfunc2 import *
 from bin import *

python/spherelib/spherelib_map.cpp

 #include "mask.h"
 #include "healpix_base.h"
 #include "fastfiller.h"
+#include "Python.h"
 
-void regress(double *tabmap, double *tabmask, int npix, int order=1)
+void regress(double *tabmap, int npix, double *tabmask, int npix2, int order=1)
 {
+  if (npix2 != npix)
+    PyErr_Format(PyExc_ValueError,"Arrays of different lengths (%d %d)",npix2, npix);
   arr<double> arrmap(tabmap, npix);
   arr<double> arrmask(tabmask, npix);
   MapExt<double> map(arrmap, RING);
@@ -25,12 +28,24 @@ void fastfill(double *tabmap, int npix,double eps=1e-6, int niter=-1, int starts
 
 }
 
-void inplace(double *invec, int n)
+void cat2mask(double *tabmask, int npix, double* theta, int ncat, double* phi, int ncat2, double radius=20 )
 {
-    int i;
-
-    for (i=0; i<n; i++)
-    {
-        invec[i] = 2*invec[i];
+  arr<double> arrmask(tabmask, npix);
+  MapExt<double> mask(arrmask, RING);
+  if (ncat2 != ncat)
+    PyErr_Format(PyExc_ValueError,"Arrays of different lengths (%d %d)",ncat2, ncat);
+  radius = radius /60.0 * degr2rad; // arcmin to degree
+  mask.fill(1.0);
+  vector<int> listpix;
+  vector<int>::iterator it;
+  for (int i=0; i<ncat; i++){
+    listpix.clear();
+    pointing p(theta[i], phi[i]); 
+    mask.query_disc(p, radius, listpix);
+    for(it = listpix.begin();it!=listpix.end();it++){
+      mask[*it]=0;
     }
+  }
+  arrmask = mask.Map();
 }
+

python/spherelib/spherelib_map.h

@@ -11,8 +11,8 @@
 #ifndef   	_SPHERELIB_MAP_H_
 # define        _SPHERELIB_MAP_H_
 
-void regress (double *mapreg, double *maskreg, int nreg, int order=1);
-void fastfill(double *mapfill, int nfill, double eps=1e-6, int niter=-1, int startscale=0);
-//void inplace (double *invec, int n);
+void regress (double *map1, int npix1, double *map2, int npix2, int order=1);
+void fastfill(double *map, int npix, double eps=1e-6, int niter=-1, int startscale=0);
+void cat2mask(double *map, int npix, double* theta, int ntheta, double* phi, int nphi, double radius=20 );
 
 #endif 	    /* !_SPHERELIB_MAP_H_ */

python/spherelib/spherelib_map.i

@@ -12,27 +12,20 @@
 %}
 
 // map regression
-%apply (double* INPLACE_ARRAY1, int DIM1) {(double* mapreg, int nreg1), (double* maskreg, int nreg2)};
+%apply (double* INPLACE_ARRAY1, int DIM1) {(double* map1, int npix1), (double* map2, int npix2)};
 %rename (regress) my_regress;
-%feature("docstring", "Remove multipole of HEALPix masked map up to order (1 or 2)") my_regress ;
-%feature("autodoc",0) my_regress;
-%inline %{
-  void my_regress (double* mapreg, int nreg1, double* maskreg, int nreg2, int order=1){
-     if (nreg1!=nreg2){
-     PyErr_Format(PyExc_ValueError, "Arrays of different lengths");
-     }
-     regress (mapreg, maskreg, nreg1, order);
-}
-%}
-%ignore regress;
+%feature("docstring", "Remove multipole of HEALPix masked map up to order (1 or 2)") regress ;
+%feature("autodoc",0) regress;
 
 // gap filling
-%apply (double* INPLACE_ARRAY1, int DIM1) {(double* mapfill, int nfill)};
-%feature("docstring", "Iteratively fills holes in an HEALPix map by putting in null and undef pixels the mean value of its neighbors") my_regress ;
+%apply (double* INPLACE_ARRAY1, int DIM1) {(double* map, int npix)};
+%feature("docstring", "Iteratively fills holes in an HEALPix map by putting in null and undef pixels the mean value of its neighbors") fastfill ;
 %feature("autodoc",0) fastfill;
 
-
-
+// cat2 mask
+%apply (double* IN_ARRAY1, int DIM1) {(double* theta, int ntheta), (double* phi, int nphi)};
+%feature("docstring", "Compute an Healpix map with null values around the sources. (input catalog should contain theta and phi coordinates (radian) as first and second column)") cat2mask ;
+%feature("autodoc",0) cat2mask;
 
 
 

python/spherelib/test.py

@@ -18,3 +18,8 @@ alm = healpy.map2alm(squeeze(map), 200)
 
 spherelib.write_alm ('alm.fits', alm)
 
+theta = array([0,0.2])
+phi = array([0.2,0])
+nside = 256
+mask = ones((nside*nside*12,))
+mask = spherelib.cat2mask(mask,theta, phi,radius=50)