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 | return prj;
|
---|
136 | }
|
---|
137 |
|
---|
138 | /* Methode */
|
---|
139 | static TArray<T> Project_Sinus(SphericalMap<T> const & map, int sx=0, int sy=0, T ump=0,
|
---|
140 | T offset=0, T scale=1, bool invy=false)
|
---|
141 | {
|
---|
142 | if ( (sx <= 0) && (sy <= 0) ) {
|
---|
143 | sy = sqrt((double)(map.NbPixels()/2));
|
---|
144 | sx = sy*2;
|
---|
145 | }
|
---|
146 | else {
|
---|
147 | if (sx <= 0) sx = 2*sy;
|
---|
148 | else sy = sx/2;
|
---|
149 | }
|
---|
150 | TArray<T> prj(sx, sy);
|
---|
151 | prj = ump; // On met tout a ump
|
---|
152 |
|
---|
153 | cout << " Sinus projection (sizeX= " << sx << " sizeY= " << sy << " )"
|
---|
154 | << " Scale=" << scale << " Offset=" << offset << endl;
|
---|
155 | r_8 xa, yd, teta,phi, facteur;
|
---|
156 | int_4 ix, jy, k;
|
---|
157 | for(jy=0; jy<sy; jy++) {
|
---|
158 | yd = (r_8)(jy+0.5)/(r_8)sy-0.5;
|
---|
159 | if (invy) teta = (0.5-yd)*M_PI;
|
---|
160 | else teta = (yd+0.5)*M_PI;
|
---|
161 | facteur=1./sin(teta);
|
---|
162 | for(ix=0; ix<sx; ix++) {
|
---|
163 | xa = (r_8)(ix+0.5)/(r_8)sx-0.5;
|
---|
164 | phi = 2*M_PI*xa*facteur+M_PI;
|
---|
165 | if ( (phi <= 2*M_PI) && (phi >= 0.) ) {
|
---|
166 | k = map.PixIndexSph(teta, phi);
|
---|
167 | prj(ix,jy,0) = offset + scale*map(k);
|
---|
168 | }
|
---|
169 | }
|
---|
170 | }
|
---|
171 | return prj;
|
---|
172 | }
|
---|
173 |
|
---|
174 | /* Methode */
|
---|
175 | static void Project(string & infile, string & outfile, int sx, int sy, T ump, T offset,
|
---|
176 | T scale, bool fgsin, bool invy, bool fgfitsin, bool fgfitsout)
|
---|
177 | {
|
---|
178 |
|
---|
179 | double deg2rad = M_PI/180.;
|
---|
180 | double minute2rad = M_PI/(180.*60.);
|
---|
181 |
|
---|
182 | SphericalMap<T> * sphp = NULL;
|
---|
183 | SphereHEALPix<T> sphh;
|
---|
184 |
|
---|
185 | PPersist * pp = NULL;
|
---|
186 |
|
---|
187 | if (fgfitsin) {
|
---|
188 | cout << "--- Reading Input FITS file " << infile << endl;
|
---|
189 | FitsInFile fii(infile);
|
---|
190 | fii >> sphh;
|
---|
191 | sphp = &sphh;
|
---|
192 | }
|
---|
193 | else {
|
---|
194 | cout << "--- Reading Input PPF file " << infile << endl;
|
---|
195 | PInPersist ppi(infile);
|
---|
196 | pp = ppi.ReadObject();
|
---|
197 | if (pp)
|
---|
198 | sphp = dynamic_cast< SphericalMap<T> * > (pp->DataObj());
|
---|
199 | else sphp = NULL;
|
---|
200 | if (sphp == NULL) {
|
---|
201 | cout << " Object in PPF file is not a SphericalMap<T> or wrong data type ! " << endl;
|
---|
202 | throw TypeMismatchExc("_PrjMap::Project(...) Error reading input PPF file !");
|
---|
203 | }
|
---|
204 | }
|
---|
205 |
|
---|
206 | SphericalMap<T> & sph = *sphp;
|
---|
207 | T min, max;
|
---|
208 | double mean, sigma;
|
---|
209 | PM_MeanSigma(sph, mean, sigma, min, max);
|
---|
210 |
|
---|
211 | cout << " Input map type= " << sph.TypeOfMap() << " PixParm (=NSide for HEALPix)= "
|
---|
212 | << sph.SizeIndex() << endl;
|
---|
213 | cout << " Input map : NbPixels= " << sph.NbPixels() << " Resolution= "
|
---|
214 | << sqrt(sph.PixSolAngle(0))/minute2rad << " Arcminutes " << endl;
|
---|
215 | cout << " map.Min= " << min << " map.Max= " << max
|
---|
216 | << " map.Mean= " << mean << " map.Sigma= " << sigma << endl;
|
---|
217 |
|
---|
218 | cout << "--- Sky map projection " << endl;
|
---|
219 | TArray<T> prj;
|
---|
220 |
|
---|
221 | if (fgsin) prj = Project_Sinus(sph, sx, sy, ump, offset, scale, invy);
|
---|
222 | else prj = Project_Molleweide(sph, sx, sy, ump, offset, scale, invy);
|
---|
223 |
|
---|
224 | cout << "--- Statistics on the projected map :" ;
|
---|
225 | prj.Show();
|
---|
226 | prj.MinMax(min, max);
|
---|
227 | MeanSigma(prj, mean, sigma);
|
---|
228 | cout << " prj.Min= " << min << " prj.Max= " << max
|
---|
229 | << " prj.Mean= " << mean << " prj.Sigma= " << sigma << endl;
|
---|
230 |
|
---|
231 | if (fgfitsout) {
|
---|
232 | cout << "--- Writing projection to Output FITS file " << outfile << endl;
|
---|
233 | FitsOutFile fio(outfile, FitsFile::clear);
|
---|
234 | fio << prj;
|
---|
235 | }
|
---|
236 | else {
|
---|
237 | cout << "--- Writing projection to Output PPF file " << outfile << endl;
|
---|
238 | POutPersist ppo(outfile);
|
---|
239 | ppo << prj;
|
---|
240 | }
|
---|
241 |
|
---|
242 | if (pp) delete pp;
|
---|
243 | }
|
---|
244 |
|
---|
245 | };
|
---|
246 |
|
---|
247 | /* Main-Program */
|
---|
248 | int main(int narg, char *arg[])
|
---|
249 | {
|
---|
250 | if ((narg > 1) && (strcmp(arg[1],"-h") == 0) ) Usage(false);
|
---|
251 |
|
---|
252 | double ump = 0.;
|
---|
253 | double off = 0.;
|
---|
254 | double scale = 1.;
|
---|
255 | int sx = 0;
|
---|
256 | int sy = 0;
|
---|
257 | string infile;
|
---|
258 | string outfile;
|
---|
259 | bool fgfitsin = false;
|
---|
260 | bool fgfitsout = false;
|
---|
261 | bool fgr4 = false;
|
---|
262 | bool fgsinus = false;
|
---|
263 | bool fginvy = false;
|
---|
264 |
|
---|
265 | cout << " prjsmap : Decoding command line options ... " << endl;
|
---|
266 |
|
---|
267 | int ko=1;
|
---|
268 | for (int k=1; k<narg; k++) {
|
---|
269 | if (strcmp(arg[k], "-sx") == 0) {
|
---|
270 | if (k == narg-1) Usage(true); // -sx est suivi d'un argument
|
---|
271 | sx = atoi(arg[k+1]); k++; // k++ pour sauter au suivant
|
---|
272 | }
|
---|
273 | else if (strcmp(arg[k], "-sy") == 0) {
|
---|
274 | if (k == narg-1) Usage(true); // -sy est suivi d'un argument
|
---|
275 | sy = atoi(arg[k+1]); k++; // k++ pour sauter au suivant
|
---|
276 | }
|
---|
277 | else if (strcmp(arg[k], "-scale") == 0) {
|
---|
278 | if (k == narg-1) Usage(true);
|
---|
279 | scale = atof(arg[k+1]); k++;
|
---|
280 | }
|
---|
281 | else if (strcmp(arg[k], "-off") == 0) {
|
---|
282 | if (k == narg-1) Usage(true);
|
---|
283 | off = atof(arg[k+1]); k++;
|
---|
284 | }
|
---|
285 | else if (strcmp(arg[k], "-ump") == 0) {
|
---|
286 | if (k == narg-1) Usage(true);
|
---|
287 | ump = atof(arg[k+1]); k++;
|
---|
288 | }
|
---|
289 | else if (strcmp(arg[k], "-sinus") == 0) {
|
---|
290 | fgsinus = true;
|
---|
291 | }
|
---|
292 | else if (strcmp(arg[k], "-invy") == 0) {
|
---|
293 | fginvy = true;
|
---|
294 | }
|
---|
295 | else if (strcmp(arg[k], "-fitsin") == 0) {
|
---|
296 | fgfitsin = true;
|
---|
297 | }
|
---|
298 | else if (strcmp(arg[k], "-fitsout") == 0) {
|
---|
299 | fgfitsout = true;
|
---|
300 | }
|
---|
301 | else if ((strcmp(arg[k], "-float") == 0) || (strcmp(arg[k], "-r4") == 0) ) {
|
---|
302 | fgr4 = true;
|
---|
303 | }
|
---|
304 |
|
---|
305 | else { ko = k; break; } // Debut des noms
|
---|
306 | }
|
---|
307 |
|
---|
308 | if ((narg-ko) < 2) Usage(true);
|
---|
309 | infile = arg[ko];
|
---|
310 | outfile = arg[ko+1];
|
---|
311 |
|
---|
312 | // Bloc try de gestion des exception
|
---|
313 | try {
|
---|
314 | InitTim();
|
---|
315 | SophyaInit();
|
---|
316 | if (fgr4) {
|
---|
317 | cout << " Projecting SphericalMap<r_4> --> Array<r_4> (float)" << endl;
|
---|
318 | _PrjMap<r_4>::Project(infile, outfile, sx, sy, ump, off, scale,
|
---|
319 | fgsinus, fginvy, fgfitsin, fgfitsout);
|
---|
320 | }
|
---|
321 | else {
|
---|
322 | cout << " Projecting SphericalMap<r_8> --> Array<r_8> (double)" << endl;
|
---|
323 | _PrjMap<r_8>::Project(infile, outfile, sx, sy, ump, off, scale,
|
---|
324 | fgsinus, fginvy, fgfitsin, fgfitsout);
|
---|
325 | }
|
---|
326 | }
|
---|
327 | catch (PThrowable & exc) { // Exceptions de SOPHYA
|
---|
328 | cerr << " prjsmap: Catched Exception " << (string)typeid(exc).name()
|
---|
329 | << " - Msg= " << exc.Msg() << endl;
|
---|
330 | }
|
---|
331 | catch (...) { // Autres Exceptions
|
---|
332 | cerr << " prjsmap: some other exception was caught ! " << endl;
|
---|
333 | }
|
---|
334 |
|
---|
335 | PrtTim("End of prjsmap ");
|
---|
336 | return(0);
|
---|
337 | }
|
---|
338 |
|
---|
339 |
|
---|
340 |
|
---|