source: Sophya/trunk/SophyaLib/Samba/localmap.cc@ 254

Last change on this file since 254 was 228, checked in by ansari, 26 years ago

Creation du module DPC/Samba Reza 13/04/99
File size: 10.9 KB
Line 
1#include "localmap.h"
2#include <iostream.h>
3#include "nbmath.h"
4//*****************************************************************************
5//++
6// Class LocalMap
7//
8// include localmap.h nbmath.h
9//
10// A local map of a region of the sky, in cartesian coordinates.
11// It has an origin in (theta0, phi0), mapped to pixel(x0, y0)
12// (x0, y0 might be outside of this local map)
13// default value of (x0, y0) is middle of the map, center of
14// pixel(nx/2, ny/2)
15//
16// la carte est consideree comme un tableau a deux indices i et j, i etant
17// indice de colonne et j indice de ligne. La carte est supposee resider
18// dans un plan tangent, dont le point de tangence est repere (x0,y0) dans
19// la carte et (theta0, phi0) sur la sphere celeste. L'extension de la
20// carte est definie par les valeurs de deux angles couverts respectivement
21// par la totalite des pixels en x de la carte et la totalite des pixels
22// en y. (SetSize()).
23// On considere un "plan de reference" : plan tangent a la sphere celeste
24// aux angles theta=Pi/2 et phi=0. Dans ce plan L'origine des coordonnees
25// est le point de tangence. L'axe Ox est la tangente parallele a
26// l'equateur, dirige vers les phi croissants, l'axe Oy est parallele
27// au meridien, dirige vers le pole nord.
28// De maniere interne a la classe une carte est definie dans ce plan de
29// reference et transportee jusqu'au point (theta0, phi0) de sorte que
30// les axes restent paralleles aux meridiens et paralleles. L'utilisateur
31// peut definir sa carte selon un repere en rotation par rapport au repere
32// de reference (par l'angle entre le parallele et l'axe Ox souhaite --
33// methode SetOrigin(...))
34
35//
36//--
37//++
38//
39// Links Parents
40//
41// PixelMap
42//
43//--
44//++
45//
46// Links Descendants
47//
48//
49//--
50//++
51// Titre Constructeurs
52//--
53//++
54
55LocalMap::LocalMap()
56
57//--
58
59{InitNul();}
60
61//++
62
63LocalMap::LocalMap(int_4 nx, int_4 ny) : mSzX_(nx), mSzY_(ny)
64
65//--
66{
67 InitNul();
68 mNPix_=nx*ny;
69 mPix_ = new r_8[mNPix_];
70}
71
72
73//++
74// Titre Destructeur
75//--
76//++
77LocalMap::~LocalMap()
78
79//--
80{
81 Clear();
82}
83
84
85
86void LocalMap::WriteSelf(POutPersist& s) const
87{
88 char strg[256];
89 if (mInfo_) {sprintf(strg, "LocalMap: NPixX=%6d NPixY=%9d HasInfo", mSzX_, mSzY_);} else {
90 sprintf(strg, "LocalMap: NPixX=%6d NPixY=%9d ", mSzX_, mSzY_);
91 }
92 s.PutLine(strg);
93 if (mInfo_) mInfo_->Write(s);
94 s.PutI4(mSzX_);
95 s.PutI4(mSzY_);
96 s.PutI4(mNPix_);
97 s.PutI4(x0_);
98 s.PutI4(y0_);
99 s.PutI4(originFlag_);
100 s.PutI4(SzFlag_);
101 s.PutR4(cos_angle_);
102 s.PutR4(sin_angle_);
103 s.PutR8(theta0_);
104 s.PutR8( phi0_);
105 s.PutR8( tgAngleX_);
106 s.PutR8(tgAngleY_);
107 s.PutR8s(mPix_, mNPix_);
108}
109void LocalMap::ReadSelf(PInPersist& s)
110{
111 Clear();
112 char strg[256];
113 s.GetLine(strg, 255);
114// Pour savoir s'il y avait un DVList Info associe
115 bool hadinfo = false;
116 if (strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0) hadinfo = true;
117 if (hadinfo) { // Lecture eventuelle du DVList Info
118 if (mInfo_ == NULL) mInfo_ = new DVList;
119 mInfo_->Read(s);
120 }
121 s.GetI4(mSzX_);
122 s.GetI4(mSzY_);
123 s.GetI4(mNPix_);
124 s.GetI4(x0_);
125 s.GetI4(y0_);
126 s.GetI4(originFlag_);
127 s.GetI4(SzFlag_);
128 s.GetR4(cos_angle_);
129 s.GetR4(sin_angle_);
130 s.GetR8(theta0_);
131 s.GetR8( phi0_);
132 s.GetR8( tgAngleX_);
133 s.GetR8(tgAngleY_);
134 mPix_ = new r_8[mNPix_+1];
135 s.GetR8s(mPix_, mNPix_);
136}
137
138//++
139// Titre Methodes
140//--
141
142
143//++
144int_4 LocalMap::NbPixels() const
145
146// Retourne le nombre de pixels du découpage
147//--
148{
149 return(mNPix_);
150}
151
152//++
153r_8& LocalMap::PixVal(int_4 k)
154
155// Retourne la valeur du contenu du pixel d'indice k
156//--
157{
158 if ( (k<0) || (k >= mNPix_) ) {
159 cout << " LocalMap::PIxVal : exceptions a mettre en place" <<endl;
160 // THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
161 THROW(rangeCheckErr);
162
163 //throw "LocalMap::PIxVal Pixel index out of range";
164 }
165 //return (*(pixel_->Data()+k));
166 return(mPix_[k]);
167}
168
169//++
170
171r_8 const& LocalMap::PixVal(int_4 k) const
172
173// Retourne la valeur du contenu du pixel d'indice k
174//--
175{
176 if ( (k<0) || (k >= mNPix_) ) {
177 cout << " LocalMap::PIxVal : exceptions a mettre en place" <<endl;
178 // THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
179
180 throw "LocalMap::PIxVal Pixel index out of range";
181 }
182 //return *(pixel_->Data()+k);
183 return(mPix_[k]);
184}
185
186//++
187int_4 LocalMap::PixIndexSph(float theta, float phi) const
188
189// Retourne l'indice du pixel vers lequel pointe une direction définie par
190// ses coordonnées sphériques
191//--
192{
193 int i,j;
194 if (!(originFlag_==1) || !(SzFlag_==1) ) {
195 cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
196 exit(0);
197 }
198 // theta et phi en coordonnees relatives (on se ramene a une situation
199 // par rapport au plan de reference)
200 float theta_aux=theta;
201 float phi_aux=phi;
202 UserToReference(theta_aux, phi_aux);
203 // coordonnees dans le plan local en unites de pixels
204 float x,y;
205 AngleProjToPix(theta_aux, phi_aux, x,y);
206
207 float xmin=-x0_-0.5;
208 float xmax=xmin+mSzX_;
209 if( (x > xmax) || (x < xmin ) ) return(-1);
210 float xcurrent=xmin;
211 for( i=0; i < mSzX_; i++ ) {
212 xcurrent += 1.;
213 if( x < xcurrent ) break;
214 }
215 float ymin=-y0_-0.5;
216 float ymax=ymin+mSzY_;
217 if( (y > ymax) || (y < ymin ) ) return(-1);
218 float ycurrent=ymin;
219 for( j=0; j < mSzY_; j++ ) {
220 ycurrent += 1.;
221 if( y < ycurrent ) break;
222 }
223
224
225 return (j*mSzX_+i);
226}
227
228//++
229void LocalMap::PixThetaPhi(int_4 k, float& theta, float& phi) const
230
231// Retourne les coordonnées (theta,phi) du milieu du pixel d'indice k
232//--
233{
234 if (!(originFlag_==1) || !(SzFlag_==1) ) {
235 cout << " LocalMap: correspondance carte-sphere non etablie" << endl;
236 exit(0);
237 }
238 int i,j;
239 Getij(k,i,j);
240 float X=float(i-x0_);
241 float Y=float(j-y0_);
242 // situation de ce pixel dans le plan de reference
243 float x=X*cos_angle_-Y*sin_angle_;
244 float y=X*sin_angle_+Y* cos_angle_;
245 // projection sur la sphere
246 PixProjToAngle(x, y,theta, phi);
247 // retour au plan utilisateur
248 ReferenceToUser(theta, phi);
249}
250
251//++
252r_8 LocalMap::PixSolAngle(int_4 k) const
253// Pixel Solid angle (steradians)
254// All the pixels have not necessarly the same size in (theta, phi)
255// because of the projection scheme which is not yet fixed.
256//
257//--
258{
259 int i,j;
260 Getij(k,i,j);
261 float X=float(i-x0_);
262 float Y=float(j-y0_);
263 float XR=X+float(i)*0.5;
264 float XL=X-float(i)*0.5;
265 float YU=Y+float(j)*0.5;
266 float YL=Y-float(j)*0.5;
267 // situation dans le plan de reference
268 float x0=XL*cos_angle_-YL*sin_angle_;
269 float y0=XL*sin_angle_+YL* cos_angle_;
270 float xa=XR*cos_angle_-YL*sin_angle_;
271 float ya=XR*sin_angle_+YL* cos_angle_;
272 float xb=XL*cos_angle_-YU*sin_angle_;
273 float yb=XL*sin_angle_+YU* cos_angle_;
274 // projection sur la sphere
275 float tet0,phi0,teta,phia,tetb,phib;
276 PixProjToAngle(x0, y0,tet0, phi0);
277 PixProjToAngle(xa, ya,teta, phia);
278 PixProjToAngle(xb, yb,tetb, phib);
279 // angle solide
280 float sol=fabs((xa-x0)*(yb-y0)-(xb-x0)*(ya-y0));
281 return r_8(sol);
282}
283
284
285//++
286void LocalMap::SetOrigin(float theta0, float phi0, float angle)
287
288//
289//--
290{
291 theta0_=theta0;
292 phi0_=phi0;
293 x0_=mSzX_/2;
294 y0_=mSzY_/2;
295 angle=angle*Pi/180.;
296 cos_angle_=cos(angle);
297 sin_angle_=sin(angle);
298 originFlag_ =1;
299}
300
301
302//++
303void LocalMap::SetOrigin(float theta0, float phi0, int_4 x0, int_4 y0, float angle)
304
305//
306//--
307{
308 theta0_=theta0;
309 phi0_=phi0;
310 x0_=x0;
311 y0_=y0;
312 angle=angle*Pi/180.;
313 cos_angle_=cos(angle);
314 sin_angle_=sin(angle);
315 originFlag_ =1;
316}
317
318//++
319void LocalMap::SetSize(float angleX, float angleY)
320
321//
322//--
323{
324 tgAngleX_=r_8(tan(angleX*Pi/360.));
325 tgAngleY_=r_8(tan(angleY*Pi/360.));
326 SzFlag_=1;
327}
328//++
329void LocalMap::Project(SphericalMap& sphere) const
330
331// Projection to spherical map
332//--
333{
334 for (int m=0; m<mNPix_; m++) {
335 float theta,phi;
336 PixThetaPhi(m,theta,phi);
337 sphere(theta,phi)=mPix_[m];
338 // cout << "theta " << theta << " phi " << phi << " valeur " << sphere(theta,phi)<< endl;
339 }
340}
341
342void LocalMap::InitNul()
343//
344// initialise à zéro certaines variables internes
345{
346originFlag_=0;
347SzFlag_=0;
348cos_angle_=1.;
349sin_angle_=0.;
350mPix_=NULL;
351}
352
353void LocalMap::Clear()
354//
355//
356{
357 if (mPix_) delete[] mPix_;
358}
359
360void LocalMap::Getij(int k, int& i, int& j) const
361{
362 i=(k+1)%mSzX_-1;
363 if (i==-1) i=mSzX_-1;
364 j=(k-i+2)/mSzX_;
365}
366
367void LocalMap::ReferenceToUser(float &theta, float &phi) const
368
369//
370{
371 if (theta > (r_4)Pi || theta < 0. || phi<0. || phi >=2* (r_4)Pi ) {
372 //cout << " LocalMap::ReferenceToUser : exceptions a mettre en place" <<endl;
373 // THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
374
375 throw "LocalMap::ReferenceToUser (theta,phi) out of range";
376 }
377 //cout << " ReferenceToUser entree, t= " << theta << " phi= " << phi << endl;
378 theta=(r_4)theta0_+theta-(r_4)Pi*0.5;
379 if (theta<0.) {
380 theta=-theta;
381 phi+=(r_4)Pi;
382 } else
383 if (theta>(r_4)Pi) {
384 theta=2.*(r_4)Pi-theta;
385 phi+=(r_4)Pi;
386 }
387
388 phi=(r_4)phi0_+phi;
389 while (phi>=2.*(r_4)Pi) phi-=2.*(r_4)Pi;
390 if (theta > (r_4)Pi || theta < 0. || phi<0. || phi >=2* (r_4)Pi ) {
391 cout << " LocalMap::ReferenceToUser : erreur bizarre dans le transfert a la carte utilisateur " << endl;
392 cout << " theta= " << theta << " phi= " << phi << endl;
393 exit(0);
394 }
395 //cout << " ReferenceToUser sortie, t= " << theta << " phi= " << phi << endl;
396}
397
398void LocalMap::UserToReference(float &theta, float &phi) const
399
400//
401{
402 if (theta > (r_4)Pi || theta < 0. || phi<0. || phi >=2* (r_4)Pi ) {
403 cout << " LocalMap::UserToReference : exceptions a mettre en place" <<endl;
404 // THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));
405
406 throw "LocalMap::UserToReference (theta,phi) out of range";
407 }
408 float phi1=phi-(r_4)phi0_;
409 if (phi1<0.) phi1+=2.*(r_4)Pi;
410
411 float theta1=theta-(r_4)theta0_+(r_4)Pi*0.5;
412 if (theta1<0.) {
413 theta=-theta1;
414 phi1+=(r_4)Pi;
415 } else
416 if (theta1>(r_4)Pi) {
417 theta=2.*(r_4)Pi-theta1;
418 phi1+=(r_4)Pi;
419 }
420 while (phi1>=2.*(r_4)Pi) phi1-=2.*(r_4)Pi;
421 phi=phi1;
422 if (theta > (r_4)Pi || theta < 0. || phi<0. || phi >=2* (r_4)Pi ) {
423 cout << " LocalMap::UserToReference : erreur bizarre dans le transfert a la carte de reference " << endl;
424 cout << " theta= " << theta << " phi= " << phi << endl;
425 exit(0);
426 }
427}
428
429void LocalMap::PixProjToAngle(float x, float y,float &theta, float &phi) const {
430
431// (x,y) representent les coordonnees en unites de pixels d'un point DANS
432// LE PLAN DE REFERENCE. On recupere (theta,phi) par rapport au repere "absolu"
433// theta=pi/2,phi=0.
434
435theta=(r_4)Pi*0.5-atan(2*y* tgAngleY_/(float) mSzY_);
436phi=atan2(2*x* tgAngleX_,(float) mSzX_);
437if( phi< 0. ) phi+= DeuxPi;
438}
439
440void LocalMap::AngleProjToPix(float theta, float phi, float& x, float& y) const {
441
442// (theta,phi) par rapport au repere "absolu" theta=pi/2,phi=0. On recupere
443// (i,j) DANS LE PLAN DE REFERENCE.
444 if (phi>=(r_4)Pi) phi-= DeuxPi;
445 // y=0.5*mSzY_*cot(theta)/tgAngleY_; $CHECK-REZA-04/99$
446 y=0.5*mSzY_/tan(theta)/tgAngleY_; // ? cot = 1/tan ?
447 x=0.5*mSzX_*tan(phi)/tgAngleX_;
448}
449
450
451
452
453
Note: See TracBrowser for help on using the repository browser.