| 1 | #include "machdefs.h" // Definitions specifiques SOPHYA
|
|---|
| 2 |
|
|---|
| 3 | #include <math.h>
|
|---|
| 4 | #include <iostream.h>
|
|---|
| 5 |
|
|---|
| 6 | #include "nbmath.h"
|
|---|
| 7 | #include "timing.h"
|
|---|
| 8 |
|
|---|
| 9 | #include "array.h"
|
|---|
| 10 | #include "skymap.h"
|
|---|
| 11 | #include "samba.h"
|
|---|
| 12 | #include "sambainit.h"
|
|---|
| 13 | #include "fitsspherehealpix.h"
|
|---|
| 14 | #include "fitstarray.h"
|
|---|
| 15 |
|
|---|
| 16 | /*!
|
|---|
| 17 | \ingroup PrgMap
|
|---|
| 18 | \file prjsmap.cc
|
|---|
| 19 | \brief \b prjsmap: Molleweide and sinus projection of sky maps
|
|---|
| 20 |
|
|---|
| 21 | \verbatim
|
|---|
| 22 |
|
|---|
| 23 | csh> prjsmap -h
|
|---|
| 24 | SOPHYA Version 1.1 Revision 0 (V_Fev2001) -- Mar 9 2001 15:45:31 cxx
|
|---|
| 25 |
|
|---|
| 26 | prjsmap : Molleweide and Sinus projection for sky maps
|
|---|
| 27 | Usage: prjsmap [-float/-r4] [-sinus] [-sx sizeX] [-sy sizeY] [-scale fact]
|
|---|
| 28 | [-off val][-ump val] [-invy] [-fitsin] [-fitsout] InFileName OutFileName
|
|---|
| 29 | -float (-r4): single precision sky map and projection (default = double)
|
|---|
| 30 | -sinus: Selects Sinus projection (default Molleweide projection)
|
|---|
| 31 | -sx sizeX: Projected map size in X
|
|---|
| 32 | -sy sizeY: Projected map size in Y
|
|---|
| 33 | default: sizeX = 2*sizeY ; sizeX*sizeY ~ map.NbPixels()/2
|
|---|
| 34 | -scale fact: scale factor applied to skymap pixels (def=1.)
|
|---|
| 35 | -off val: offset value for projected map (def=0.)
|
|---|
| 36 | prj(ix, jy) = offset + scale*skymap(theta, phi)
|
|---|
| 37 | -ump val: Unmapped pixel value in projected map (def=0.)
|
|---|
| 38 | -invy: reverses the projection direction along Y (Theta)
|
|---|
| 39 | -fitsout: Select the FITS format for the output map (default PPF format)
|
|---|
| 40 | -fitsin : Select the FITS format for the input vector (default PPF format)
|
|---|
| 41 | InFileName : Input file name (HEALPix map)
|
|---|
| 42 | OutFileName : Output file name (the C_l vector)
|
|---|
| 43 |
|
|---|
| 44 | \endverbatim
|
|---|
| 45 | */
|
|---|
| 46 |
|
|---|
| 47 |
|
|---|
| 48 | /* Nouvelle-Fonction */
|
|---|
| 49 | void Usage(bool fgerr)
|
|---|
| 50 | {
|
|---|
| 51 | cout << endl;
|
|---|
| 52 | if (fgerr) {
|
|---|
| 53 | cout << " prjsmap : Argument Error ! prjsmap -h for usage " << endl;
|
|---|
| 54 | exit(1);
|
|---|
| 55 | }
|
|---|
| 56 | else {
|
|---|
| 57 | cout << " prjsmap : Molleweide and Sinus projection for sky maps \n"
|
|---|
| 58 | << " Usage: prjsmap [-float/-r4] [-sinus] [-sx sizeX] [-sy sizeY] [-scale fact] \n"
|
|---|
| 59 | << " [-off val][-ump val] [-invy] [-fitsin] [-fitsout] InFileName OutFileName\n"
|
|---|
| 60 | << " -float (-r4): single precision sky map and projection (default = double) \n"
|
|---|
| 61 | << " -sinus: Selects Sinus projection (default Molleweide projection) \n"
|
|---|
| 62 | << " -sx sizeX: Projected map size in X \n"
|
|---|
| 63 | << " -sy sizeY: Projected map size in Y \n"
|
|---|
| 64 | << " default: sizeX = 2*sizeY ; sizeX*sizeY ~ map.NbPixels()/2 \n"
|
|---|
| 65 | << " -scale fact: scale factor applied to skymap pixels (def=1.) \n"
|
|---|
| 66 | << " -off val: offset value for projected map (def=0.) \n"
|
|---|
| 67 | << " prj(ix, jy) = offset + scale*skymap(theta, phi) \n"
|
|---|
| 68 | << " -ump val: Unmapped pixel value in projected map (def=0.) \n"
|
|---|
| 69 | << " -invy: reverses the projection direction along Y (Theta) \n"
|
|---|
| 70 | << " -fitsout: Select the FITS format for the output map (default PPF format) \n"
|
|---|
| 71 | << " -fitsin : Select the FITS format for the input vector (default PPF format) \n"
|
|---|
| 72 | << " InFileName : Input file name (HEALPix map) \n"
|
|---|
| 73 | << " OutFileName : Output file name (the C_l vector) \n" << endl;
|
|---|
| 74 | exit(0);
|
|---|
| 75 | }
|
|---|
| 76 | }
|
|---|
| 77 |
|
|---|
| 78 | /* Nouvelle-Fonction */
|
|---|
| 79 | template <class T>
|
|---|
| 80 | class _PrjMap {
|
|---|
| 81 | public :
|
|---|
| 82 | /* Methode */
|
|---|
| 83 | static void PM_MeanSigma(PixelMap<T> const & map, double& gmoy, double& gsig, T& min, T& max)
|
|---|
| 84 | {
|
|---|
| 85 | min = max = map(0);
|
|---|
| 86 | gmoy=0.;
|
|---|
| 87 | gsig = 0.;
|
|---|
| 88 | double valok;
|
|---|
| 89 | for(int k=0; k<map.NbPixels(); k++) {
|
|---|
| 90 | valok = map(k);
|
|---|
| 91 | gmoy += valok; gsig += valok*valok;
|
|---|
| 92 | if (map(k)<min) min = map(k);
|
|---|
| 93 | else if (map(k)>max) max = map(k);
|
|---|
| 94 | }
|
|---|
| 95 | gmoy /= (double)map.NbPixels();
|
|---|
| 96 | gsig = gsig/(double)map.NbPixels() - gmoy*gmoy;
|
|---|
| 97 | if (gsig >= 0.) gsig = sqrt(gsig);
|
|---|
| 98 | return;
|
|---|
| 99 | }
|
|---|
| 100 |
|
|---|
| 101 | /* Methode */
|
|---|
| 102 | static TArray<T> Project_Molleweide(SphericalMap<T> const & map, int sx=0, int sy=0, T ump=0,
|
|---|
| 103 | T offset=0, T scale=1, bool invy=false)
|
|---|
| 104 | {
|
|---|
| 105 | if ( (sx <= 0) && (sy <= 0) ) {
|
|---|
| 106 | sy = sqrt((double)(map.NbPixels()/2));
|
|---|
| 107 | sx = sy*2;
|
|---|
| 108 | }
|
|---|
| 109 | else {
|
|---|
| 110 | if (sx <= 0) sx = 2*sy;
|
|---|
| 111 | else sy = sx/2;
|
|---|
| 112 | }
|
|---|
| 113 | TArray<T> prj(sx, sy);
|
|---|
| 114 | prj = ump; // On met tout a ump
|
|---|
| 115 |
|
|---|
| 116 | cout << " Molleweide projection (sizeX= " << sx << " sizeY= " << sy << " )"
|
|---|
| 117 | << " Scale=" << scale << " Offset=" << offset << endl;
|
|---|
| 118 | r_8 xa, yd, teta,phi, facteur;
|
|---|
| 119 | int_4 ix, jy, k;
|
|---|
| 120 | for(jy=0; jy<sy; jy++) {
|
|---|
| 121 | yd = (r_8)(jy+0.5)/(r_8)sy-0.5;
|
|---|
| 122 | if (invy) teta = (0.5-yd)*M_PI;
|
|---|
| 123 | else teta = (yd+0.5)*M_PI;
|
|---|
| 124 | facteur=2.*M_PI/sin(acos((double)yd*2));
|
|---|
| 125 | for(ix=0; ix<sx; ix++) {
|
|---|
| 126 | xa = (r_8)(ix+0.5)/(r_8)sx-0.5;
|
|---|
| 127 | phi = xa*facteur+M_PI;
|
|---|
| 128 | if ( (phi <= 2*M_PI) && (phi >= 0.) ) {
|
|---|
| 129 | k = map.PixIndexSph(teta, phi);
|
|---|
| 130 | prj(ix,jy,0) = offset + scale*map(k);
|
|---|
| 131 | }
|
|---|
| 132 | }
|
|---|
| 133 | }
|
|---|
| 134 |
|
|---|
| 135 | prj.SetTemp(true);
|
|---|
| 136 | return prj;
|
|---|
| 137 | }
|
|---|
| 138 |
|
|---|
| 139 | /* Methode */
|
|---|
| 140 | static TArray<T> Project_Sinus(SphericalMap<T> const & map, int sx=0, int sy=0, T ump=0,
|
|---|
| 141 | T offset=0, T scale=1, bool invy=false)
|
|---|
| 142 | {
|
|---|
| 143 | if ( (sx <= 0) && (sy <= 0) ) {
|
|---|
| 144 | sy = sqrt((double)(map.NbPixels()/2));
|
|---|
| 145 | sx = sy*2;
|
|---|
| 146 | }
|
|---|
| 147 | else {
|
|---|
| 148 | if (sx <= 0) sx = 2*sy;
|
|---|
| 149 | else sy = sx/2;
|
|---|
| 150 | }
|
|---|
| 151 | TArray<T> prj(sx, sy);
|
|---|
| 152 | prj = ump; // On met tout a ump
|
|---|
| 153 |
|
|---|
| 154 | cout << " Sinus projection (sizeX= " << sx << " sizeY= " << sy << " )"
|
|---|
| 155 | << " Scale=" << scale << " Offset=" << offset << endl;
|
|---|
| 156 | r_8 xa, yd, teta,phi, facteur;
|
|---|
| 157 | int_4 ix, jy, k;
|
|---|
| 158 | for(jy=0; jy<sy; jy++) {
|
|---|
| 159 | yd = (r_8)(jy+0.5)/(r_8)sy-0.5;
|
|---|
| 160 | if (invy) teta = (0.5-yd)*M_PI;
|
|---|
| 161 | else teta = (yd+0.5)*M_PI;
|
|---|
| 162 | facteur=1./sin(teta);
|
|---|
| 163 | for(ix=0; ix<sx; ix++) {
|
|---|
| 164 | xa = (r_8)(ix+0.5)/(r_8)sx-0.5;
|
|---|
| 165 | phi = 2*M_PI*xa*facteur+M_PI;
|
|---|
| 166 | if ( (phi <= 2*M_PI) && (phi >= 0.) ) {
|
|---|
| 167 | k = map.PixIndexSph(teta, phi);
|
|---|
| 168 | prj(ix,jy,0) = offset + scale*map(k);
|
|---|
| 169 | }
|
|---|
| 170 | }
|
|---|
| 171 | }
|
|---|
| 172 | prj.SetTemp(true);
|
|---|
| 173 | return prj;
|
|---|
| 174 | }
|
|---|
| 175 |
|
|---|
| 176 | /* Methode */
|
|---|
| 177 | static void Project(string & infile, string & outfile, int sx, int sy, T ump, T offset,
|
|---|
| 178 | T scale, bool fgsin, bool invy, bool fgfitsin, bool fgfitsout)
|
|---|
| 179 | {
|
|---|
| 180 |
|
|---|
| 181 | double deg2rad = M_PI/180.;
|
|---|
| 182 | double minute2rad = M_PI/(180.*60.);
|
|---|
| 183 |
|
|---|
| 184 | SphericalMap<T> * sphp = NULL;
|
|---|
| 185 | SphereHEALPix<T> sphh;
|
|---|
| 186 |
|
|---|
| 187 | PPersist * pp = NULL;
|
|---|
| 188 |
|
|---|
| 189 | if (fgfitsin) {
|
|---|
| 190 | cout << "--- Reading Input FITS file " << infile << endl;
|
|---|
| 191 | FitsInFile fii(infile);
|
|---|
| 192 | fii >> sphh;
|
|---|
| 193 | sphp = &sphh;
|
|---|
| 194 | }
|
|---|
| 195 | else {
|
|---|
| 196 | cout << "--- Reading Input PPF file " << infile << endl;
|
|---|
| 197 | PInPersist ppi(infile);
|
|---|
| 198 | pp = ppi.ReadObject();
|
|---|
| 199 | if (pp)
|
|---|
| 200 | sphp = dynamic_cast< SphericalMap<T> * > (pp->DataObj());
|
|---|
| 201 | else sphp = NULL;
|
|---|
| 202 | if (sphp == NULL) {
|
|---|
| 203 | cout << " Object in PPF file is not a SphericalMap<T> or wrong data type ! " << endl;
|
|---|
| 204 | throw TypeMismatchExc("_PrjMap::Project(...) Error reading input PPF file !");
|
|---|
| 205 | }
|
|---|
| 206 | }
|
|---|
| 207 |
|
|---|
| 208 | SphericalMap<T> & sph = *sphp;
|
|---|
| 209 | T min, max;
|
|---|
| 210 | double mean, sigma;
|
|---|
| 211 | PM_MeanSigma(sph, mean, sigma, min, max);
|
|---|
| 212 |
|
|---|
| 213 | cout << " Input map type= " << sph.TypeOfMap() << " PixParm (=NSide for HEALPix)= "
|
|---|
| 214 | << sph.SizeIndex() << endl;
|
|---|
| 215 | cout << " Input map : NbPixels= " << sph.NbPixels() << " Resolution= "
|
|---|
| 216 | << sqrt(sph.PixSolAngle(0))/minute2rad << " Arcminutes " << endl;
|
|---|
| 217 | cout << " map.Min= " << min << " map.Max= " << max
|
|---|
| 218 | << " map.Mean= " << mean << " map.Sigma= " << sigma << endl;
|
|---|
| 219 |
|
|---|
| 220 | cout << "--- Sky map projection " << endl;
|
|---|
| 221 | TArray<T> prj;
|
|---|
| 222 |
|
|---|
| 223 | if (fgsin) prj = Project_Sinus(sph, sx, sy, ump, offset, scale, invy);
|
|---|
| 224 | else prj = Project_Molleweide(sph, sx, sy, ump, offset, scale, invy);
|
|---|
| 225 |
|
|---|
| 226 | cout << "--- Statistics on the projected map :" ;
|
|---|
| 227 | prj.Show();
|
|---|
| 228 | prj.MinMax(min, max);
|
|---|
| 229 | MeanSigma(prj, mean, sigma);
|
|---|
| 230 | cout << " prj.Min= " << min << " prj.Max= " << max
|
|---|
| 231 | << " prj.Mean= " << mean << " prj.Sigma= " << sigma << endl;
|
|---|
| 232 |
|
|---|
| 233 | if (fgfitsout) {
|
|---|
| 234 | cout << "--- Writing projection to Output FITS file " << outfile << endl;
|
|---|
| 235 | FitsOutFile fio(outfile);
|
|---|
| 236 | fio << prj;
|
|---|
| 237 | }
|
|---|
| 238 | else {
|
|---|
| 239 | cout << "--- Writing projection to Output PPF file " << outfile << endl;
|
|---|
| 240 | POutPersist ppo(outfile);
|
|---|
| 241 | ppo << prj;
|
|---|
| 242 | }
|
|---|
| 243 |
|
|---|
| 244 | if (pp) delete pp;
|
|---|
| 245 | }
|
|---|
| 246 |
|
|---|
| 247 | };
|
|---|
| 248 |
|
|---|
| 249 | /* Main-Program */
|
|---|
| 250 | int main(int narg, char *arg[])
|
|---|
| 251 | {
|
|---|
| 252 | if ((narg > 1) && (strcmp(arg[1],"-h") == 0) ) Usage(false);
|
|---|
| 253 |
|
|---|
| 254 | double ump = 0.;
|
|---|
| 255 | double off = 0.;
|
|---|
| 256 | double scale = 1.;
|
|---|
| 257 | int sx = 0;
|
|---|
| 258 | int sy = 0;
|
|---|
| 259 | string infile;
|
|---|
| 260 | string outfile;
|
|---|
| 261 | bool fgfitsin = false;
|
|---|
| 262 | bool fgfitsout = false;
|
|---|
| 263 | bool fgr4 = false;
|
|---|
| 264 | bool fgsinus = false;
|
|---|
| 265 | bool fginvy = false;
|
|---|
| 266 |
|
|---|
| 267 | cout << " prjsmap : Decoding command line options ... " << endl;
|
|---|
| 268 |
|
|---|
| 269 | int ko=1;
|
|---|
| 270 | for (int k=1; k<narg; k++) {
|
|---|
| 271 | if (strcmp(arg[k], "-sx") == 0) {
|
|---|
| 272 | if (k == narg-1) Usage(true); // -sx est suivi d'un argument
|
|---|
| 273 | sx = atoi(arg[k+1]); k++; // k++ pour sauter au suivant
|
|---|
| 274 | }
|
|---|
| 275 | if (strcmp(arg[k], "-sy") == 0) {
|
|---|
| 276 | if (k == narg-1) Usage(true); // -sy est suivi d'un argument
|
|---|
| 277 | sy = atoi(arg[k+1]); k++; // k++ pour sauter au suivant
|
|---|
| 278 | }
|
|---|
| 279 | else if (strcmp(arg[k], "-scale") == 0) {
|
|---|
| 280 | if (k == narg-1) Usage(true);
|
|---|
| 281 | scale = atof(arg[k+1]); k++;
|
|---|
| 282 | }
|
|---|
| 283 | else if (strcmp(arg[k], "-off") == 0) {
|
|---|
| 284 | if (k == narg-1) Usage(true);
|
|---|
| 285 | off = atof(arg[k+1]); k++;
|
|---|
| 286 | }
|
|---|
| 287 | else if (strcmp(arg[k], "-ump") == 0) {
|
|---|
| 288 | if (k == narg-1) Usage(true);
|
|---|
| 289 | ump = atof(arg[k+1]); k++;
|
|---|
| 290 | }
|
|---|
| 291 | else if (strcmp(arg[k], "-sinus") == 0) {
|
|---|
| 292 | fgsinus = true;
|
|---|
| 293 | }
|
|---|
| 294 | else if (strcmp(arg[k], "-invy") == 0) {
|
|---|
| 295 | fginvy = true;
|
|---|
| 296 | }
|
|---|
| 297 | else if (strcmp(arg[k], "-fitsin") == 0) {
|
|---|
| 298 | fgfitsin = true;
|
|---|
| 299 | }
|
|---|
| 300 | else if (strcmp(arg[k], "-fitsout") == 0) {
|
|---|
| 301 | fgfitsout = true;
|
|---|
| 302 | }
|
|---|
| 303 | else if ((strcmp(arg[k], "-float") == 0) || (strcmp(arg[k], "-r4") == 0) ) {
|
|---|
| 304 | fgr4 = true;
|
|---|
| 305 | }
|
|---|
| 306 |
|
|---|
| 307 | else { ko = k; break; } // Debut des noms
|
|---|
| 308 | }
|
|---|
| 309 |
|
|---|
| 310 | if ((narg-ko) < 2) Usage(true);
|
|---|
| 311 | infile = arg[ko];
|
|---|
| 312 | outfile = arg[ko+1];
|
|---|
| 313 |
|
|---|
| 314 | // Bloc try de gestion des exception
|
|---|
| 315 | try {
|
|---|
| 316 | InitTim();
|
|---|
| 317 | SophyaInit();
|
|---|
| 318 | if (fgr4) {
|
|---|
| 319 | cout << " Projecting SphericalMap<r_4> --> Array<r_4> (float)" << endl;
|
|---|
| 320 | _PrjMap<r_4>::Project(infile, outfile, sx, sy, ump, off, scale,
|
|---|
| 321 | fgsinus, fginvy, fgfitsin, fgfitsout);
|
|---|
| 322 | }
|
|---|
| 323 | else {
|
|---|
| 324 | cout << " Projecting SphericalMap<r_8> --> Array<r_8> (double)" << endl;
|
|---|
| 325 | _PrjMap<r_8>::Project(infile, outfile, sx, sy, ump, off, scale,
|
|---|
| 326 | fgsinus, fginvy, fgfitsin, fgfitsout);
|
|---|
| 327 | }
|
|---|
| 328 | }
|
|---|
| 329 | catch (PThrowable & exc) { // Exceptions de SOPHYA
|
|---|
| 330 | cerr << " prjsmap: Catched Exception " << (string)typeid(exc).name()
|
|---|
| 331 | << " - Msg= " << exc.Msg() << endl;
|
|---|
| 332 | }
|
|---|
| 333 | catch (...) { // Autres Exceptions
|
|---|
| 334 | cerr << " prjsmap: some other exception was caught ! " << endl;
|
|---|
| 335 | }
|
|---|
| 336 |
|
|---|
| 337 | PrtTim("End of prjsmap ");
|
|---|
| 338 | return(0);
|
|---|
| 339 | }
|
|---|
| 340 |
|
|---|
| 341 |
|
|---|
| 342 |
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
