Commit f0461cda authored by Maude Le Jeune's avatar Maude Le Jeune
Browse files

adapt to new healpix

parent 1aa27558
This diff is collapsed.
......@@ -139,7 +139,7 @@ int main(int argc, char ** argv){
else apodize(distmap, mask, radius, pixbool,inside);
//writing results
outfits.create (outfile);
write_Healpix_map_to_fits (outfits,mask,FITSUTIL<double>::DTYPE);
write_Healpix_map_to_fits (outfits,mask,planckType<double>());
}
// save distance map
......@@ -148,7 +148,7 @@ int main(int argc, char ** argv){
// if(distmap[i] > radius) distmap[i] = Healpix_undef;
fitshandle outfits2;
outfits2.create (filedistmap);
write_Healpix_map_to_fits (outfits2,distmap,FITSUTIL<double>::DTYPE);
write_Healpix_map_to_fits (outfits2,distmap,planckType<double>());
// save options of distance computation in fits
char* keyname = "minpix";
......
......@@ -20,6 +20,7 @@
* You should have received a copy of the GNU General Public License
* along with this program; if not, see http://www.gnu.org/licenses/gpl.html
*/
#include "lsconstants.h"
#include <gsl/gsl_blas.h>
#include <gsl/gsl_blas_types.h>
#include <gsl/gsl_rng.h>
......
......@@ -30,7 +30,6 @@
#include "alm.h"
#include "healpix_base.h"
#include "fitshandle.h"
#include "alm_map_tools.h"
#include "healpix_data_io.h"
#include "AlmExt.h"
......
......@@ -28,6 +28,9 @@
#include "alm.h"
#include "trafos.h"
#include "alm_powspec_tools.h"
#include <stdio.h>
#include "lsconstants.h"
#include <cstdio>
using namespace std;
......@@ -165,7 +168,7 @@ class AlmExt : public Alm< xcomplex<T> >
mapU.SetUndefTo0();
double avg=mapT.average();
mapT.add(-avg);
mapT.Add(-avg);
if (mapT.Scheme()==NEST) mapT.swap_scheme();
if (mapQ.Scheme()==NEST) mapQ.swap_scheme();
......@@ -174,7 +177,7 @@ class AlmExt : public Alm< xcomplex<T> >
map2alm_pol_iter (mapT,mapQ,mapU,(*this),AlmE,AlmB,N_iter,weight_T);
(*this)(0,0) += avg*sqrt(fourpi);
mapT.add(avg);
mapT.Add(avg);
CD_returnOK;
}
......@@ -194,7 +197,7 @@ class AlmExt : public Alm< xcomplex<T> >
case 10: return Trafo(1950,1950,Ecliptic,Equatorial);
case 11: return Trafo(1950,1950,Ecliptic,Galactic);
case 12: return Trafo(1950,1950,Galactic,Ecliptic);
default: throw Message_error ("Unsupported transformation");
default: planck_fail("Unsupported transformation "+dataToString(num));
}
}
......@@ -245,7 +248,7 @@ class AlmExt : public Alm< xcomplex<T> >
Map.SetUndefTo0();
// map minus mean like anafast
double avg = Map.average();
Map.add(-avg);
Map.Add(-avg);
Healpix_Ordering_Scheme s = Map.Scheme();
if (s==NEST)
Map.swap_scheme();
......@@ -253,7 +256,7 @@ class AlmExt : public Alm< xcomplex<T> >
// like anafast
(*this)(0, 0) += avg*sqrt(fourpi);
// Add avg don't change map
Map.add(avg);
Map.Add(avg);
CD_returnOK;
}
......
......@@ -187,66 +187,179 @@ class CovMat
CD_returnOK;
}
///////////////////////////////////////////////// /////////////////////////////////////////////////
/// constructor : /// constructor :
/// \param Maps : collection Healpix maps
/// \param NbMap : number of maps
/// \param fwhm : Gaussian filter
/// \param lmax : maximum order of analysis
/////////////////////////////////////////////////
Status FromMap (MapExt<T> *Maps, int NbMap, double fwhm, int lmax, int nside)
{
// check number of maps
if (NbMap>0)
m_NbIn = NbMap;
else
CD_return(12);
// check nside and ordering compatibility
for (int m1=1; m1<NbMap; m1++)
CD_iCheck( Maps[m1].Conformable(Maps[0]) );
m_NbBin = 12*nside*nside;
// set size of covmat
m_NbCross = m_NbIn*(m_NbIn+1)/2;
m_CovMat.alloc (m_NbCross, m_NbBin);
m_NbMode.alloc (m_NbBin);
m_CovMat.fill(0);
arr<T> Gaussbeam (lmax+1);
GaussBeam (Gaussbeam, fwhm);
int i = 0;
for (int m1=0; m1<NbMap; m1++)
for (int m2=0; m2<=m1; m2++)
{
MapExt<T> tmp (Maps[m1].Map(), Maps[m1].Scheme());
tmp *= Maps[m2];
AlmExt<T> alm (lmax, lmax);
alm.FromMap(tmp);
alm.ScaleL (Gaussbeam);
tmp.FromAlm(alm , nside, Maps[m1].Scheme());
for (int pix=0; pix<m_NbBin ; pix++)
m_CovMat[i][pix] = tmp[pix];
i++;
}
int w = (int)floor( (fwhm*(Maps[0].Nside()/2)*pi) / (2*60*180));
for (int pix=0; pix<m_NbBin ; pix++)
m_NbMode[pix] = w*w;
/////////////////////////////////////////////////
/// constructor :
/// \param Maps : collection Healpix maps
/// \param NbMap : number of maps
/// \param Mask : Healpix map which sets spatial zones on which average is performed
/////////////////////////////////////////////////
template <typename T2>Status FromMapMinusMean (MapExt<T> *Maps, long NbMap, MapExt<T2> &Mask)
{
// check number of maps
if (NbMap>0)
m_NbIn = NbMap;
else
CD_return(12);
// get number of zones
T2 min , max;
Mask.minmax(min,max);
m_NbBin = (long) max+1;
// set size of covmat
m_NbCross = m_NbIn*(m_NbIn+1)/2;
m_CovMat.alloc (m_NbCross, m_NbBin);
m_NbMode.alloc (m_NbBin);
m_CovMat.fill(0);
arr2<T> Vmean (m_NbIn, m_NbBin);
Vmean.fill(0);
// check nside and ordering compatibility
for (int m1=0; m1<NbMap; m1++)
CD_iCheck( Maps[m1].Conformable(Mask) );
for (int m1=0; m1<NbMap; m1++)
{
#pragma omp parallel
{
// cross term index
T tempwei[m_NbBin];
T tempval[m_NbBin];
for (long z=0 ; z<m_NbBin; z++){
tempwei[z] = 0;
tempval[z] = 0;
}
#pragma omp for
for (long pix=0; pix<Mask.Npix() ; pix++)
{
tempval[(int)Mask[pix]] += (Maps[m1])[pix];
tempwei[(int)Mask[pix]] += 1;
}
#pragma omp critical
{
//i = (int)((m1+1)*(m1+2)/2)-1;
for (long z=0 ; z<m_NbBin; z++)
{
m_NbMode[z] += tempwei[z];
Vmean[m1][z] += tempval[z];
}
}
}
for (long z=0 ; z<m_NbBin; z++)
Vmean[m1][z] /= m_NbMode[z]; // weight stats with number of pixel
}
int i = 0;
for (int m1=0; m1<NbMap; m1++)
for (int m2=0; m2<=m1; m2++)
{
m_NbMode.fill(0);
#pragma omp parallel
{
// cross term index
T tempwei[m_NbBin];
T tempval[m_NbBin];
for (long z=0 ; z<m_NbBin; z++){
tempwei[z] = 0;
tempval[z] = 0;
}
/* arr<T> tempwei (m_NbBin); */
/* arr<T> tempval (m_NbBin); // sum pixel by pixel product for each cross term */
/* tempwei.fill(0); */
/* tempval.fill(0); */
#pragma omp for
for (long pix=0; pix<Mask.Npix() ; pix++)
{
tempval[(int)Mask[pix]] += ((Maps[m1])[pix]-Vmean[m1][(int)Mask[pix]]) * ((Maps[m2])[pix]-Vmean[m2][(int)Mask[pix]]);// sum product
tempwei[(int)Mask[pix]] += 1;
}
#pragma omp critical
{
//i = (int)((m1+1)*(m1+2)/2)-1;
for (long z=0 ; z<m_NbBin; z++)
{
m_NbMode[z] += tempwei[z];
m_CovMat[i][z] += tempval[z];
}
}
}
for (long z=0 ; z<m_NbBin; z++)
m_CovMat[i][z] /= m_NbMode[z]; // weight stats with number of pixel
i++;
}
CD_returnOK;
}
/////////////////////////////////////////////////
/// constructor : /// constructor :
/// \param Maps : collection Healpix maps
/// \param NbMap : number of maps
/// \param fwhm : Gaussian filter
/// \param lmax : maximum order of analysis
/////////////////////////////////////////////////
Status FromMap (MapExt<T> *Maps, int NbMap, double fwhm, int lmax, int nside)
{
// check number of maps
if (NbMap>0)
m_NbIn = NbMap;
else
CD_return(12);
CD_returnOK;
}
// check nside and ordering compatibility
for (int m1=1; m1<NbMap; m1++)
CD_iCheck( Maps[m1].Conformable(Maps[0]) );
m_NbBin = 12*nside*nside;
// set size of covmat
m_NbCross = m_NbIn*(m_NbIn+1)/2;
m_CovMat.alloc (m_NbCross, m_NbBin);
m_NbMode.alloc (m_NbBin);
m_CovMat.fill(0);
arr<T> Gaussbeam (lmax+1);
GaussBeam (Gaussbeam, fwhm);
int i = 0;
for (int m1=0; m1<NbMap; m1++)
for (int m2=0; m2<=m1; m2++)
{
MapExt<T> tmp (Maps[m1].Map(), Maps[m1].Scheme());
tmp *= Maps[m2];
AlmExt<T> alm (lmax, lmax);
alm.FromMap(tmp);
alm.ScaleL (Gaussbeam);
tmp.FromAlm(alm , nside, Maps[m1].Scheme());
for (int pix=0; pix<m_NbBin ; pix++)
m_CovMat[i][pix] = tmp[pix];
i++;
}
int w = (int)floor( (fwhm*(Maps[0].Nside()/2)*pi) / (2*60*180));
for (int pix=0; pix<m_NbBin ; pix++)
m_NbMode[pix] = w*w;
CD_returnOK;
}
/////////////////////////////////////////////////
/// constructor :
/// \param Alms : collection Healpix alms
/// \param NbMap : number of maps
/////////////////////////////////////////////////
Status FromAlm (AlmExt<T> *Alms, long NbMap, double fsky=1.0)
{
// check number of maps
if (NbMap>0)
Status FromAlm (AlmExt<T> *Alms, long NbMap, double fsky=1.0)
{
// check number of maps
if (NbMap>0)
m_NbIn = NbMap;
else
CD_return(12);
......
......@@ -23,7 +23,7 @@
#define HME_DBG
#include "alm_map_tools.h"
#include "healpix_base.h"
#include "healpix_map.h"
#include "CD_CommonDefinition.h"
......
......@@ -30,6 +30,7 @@
#include <omp.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sort_double.h>
#include <stdio.h>
using namespace std;
......
......@@ -30,7 +30,6 @@
#include "alm.h"
#include "healpix_base.h"
#include "fitshandle.h"
#include "alm_map_tools.h"
#include "healpix_data_io.h"
#include "powspec.h"
#include "AlmExt.h"
......
......@@ -33,9 +33,10 @@
#include "alm.h"
#include "healpix_base.h"
#include "fitshandle.h"
#include "alm_map_tools.h"
//#include "alm_map_tools.h"
#include "healpix_data_io.h"
#include "powspec.h"
#include <stdio.h>
using namespace std;
......
......@@ -97,8 +97,8 @@ void AtlasManager::setMask(Healpix_Map<double> &mask)
double max;
for (int i=0; i<centers.size(); i++)
{
pointing pc;
vec2pnt(centers[i],pc);
pointing pc(centers[i]);
//vec2pnt(centers[i],pc);
sqmap.alloc(mapSize,mapSize,true);
sqmap.gnomonicProjection(mask, pc, mapRes, mapSize);
sqmap.minMax(min, max);
......@@ -109,15 +109,15 @@ void AtlasManager::setMask(Healpix_Map<double> &mask)
//---------------------------GETTERS------------------------------------
void AtlasManager::writeApodWin(){
ostringstream outfilename;
outfilename << atlasName << "mask.fits";
apodWin.toFits(outfilename.str());
}
// void AtlasManager::writeApodWin(){
// ostringstream outfilename;
// outfilename << atlasName << "mask.fits";
// apodWin.toFits(outfilename.str());
// }
void AtlasManager::getMap(int index, TangentPlane & sqmap, bool dftReady){
pointing pc;
vec2pnt(centers[index],pc);
pointing pc(centers[index]);
//vec2pnt(centers[index],pc);
sqmap.alloc(mapSize,mapSize,dftReady);
sqmap.gnomonicProjection(*sphere, pc, mapRes, mapSize);
//sqmap *= apodWin;
......@@ -143,22 +143,22 @@ void AtlasManager::whichMap(pointing pos, vector<int> & mI){
//-------------------------High Level Functions-----------------------
void AtlasManager::sphere2atlas(){
TangentPlane sqmap(mapSize,mapSize);
TGA_Image im;
for(unsigned int i = 0; i<centers.size(); i++){
pointing pc;
vec2pnt(centers[i],pc);
sqmap.gnomonicProjection(*sphere, pc, mapRes, mapSize);
sqmap *= apodWin;
ostringstream outfilename;
outfilename << atlasName << i << ".fits";
//outfilename << atlasName << i << ".tga";
sqmap.toFits(outfilename.str());
//sqmap.toTGA(im);
//im.write(outfilename.str());
}
}
// void AtlasManager::sphere2atlas(){
// TangentPlane sqmap(mapSize,mapSize);
// TGA_Image im;
// for(unsigned int i = 0; i<centers.size(); i++){
// pointing pc;
// vec2pnt(centers[i],pc);
// sqmap.gnomonicProjection(*sphere, pc, mapRes, mapSize);
// sqmap *= apodWin;
// ostringstream outfilename;
// outfilename << atlasName << i << ".fits";
// //outfilename << atlasName << i << ".tga";
// sqmap.toFits(outfilename.str());
// //sqmap.toTGA(im);
// //im.write(outfilename.str());
// }
// }
void AtlasManager::atlas2sphere(){
TangentPlane sqmap;
......@@ -174,7 +174,7 @@ void AtlasManager::atlas2sphere(){
for(int i = 0; i<centers.size(); i++){
ostringstream infilename;
infilename << atlasName << i << ".fits";
sqmap.fromFits(infilename.str());
//sqmap.fromFits(infilename.str());
reprojMap(sqmap, i, *sphere, weight,tv,tw,tn);
}
for(int ipix = 0; ipix<sphere->Npix(); ipix++){
......@@ -192,8 +192,8 @@ void AtlasManager::coverage(Healpix_Map<double> &hmap, double radius, bool arcmi
if (arcmin) radius *= degr2rad/60.0;
vector< int > listpix;
for(int i = 0; i<centers.size(); i++){
pointing pc;
vec2pnt(centers[i],pc);
pointing pc(centers[i]);
//vec2pnt(centers[i],pc);
hmap.query_disc (pc, radius, listpix);
vector<int>::iterator pix;
for(pix = listpix.begin(); pix != listpix.end(); pix++)
......@@ -222,8 +222,8 @@ void AtlasManager::reprojMap(TangentPlane & sqmap, int mapIndex,
Healpix_Map<double> &tempv,
Healpix_Map<double> &tempw,
Healpix_Map<unsigned char> &tempn){
pointing pntc;
vec2pnt(centers[mapIndex], pntc);
pointing pntc(centers[mapIndex]);
//vec2pnt(centers[mapIndex], pntc);
double lon0rad = pntc.phi;
double lat0rad = pi/2-pntc.theta;
......
......@@ -25,7 +25,7 @@
#ifndef ATLASMANAGER_H
#define ATLASMANAGER_H
#include "lsconstants.h"
#include <stdlib.h>
#include <math.h>
#include <fitshandle.h>
......@@ -220,8 +220,8 @@ public :
* \return The spherical coordinates of the specified pixel.
*/
pointing ixy2thetaphi(int mapIndex, int x, int y){
pointing pntc;
vec2pnt(centers[mapIndex], pntc);
pointing pntc (centers[mapIndex]);
//vec2pnt(centers[mapIndex], pntc);
double lon0rad = pntc.phi;
double lat0rad = pi/2-pntc.theta;
......@@ -245,8 +245,8 @@ public :
* \param y Y-coordinate on the map in pixels.
*/
void ithetaphi2xy(int mapIndex, pointing p, int & x, int & y){
pointing pntc;
vec2pnt(centers[mapIndex], pntc);
pointing pntc(centers[mapIndex]);;
//vec2pnt(centers[mapIndex], pntc);
double lon0rad = pntc.phi;
double lat0rad = pi/2-pntc.theta;
......
......@@ -199,11 +199,11 @@ void Catalog::extractSZ(double alpha, double beta, double gamma,Healpix_Map<doub
profileTemplate.fill(0.0);
profileTemplate.setWCS(sky->getRes(),sky->getRes(), 0,0,xsize/2,xsize/2);
profileTemplate.addXMMProfile(radius[idx],alpha,beta,gamma,1,xsize/2,xsize/2);
TGA_Image im;
if(show){
profileTemplate.toTGA(im);
im.write("profiletemplates.tga");
}
//TGA_Image im;
//if(show){
//profileTemplate.toTGA(im);
//im.write("profiletemplates.tga");
//}
profileTemplate.recenter(xsize/2,xsize/2);
valarray<double> noisespec;
......@@ -270,11 +270,11 @@ void Catalog::extractCatalog(double beamfwhm, Healpix_Map<double> & hmap,
double radpix = beamfwhm / sky->getRes();
beamTemplate.addGaussianBeam(beamfwhm,1,xsize/2,xsize/2);
beamTemplate.recenter(xsize/2,xsize/2);
TGA_Image im;
if(show){
beamTemplate.toTGA(im);
im.write("beamtemplates.tga");
}
//TGA_Image im;
//if(show){
// beamTemplate.toTGA(im);
// im.write("beamtemplates.tga");
//}
valarray<double> noisespec;
if(sigmanoise != (double)0){
......@@ -283,12 +283,12 @@ void Catalog::extractCatalog(double beamfwhm, Healpix_Map<double> & hmap,
}
for(int i = start; i < stop; i++){
sky->getMap(i, sqmap, true);
if(show){
ostringstream outfilename;
outfilename << "map_" << i << ".tga";
sqmap.toTGA(im);
im.write(outfilename.str());
}
//if(show){
// ostringstream outfilename;
// outfilename << "map_" << i << ".tga";
// sqmap.toTGA(im);
// im.write(outfilename.str());
//}
if(apodize)
sqmap *= (*apwin);
......@@ -304,12 +304,12 @@ void Catalog::extractCatalog(double beamfwhm, Healpix_Map<double> & hmap,
sqmap.matchedFilter(beamTemplate, noisespec);
else
sqmap.matchedFilter(beamTemplate);
if(show){
ostringstream outfilename;
outfilename << "filtered_map_" << i << ".tga";
sqmap.toTGA(im);
im.write(outfilename.str());
}
// if(show){
// ostringstream outfilename;
// outfilename << "filtered_map_" << i << ".tga";
// sqmap.toTGA(im);
// im.write(outfilename.str());
// }
sqmap.Zero2Undef();
sigma = sqmap.detectSources(xpos, ypos, flux, nsigma, beamfwhm, 1);
for(int is = 0 ; is < flux.size(); is++){
......@@ -320,15 +320,15 @@ void Catalog::extractCatalog(double beamfwhm, Healpix_Map<double> & hmap,
//pos = hmap.pix2ang(pix);
sources.push_back(Source(i, xpos[is], ypos[is], flux[is],
beamfwhm, pos, pix,sigma));
if(show) drawCircle(im, xpos[is], xsize-ypos[is],radpix);
//if(show) drawCircle(im, xpos[is], xsize-ypos[is],radpix);
}
}
if(show){
ostringstream outfilename;
outfilename << "detection_map_" << i << ".tga";
//sqmap.toTGA(im);
im.write(outfilename.str());
}
// if(show){
// ostringstream outfilename;
// outfilename << "detection_map_" << i << ".tga";
// //sqmap.toTGA(im);
// im.write(outfilename.str());
// }
}
sources.sort();
sources.unique();
......@@ -337,7 +337,7 @@ void Catalog::extractCatalog(double beamfwhm, Healpix_Map<double> & hmap,
void Catalog::medianDetection(double beamfwhm, Healpix_Map<double> & hmap, bool show){
int xsize = sky->getSize();
TangentPlane sqmap(xsize, xsize);
TGA_Image im;
//TGA_Image im;
double sigma;
double radpix = beamfwhm/ sky->getRes();
for(int i = 0; i < sky->nMaps(); i++){
......@@ -346,9 +346,9 @@ void Catalog::medianDetection(double beamfwhm, Healpix_Map<double> & hmap, bool
vector<double> flux;
sky->getMap(i, sqmap, false);
sqmap.medianFilter(2*radpix);
if(show){
sqmap.toTGA(im);
}
// if(show){
// sqmap.toTGA(im);
// }
sigma = sqmap.detectSources(xpos, ypos, flux, nsigma, beamfwhm, 1);
for(int is = 0 ; is < flux.size(); is++){
pointing pos;
......@@ -356,13 +356,13 @@ void Catalog::medianDetection(double beamfwhm, Healpix_Map<double> & hmap, bool