source: Sophya/trunk/SophyaProg/Tests/sph2lm.cc@ 2025

Last change on this file since 2025 was 2005, checked in by ansari, 23 years ago

Ajout programme sph2lm.cc , projection/calcul FFT healpix->cartes locales - Reza 15/5/2002
File size: 17.9 KB
Line 
1#include "tarray.h"
2#include "cimage.h"
3#include "skymap.h"
4#include "mapoperation.h"
5#include "sambainit.h"
6#include "fftwserver.h"
7
8#include <typeinfo>
9
10
11//--------------------------------------------------------------------------
12// Classe pour projeter localement (dans un tableau indexe i,j)
13// des coordonnees spheriques
14
15class MtxLM {
16public:
17 MtxLM(double latdeg, double longdeg,
18 double deltheta_min, double delphi_min,
19 int_4 npixtheta, int_4 npixphi);
20 MtxLM(TMatrix<r_8> & mtx, bool uselat=true);
21
22 virtual ~MtxLM();
23
24 void Init(double latdeg, double longdeg,
25 double deltheta_min, double delphi_min,
26 int_4 npixtheta, int_4 npixphi);
27
28 // Conversion d'index de pixel (ip,jt) (ip: Phi/X/Column, jt: Theta/Y/Row )
29 // en Theta,Phi
30 void ip_jt2tp(int_4 ip, int_4 jt, double& theta, double& phi);
31 // Conversion de (Theta,Phi) en index de pixel (ip,jt)
32 // (ip: Phi/X/Column, jt: Theta/Y/Row )
33 void tp2ipjt(double theta, double phi, int_4& ip, int_4& jt);
34
35 double _latdeg, _longdeg;
36 double _deltheta_min, _delphi_min;
37 int_4 _npixtheta, _npixphi;
38 double _theta0, _phi0;
39 double _thetaC, _phiC;
40 double _deltatheta, _deltaphi;
41 static double _deuxpi;
42};
43
44double MtxLM::_deuxpi = 2.*M_PI;
45
46MtxLM::MtxLM(double latdeg, double longdeg,
47 double deltheta_min, double delphi_min,
48 int_4 npixtheta, int_4 npixphi)
49{
50 Init(latdeg, longdeg, deltheta_min, delphi_min, npixtheta, npixphi );
51}
52
53//--------------------------------------------------------------------------
54
55MtxLM::MtxLM(TMatrix<r_8> & mtx, bool uselat)
56{
57 double latdeg = (uselat) ? (double)mtx.Info()["Latitude_C"] : 0.;
58 double longdeg = mtx.Info()["Longitude_C"];
59 double deltheta_min = mtx.Info()["LatPixSize"];
60 double delphi_min = mtx.Info()["LongPixSize"];
61 int_4 npixtheta = mtx.NRows();
62 int_4 npixphi = mtx.NCols();
63 Init(latdeg, longdeg, deltheta_min, delphi_min, npixtheta, npixphi);
64}
65
66MtxLM::~MtxLM()
67{
68}
69
70void MtxLM::Init(double latdeg, double longdeg,
71 double deltheta_min, double delphi_min,
72 int_4 npixtheta, int_4 npixphi)
73{
74 _latdeg = latdeg;
75 _longdeg = longdeg;
76 _deltheta_min = deltheta_min;
77 _delphi_min = delphi_min;
78 _npixtheta = npixtheta;
79 _npixphi = npixphi;
80 _thetaC = (90.-latdeg)*M_PI/180.;
81 _phiC = longdeg*M_PI/180.;
82 _deltatheta = deltheta_min*M_PI/(180.*60.);
83 _deltaphi = delphi_min*M_PI/(180.*60.)/sin(_thetaC);
84 _theta0 = _thetaC-0.5*(double)_npixtheta*_deltatheta;
85 if ((_thetaC + 0.5*(double)_npixtheta*_deltatheta) > M_PI-1.e-6)
86 throw RangeCheckError("MtxLM::Init() ThetaMax out of range (> M_PI-1.e-6)");
87 if ((_thetaC - 0.5*(double)_npixtheta*_deltatheta) < 1.e-6)
88 throw RangeCheckError("MtxLM::Init() ThetaMin out of range (< 1.e-6)");
89
90 double deltaphitot = delphi_min*M_PI/(180.*60.)/sin(_thetaC + 0.5*(double)_npixtheta*_deltatheta);
91 if (deltaphitot*_npixphi > _deuxpi)
92 throw RangeCheckError("MtxLM::Init() Wrapping will occur for Phi at ThetaMax ! ");
93 deltaphitot = delphi_min*M_PI/(180.*60.)/sin(_thetaC - 0.5*(double)_npixtheta*_deltatheta);
94 if (deltaphitot*_npixphi > _deuxpi)
95 throw RangeCheckError("MtxLM::Init() Wrapping will occur for Phi at ThetaMin ! ");
96
97 _phi0 = _phiC-0.5*(double)_npixphi*_deltaphi;
98 return;
99}
100
101void MtxLM::ip_jt2tp(int_4 ip, int_4 jt, double& theta, double& phi)
102{
103 theta = jt*_deltatheta+_theta0;
104 if ( (theta<0.) || (theta > M_PI) )
105 throw RangeCheckError(" MtxLM::ip_jt2tp - Theta out of range !");
106 double delphi = _delphi_min*M_PI/(180.*60.)/sin(theta);
107 double phi0 = _phiC-0.5*(double)_npixphi*delphi;
108 // phi = ip*_deltaphi+_phi0;
109 phi = ip*delphi+phi0;
110
111 while (phi < 0.) phi += _deuxpi;
112 while (phi >=_deuxpi ) phi -= _deuxpi;
113
114 if ( (phi<0.) || (phi > _deuxpi) )
115 throw RangeCheckError(" MtxLM::ip_jt2tp - Phi out of range !");
116 return;
117}
118
119void MtxLM::tp2ipjt(double theta, double phi, int_4& ip, int_4& jt)
120{
121 if ( (theta<0.) || (theta > M_PI) )
122 throw RangeCheckError(" MtxLM::tp2ipjt - Theta out of range !");
123 if ( (phi<0.) || (phi > _deuxpi) )
124 throw RangeCheckError(" MtxLM::tp2ipjt - Phi out of range !");
125
126 jt = (theta-_theta0)/_deltatheta;
127 double delphi = _delphi_min*M_PI/(180.*60.)/sin(theta);
128 double phi0 = _phiC-0.5*(double)_npixphi*delphi;
129 if (phi0 < 0.) phi0 += 2*M_PI;
130 // ip = (phi-_phi0)/_deltaphi;
131 ip = (phi-phi0)/delphi;
132 if ((ip < 0) || (ip >= _npixphi) )
133 throw RangeCheckError(" MtxLM::ip_jt2tp - ip out of range !");
134 if ((jt < 0) || (jt >= _npixtheta) )
135 throw RangeCheckError(" MtxLM::ip_jt2tp - jt out of range !");
136 return;
137}
138
139//--------------------------------------------------------------------------
140// Fonction de filtrage de carte de type Matrix
141//--------------------------------------------------------------------------
142template <class T>
143void FilterMtxMap(TMatrix<T> & mtx, TMatrix<T> & mxout, int filw=1)
144{
145 if (filw < 1) return;
146 mxout = mtx;
147 sa_size_t r,c;
148 sa_size_t fw = 2*filw+1;
149 // Filtering
150 // Filling holes in Matrix
151 for(r=filw; r<mtx.NRows()-filw; r++)
152 for(c=filw; c<mtx.NCols()-filw; c++)
153 mxout(r,c) = mtx(Range(r-filw, 0, fw), Range(c-filw, 0, fw)).Sum()/fw*fw;
154
155 return;
156}
157
158//--------------------------------------------------------------------------
159// Projection depuis une carte spherique dans un local map
160//--------------------------------------------------------------------------
161
162template <class T>
163void Sph2LocalMap(SphericalMap<T> & in, LocalMap<T> & out)
164{
165 out.SetPixels(0);
166 int kout;
167 double teta,phi;
168
169 // Projecting to localMap
170 for(kout=0; kout<out.NbPixels(); kout++) {
171 out.PixThetaPhi(kout, teta, phi);
172 int pixel = in.PixIndexSph(teta,phi);
173 out(kout) = in.PixVal(pixel);
174 }
175
176 return;
177}
178
179//--------------------------------------------------------------------------
180// Projection depuis une matrice ds carte spherique
181//--------------------------------------------------------------------------
182template <class T>
183void ProjectMatrix2Sph(TMatrix<T> & mtx, SphericalMap<T> & out, bool userot=false)
184{
185 if (userot) {
186 double phi = (double)mtx.Info()["Longitude_C"]+90.;
187 double psi = -((double)mtx.Info()["Latitude_C"])*M_PI/180.;
188 while (phi >= 360.) phi -= 360.;
189 phi *= (M_PI/180.);
190 Vector3d omega(M_PI/2., phi);
191 ProjMatrix2RotateSph(mtx, out, omega, psi, false);
192 }
193 else ProjMatrix2Sph(mtx, out);
194 return;
195}
196
197//--------------------------------------------------------------------------
198// Projection depuis une matrice ds carte spherique sans rotation explicite
199//--------------------------------------------------------------------------
200template <class T>
201void ProjMatrix2Sph(TMatrix<T> & mtx, SphericalMap<T> & out)
202{
203 int kout;
204 double teta,phi;
205
206 // Projecting to Matrix
207 sa_size_t i,j,r,c;
208 MtxLM mtxlm(mtx);
209 for(r=0; r<mtx.NRows(); r++)
210 for(c=0; c<mtx.NCols(); c++) {
211 mtxlm.ip_jt2tp(c, r, teta, phi);
212 out(teta, phi) = mtx(r,c);
213 }
214 return;
215}
216
217//--------------------------------------------------------------------------
218// Projection depuis une matrice ds carte spherique avec rotation explicite
219//--------------------------------------------------------------------------
220template <class T>
221void ProjMatrix2RotateSph(TMatrix<T> & mtx, SphericalMap<T> & out,
222 Vector3d& omega, double phirot, bool uselat=true)
223{
224 double teta,phi;
225
226 // Projecting to Matrix
227 sa_size_t r,c;
228 MtxLM mtxlm(mtx, uselat);
229 for(r=0; r<mtx.NRows(); r++)
230 for(c=0; c<mtx.NCols(); c++) {
231 mtxlm.ip_jt2tp(c, r, teta, phi);
232 Vector3d vv1(teta, phi);
233 Vector3d vv2 = vv1.Rotate(omega, phirot);
234 out(vv2.Theta(), vv2.Phi()) = mtx(r,c);
235 }
236 return;
237}
238//--------------------------------------------------------------------------
239// Projection depuis carte spherique ds une matrice
240//--------------------------------------------------------------------------
241template <class T>
242void Sphere2Matrix(SphericalMap<T> & in, TMatrix<T> & mtx, bool userot=false)
243{
244 if (userot) {
245 double phi = (double)mtx.Info()["Longitude_C"]+90.;
246 double psi = -((double)mtx.Info()["Latitude_C"])*M_PI/180.;
247 while (phi >= 360.) phi -= 360.;
248 phi *= (M_PI/180.);
249 Vector3d omega(M_PI/2., phi);
250 RotateSph2Matrix(in, mtx, omega, psi, false);
251 }
252 else Sph2Matrix(in, mtx);
253 return;
254
255}
256
257//--------------------------------------------------------------------------
258// Projection depuis carte spherique ds une matrice sans rotation explicite
259//--------------------------------------------------------------------------
260template <class T>
261void Sph2Matrix(SphericalMap<T> & in, TMatrix<T> & mtx)
262{
263 int kout;
264 double teta,phi;
265
266 // Projecting to Matrix
267 sa_size_t i,j,r,c;
268 MtxLM mtxlm(mtx);
269 for(r=0; r<mtx.NRows(); r++)
270 for(c=0; c<mtx.NCols(); c++) {
271 mtxlm.ip_jt2tp(c, r, teta, phi);
272 mtx(r, c) = in.PixVal(in.PixIndexSph(teta,phi));
273 }
274
275 // Filling holes in Matrix
276 int_4 nholes = 0;
277 for(r=2; r<mtx.NRows()-2; r++)
278 for(c=2; c<mtx.NCols()-2; c++) {
279 if (fabs(mtx(r,c)) < 1.e-31) {
280 nholes++;
281 mtx(r,c) = mtx(Range(r-2, 0, 5), Range(c-2, 0, 5)).Sum()/24.;
282 }
283 }
284 cout << " Sph2Matrix()/Info: NbHoles= " << nholes << " filled " << endl;
285 return;
286}
287
288//--------------------------------------------------------------------------
289// Projection depuis carte spherique ds une matrice, avec une rotation
290// d'angle phi autour du vecteur omega
291//--------------------------------------------------------------------------
292template <class T>
293void RotateSph2Matrix(SphericalMap<T> & in, TMatrix<T> & mtx,
294 Vector3d& omega,double phirot, bool uselat=true)
295{
296 double teta,phi;
297
298 // Projecting to Matrix
299 sa_size_t r,c;
300 MtxLM mtxlm(mtx,uselat);
301 for(r=0; r<mtx.NRows(); r++)
302 for(c=0; c<mtx.NCols(); c++) {
303 mtxlm.ip_jt2tp(c, r, teta, phi);
304 Vector3d vv1(teta, phi);
305 Vector3d vv2 = vv1.Rotate(omega, phirot);
306 mtx(r, c) = in.PixVal(in.PixIndexSph(vv2.Theta(), vv2.Phi()));
307 }
308
309 // Filling holes in Matrix
310 int_4 nholes = 0;
311 for(r=2; r<mtx.NRows()-2; r++)
312 for(c=2; c<mtx.NCols()-2; c++) {
313 if (fabs(mtx(r,c)) < 1.e-31) {
314 nholes++;
315 mtx(r,c) = mtx(Range(r-2, 0, 5), Range(c-2, 0, 5)).Sum()/24.;
316 }
317 }
318 cout << " RotateSph2Matrix()/Info: NbHoles= " << nholes << " filled " << endl;
319 return;
320}
321
322
323//--------------------------------------------------------------------------
324//--------------------------------------------------------------------------
325
326template <class T>
327void CheckRotatedProjection(SphericalMap<T> & in);
328
329//--------------------------------------------------------------------------
330// ************** Le main *****************
331//--------------------------------------------------------------------------
332
333int main(int narg, char* arg[])
334{
335 #define NLM 9
336 // Tableau des theta,phi des cartes partielles
337 // latitude_map : Latitude en degre
338 // longitude_map : longitude en degre
339 // nx_map : Nb de pixels selon la longitude (Phi)
340 // ny_map : Nb de pixels selon la latitude (Theta)
341 // resolx_map : Resolution cartes selon la longitude (phi)
342 // resoly_map : Resolution cartes selon la latitue (theta)
343 // double latitude_map[NLM] = {0., -15., -15., 65., 65., 65., 65.};
344 double latitude_map[NLM] = {0., -15., -15., -28., 50., 60., 50., 75., 75.};
345 // double thetadeg_map[NLM] = {90., 105., 105., 25., 25., 25., 25.};
346 // double thetadeg_map[NLM] = {90., 105., 105., 115., 40., 30., 40., 15., 15.};
347 double thetadeg_map[NLM]; // = 90.-latitude_map[lnm]
348 double longitude_map[NLM] = {95., 100., 120., 130., 82., 90., 200., 140., 200.};
349 int nx_map[NLM] = {512, 512, 512, 512, 512, 512, 512, 512, 512};
350 int ny_map[NLM] = {512, 512, 512, 512, 512, 512, 512, 512, 512};
351 double resolx_map[NLM] = {2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5};
352 double resoly_map[NLM] = {2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5};
353
354
355 if ( (narg > 1) && (strcmp(arg[1], "-h") == 0) ) {
356 cout << " Usage: sph2lm InHealPixPPFName OutPPFName [-userot/-tstrot] [Filter1/2W = 0] "
357 << endl;
358 exit(0);
359 }
360 if (narg < 3) {
361 cout << " sph2lm/Erreur argument - sph2lm -h for usage" << endl;
362 exit(0);
363 }
364
365 bool userot = false;
366 bool tstrot = false;
367 if (narg > 3) {
368 if (strcmp(arg[3],"-userot") == 0) {
369 userot = true;
370 cout << " sph2lm: Explicit rotation will be used -> userot=true " << endl;
371 }
372 if (strcmp(arg[3],"-tstrot") == 0) {
373 userot = true;
374 cout << " sph2lm: Test rotated projection -> tstrot=true / CheckRotatedProjection() " << endl;
375 }
376 }
377 int filw = 0;
378 if (narg > 4) filw = atoi(arg[3]);
379
380 cout << "\n >>>> Starting sph2lm: arg[1]= " << arg[1]
381 << " arg[2]= " << arg[2] << " FW=" << filw << endl;
382
383 // We handle exception at the high level
384 try {
385 // This macro initialize the library
386 // static objects handle this - However, not all loader call
387 // the constructor for static objects
388 SophyaInit();
389
390 cout << " --- Reading input SphericalMap from file " << arg[1] << endl;
391 PInPersist inppf(arg[1]);
392 PPersist * pps = inppf.ReadObject();
393 AnyDataObj * obj = pps->DataObj();
394 cout << " Object type read from input PPF file : "
395 << typeid(*obj).name() << endl;
396 SphereHEALPix<r_8> * ps8 = dynamic_cast< SphereHEALPix<r_8>* > (obj);
397 SphereHEALPix<r_8> sphin;
398 if (ps8 != NULL) sphin.Share(*ps8);
399 else {
400 SphereHEALPix<r_4> * ps4 = dynamic_cast< SphereHEALPix<r_4>* > (obj);
401 if (ps4 != NULL) {
402 cout << " -- Converting SphereHEALPix<r_4> to SphereHEALPix<r_8> " << endl;
403 sphin.Resize(ps4->SizeIndex());
404 for(int kp=0; kp<sphin.NbPixels(); kp++) sphin(kp) = (*ps4)(kp);
405 }
406 else {
407 cout << " sph2lm/Error : Unsupported input object ! " << endl;
408 throw TypeMismatchExc("sph2lm/Error wrong input object type ");
409 }
410 }
411
412 // inppf.GetObject(sphin);
413 cout << " Sphere read: SizeIndex()= " << sphin.SizeIndex()
414 << " NbPixels()= " << sphin.NbPixels() << endl;
415
416 if (tstrot) {
417 CheckRotatedProjection(sphin);
418 cout << " sph2lm: CheckRotatedProjection() done --> exit() " << endl;
419 return(0);
420 }
421 cout << " --- Opening output PPF file " << arg[1] << endl;
422 POutPersist outppf(arg[2]);
423
424 FFTWServer FFTServ;
425 SphereHEALPix<r_8> sphck(128);
426 SphereHEALPix<r_8> sphckmtx(128);
427
428 int nlm;
429 for(nlm=0; nlm<NLM; nlm++) {
430 thetadeg_map[nlm] = 90.-latitude_map[nlm];
431 cout << "\n ------ LocalMap No " << nlm+1 << " Latitude=" << latitude_map[nlm]
432 << " Longitude= " << longitude_map[nlm] << endl;
433
434 LocalMap<r_8> lm(nx_map[nlm], ny_map[nlm]);
435 TMatrix<r_8> mtx(nx_map[nlm], ny_map[nlm]);
436
437 lm.SetOrigin(thetadeg_map[nlm], longitude_map[nlm]);
438 lm.SetSize(nx_map[nlm]*resolx_map[nlm]/60.,
439 ny_map[nlm]*resoly_map[nlm]/60.);
440
441 lm.SetPixels((nlm+1)*100.);
442 mtx.Info()["Latitude_C"] = latitude_map[nlm];
443 mtx.Info()["Longitude_C"] = longitude_map[nlm];
444 mtx.Info()["LatPixSize"] = resoly_map[nlm];
445 mtx.Info()["LongPixSize"] = resolx_map[nlm];
446 mtx = (nlm+1)*100.;
447
448 cout << " Doing lm.Project(sphck) ... " << endl;
449 lm.Project(sphck);
450 cout << " Doing mtx.Project(sphck) ... " << endl;
451 ProjectMatrix2Sph(mtx, sphckmtx, userot);
452
453 cout << " Doing Sph2LocalMap(sphin, lm) ... " << endl;
454 Sph2LocalMap(sphin, lm);
455 cout << " Doing Sphere2Matrix(sphin, lm) ... " << endl;
456 Sphere2Matrix(sphin, mtx, userot);
457 MuTyV num = nlm+1;
458 string name = "lm" + (string)num;
459 cout << " Writing local map and mtx to OutPPF " << endl;
460 outppf.PutObject(lm, name);
461 name = "mtx" + (string)num;
462 outppf.PutObject(mtx, name);
463
464 cout << " Calling FFTServ.FFTForward(mtx, fftmtx) " << endl;
465 TMatrix< complex<r_8> > fftmtx;
466 FFTServ.FFTForward(mtx, fftmtx);
467
468 name = "fftmtx" + (string)num;
469 cout << " Writing fftmtx to OutPPF " << endl;
470 outppf.PutObject(fftmtx, name);
471
472 if (filw < 1) continue;
473 TMatrix<r_8> fmtx;
474 cout << " Filtering mtx FilterW= " << filw << endl;
475 FilterMtxMap(mtx, fmtx, filw);
476 cout << " Calling FFTServ.FFTForward(fmtx, fftfmtx) " << endl;
477 TMatrix< complex<r_8> > fftfmtx;
478 FFTServ.FFTForward(fmtx, fftfmtx);
479
480 cout << " Writing fmtx and fftfmtx to OutPPF " << endl;
481 name = "fmtx" + (string)num;
482 outppf.PutObject(fmtx, name);
483 name = "fftfmtx" + (string)num;
484 outppf.PutObject(fftfmtx, name);
485 }
486
487 // cout << " Doing Sph2Sph(sphin, lm) ... " << endl;
488 // Sph2Sph(sphin, lm);
489 cout << "\n --- Writing sphck and sphckmtxsphericalmap to OutPPF " << endl;
490 outppf.PutObject(sphck, "sphck");
491 outppf.PutObject(sphckmtx, "sphckmtx");
492
493 }
494 catch (PThrowable & exc) {
495 cerr << " Catched Exception " << (string)typeid(exc).name()
496 << " - Msg= " << exc.Msg() << endl;
497 }
498 catch (...) {
499 cerr << " some other exception was caught ! " << endl;
500 }
501
502 cout << " >>>> End of sph2lm <<<< " << endl;
503}
504
505//--------------------------------------------------------------------------
506// Fonction pour la verification des projection avec rotation
507//--------------------------------------------------------------------------
508
509template <class T>
510void CheckRotatedProjection(SphericalMap<T> & in)
511{
512
513 cout << " >>> CheckRotatedProjection() : Testing rotated projection with matrices " << endl;
514 SphereHEALPix<r_8> sphckrot(128);
515 double latdeg = 0.;
516 double longdeg = 90.;
517 int nr = 512;
518 int nc = 512;
519 double resol = 2.5;
520
521 TMatrix<r_8> mtx(nr, nc);
522 TMatrix<r_8> mtxA, mtxB, mtxC;
523
524 mtx.Info()["Latitude_C"] = latdeg;
525 mtx.Info()["Longitude_C"] = longdeg;
526 mtx.Info()["LatPixSize"] = resol;
527 mtx.Info()["LongPixSize"] = resol;
528
529 mtx = 100.;
530 ProjMatrix2Sph(mtx, sphckrot);
531 Sph2Matrix(in, mtx);
532 mtxA = mtx;
533 Vector3d omega(M_PI/2., M_PI);
534 mtx = 200.;
535 ProjMatrix2RotateSph(mtx, sphckrot, omega, -M_PI/2.);
536 RotateSph2Matrix(in, mtx, omega, -M_PI/2.);
537 mtxB = mtx;
538 mtx = 300.;
539 ProjMatrix2RotateSph(mtx, sphckrot, omega, -M_PI/4.);
540 RotateSph2Matrix(in, mtx, omega, -M_PI/4.);
541 mtxC = mtx;
542
543 string ppfname = "rotproj.ppf";
544 cout << " -- Writing mtxA,B,C & sphckrot to OutPPF " << ppfname << endl;
545 POutPersist outppf(ppfname);
546
547 outppf.PutObject(mtxA, "mtxA");
548 outppf.PutObject(mtxB, "mtxB");
549 outppf.PutObject(mtxC, "mtxC");
550
551 outppf.PutObject(sphckrot, "sphckrot");
552
553 cout << " >>> End of CheckRotatedProjection() " << endl;
554
555 return;
556}
Note: See TracBrowser for help on using the repository browser.