Commit 88a5aba4 authored by Maude Le Jeune's avatar Maude Le Jeune
Browse files

+ ps gauss fit

parent e4dbd10b
......@@ -152,6 +152,14 @@ void Catalog::input(arr2<double> catarr){
}
}
void Catalog::input(double* tabcat, int ns, int ncol){
for (int i=0; i<ns; i++){
Source s;
s.input(tabcat[i*ncol+0], tabcat[i*ncol+1], tabcat[i*ncol+2], tabcat[i*ncol+3] , tabcat[i*ncol+4] );
sources.push_back(s);
}
}
void Catalog::input(istream &is){
char comment[1000];
......
......@@ -314,6 +314,7 @@ public:
void input(istream &is);
void input(arr2<double> catarr);
void input(double* tabcat, int ns, int ncol);
private:
double nsigma;
AtlasManager * sky;
......
......@@ -23,7 +23,7 @@
"""
from map import *
from spherelib_map import regress, fastfill, cat2mask, apodize_mask, _psdetection
from spherelib_map import regress, fastfill, _cat2mask, apodize_mask, _psdetection, _psgaussfit
from spherelib_alm import _alm2cov
from spherelib_bin import _make_prolate
from fitsfunc2 import *
......
......@@ -191,3 +191,9 @@ def glat2theta(glat):
def theta2glat(theta):
rad2deg = 180.0/pi
return 90- (theta*rad2deg)
def deg2rad(x):
return x*pi/180.0
def rad2deg(x):
return x*180.0/pi
......@@ -156,3 +156,48 @@ def ps_detection (map, fwhm_arcmin, nsigma=3):
lcat = len(cat)
cat = cat.reshape(len(cat)/5, 5)
return cat
def ps_gaussfit (map, smap, cat, fwhm_arcmin, nsigma=1, latcut=0, maxxi2=10):
""" Perform gaussian fit.
Parameters
----------
map: healpix map
smap: healpix map of the background
cat : ps catalog array like (nsource, 5) with columns corresponding to (pixel, glon, glat, Ampl, sigmaAmpl)
fwhm_arcmin: map resolution
nsigma: fit radius
latcut: cut on latitude (degree)
maxxi2: maximum x2 value
Returns
-------
array-like (nsource, 7) with columns corresponding to (theta, phi, Ampl, s_theta, s_phi, s_flux, xi2/dof)
"""
nsource = shape(cat)[0]
Bcat = zeros((nsource, 7))
Bcat[:,0:5] = cat
spherelib_map._psgaussfit(map, smap, Bcat, fwhm_arcmin, nsigma, latcut, maxxi2)
return Bcat
def cat2mask (mask, theta, phi, fwhm, snr=None):
""" Compute an Healpix map with null values around the
sources.
Parameters
----------
mask: array-like, healpix map
theta: array_like, theta coordinate of the sources
phi: array_like, phi coordinate of the sources
fwhm: float, radius or fwhm of the gaussian shape
theta: array_like, snr of the sources. If None, use fwhm as constant radius.
Returns
-------
array-like, healpix map
"""
if snr is None:
snr = zeros((0))
return spherelib_map._cat2mask (mask, theta, phi, snr, fwhm)
......@@ -30,6 +30,8 @@
#include <iterator>
#include "atlasmanager.h"
#include "catalog.h"
#include "math.h"
using namespace std;
void regress(double *tabmap, int npix, double *tabmask, int npix2, int order=1)
......@@ -57,19 +59,29 @@ void fastfill(double *tabmap, int npix,double eps=1e-6, int niter=-1, int starts
arrmap = map.Map();
}
void cat2mask(double *tabmask, int npix, double* theta, int ncat, double* phi, int ncat2, double radius=20 )
void _cat2mask(double *tabmask, int npix, double* theta, int ncat, double* phi, int ncat2, double* snr, int ncat3, double fwhm=20 )
{
arr<double> arrmask(tabmask, npix);
MapExt<double> mask(arrmask, RING);
bool cst = 0;
if (ncat2 != ncat)
PyErr_Format(PyExc_ValueError,"Arrays of different lengths (%d %d)",ncat2, ncat);
radius = radius /60.0 * degr2rad; // arcmin to degree
if (ncat3 ==0)
cst = 1;
else if (ncat3 != ncat2)
PyErr_Format(PyExc_ValueError,"Arrays of different lengths (%d %d)",ncat2, ncat3);
fwhm = fwhm /60.0 * degr2rad; // arcmin to degree
mask.fill(1.0);
vector<int> listpix;
vector<int>::iterator it;
double radius;
for (int i=0; i<ncat; i++){
listpix.clear();
pointing p(theta[i], phi[i]);
if (cst==1)
radius = fwhm;
else
radius = fwhm/2.0*sqrt(log(snr[i])/log(2));
mask.query_disc(p, radius, listpix);
for(it = listpix.begin();it!=listpix.end();it++){
mask[*it]=0;
......@@ -341,7 +353,71 @@ void _psdetection(double* tabmap, int npix, double** tabcat, int* ncat, double f
*tabcat = temp;
*ncat = ns*ncol;
}
void _psgaussfit (double* tabmap, int npix, double* sigmamap, int npix2, double* tabcat, int ns,int ncol, double fwhm_arcmin,double nfwhm ,double latcut,double maxxi2 ) {
if (npix != npix2)
PyErr_Format(PyExc_ValueError,"Arrays of different lengths (%d %d)",npix, npix2);
if (ncol != 7)
PyErr_Format(PyExc_ValueError,"Catalog arrays should have 7 columns (%d)",ncol);
arr<double> arrmap(tabmap, npix);
MapExt<double> hmap(arrmap, RING);
// -------- GAUSSIAN FIT
bool sigmap = true;
arr<double> arrmap2(sigmamap, npix);
Healpix_Map <double> smap(arrmap2, RING);
sigmap = true;
Catalog cat(NULL,0); // todo load tabcat in cat
//# pix, glon, glat, A, sigmaA
//# theta phi A sigmaTheta sigmaPhi sigmaA chi2/dof
cat.input(tabcat, ns, ncol);
list<Source>::iterator s;
double guess[3];
double err[3];
double chisqrdof;
double fwhm = fwhm_arcmin * degr2rad /60.;
double radius = nfwhm;
int i=0;
for(s = cat.begin(); s != cat.end(); s++){
pointing p = (*s).getPointing();
guess[1] = p.theta;
guess[2] = p.phi;
guess[0] = (*s).getAmplitude();
if (abs(90-(p.theta*180.0/pi)) < latcut)
chisqrdof = INFINITY;
else
healpixGaussianFit(hmap, fwhm, guess, err, &chisqrdof,smap,radius);
// fit is ok, set to new value
if (chisqrdof < maxxi2 && guess[0]>0){
addGaussianPatternOnSky(hmap, pointing(guess[1],guess[2]), fwhm, guess[0], err[0]/10, true);
tabcat[i*ncol+0] = guess[1];
tabcat[i*ncol+1] = guess[2];
tabcat[i*ncol+2] = guess[0];
tabcat[i*ncol+3] = err[1];
tabcat[i*ncol+4] = err[2];
tabcat[i*ncol+5] = err[0];
tabcat[i*ncol+6] = chisqrdof;
}
// fit is bad, keep old values
else {
tabcat[i*ncol+0] = p.theta;
tabcat[i*ncol+1] = p.phi;
tabcat[i*ncol+2] = (*s).getAmplitude();
tabcat[i*ncol+3] = 0;
tabcat[i*ncol+4] = 0;
tabcat[i*ncol+5] = (*s).getSigma();
tabcat[i*ncol+6] = chisqrdof;
}
i++;
}
arrmap = hmap.Map();
arrmap2.dealloc();
}
......@@ -30,7 +30,7 @@
using namespace std;
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 );
void _cat2mask(double *map, int npix, double* theta, int ntheta, double* phi, int nphi, double *snr, int nsnr, double fwhm=20 );
void apodize_mask(double *map, int npix, double radius=20,bool inside=false ,bool pixbool=0,bool spline=false );
void _mask2mll(double *map, int npix, double* mll, int nell1, int nell2);
void _make_kp_zone(double *map, int npix,double* covers, int nz, double* mapref, int npixref );
......@@ -38,5 +38,5 @@ void _make_grid_zone(double *map, int npix , int nside_in);
void _map2cov(double* map, int nin, int npix, double* mask, int npix2, double* cov, int ncross, int nz, double* nmode, int nz2, int reg=0);
void _combine_zone(double* map, int nin, int npix, double* mask, int npix2, double* cmap, int npix3, double* W, int nin2, int nz);
void _psdetection(double* map, int npix, double** cat, int* ncat, double fwhm_arcmin, double nsigma, int nside, double resolution, char* centers);
void _psgaussfit (double* cmap, int npix3, double* smap, int npix2, double* cat, int ns, int ncol,double fwhm_arcmin,double nfwhm ,double latcut,double maxxi2 ) ;
#endif /* !_SPHERELIB_MAP_H_ */
......@@ -42,9 +42,9 @@
%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;
%apply (double* IN_ARRAY1, int DIM1) {(double* theta, int ntheta), (double* phi, int nphi), (double* snr, int nsnr)};
%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;
// apodize mask
%feature("docstring", "Convert a binary mask (0 or 1) of arbitrary shape into smooth mask. ") apodize_mask ;
......@@ -86,6 +86,12 @@
%feature("docstring", "Perform point source detection using match filter.") _psdetection ;
%feature("autodoc",0) _psdetection;
//psgaussfit
%apply (double* INPLACE_ARRAY1, int DIM1) {(double* cmap, int npix3)};
%apply (double* IN_ARRAY1, int DIM1) {(double* smap, int npix2)};
%apply (double* INPLACE_ARRAY2, int DIM1, int DIM2) {(double* cat, int ns, int ncol)};
%feature("docstring", "Perform gaussian fit.") _psgaussfit ;
%feature("autodoc",0) _psgaussfit;
%include "spherelib_map.h"
......
......@@ -39,7 +39,7 @@ def configure( conf ):
conf.env.LIBPATH_HEALPIX = [op.join(hp_base_dir,'lib')]
conf.env.LIB_HEALPIX = ['cfitsio']
conf.env.STLIB_HEALPIX = ['healpix_cxx', 'cxxsupport', 'psht', 'c_utils', 'fftpack']
conf.env.DEFINES_HEALPIX = ['HEALPIXDATA="'+op.join(hp_base_dir,'../data"')]
conf.env.DEFINES_HEALPIX = ['HEALPIXDATA="'+op.abspath(op.join(hp_base_dir,'../data"'))]
conf.check_cxx(stlib='healpix_cxx', use='HEALPIX')
conf.check_cxx(stlib='cxxsupport', use='HEALPIX')
conf.check_cxx(stlib='psht', use='HEALPIX')
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment