Changeset 682 in Sophya for trunk/SophyaLib/Samba/localmap.cc
- Timestamp:
- Dec 10, 1999, 5:56:03 PM (26 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/SophyaLib/Samba/localmap.cc
r604 r682 6 6 #include <string.h> 7 7 #include <iostream.h> 8 9 10 //***************************************************************************** 11 //++ 12 // Class LocalMap 13 // 14 // include localmap.h nbmath.h 15 // 16 // A local map of a region of the sky, in cartesian coordinates. 17 // It has an origin in (theta0, phi0), mapped to pixel(x0, y0) 18 // (x0, y0 might be outside of this local map) 19 // default value of (x0, y0) is middle of the map, center of 20 // pixel(nx/2, ny/2) 21 // 22 // A local map is a 2 dimensional array, with i as column index and j 23 // as row index. The map is supposed to lie on a plan tangent to the 24 // celestial sphere in a point whose coordinates are (x0,y0) on the local 25 // map and (theta0, phi0) on the sphere. The range of the map is defined 26 // by two values of angles covered respectively by all the pixels in 27 // x direction and all the pixels in y direction (SetSize()). 28 // 29 // A "reference plane" is considered : this plane is tangent to the 30 // celestial sphere in a point with angles theta=Pi/2 and phi=0. This 31 // point is the origine of coordinates is of the reference plane. The 32 // x-axis is the tangent parallel to the equatorial line and oriented 33 // toward the increasing phi's ; the y-axis is parallel to the meridian 34 // line and oriented toward the north pole. 35 // 36 // Internally, a map is first defined within this reference plane and 37 // tranported until the point (theta0, phi0) in such a way that both 38 // axes are kept parallel to meridian and parallel lines of the sphere. 39 // The user can define its own map with axes rotated with respect to 40 // reference axes (this rotation is characterized by angle between 41 // the local parallel line and the wanted x-axis-- see method 42 // SetOrigin(...)) 43 // 44 // 45 46 // 47 //-- 48 //++ 49 // 50 // Links Parents 51 // 52 // PixelMap 53 // 54 //-- 55 //++ 56 // 57 // Links Childs 58 // 59 // 60 //-- 61 //++ 62 // Titre Constructors 63 //-- 64 //++ 65 template<class T> 66 LocalMap<T>::LocalMap() 67 // 68 // 69 //-- 70 { 71 InitNul(); 72 pixels_.Reset(); 73 } 74 75 //++ 76 template<class T> 77 LocalMap<T>::LocalMap(int_4 nx, int_4 ny) : nSzX_(nx), nSzY_(ny) 78 // 79 // 80 //-- 81 { 82 InitNul(); 83 nPix_= nx*ny; 84 pixels_.ReSize(nPix_); 85 pixels_.Reset(); 86 } 87 88 //++ 89 template<class T> 90 LocalMap<T>::LocalMap(const LocalMap<T>& lm, bool share) 91 : pixels_(lm.pixels_, share) 92 93 // 94 // copy constructor 95 //-- 96 { 97 cout<<" LocalMap:: Appel du constructeur de recopie " << endl; 98 99 if(lm.mInfo_) mInfo_= new DVList(*lm.mInfo_); 100 nSzX_= lm.nSzX_; 101 nSzY_= lm.nSzY_; 102 nPix_= lm.nPix_; 103 originFlag_= lm.originFlag_; 104 extensFlag_= lm.extensFlag_; 105 x0_= lm.x0_; 106 y0_= lm.y0_; 107 theta0_= lm.theta0_; 108 phi0_= lm.phi0_; 109 angle_= lm.angle_; 110 cos_angle_= lm.cos_angle_; 111 sin_angle_= lm.sin_angle_; 112 angleX_= lm.angleX_; 113 angleY_= lm.angleY_; 114 tgAngleX_= lm.tgAngleX_; 115 tgAngleY_= lm.tgAngleY_; 116 } 117 118 //++ 119 // Titre Destructor 120 //-- 121 //++ 122 template<class T> 123 LocalMap<T>::~LocalMap() 124 // 125 //-- 126 { 127 InitNul(); 128 } 129 130 131 132 //++ 133 // Titre Public Methods 134 //-- 135 136 //++ 137 template<class T> 138 void LocalMap<T>::ReSize(int_4 nx, int_4 ny) 139 // 140 // Resize storage area for pixels 141 //-- 142 { 143 InitNul(); 144 nSzX_ = nx; 145 nSzY_ = ny; 146 nPix_= nx*ny; 147 pixels_.ReSize(nPix_); 148 pixels_.Reset(); 149 } 150 151 //++ 152 template<class T> 153 int_4 LocalMap<T>::NbPixels() const 154 // 155 // Return number of pixels 156 //-- 157 { 158 return(nPix_); 159 } 160 161 //++ 162 template<class T> 163 T& LocalMap<T>::PixVal(int_4 k) 164 // 165 // Return value of pixel with index k 166 //-- 167 { 168 if((k < 0) || (k >= nPix_)) 169 { 170 cout << " LocalMap::PIxVal : exceptions a mettre en place" <<endl; 171 // THROW(out_of_range("LocalMap::PIxVal Pixel index out of range")); 172 THROW(rangeCheckErr); 173 //throw "LocalMap::PIxVal Pixel index out of range"; 174 } 175 return(pixels_(k)); 176 } 177 178 //++ 179 180 template<class T> 181 T const& LocalMap<T>::PixVal(int_4 k) const 182 // 183 // const version of previous method 184 //-- 185 { 186 if((k < 0) || (k >= nPix_)) 187 { 188 cout << " LocalMap::PIxVal : exceptions a mettre en place" <<endl; 189 // THROW(out_of_range("LocalMap::PIxVal Pixel index out of range")); 190 191 throw "LocalMap::PIxVal Pixel index out of range"; 192 } 193 return *(pixels_.Data()+k); 194 } 195 196 template<class T> 197 bool LocalMap<T>::ContainsSph(double theta, double phi) const 198 { 199 // $CHECK$ Reza 11/11/99 - A modifier 200 return(true); 201 } 202 203 //++ 204 template<class T> 205 int_4 LocalMap<T>::PixIndexSph(double theta,double phi) const 206 // 207 // Return index of the pixel with spherical coordinates (theta,phi) 208 //-- 209 { 210 int_4 i,j; 211 if(!(originFlag_) || !(extensFlag_)) 212 { 213 cout << " LocalMap: correspondance carte-sphere non etablie" << endl; 214 exit(0); 215 } 216 217 // theta et phi en coordonnees relatives (on se ramene a une situation par rapport au plan de reference) 218 double theta_aux= theta; 219 double phi_aux = phi; 220 UserToReference(theta_aux, phi_aux); 221 222 // coordonnees dans le plan local en unites de pixels 223 double x,y; 224 AngleProjToPix(theta_aux,phi_aux, x, y); 225 226 double xmin= -x0_-0.5; 227 double xmax= xmin+nSzX_; 228 if((x > xmax) || (x < xmin)) return(-1); 229 double xcurrent= xmin; 230 for(i = 0; i < nSzX_; i++ ) 231 { 232 xcurrent += 1.; 233 if( x < xcurrent ) break; 234 } 235 double ymin= -y0_-0.5; 236 double ymax= ymin+nSzY_; 237 if((y > ymax) || (y < ymin)) return(-1); 238 double ycurrent= ymin; 239 for(j = 0; j < nSzY_; j++ ) 240 { 241 ycurrent += 1.; 242 if( y < ycurrent ) break; 243 } 244 return (j*nSzX_+i); 245 } 246 247 //++ 248 template<class T> 249 void LocalMap<T>::PixThetaPhi(int_4 k,double& theta,double& phi) const 250 // 251 // Return (theta, phi) coordinates of pixel with index k 252 //-- 253 { 254 if(!(originFlag_) || !(extensFlag_)) 255 { 256 cout << " LocalMap: correspondance carte-sphere non etablie" << endl; 257 exit(0); 258 } 259 260 int_4 i,j; 261 Getij(k,i,j); 262 263 double X= double(i-x0_); 264 double Y= double(j-y0_); 265 // situation de ce pixel dans le plan de reference 266 double x= X*cos_angle_-Y*sin_angle_; 267 double y= X*sin_angle_+Y* cos_angle_; 268 // projection sur la sphere 269 PixProjToAngle(x, y, theta, phi); 270 // retour au plan utilisateur 271 ReferenceToUser(theta, phi); 272 } 273 274 template <class T> 275 T LocalMap<T>::SetPixels(T v) 276 { 277 pixels_.Reset(v); 278 return(v); 279 } 280 281 //++ 282 template<class T> 283 double LocalMap<T>::PixSolAngle(int_4 k) const 284 // 285 // Pixel Solid angle (steradians) 286 // All the pixels have not necessarly the same size in (theta, phi) 287 // because of the projection scheme which is not yet fixed. 288 //-- 289 { 290 int_4 i,j; 291 Getij(k,i,j); 292 double X= double(i-x0_); 293 double Y= double(j-y0_); 294 double XR= X+double(i)*0.5; 295 double XL= X-double(i)*0.5; 296 double YU= Y+double(j)*0.5; 297 double YL= Y-double(j)*0.5; 298 299 // situation dans le plan de reference 300 double x0= XL*cos_angle_-YL*sin_angle_; 301 double y0= XL*sin_angle_+YL*cos_angle_; 302 double xa= XR*cos_angle_-YL*sin_angle_; 303 double ya= XR*sin_angle_+YL*cos_angle_; 304 double xb= XL*cos_angle_-YU*sin_angle_; 305 double yb= XL*sin_angle_+YU*cos_angle_; 306 307 // projection sur la sphere 308 double thet0,phi0,theta,phia,thetb,phib; 309 PixProjToAngle(x0, y0, thet0, phi0); 310 PixProjToAngle(xa, ya, theta, phia); 311 PixProjToAngle(xb, yb, thetb, phib); 312 313 // angle solide 314 double sol= fabs((xa-x0)*(yb-y0)-(xb-x0)*(ya-y0)); 315 return sol; 316 } 317 318 //++ 319 template<class T> 320 void LocalMap<T>::SetOrigin(double theta0,double phi0,double angle) 321 // 322 // set the referential of the map (angles in degrees) 323 // (default x0=siz_x/2, y0=siz_y/2) 324 //-- 325 { 326 theta0_= theta0; 327 phi0_ = phi0; 328 angle_ = angle; 329 x0_= nSzX_/2; 330 y0_= nSzY_/2; 331 cos_angle_= cos(angle*Pi/180.); 332 sin_angle_= sin(angle*Pi/180.); 333 originFlag_= true; 334 cout << " LocalMap:: set origin 1 done" << endl; 335 } 336 337 //++ 338 template<class T> 339 void LocalMap<T>::SetOrigin(double theta0,double phi0,int_4 x0,int_4 y0,double angle) 340 // 341 // set the referential of the map (angles in degrees) 342 //-- 343 { 344 theta0_= theta0; 345 phi0_ = phi0; 346 angle_ = angle; 347 x0_= x0; 348 y0_= y0; 349 cos_angle_= cos(angle*Pi/180.); 350 sin_angle_= sin(angle*Pi/180.); 351 originFlag_= true; 352 cout << " LocalMap:: set origin 2 done" << endl; 353 } 354 355 //++ 356 template<class T> 357 void LocalMap<T>::SetSize(double angleX,double angleY) 358 // 359 // angle range of tthe map (angles in degrees) 360 //-- 361 { 362 angleX_= angleX; 363 angleY_= angleY; 364 365 // tangente de la moitie de l'ouverture angulaire totale 366 tgAngleX_= tan(0.5*angleX_*Pi/180.); 367 tgAngleY_= tan(0.5*angleY_*Pi/180.); 368 369 extensFlag_= true; 370 cout << " LocalMap:: set extension done" << endl; 371 } 372 373 //++ 374 template<class T> 375 void LocalMap<T>::Project(SphericalMap<T>& sphere) const 376 // 377 // Projection to a spherical map 378 //-- 379 { 380 for(int m = 0; m < nPix_; m++) 381 { 382 double theta,phi; 383 PixThetaPhi(m,theta,phi); 384 sphere(theta,phi)= pixels_(m); 385 // cout << "theta " << theta << " phi " << phi << " valeur " << sphere(theta,phi)<< endl; 386 } 387 } 388 // Titre Private Methods 389 //++ 390 template<class T> 391 void LocalMap<T>::InitNul() 392 // 393 // set some attributes to zero 394 //-- 395 { 396 originFlag_= false; 397 extensFlag_= false; 398 cos_angle_= 1.0; 399 sin_angle_= 0.0; 400 // pixels_.Reset(); Pas de reset par InitNul (en cas de share) - Reza 20/11/99 $CHECK$ 401 } 402 403 //++ 404 template<class T> 405 void LocalMap<T>::Getij(int_4 k,int_4& i,int_4& j) const 406 // 407 // Return 2 indices corresponding to the pixel number k 408 //-- 409 { 410 i= (k+1)%nSzX_-1; 411 if(i == -1) i= nSzX_-1; 412 j= (k-i+2)/nSzX_; 413 } 414 415 //++ 416 template<class T> 417 void LocalMap<T>::ReferenceToUser(double& theta,double& phi) const 418 // 419 // Transform a pair of coordinates (theta, phi) given in 420 // reference coordinates into map coordinates 421 //-- 422 { 423 if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi) 424 { 425 throw "LocalMap::ReferenceToUser (theta,phi) out of range"; 426 } 427 428 theta= theta0_*Pi/180.+theta-Pi*0.5; 429 if(theta < 0.) 430 { 431 theta= -theta; 432 phi += Pi; 433 } 434 else 435 { 436 if(theta > Pi) 437 { 438 theta= 2.*Pi-theta; 439 phi += Pi; 440 } 441 } 442 443 phi= phi0_*Pi/180.+phi; 444 while(phi >= 2.*Pi) phi-= 2.*Pi; 445 446 if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi) 447 { 448 cout << " LocalMap::ReferenceToUser : erreur bizarre dans le transfert a la carte utilisateur " << endl; 449 cout << " theta= " << theta << " phi= " << phi << endl; 450 exit(0); 451 } 452 } 453 454 //++ 455 template<class T> 456 void LocalMap<T>::UserToReference(double& theta,double& phi) const 457 // 458 // Transform a pair of coordinates (theta, phi) given in 459 // map coordinates into reference coordinates 460 //-- 461 { 462 if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi) 463 { 464 cout<<" LocalMap::UserToReference: exceptions a mettre en place" <<endl; 465 // THROW(out_of_range("LocalMap::PIxVal Pixel index out of range")); 466 throw "LocalMap::UserToReference (theta,phi) out of range"; 467 } 468 469 double phi1= phi-phi0_*Pi/180.; 470 if(phi1 < 0.) phi1+= 2.*Pi; 471 472 double theta1= theta-theta0_*Pi/180.+Pi*0.5; 473 if(theta1 < 0.) 474 { 475 theta= -theta1; 476 phi1+= Pi; 477 } 478 else 479 { 480 if(theta1 > Pi) 481 { 482 theta= 2.*Pi-theta1; 483 phi1+= Pi; 484 } 485 } 486 487 while(phi1 >= 2.*Pi) phi1-= 2.*Pi; 488 phi= phi1; 489 if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi) 490 { 491 cout << " LocalMap::UserToReference : erreur bizarre dans le transfert a la carte de reference " << endl; 492 cout << " theta= " << theta << " phi= " << phi << endl; 493 exit(0); 494 } 495 } 496 497 //++ 498 template<class T> 499 void LocalMap<T>::PixProjToAngle(double x,double y,double& theta,double& phi) const 500 // 501 // 502 // Given coordinates in pixel units in the REFERENCE PLANE, return 503 // (theta, phi) in "absolute" referential theta=pi/2 ,phi=0. 504 //-- 505 { 506 theta= Pi*0.5-atan(2.*y*tgAngleY_/(double)nSzY_); 507 phi= atan2(2.*x*tgAngleX_,(double)nSzX_); 508 if(phi < 0.) phi += DeuxPi; 509 } 510 511 //++ 512 template<class T> 513 void LocalMap<T>::AngleProjToPix(double theta,double phi,double& x,double& y) const 514 // 515 // Given coordinates (theta, phi) in "absolute" referential 516 // theta=pi/2 ,phi=0 return pixel indices (i,j) in the REFERENCE PLANE. 517 //-- 518 { 519 if(phi >= Pi) phi-= DeuxPi; 520 // y=0.5*mSzY_*cot(theta)/tgAngleY_; $CHECK-REZA-04/99$ 521 y= 0.5*nSzY_/tan(theta)/tgAngleY_; // ? cot = 1/tan ? 522 x= 0.5*nSzX_*tan(phi)/tgAngleX_; 523 } 524 525 template<class T> 526 void LocalMap<T>::print(ostream& os) const 527 { 528 os<<" SzX= "<<nSzX_<<", SzY= "<<nSzY_<<", NPix= "<<nPix_<<endl; 529 if(LocalMap_isDone()) 530 { 531 os<<" theta0= "<<theta0_<<", phi0= "<<phi0_<<", angle= "<<angle_<<endl; 532 os<<" x0= "<<x0_<<", y0= "<<y0_<<endl; 533 os<<" cos= "<<cos_angle_<<", & sin= "<<sin_angle_<<endl; 534 os<<" angleX= "<<angleX_<<", angleY= "<<angleY_<<endl; 535 os<<" tg(angleX)= "<<tgAngleX_<<", tg(angleY)= "<<tgAngleY_<<endl; 536 } 537 538 os << " contenu de pixels : "; 539 for(int i=0; i < nPix_; i++) 540 { 541 if(i%5 == 0) os << endl; 542 os << pixels_(i) <<", "; 543 } 544 os << endl; 545 } 546 //++ 547 // Titre class FIO_LocalMap 548 // Delegated objects for persitance management 549 //-- 550 551 //******************************************************************* 552 // class FIO_LocalMap<T> 553 // Les objets delegues pour la gestion de persistance 554 //******************************************************************* 555 556 //++ 557 template <class T> 558 FIO_LocalMap<T>::FIO_LocalMap() 559 // 560 //-- 561 { 562 dobj= new LocalMap<T>; 563 ownobj= true; 564 } 565 //++ 566 template <class T> 567 FIO_LocalMap<T>::FIO_LocalMap(string const& filename) 568 // 569 //-- 570 { 571 dobj= new LocalMap<T>; 572 dobj->DataBlock().SetTemp(true); 573 ownobj= true; 574 Read(filename); 575 } 576 577 //++ 578 template <class T> 579 FIO_LocalMap<T>::FIO_LocalMap(const LocalMap<T>& obj) 580 // 581 //-- 582 { 583 dobj= new LocalMap<T>(obj, true); 584 dobj->DataBlock().SetTemp(true); 585 ownobj= true; 586 } 587 588 template <class T> 589 FIO_LocalMap<T>::FIO_LocalMap(LocalMap<T>* obj) 590 { 591 dobj= obj; 592 ownobj= false; 593 } 594 595 //++ 596 template <class T> 597 FIO_LocalMap<T>::~FIO_LocalMap() 598 // 599 //-- 600 { 601 if (ownobj && dobj) delete dobj; 602 } 603 //++ 604 template <class T> 605 AnyDataObj* FIO_LocalMap<T>::DataObj() 606 // 607 //-- 608 { 609 return(dobj); 610 } 611 612 //++ 613 template <class T> 614 void FIO_LocalMap<T>::ReadSelf(PInPersist& is) 615 // 616 //-- 617 { 618 619 if(dobj == NULL) 620 { 621 dobj= new LocalMap<T>; 622 dobj->DataBlock().SetTemp(true); 623 ownobj= true; 624 } 625 626 // Pour savoir s'il y avait un DVList Info associe 627 char strg[256]; 628 is.GetLine(strg, 255); 629 bool hadinfo= false; 630 if(strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0) hadinfo= true; 631 if(hadinfo) 632 { // Lecture eventuelle du DVList Info 633 is >> dobj->Info(); 634 } 635 636 int_4 nSzX; 637 is.GetI4(nSzX); 638 dobj->setSize_x(nSzX); 639 640 int_4 nSzY; 641 is.GetI4(nSzY); 642 dobj->setSize_y(nSzY); 643 644 int_4 nPix; 645 is.GetI4(nPix); 646 dobj->setNbPixels(nPix); 647 648 string ss("local mapping is done"); 649 string sso; 650 is.GetStr(sso); 651 if(sso == ss) 652 { 653 cout<<" ReadSelf:: local mapping"<<endl; 654 int_4 x0, y0; 655 double theta, phi, angle; 656 is.GetI4(x0); 657 is.GetI4(y0); 658 is.GetR8(theta); 659 is.GetR8(phi); 660 is.GetR8(angle); 661 dobj->SetOrigin(theta, phi, x0, y0, angle); 662 663 double angleX, angleY; 664 is.GetR8(angleX); 665 is.GetR8(angleY); 666 dobj->SetSize(angleX, angleY); 667 } 668 669 // On lit le DataBlock; 670 is >> dobj->DataBlock(); 671 } 672 673 //++ 674 template <class T> 675 void FIO_LocalMap<T>::WriteSelf(POutPersist& os) const 676 // 677 //-- 678 { 679 if(dobj == NULL) 680 { 681 cout << " FIO_LocalMap::WriteSelf:: dobj= null " << endl; 682 return; 683 } 684 685 char strg[256]; 686 int_4 nSzX= dobj->Size_x(); 687 int_4 nSzY= dobj->Size_y(); 688 int_4 nPix= dobj->NbPixels(); 689 690 if(dobj->ptrInfo()) 691 { 692 sprintf(strg,"LocalMap: NPixX=%6d NPixY=%9d HasInfo",nSzX,nSzY); 693 os.PutLine(strg); 694 os << dobj->Info(); 695 } 696 else 697 { 698 sprintf(strg,"LocalMap: NPixX=%6d NPixY=%9d ",nSzX,nSzY); 699 os.PutLine(strg); 700 } 701 702 os.PutI4(nSzX); 703 os.PutI4(nSzY); 704 os.PutI4(nPix); 705 706 if(dobj->LocalMap_isDone()) 707 { 708 string ss("local mapping is done"); 709 os.PutStr(ss); 710 int_4 x0, y0; 711 double theta, phi, angle; 712 dobj->Origin(theta, phi, x0, y0, angle); 713 os.PutI4(x0); 714 os.PutI4(y0); 715 os.PutR8(theta); 716 os.PutR8(phi); 717 os.PutR8(angle); 718 719 double angleX, angleY; 720 dobj->Aperture(angleX, angleY); 721 os.PutR8(angleX); 722 os.PutR8(angleY); 723 } 724 else 725 { 726 string ss("no local mapping"); 727 os.PutStr(ss); 728 } 729 730 // On ecrit le dataBlock 731 os << dobj->DataBlock(); 732 } 8 733 9 734 #ifdef __CXX_PRAGMA_TEMPLATES__ … … 27 752 template class FIO_LocalMap< complex<double> >; 28 753 #endif 29 30 //*****************************************************************************31 //++32 // Class LocalMap33 //34 // include localmap.h nbmath.h35 //36 // A local map of a region of the sky, in cartesian coordinates.37 // It has an origin in (theta0, phi0), mapped to pixel(x0, y0)38 // (x0, y0 might be outside of this local map)39 // default value of (x0, y0) is middle of the map, center of40 // pixel(nx/2, ny/2)41 //42 // A local map is a 2 dimensional array, with i as column index and j43 // as row index. The map is supposed to lie on a plan tangent to the44 // celestial sphere in a point whose coordinates are (x0,y0) on the local45 // map and (theta0, phi0) on the sphere. The range of the map is defined46 // by two values of angles covered respectively by all the pixels in47 // x direction and all the pixels in y direction (SetSize()).48 //49 // A "reference plane" is considered : this plane is tangent to the50 // celestial sphere in a point with angles theta=Pi/2 and phi=0. This51 // point is the origine of coordinates is of the reference plane. The52 // x-axis is the tangent parallel to the equatorial line and oriented53 // toward the increasing phi's ; the y-axis is parallel to the meridian54 // line and oriented toward the north pole.55 //56 // Internally, a map is first defined within this reference plane and57 // tranported until the point (theta0, phi0) in such a way that both58 // axes are kept parallel to meridian and parallel lines of the sphere.59 // The user can define its own map with axes rotated with respect to60 // reference axes (this rotation is characterized by angle between61 // the local parallel line and the wanted x-axis-- see method62 // SetOrigin(...))63 //64 //65 66 //67 //--68 //++69 //70 // Links Parents71 //72 // PixelMap73 //74 //--75 //++76 //77 // Links Childs78 //79 //80 //--81 //++82 // Titre Constructors83 //--84 //++85 template<class T>86 LocalMap<T>::LocalMap()87 //88 //89 //--90 {91 InitNul();92 pixels_.Reset();93 }94 95 //++96 template<class T>97 LocalMap<T>::LocalMap(int nx, int ny) : nSzX_(nx), nSzY_(ny)98 //99 //100 //--101 {102 InitNul();103 nPix_= nx*ny;104 pixels_.ReSize(nPix_);105 pixels_.Reset();106 }107 108 //++109 template<class T>110 LocalMap<T>::LocalMap(const LocalMap<T>& lm, bool share)111 : pixels_(lm.pixels_, share)112 113 //114 // copy constructor115 //--116 {117 cout<<" LocalMap:: Appel du constructeur de recopie " << endl;118 119 if(lm.mInfo_) mInfo_= new DVList(*lm.mInfo_);120 nSzX_= lm.nSzX_;121 nSzY_= lm.nSzY_;122 nPix_= lm.nPix_;123 originFlag_= lm.originFlag_;124 extensFlag_= lm.extensFlag_;125 x0_= lm.x0_;126 y0_= lm.y0_;127 theta0_= lm.theta0_;128 phi0_= lm.phi0_;129 angle_= lm.angle_;130 cos_angle_= lm.cos_angle_;131 sin_angle_= lm.sin_angle_;132 angleX_= lm.angleX_;133 angleY_= lm.angleY_;134 tgAngleX_= lm.tgAngleX_;135 tgAngleY_= lm.tgAngleY_;136 }137 138 //++139 // Titre Destructor140 //--141 //++142 template<class T>143 LocalMap<T>::~LocalMap()144 //145 //--146 {147 InitNul();148 }149 150 151 152 //++153 // Titre Public Methods154 //--155 156 //++157 template<class T>158 void LocalMap<T>::ReSize(int nx, int ny)159 //160 // Resize storage area for pixels161 //--162 {163 InitNul();164 nSzX_ = nx;165 nSzY_ = ny;166 nPix_= nx*ny;167 pixels_.ReSize(nPix_);168 pixels_.Reset();169 }170 171 //++172 template<class T>173 int LocalMap<T>::NbPixels() const174 //175 // Return number of pixels176 //--177 {178 return(nPix_);179 }180 181 //++182 template<class T>183 T& LocalMap<T>::PixVal(int k)184 //185 // Return value of pixel with index k186 //--187 {188 if((k < 0) || (k >= nPix_))189 {190 cout << " LocalMap::PIxVal : exceptions a mettre en place" <<endl;191 // THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));192 THROW(rangeCheckErr);193 //throw "LocalMap::PIxVal Pixel index out of range";194 }195 return(pixels_(k));196 }197 198 //++199 200 template<class T>201 T const& LocalMap<T>::PixVal(int k) const202 //203 // const version of previous method204 //--205 {206 if((k < 0) || (k >= nPix_))207 {208 cout << " LocalMap::PIxVal : exceptions a mettre en place" <<endl;209 // THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));210 211 throw "LocalMap::PIxVal Pixel index out of range";212 }213 return *(pixels_.Data()+k);214 }215 216 template<class T>217 bool LocalMap<T>::ContainsSph(double theta, double phi) const218 {219 // $CHECK$ Reza 11/11/99 - A modifier220 return(true);221 }222 223 //++224 template<class T>225 int LocalMap<T>::PixIndexSph(double theta,double phi) const226 //227 // Return index of the pixel with spherical coordinates (theta,phi)228 //--229 {230 int i,j;231 if(!(originFlag_) || !(extensFlag_))232 {233 cout << " LocalMap: correspondance carte-sphere non etablie" << endl;234 exit(0);235 }236 237 // theta et phi en coordonnees relatives (on se ramene a une situation par rapport au plan de reference)238 double theta_aux= theta;239 double phi_aux = phi;240 UserToReference(theta_aux, phi_aux);241 242 // coordonnees dans le plan local en unites de pixels243 double x,y;244 AngleProjToPix(theta_aux,phi_aux, x, y);245 246 double xmin= -x0_-0.5;247 double xmax= xmin+nSzX_;248 if((x > xmax) || (x < xmin)) return(-1);249 double xcurrent= xmin;250 for(i = 0; i < nSzX_; i++ )251 {252 xcurrent += 1.;253 if( x < xcurrent ) break;254 }255 double ymin= -y0_-0.5;256 double ymax= ymin+nSzY_;257 if((y > ymax) || (y < ymin)) return(-1);258 double ycurrent= ymin;259 for(j = 0; j < nSzY_; j++ )260 {261 ycurrent += 1.;262 if( y < ycurrent ) break;263 }264 return (j*nSzX_+i);265 }266 267 //++268 template<class T>269 void LocalMap<T>::PixThetaPhi(int k,double& theta,double& phi) const270 //271 // Return (theta, phi) coordinates of pixel with index k272 //--273 {274 if(!(originFlag_) || !(extensFlag_))275 {276 cout << " LocalMap: correspondance carte-sphere non etablie" << endl;277 exit(0);278 }279 280 int i,j;281 Getij(k,i,j);282 283 double X= double(i-x0_);284 double Y= double(j-y0_);285 // situation de ce pixel dans le plan de reference286 double x= X*cos_angle_-Y*sin_angle_;287 double y= X*sin_angle_+Y* cos_angle_;288 // projection sur la sphere289 PixProjToAngle(x, y, theta, phi);290 // retour au plan utilisateur291 ReferenceToUser(theta, phi);292 }293 294 template <class T>295 T LocalMap<T>::SetPixels(T v)296 {297 pixels_.Reset(v);298 return(v);299 }300 301 //++302 template<class T>303 double LocalMap<T>::PixSolAngle(int k) const304 //305 // Pixel Solid angle (steradians)306 // All the pixels have not necessarly the same size in (theta, phi)307 // because of the projection scheme which is not yet fixed.308 //--309 {310 int i,j;311 Getij(k,i,j);312 double X= double(i-x0_);313 double Y= double(j-y0_);314 double XR= X+double(i)*0.5;315 double XL= X-double(i)*0.5;316 double YU= Y+double(j)*0.5;317 double YL= Y-double(j)*0.5;318 319 // situation dans le plan de reference320 double x0= XL*cos_angle_-YL*sin_angle_;321 double y0= XL*sin_angle_+YL*cos_angle_;322 double xa= XR*cos_angle_-YL*sin_angle_;323 double ya= XR*sin_angle_+YL*cos_angle_;324 double xb= XL*cos_angle_-YU*sin_angle_;325 double yb= XL*sin_angle_+YU*cos_angle_;326 327 // projection sur la sphere328 double thet0,phi0,theta,phia,thetb,phib;329 PixProjToAngle(x0, y0, thet0, phi0);330 PixProjToAngle(xa, ya, theta, phia);331 PixProjToAngle(xb, yb, thetb, phib);332 333 // angle solide334 double sol= fabs((xa-x0)*(yb-y0)-(xb-x0)*(ya-y0));335 return sol;336 }337 338 //++339 template<class T>340 void LocalMap<T>::SetOrigin(double theta0,double phi0,double angle)341 //342 // set the referential of the map (angles in degrees)343 // (default x0=siz_x/2, y0=siz_y/2)344 //--345 {346 theta0_= theta0;347 phi0_ = phi0;348 angle_ = angle;349 x0_= nSzX_/2;350 y0_= nSzY_/2;351 cos_angle_= cos(angle*Pi/180.);352 sin_angle_= sin(angle*Pi/180.);353 originFlag_= true;354 cout << " LocalMap:: set origin 1 done" << endl;355 }356 357 //++358 template<class T>359 void LocalMap<T>::SetOrigin(double theta0,double phi0,int x0,int y0,double angle)360 //361 // set the referential of the map (angles in degrees)362 //--363 {364 theta0_= theta0;365 phi0_ = phi0;366 angle_ = angle;367 x0_= x0;368 y0_= y0;369 cos_angle_= cos(angle*Pi/180.);370 sin_angle_= sin(angle*Pi/180.);371 originFlag_= true;372 cout << " LocalMap:: set origin 2 done" << endl;373 }374 375 //++376 template<class T>377 void LocalMap<T>::SetSize(double angleX,double angleY)378 //379 // angle range of tthe map (angles in degrees)380 //--381 {382 angleX_= angleX;383 angleY_= angleY;384 385 // tangente de la moitie de l'ouverture angulaire totale386 tgAngleX_= tan(0.5*angleX_*Pi/180.);387 tgAngleY_= tan(0.5*angleY_*Pi/180.);388 389 extensFlag_= true;390 cout << " LocalMap:: set extension done" << endl;391 }392 393 //++394 template<class T>395 void LocalMap<T>::Project(SphericalMap<T>& sphere) const396 //397 // Projection to a spherical map398 //--399 {400 for(int m = 0; m < nPix_; m++)401 {402 double theta,phi;403 PixThetaPhi(m,theta,phi);404 sphere(theta,phi)= pixels_(m);405 // cout << "theta " << theta << " phi " << phi << " valeur " << sphere(theta,phi)<< endl;406 }407 }408 // Titre Private Methods409 //++410 template<class T>411 void LocalMap<T>::InitNul()412 //413 // set some attributes to zero414 //--415 {416 originFlag_= false;417 extensFlag_= false;418 cos_angle_= 1.0;419 sin_angle_= 0.0;420 // pixels_.Reset(); Pas de reset par InitNul (en cas de share) - Reza 20/11/99 $CHECK$421 }422 423 //++424 template<class T>425 void LocalMap<T>::Getij(int k,int& i,int& j) const426 //427 // Return 2 indices corresponding to the pixel number k428 //--429 {430 i= (k+1)%nSzX_-1;431 if(i == -1) i= nSzX_-1;432 j= (k-i+2)/nSzX_;433 }434 435 //++436 template<class T>437 void LocalMap<T>::ReferenceToUser(double& theta,double& phi) const438 //439 // Transform a pair of coordinates (theta, phi) given in440 // reference coordinates into map coordinates441 //--442 {443 if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi)444 {445 throw "LocalMap::ReferenceToUser (theta,phi) out of range";446 }447 448 theta= theta0_*Pi/180.+theta-Pi*0.5;449 if(theta < 0.)450 {451 theta= -theta;452 phi += Pi;453 }454 else455 {456 if(theta > Pi)457 {458 theta= 2.*Pi-theta;459 phi += Pi;460 }461 }462 463 phi= phi0_*Pi/180.+phi;464 while(phi >= 2.*Pi) phi-= 2.*Pi;465 466 if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi)467 {468 cout << " LocalMap::ReferenceToUser : erreur bizarre dans le transfert a la carte utilisateur " << endl;469 cout << " theta= " << theta << " phi= " << phi << endl;470 exit(0);471 }472 }473 474 //++475 template<class T>476 void LocalMap<T>::UserToReference(double& theta,double& phi) const477 //478 // Transform a pair of coordinates (theta, phi) given in479 // map coordinates into reference coordinates480 //--481 {482 if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi)483 {484 cout<<" LocalMap::UserToReference: exceptions a mettre en place" <<endl;485 // THROW(out_of_range("LocalMap::PIxVal Pixel index out of range"));486 throw "LocalMap::UserToReference (theta,phi) out of range";487 }488 489 double phi1= phi-phi0_*Pi/180.;490 if(phi1 < 0.) phi1+= 2.*Pi;491 492 double theta1= theta-theta0_*Pi/180.+Pi*0.5;493 if(theta1 < 0.)494 {495 theta= -theta1;496 phi1+= Pi;497 }498 else499 {500 if(theta1 > Pi)501 {502 theta= 2.*Pi-theta1;503 phi1+= Pi;504 }505 }506 507 while(phi1 >= 2.*Pi) phi1-= 2.*Pi;508 phi= phi1;509 if(theta > Pi || theta < 0. || phi < 0. || phi >= 2*Pi)510 {511 cout << " LocalMap::UserToReference : erreur bizarre dans le transfert a la carte de reference " << endl;512 cout << " theta= " << theta << " phi= " << phi << endl;513 exit(0);514 }515 }516 517 //++518 template<class T>519 void LocalMap<T>::PixProjToAngle(double x,double y,double& theta,double& phi) const520 //521 //522 // Given coordinates in pixel units in the REFERENCE PLANE, return523 // (theta, phi) in "absolute" referential theta=pi/2 ,phi=0.524 //--525 {526 theta= Pi*0.5-atan(2.*y*tgAngleY_/(double)nSzY_);527 phi= atan2(2.*x*tgAngleX_,(double)nSzX_);528 if(phi < 0.) phi += DeuxPi;529 }530 531 //++532 template<class T>533 void LocalMap<T>::AngleProjToPix(double theta,double phi,double& x,double& y) const534 //535 // Given coordinates (theta, phi) in "absolute" referential536 // theta=pi/2 ,phi=0 return pixel indices (i,j) in the REFERENCE PLANE.537 //--538 {539 if(phi >= Pi) phi-= DeuxPi;540 // y=0.5*mSzY_*cot(theta)/tgAngleY_; $CHECK-REZA-04/99$541 y= 0.5*nSzY_/tan(theta)/tgAngleY_; // ? cot = 1/tan ?542 x= 0.5*nSzX_*tan(phi)/tgAngleX_;543 }544 545 template<class T>546 void LocalMap<T>::print(ostream& os) const547 {548 os<<" SzX= "<<nSzX_<<", SzY= "<<nSzY_<<", NPix= "<<nPix_<<endl;549 if(LocalMap_isDone())550 {551 os<<" theta0= "<<theta0_<<", phi0= "<<phi0_<<", angle= "<<angle_<<endl;552 os<<" x0= "<<x0_<<", y0= "<<y0_<<endl;553 os<<" cos= "<<cos_angle_<<", & sin= "<<sin_angle_<<endl;554 os<<" angleX= "<<angleX_<<", angleY= "<<angleY_<<endl;555 os<<" tg(angleX)= "<<tgAngleX_<<", tg(angleY)= "<<tgAngleY_<<endl;556 }557 558 os << " contenu de pixels : ";559 for(int i=0; i < nPix_; i++)560 {561 if(i%5 == 0) os << endl;562 os << pixels_(i) <<", ";563 }564 os << endl;565 }566 //++567 // Titre class FIO_LocalMap568 // Delegated objects for persitance management569 //--570 571 //*******************************************************************572 // class FIO_LocalMap<T>573 // Les objets delegues pour la gestion de persistance574 //*******************************************************************575 576 //++577 template <class T>578 FIO_LocalMap<T>::FIO_LocalMap()579 //580 //--581 {582 dobj= new LocalMap<T>;583 ownobj= true;584 }585 //++586 template <class T>587 FIO_LocalMap<T>::FIO_LocalMap(string const& filename)588 //589 //--590 {591 dobj= new LocalMap<T>;592 dobj->DataBlock().SetTemp(true);593 ownobj= true;594 Read(filename);595 }596 597 //++598 template <class T>599 FIO_LocalMap<T>::FIO_LocalMap(const LocalMap<T>& obj)600 //601 //--602 {603 dobj= new LocalMap<T>(obj, true);604 dobj->DataBlock().SetTemp(true);605 ownobj= true;606 }607 608 template <class T>609 FIO_LocalMap<T>::FIO_LocalMap(LocalMap<T>* obj)610 {611 dobj= obj;612 ownobj= false;613 }614 615 //++616 template <class T>617 FIO_LocalMap<T>::~FIO_LocalMap()618 //619 //--620 {621 if (ownobj && dobj) delete dobj;622 }623 //++624 template <class T>625 AnyDataObj* FIO_LocalMap<T>::DataObj()626 //627 //--628 {629 return(dobj);630 }631 632 //++633 template <class T>634 void FIO_LocalMap<T>::ReadSelf(PInPersist& is)635 //636 //--637 {638 639 if(dobj == NULL)640 {641 dobj= new LocalMap<T>;642 dobj->DataBlock().SetTemp(true);643 ownobj= true;644 }645 646 // Pour savoir s'il y avait un DVList Info associe647 char strg[256];648 is.GetLine(strg, 255);649 bool hadinfo= false;650 if(strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0) hadinfo= true;651 if(hadinfo)652 { // Lecture eventuelle du DVList Info653 is >> dobj->Info();654 }655 656 int nSzX;657 is.GetI4(nSzX);658 dobj->setSize_x(nSzX);659 660 int nSzY;661 is.GetI4(nSzY);662 dobj->setSize_y(nSzY);663 664 int nPix;665 is.GetI4(nPix);666 dobj->setNbPixels(nPix);667 668 string ss("local mapping is done");669 string sso;670 is.GetStr(sso);671 if(sso == ss)672 {673 cout<<" ReadSelf:: local mapping"<<endl;674 int x0, y0;675 double theta, phi, angle;676 is.GetI4(x0);677 is.GetI4(y0);678 is.GetR8(theta);679 is.GetR8(phi);680 is.GetR8(angle);681 dobj->SetOrigin(theta, phi, x0, y0, angle);682 683 double angleX, angleY;684 is.GetR8(angleX);685 is.GetR8(angleY);686 dobj->SetSize(angleX, angleY);687 }688 689 // On lit le DataBlock;690 is >> dobj->DataBlock();691 }692 693 //++694 template <class T>695 void FIO_LocalMap<T>::WriteSelf(POutPersist& os) const696 //697 //--698 {699 if(dobj == NULL)700 {701 cout << " FIO_LocalMap::WriteSelf:: dobj= null " << endl;702 return;703 }704 705 char strg[256];706 int nSzX= dobj->Size_x();707 int nSzY= dobj->Size_y();708 int nPix= dobj->NbPixels();709 710 if(dobj->ptrInfo())711 {712 sprintf(strg,"LocalMap: NPixX=%6d NPixY=%9d HasInfo",nSzX,nSzY);713 os.PutLine(strg);714 os << dobj->Info();715 }716 else717 {718 sprintf(strg,"LocalMap: NPixX=%6d NPixY=%9d ",nSzX,nSzY);719 os.PutLine(strg);720 }721 722 os.PutI4(nSzX);723 os.PutI4(nSzY);724 os.PutI4(nPix);725 726 if(dobj->LocalMap_isDone())727 {728 string ss("local mapping is done");729 os.PutStr(ss);730 int x0, y0;731 double theta, phi, angle;732 dobj->Origin(theta, phi, x0, y0, angle);733 os.PutI4(x0);734 os.PutI4(y0);735 os.PutR8(theta);736 os.PutR8(phi);737 os.PutR8(angle);738 739 double angleX, angleY;740 dobj->Aperture(angleX, angleY);741 os.PutR8(angleX);742 os.PutR8(angleY);743 }744 else745 {746 string ss("no local mapping");747 os.PutStr(ss);748 }749 750 // On ecrit le dataBlock751 os << dobj->DataBlock();752 }753
Note:
See TracChangeset
for help on using the changeset viewer.