source: Sophya/trunk/SophyaExt/FitsIOServer/fitsspherehealpix.cc@ 3072

Last change on this file since 3072 was 3047, checked in by ansari, 19 years ago

Ajout FitsInOutFile::SkipEmptyFirstHDU() , positionnement sur HDU 2 si HDU 1 vide ds operateur >> lisant des tables + autres petites corrections , Reza 11/08/2006
File size: 8.7 KB
Line 
1#include "sopnamsp.h"
2#include "pexceptions.h"
3#include "fitsspherehealpix.h"
4#include "tarray.h"
5#include "fitstarray.h"
6
7
8
9///////////////////////////////////////////////////////////////////////
10// ----objets delegues pour la gestion de persistance I/O format fits--
11// SphereHealpix
12////////////////////////////////////////////////////////////////////
13
14/*!
15 \class SOPHYA::FITS_SphereHEALPix
16 \ingroup FitsIOServer
17 FITS format I/O handler for SOPHYA::SphereHEALPix objects.
18*/
19
20
21template <class T>
22FITS_SphereHEALPix<T>::FITS_SphereHEALPix()
23{
24 dobj_= new SphereHEALPix<T>;
25 ownobj_= true;
26}
27
28template <class T>
29FITS_SphereHEALPix<T>::FITS_SphereHEALPix(char inputfile[],int hdunum)
30{
31 dobj_= new SphereHEALPix<T>;
32 ownobj_= true;
33 Read(inputfile,hdunum);
34}
35
36
37template <class T>
38FITS_SphereHEALPix<T>::FITS_SphereHEALPix(const SphereHEALPix<T>& obj)
39{
40 dobj_= new SphereHEALPix<T>(obj);
41 ownobj_= true;
42}
43
44template <class T>
45FITS_SphereHEALPix<T>::FITS_SphereHEALPix(SphereHEALPix<T>* obj)
46{
47 dobj_= obj;
48 ownobj_= false;
49}
50
51template <class T>
52FITS_SphereHEALPix<T>::~FITS_SphereHEALPix()
53{
54 if (ownobj_ && dobj_) delete dobj_;
55}
56
57template <class T>
58AnyDataObj* FITS_SphereHEALPix<T>::DataObj()
59{
60 return(dobj_);
61}
62
63template <class T>
64void FITS_SphereHEALPix<T>::SetDataObj(AnyDataObj & o)
65{
66 SphereHEALPix<T> * po = dynamic_cast< SphereHEALPix<T> * >(&o);
67 if (po == NULL) return;
68 if (ownobj_ && dobj_) delete dobj_;
69 dobj_ = po;
70 ownobj_ = false;
71}
72
73template <class T>
74int FITS_SphereHEALPix<T>::CheckHandling(AnyDataObj & o)
75{
76 if ( typeid(o) == typeid( SphereHEALPix< T > ) ) return 2;
77 SphereHEALPix<T> * po = dynamic_cast< SphereHEALPix<T> * >(&o);
78 if (po == NULL) return 0;
79 else return 1;
80}
81
82template <class T>
83int FITS_SphereHEALPix<T>::CheckReadability(FitsInOutFile& is)
84{
85 if (is.CurrentHDUType() != BINARY_TBL ) return 0;
86 string key;
87 key = "PIXTYPE";
88 if ( is.KeyValue(key) != "HEALPIX") return 0;
89
90 bool nosk = false;
91 key = "ORDERING";
92 is.KeyValue(key, nosk);
93 bool nosk2 = false;
94 key = "NSIDE";
95 is.KeyValue(key, nosk2);
96 if (nosk || nosk2) return 0;
97
98 vector<string> colnames;
99 vector<int> coltypes;
100 vector<long> repcnt;
101 vector<long> width;
102 long ncols = is.GetColInfo(colnames, coltypes, repcnt, width);
103 if (ncols < 1) return 0;
104 T x = 0;
105 if (coltypes[0] == FitsTypes::DataType(x)) return 2 ;
106 else return 1;
107}
108
109template <class T>
110FitsHandlerInterface* FITS_SphereHEALPix<T>::Clone()
111{
112 return new FITS_SphereHEALPix<T>() ;
113}
114
115
116template <class T>
117void FITS_SphereHEALPix<T>::WriteToFits(FitsOutFile& os)
118{
119 if(dobj_ == NULL)
120 {
121 cout << " WriteToFits:: dobj_= null " << endl;
122 return;
123 }
124
125 DVList dvl( dobj_->Info() );
126
127 SphereCoordSys* cs= dobj_->GetCoordSys();
128 string description= cs->description();
129 dvl["COORDSYS"]= description;
130
131 dvl["PIXTYPE"] = "HEALPIX";
132 dvl.SetComment("PIXTYPE", "HEALPIX Pixelization");
133 dvl["ORDERING"]= dobj_->TypeOfMap();
134 dvl.SetComment("ORDERING", "Pixel ordering scheme, either RING or NESTED");
135
136 int nSide= dobj_->SizeIndex();
137 dvl["NSIDE"]= (int_4) nSide;
138 dvl.SetComment("NSIDE","Resolution parameter for HEALPIX" );
139
140 int nPix= dobj_->NbPixels();
141 dvl["FIRSTPIX"]= (int_4) 0;
142 dvl.SetComment("FIRSTPIX", "First pixel # (0 based)");
143 dvl["LASTPIX"]= (int_4) nPix-1;
144 dvl.SetComment("LASTPIX", "Last pixel # (0 based)");
145 dvl["Content"]= "SphereHEALPix";
146 dvl.SetComment("Content", "name of SOPHYA object");
147 dvl["SOPCLSNM"]= "SOPHYA::SphereHEALPix<T>";
148 dvl.SetComment("SOPCLSNM", "SOPHYA class name");
149
150 // On ecrit les dataBlocks
151 vector<string> Noms;
152 Noms.push_back(dvl.GetS("Content"));
153 // string extname("SIMULATION");
154 string extname = os.NextExtensionName();
155
156 string Type;
157 if (typeid(T) == typeid(r_8) ) Type+='D';
158 else if (typeid(T) == typeid(r_4) ) Type+='E';
159 else if (typeid(T) == typeid(int_4) ) Type+='J';
160 else
161 {
162 cout << "FITS_SphereHEALPix<T>::WriteToFits()- Unsupported data type ("
163 << typeid(T).name() << ") throwing IOExc " << endl;
164 throw IOExc("FITS_SphereHEALPix<T>::WriteToFits() unsupported data type");
165 }
166 vector<int> dummy;
167 os.makeHeaderBntblOnFits(Type, Noms, nPix, 1, &dvl, extname, dummy);
168 os.PutColToFits(0, nPix, dobj_->pixels_.Data());
169}
170
171template <class T>
172void FITS_SphereHEALPix<T>::ReadFromFits(FitsInFile& is)
173{
174 if (is.CurrentHDUType() != BINARY_TBL )
175 throw FitsIOException("FITS_SphereHEALPix<T>::ReadFromFits Not a binary table HDU");
176
177 if(dobj_ == NULL) {
178 dobj_= new SphereHEALPix<T>;
179 ownobj_= true;
180 }
181 int nbcols, nbentries;
182 //
183
184 nbcols = is.NbColsFromFits();
185 if (nbcols != 1)
186 cout << "FITS_SphereHEALPix<T>::ReadFromFits()/Warning Fits::NbCols="
187 << nbcols << " not equal 1" << endl;
188
189 DVList dvl=is.DVListFromFits();
190 nbentries = is.NentriesFromFits(0);
191 int lastpix=dvl.GetI("LASTPIX");
192 if (lastpix>0)
193 {
194 if (nbentries!=lastpix+1)
195 {
196 nbentries=lastpix+1;
197 }
198 }
199 // Let's Read the SphereCoordSys object
200 int id= SphereCoordSys::NEUTRAL;
201 string description= "NEUTRAL SphereCoordSystem";
202 string coordsys= dvl.GetS("COORDSYS");
203 // $CHECK$ - Reza 2/05/2000 - compare(size_t, size_t,...) passe pas sur g++
204 // if(coordsys.compare(0,7,"ROTATION",0,7) == 0) $CHECK$
205 if(coordsys.substr(0,8) == "ROTATION")
206 {
207 id= SphereCoordSys::ROTATION;
208 description= "ROTATION SphereCoordSystem";
209 }
210 // else if(coordsys.compare(0,4,"OTHER",0,4) == 0) $CHECK$ Reza 2/05/2000
211 else if(coordsys.substr(0,5) == "OTHER" )
212 {
213 id= SphereCoordSys::OTHER;
214 description= "OTHER SphereCoordSystem";
215 }
216 dobj_->SetCoordSys(new SphereCoordSys(id,description));
217
218 string ordering= dvl.GetS("ORDERING");
219 if ( (ordering.substr(0,4) != "RING" ) && (ordering.substr(0,6) != "NESTED" ) ) {
220 cerr << "FITS_SphereHEALPix/Error Not supported ORDERING= " << ordering
221 << " != RING/NESTED " << endl;
222 throw FitsIOException("FITS_SphereHEALPix/Error ordering != RING/NESTED ");
223 }
224
225 bool fgring = true;
226 if (ordering.substr(0,4) != "RING" ) fgring = false;
227
228 int nside= dvl.GetI("NSIDE");
229 int nPix = 0;
230 if(nside <= 0)
231 {
232 // throw IOExc("FITS_SphereHEALPix:: No resol parameter nside");
233 if (lastpix>0)
234 {
235 nside = sqrt((double)(nbentries/12));
236 }
237 else
238 {
239 // dvl.Print();
240 // int naxis2 = dvl.GetI("NAXIS2");
241 // if (naxis2 > 0 )
242 // {
243 nside = sqrt((double)(nbentries/12));
244 // }
245 // else
246 // throw IOExc("FITS_SphereHEALPix:: No NSIDE nor LASTPIX nor NAXIS2 on file");
247 }
248 }
249 nPix = 12*nside*nside;
250 if (nPix != nbentries)
251 {
252 cout << "WARNING: le nombre de pixels relu " << nbentries << " est incoherent avec la pixelisation Healpix qui donne nPix= " << nPix << endl;
253 }
254 double Omega= 4.0*Pi/nPix;
255 dobj_->setParameters(nside,nPix,Omega,fgring);
256 // On lit les DataBlocks;
257 //
258 dobj_->pixels_.ReSize(nPix);
259 is.GetSingleColumn(dobj_->pixels_.Data(),nPix);
260 //
261
262 // on effectue le decoupage en tranches
263 dobj_->SetThetaSlices();
264 //
265 dobj_->Info() = is.DVListFromFits();
266}
267
268template <class T>
269void FITS_SphereHEALPix<T>::Mollweide_picture_projection(char filename[])
270{
271 int ni = 300;
272 int nj = 600;
273 TMatrix<float> map(ni, nj);
274
275 int i,j;
276 for(i=0; i < ni; i++)
277 {
278 double yd = (i+0.5)/ni-0.5;
279 double facteur=2.*Pi/sin(acos(yd*2));
280 double theta = (0.5-yd)*Pi;
281 for(j=0; j < nj; j++)
282 {
283 double xa = (j+0.5)/nj-0.5;
284 double phi = xa*facteur+Pi;
285 float th=float(theta);
286 float ff=float(phi);
287 if (phi<2*Pi && phi>= 0)
288 {
289 map(i,j) = (float)dobj_->PixValSph(th, ff);
290 }
291 }
292 }
293 FitsOutFile fout(filename);
294 fout.firstImageOnPrimaryHeader(true);
295 FITS_TArray<r_4> fta(map);
296 fta.Write(fout);
297
298}
299template <class T>
300void FITS_SphereHEALPix<T>::sinus_picture_projection(char filename[])
301{
302 int ni = 300;
303 int nj = 600;
304 TMatrix<float> map(ni, nj);
305
306 int i,j;
307 for(i=0; i < ni; i++)
308 {
309 double yd = (i+0.5)/ni-0.5;
310 double theta = (0.5-yd)*Pi;
311 double facteur=1./sin(theta);
312 for(j=0; j < nj; j++)
313 {
314 double xa = (j+0.5)/nj-0.5;
315 double phi = 2.*Pi*xa*facteur+Pi;
316 float th=float(theta);
317 float ff=float(phi);
318 if (phi<2*Pi && phi>= 0)
319 {
320 map(i,j) = (float)dobj_->PixValSph(th, ff);
321 }
322 }
323 }
324 FitsOutFile fout(filename);
325 fout.firstImageOnPrimaryHeader(true);
326 FITS_TArray<r_4> fta(map);
327 fta.Write(fout);
328
329
330}
331
332
333#ifdef __CXX_PRAGMA_TEMPLATES__
334#pragma define_template FITS_SphereHEALPix<r_8>
335#pragma define_template FITS_SphereHEALPix<r_4>
336#pragma define_template FITS_SphereHEALPix<int_4>
337#endif
338#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
339namespace SOPHYA {
340template class FITS_SphereHEALPix<r_8>;
341template class FITS_SphereHEALPix<r_4>;
342template class FITS_SphereHEALPix<int_4>;
343}
344#endif
Note: See TracBrowser for help on using the repository browser.