source: Sophya/trunk/SophyaLib/SkyMap/sphereecp.cc@ 3636

Last change on this file since 3636 was 3515, checked in by ansari, 17 years ago

correction bug debordement memoire ds SphereECP lors de synfast par sphericaltransform, Reza 29/08/2008
File size: 13.4 KB
Line 
1#include "sopnamsp.h"
2#include "sphereecp.h"
3#include <math.h>
4
5/*!
6 \class SOPHYA::SphereECP
7 \ingroup SkyMap
8 \brief Spherical maps in Equatorial Cylindrical Projection
9
10 Class implementing spherical maps, with a pixelisation
11 corresponding to the an Equatorial Cylindrical Projection (ECP).
12 Pixel size varie from equator to poles. The sphere is divided
13 into a number of slices along the parallels. Each slice is subdivided
14 into the same number of pixels.
15 This class provides the possibility to have partial coverage.
16 Can be used for the following types: <br>
17 <tt> int_4 , r_4 , r_8 , complex<r_4> , complex<r_8> </tt>
18*/
19
20namespace SOPHYA {
21
22//! Default constructor
23template <class T>
24SphereECP<T>::SphereECP()
25 : _pixels()
26{
27 _partial = false;
28 _theta1 = 0.;
29 _theta2 = M_PI;
30 _phi1 = 0.;
31 _phi2 = 2*M_PI;
32 _outofmapidx = -99;
33 _outofmapnphi = 0;
34 _outofmapntet = 0;
35 _outofmappix = 0;
36 _outofmapval = 0;
37 _dtheta = 0.;
38 _dphi = 0.;
39}
40
41/*!
42 \brief Constructor with the pixelisation parameter m
43
44 Complete coverage.
45 The sphere is divided into \b m slices in theta (0..pi). each slice is divided
46 into \b 2m pixels (along phi)
47*/
48template <class T>
49SphereECP<T>::SphereECP(int m)
50 : _pixels(2*m,m,0,0,0,true)
51{
52 _partial = false;
53 _theta1 = 0.;
54 _theta2 = M_PI;
55 _phi1 = 0.;
56 _phi2 = 2*M_PI;
57 _outofmapidx = -99;
58 _outofmapnphi = 0;
59 _outofmapntet = 0;
60 _outofmappix = 0;
61 _outofmapval = 0;
62 _dtheta = _dphi = M_PI/m;
63}
64
65/*!
66 Complete coverage.
67 \b ntet slices in theta (0..pi) and \b nphi pixels (along phi) for
68 each slice in theta
69*/
70template <class T>
71SphereECP<T>::SphereECP(int ntet, int nphi)
72 : _pixels(nphi,ntet,0,0,0,true)
73{
74 _partial = false;
75 _theta1 = 0.;
76 _theta2 = M_PI;
77 _phi1 = 0.;
78 _phi2 = 2*M_PI;
79 _outofmapidx = -99;
80 _outofmapnphi = 0;
81 _outofmapntet = 0;
82 _outofmappix = 0;
83 _outofmapval = 0;
84 _dtheta = M_PI/ntet;
85 _dphi = 2*M_PI/nphi;
86}
87
88/*!
89 Partial coverage :
90 \b ntet slices in theta , in the range (tet1 .. tet2) and \b nphi pixels (along phi) for
91 each slice in theta in the range (phi1 .. phi2)
92*/
93template <class T>
94SphereECP<T>::SphereECP(r_8 tet1, r_8 tet2, int ntet, r_8 phi1, r_8 phi2, int nphi)
95 : _pixels(nphi,ntet,0,0,0,true)
96{
97 if( (tet1 > M_PI) || (tet1 < 0.) || (tet2 > M_PI) || (tet2 < 0.) ||
98 (phi1 < 0.) || (phi1 > 2*M_PI) || (phi1 < 0.) || (phi1 > 2*M_PI) ||
99 (tet2 <= tet1) || (phi2 <= phi1) )
100 throw ParmError("SphereECP::SphereECP() Bad range for theta/phi limits");
101
102 _partial = true;
103 _theta1 = tet1;
104 _theta2 = tet2;
105 _phi1 = phi1;
106 _phi2 = phi2;
107 _outofmapidx = _pixels.Size()+659;
108 // Reza, 23 sep 2004
109 // pour des cartes partielles, et des pixels sur la sphere, mais
110 // en dehors de la couverture, on code le numero de tranche en theta,
111 // compte a partir de theta=0,phi=0 + _outofmapidx
112 // theta -> idx = _outofmapidx + theta/_dtheta*_outofmapnphi + _outofmapntet
113 _outofmapntet = 0;
114 _outofmappix = 0;
115 _outofmapval = 0;
116 _dtheta = (_theta2-_theta1)/ntet;
117 _dphi = (_phi2-_phi1)/nphi;
118
119 _outofmapntet = M_PI/_dtheta;
120 _outofmapnphi = 2*M_PI/_dphi;
121}
122
123//! Copy constructor - shares the pixel data if \b share=true
124template <class T>
125SphereECP<T>::SphereECP(const SphereECP<T>& a, bool share)
126 : _pixels(a._pixels, share)
127{
128 _partial = a._partial;
129 _theta1 = a._theta1;
130 _theta2 = a._theta2;
131 _phi1 = a._phi1;
132 _phi2 = a._phi2;
133 _outofmapidx = a._outofmapidx;
134 _outofmapnphi = a._outofmapnphi;
135 _outofmapntet = a._outofmapntet;
136 _outofmappix = a._outofmappix;
137 _outofmapval = a._outofmapval;
138 _dtheta = a._dtheta;
139 _dphi = a._dphi;
140}
141
142//! Copy constructor - shares the pixel data
143template <class T>
144SphereECP<T>::SphereECP(const SphereECP<T>& a)
145 : _pixels(a._pixels)
146{
147 _partial = a._partial;
148 _theta1 = a._theta1;
149 _theta2 = a._theta2;
150 _phi1 = a._phi1;
151 _phi2 = a._phi2;
152 _outofmapidx = a._outofmapidx;
153 _outofmapnphi = a._outofmapnphi;
154 _outofmapntet = a._outofmapntet;
155 _outofmappix = a._outofmappix;
156 _outofmapval = a._outofmapval;
157 _dtheta = a._dtheta;
158 _dphi = a._dphi;
159}
160
161template <class T>
162SphereECP<T>::~SphereECP()
163{
164}
165
166/*!
167 \brief Extract a partial map from a full ECP skymap
168 Theta and Phi limits are adjusted automatically to be correspond to the source map pixel boundaries
169*/
170template <class T>
171SphereECP<T> SphereECP<T>::ExtractPartial(r_8 tet1, r_8 tet2, r_8 phi1, r_8 phi2)
172{
173 if (IsPartial())
174 throw ParmError("SphereECP::ExtractPartial() source map NOT a full sky map");
175 if( (tet1 > M_PI) || (tet1 < 0.) || (tet2 > M_PI) || (tet2 < 0.) ||
176 (phi1 < 0.) || (phi1 > 2*M_PI) || (phi1 < 0.) || (phi1 > 2*M_PI) ||
177 (tet2 <= tet1) || (phi2 <= phi1) )
178 throw ParmError("SphereECP::ExtractPartial() Bad range for theta/phi limits");
179
180// correction/calcul des limites en theta et nombre de rangees en theta
181 sa_size_t ntet = (sa_size_t)((tet2-tet1)/_dtheta + 0.5);
182 sa_size_t offt = (sa_size_t)(tet1/_dtheta + 1e-3);
183 tet1 = _dtheta*offt;
184 tet2 = tet1+ntet*_dtheta;
185// correction/calcul des limites en phi et nombre de colonnes en phi
186 sa_size_t nphi = (sa_size_t)((phi2-phi1)/_dphi + 0.5);
187 sa_size_t offp = (sa_size_t)(phi1/_dphi + 1e-3);
188 phi1 = _dphi*offp;
189 phi2 = phi1+nphi*_dphi;
190
191// On cree une carte partielle avec les bons parametres
192 SphereECP<T> pmap(tet1, tet2, ntet, phi1, phi2, nphi);
193// On recopie les valeurs de pixels de la zone correspondante
194 pmap.GetPixelArray() = _pixels.SubArray(Range(offp, 0, nphi, 0), Range(offt, 0, ntet, 0),
195 Range::first(), Range::first(), Range::first());
196 return pmap;
197}
198
199template <class T>
200string SphereECP<T>::TypeOfMap() const
201{
202 return string("ECP");
203}
204
205template <class T>
206void SphereECP<T>::SetTemp(bool temp) const
207{
208 _pixels.SetTemp(temp);
209}
210
211template <class T>
212int_4 SphereECP<T>::NbPixels() const
213{
214 return _pixels.Size();
215}
216
217template <class T>
218T& SphereECP<T>::PixVal(int_4 k)
219{
220 // On pourrait etre plus strict en demandant aussi k<_outofmapidx+_outofmapnphi*_outofmapntet (Reza 23/09/04)
221 if ((_partial) && (k>=_outofmapidx)) return _outofmappix;
222 if((k < 0) || (k >= _pixels.Size()) )
223 throw RangeCheckError("SphereECP::PIxVal() Pixel index out of range");
224 return _pixels.DataBlock()(k);
225}
226
227template <class T>
228T const& SphereECP<T>::PixVal(int_4 k) const
229{
230 // On pourrait etre plus strict en demandant aussi k<_outofmapidx+_outofmapnphi*_outofmapntet (Reza 23/09/04)
231 if ((_partial) && (k>=_outofmapidx)) return _outofmapval;
232 if((k < 0) || (k >= _pixels.Size()) )
233 throw RangeCheckError("SphereECP::PIxVal() Pixel index out of range");
234 // return _pixels.DataBlock()(k);
235 return _pixels.Data()[k];
236}
237
238template <class T>
239bool SphereECP<T>::ContainsSph(double theta, double phi) const
240{
241 if (_partial) {
242 if ( (theta >= _theta1) && (theta <= _theta2) &&
243 (phi >= _phi1) && (phi <= _phi2) ) return true;
244 else return false;
245 }
246 else return false;
247}
248
249template <class T>
250int_4 SphereECP<T>::PixIndexSph(double theta, double phi) const
251{
252 if((theta > M_PI) || (theta < 0.)) return(-1);
253 if((phi < 0.) || (phi > 2*M_PI)) return(-1);
254 if (_partial) {
255 if ( (theta < _theta1) || (theta > _theta2) ||
256 (phi < _phi1) || (phi > _phi2) ) return _outofmapidx+theta/_dtheta*_outofmapnphi+phi/_dphi;
257 }
258 int_4 it = (theta-_theta1)/_dtheta;
259 int_4 jp = (phi-_phi1)/_dphi;
260 return (it*_pixels.SizeX()+jp);
261}
262
263template <class T>
264void SphereECP<T>::PixThetaPhi(int_4 k, double& theta, double& phi) const
265{
266 if ((_partial) && (k>=_outofmapidx)) {
267 theta = (((k-_outofmapidx) / _outofmapnphi)+0.5)*_dtheta;
268 phi = (((k-_outofmapidx) % _outofmapnphi)+0.5)*_dphi;;
269 return;
270 }
271 if((k < 0) || (k >= _pixels.Size()) )
272 throw RangeCheckError("SphereECP::PixThetaPhi() Pixel index out of range");
273
274 int_4 it = k / _pixels.SizeX();
275 int_4 jp = k % _pixels.SizeX();
276 theta = _theta1+(it+0.5)*_dtheta;
277 phi = _phi1+(jp+0.5)*_dphi;
278 return;
279}
280
281template <class T>
282T SphereECP<T>::SetPixels(T v)
283{
284 _pixels = v;
285 return v;
286}
287
288template <class T>
289double SphereECP<T>::PixSolAngle(int_4 k) const
290{
291 if ((_partial) && (k>=_outofmapidx)) {
292 double theta = (((k-_outofmapidx) / _outofmapnphi)+0.5)*_dtheta;
293 return (_dtheta*_dphi*sin(theta));
294 }
295 if((k < 0) || (k >= _pixels.Size()) )
296 throw RangeCheckError("SphereECP::PixSolAngle(int_4 k) Pixel index out of range");
297
298 int_4 it = k / _pixels.SizeX();
299 double theta = _theta1+it*_dtheta;
300 return (_dtheta*_dphi*sin(theta));
301}
302
303template <class T>
304int_4 SphereECP<T>::SizeIndex() const
305{
306 return _pixels.SizeY();
307}
308
309template <class T>
310void SphereECP<T>::Resize(int_4 m)
311{
312 if ( (m <= 0) || ( m == _pixels.SizeY()) ) {
313 cout << " SphereECP<T>::Resize(int_4 " << m << ") m<0 ou m=NTheta - Ne fait rien " << endl;
314 return;
315 }
316 int mphi = m;
317 if (_pixels.Size() > 0) mphi = m*_pixels.SizeX()/_pixels.SizeY();
318 cout << " SphereECP<T>::Resize(" << m
319 << ") -> _pixels.Resize(" << mphi << "," << m << ")" << endl;
320 sa_size_t sz[5] = {0,0,0,0,0};
321 sz[0] = mphi; sz[1] = m;
322 _pixels.ReSize(2,sz);
323 _outofmapidx = _pixels.Size()+659;
324 _dtheta = (_theta2-_theta1)/m;
325 _dphi = (_phi2-_phi1)/mphi;
326
327 _outofmapntet = M_PI/_dtheta;
328 _outofmapnphi = 2*M_PI/_dphi;
329 return;
330}
331
332template <class T>
333uint_4 SphereECP<T>::NbThetaSlices() const
334{
335 return _pixels.SizeY();
336}
337
338template <class T>
339r_8 SphereECP<T>::ThetaOfSlice(int_4 index) const
340{
341 if( (index < 0) || (index >= _pixels.SizeY()) )
342 throw RangeCheckError("SphereECP::ThetaOfSlice() theta index out of range");
343 return (_theta1 + (index+0.5)*_dtheta);
344
345}
346
347template <class T>
348void SphereECP<T>::GetThetaSlice(int_4 index,r_8& theta,
349 TVector<r_8>& phi, TVector<T>& value) const
350{
351 if( (index < 0) || (index >= _pixels.SizeY()) )
352 throw RangeCheckError("SphereECP::GetThetaSlice() index out of range");
353
354 theta = _theta1 + (index+0.5)*_dtheta;
355 int_4 nphit = 2.*M_PI/_dphi;
356 phi.ReSize(nphit);
357 value.ReSize(nphit);
358 int_4 ioff = _phi1/_dphi;
359 int_4 i;
360 for(i=0; i<ioff; i++) {
361 phi(i) = _phi1 + (i-ioff+0.5)*_dphi;
362 value(i) = _outofmapval;
363 }
364 for(i=ioff; i<ioff+_pixels.SizeX(); i++) {
365 phi(i) = _phi1 + (i-ioff+0.5)*_dphi;
366 value(i) = _pixels(i-ioff,index,0);
367 }
368 for(i=ioff+_pixels.SizeX(); i<nphit; i++) {
369 phi(i) = _phi1 + (i-ioff+0.5)*_dphi;
370 value(i) = _outofmapval;
371 }
372 return;
373}
374
375template <class T>
376void SphereECP<T>::GetThetaSlice(int_4 index, r_8& theta, r_8& phi0,
377 TVector<int_4>& pixidx, TVector<T>& value) const
378{
379 if( (index < 0) || (index >= _pixels.SizeY()) )
380 throw RangeCheckError("SphereECP::GetThetaSlice() index out of range");
381
382 theta = _theta1 + (index+0.5)*_dtheta;
383 int_4 ioff = _phi1/_dphi;
384 phi0 = _phi1 - (ioff-0.5)*_dphi;
385
386 int_4 nphit = 2.*M_PI/_dphi;
387 pixidx.ReSize(nphit);
388 value.ReSize(nphit);
389
390 int_4 i;
391 for(i=0; i<ioff; i++) {
392 pixidx(i) = _outofmapidx+theta/_dtheta*_outofmapnphi+phi0/_dphi+i;
393 value(i) = _outofmapval;
394 }
395 for(i=ioff; i<ioff+_pixels.SizeX(); i++) {
396 pixidx(i) = index*_pixels.SizeX()+(i-ioff);
397 value(i) = _pixels(i-ioff,index,0);
398 }
399 for(i=ioff+_pixels.SizeX(); i<nphit; i++) {
400 pixidx(i) = _outofmapidx+theta/_dtheta*_outofmapnphi+phi0/_dphi+i;
401 value(i) = _outofmapval;
402 }
403 return;
404
405}
406
407template <class T>
408T* SphereECP<T>::GetThetaSliceDataPtr(int_4 index)
409{
410 if( (index < 0) || (index >= _pixels.SizeY()) )
411 throw RangeCheckError("SphereECP::GetThetaSliceDataPtr() index out of range");
412 int_4 ioff = _phi1/_dphi;
413 int_4 nphit = 2.*M_PI/_dphi;
414// Correction bug signale par cmv le 27/08/2008, debordement memoire lors
415// de synthese par SphericalTransformServer ...
416 if ( (ioff != 0) || (nphit!=_pixels.SizeX()) ) return NULL;
417 return _pixels.Data() + index*_pixels.SizeX();
418}
419
420template <class T>
421void SphereECP<T>::Show(ostream& os) const
422{
423 PixelMap<T>::Show(os);
424 if (_partial)
425 os << "SphereECP<T>::Print() - Partial ECP Map NPhi=" << _pixels.SizeX()
426 << " x NTheta= " << _pixels.SizeY() << " SolidAngle=" << PixSolAngle() << "\n"
427 << "ThetaLimits= " << _theta1 << " .. " << _theta2
428 << " PhiLimits= " << _phi1 << " .. " << _phi2 << endl;
429 else
430 os << "SphereECP<T>::Print() - Full ECP Map NPhi=" << _pixels.SizeX()
431 << " x NTheta= " << _pixels.SizeY() << " SolidAngle=" << PixSolAngle() << endl;
432}
433
434template <class T>
435void SphereECP<T>::Print(ostream& os) const
436{
437 Show(os);
438 os << _pixels;
439}
440
441template <class T>
442SphereECP<T>& SphereECP<T>::Set(const SphereECP<T>& a)
443{
444 _pixels.Set(a._pixels);
445 _partial = a._partial;
446 _theta1 = a._theta1;
447 _theta2 = a._theta2;
448 _phi1 = a._phi1;
449 _phi2 = a._phi2;
450 _outofmapidx = a._outofmapidx;
451 _outofmapnphi = a._outofmapnphi;
452 _outofmapntet = a._outofmapntet;
453 _outofmappix = a._outofmappix;
454 _outofmapval = a._outofmapval;
455 _dtheta = a._dtheta;
456 _dphi = a._dphi;
457 return *this ;
458}
459
460template <class T>
461SphereECP<T>& SphereECP<T>::SetCst(T x)
462{
463 _pixels.SetCst(x);
464 return *this ;
465}
466
467template <class T>
468SphereECP<T>& SphereECP<T>::AddCst(T x)
469{
470 _pixels += x;
471 return *this ;
472}
473
474template <class T>
475SphereECP<T>& SphereECP<T>::MulCst(T x)
476{
477 _pixels *= x;
478 return *this ;
479}
480
481
482
483}// Fin du namespace
484
485
486#ifdef __CXX_PRAGMA_TEMPLATES__
487#pragma define_template SphereECP<int_4>
488#pragma define_template SphereECP<r_4>
489#pragma define_template SphereECP<r_8>
490#pragma define_template SphereECP< complex<r_4> >
491#pragma define_template SphereECP< complex<r_8> >
492#endif
493#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
494namespace SOPHYA {
495template class SphereECP<int_4>;
496template class SphereECP<r_4>;
497template class SphereECP<r_8>;
498template class SphereECP< complex<r_4> >;
499template class SphereECP< complex<r_8> >;
500}
501#endif
502
Note: See TracBrowser for help on using the repository browser.