source: trunk/FACT++/src/fitsdump.cc@ 12877

Last change on this file since 12877 was 12877, checked in by tbretz, 13 years ago
Added force option.
File size: 46.2 KB
Line 
1//****************************************************************
2/** @class FitsDumper
3
4 @brief Dumps contents of fits tables to stdout or a file
5
6 */
7 //****************************************************************
8#include "Configuration.h"
9
10#include <map>
11#include <fstream>
12
13#include "externals/fits.h"
14
15//#define PLOTTING_PLEASE
16
17#ifdef PLOTTING_PLEASE
18#include <QPen>
19#include <QtGui>
20#include <QApplication>
21
22#include <qwt_plot.h>
23#include <qwt_plot_grid.h>
24#include <qwt_plot_curve.h>
25#include <qwt_plot_zoomer.h>
26#include <qwt_legend.h>
27#include <qwt_scale_draw.h>
28
29#endif
30#include "Time.h"
31
32using namespace std;
33
34
35#ifdef PLOTTING_PLEASE
36class TimeScaleDraw: public QwtScaleDraw
37{
38public:
39 virtual QwtText label(double v) const
40 {
41 Time t(v);
42 string time = t.GetAsStr("%H:%M:%S%F");
43 while (time[time.size()-1] == '0' && time.size() > 2)
44 {
45 time = time.substr(0, time.size()-1);
46 }
47 return QwtText(time.c_str());
48 }
49};
50#endif
51
52class FitsDumper
53{
54public:
55 FitsDumper();
56 ~FitsDumper();
57
58private:
59 fits* fFile;
60 bool fDotsPlease;
61 bool fNoZeroPlease;
62 string fFilename;
63
64 fits::Table::Columns fColMap;
65 fits::Table::Keys fKeyMap;
66
67 // Convert CCfits::ValueType into a human readable string
68 string ValueTypeToStr(char type) const;
69
70 // Convert CCfits::ValueType into a number of associated bytes
71 int ValueTypeToSize(char type) const;
72
73 /// Calculate the buffer size required to read a row of the fits table, as well as the offsets to each column
74// vector<int> CalculateOffsets() const;
75
76 template<class T>
77 T PtrToValue(const unsigned char* &ptr) const;
78// template<class T>
79// double PtrToDouble(const unsigned char *ptr) const;
80// double PtrToDouble(const unsigned char *ptr, CCfits::ValueType type) const;
81
82 /// Write a single row of the selected data
83 // int WriteRow(ostream &, const vector<MyColumn*> &, const vector<int> &, unsigned char *, const vector<pair<int, int> >&) const;
84
85 bool OpenFile(const string &, bool); /// Open a file
86 bool OpenTable(const string &); /// Open a table
87
88 /// Lists all columns of an open file
89 void List();
90 void ListHeader(const string& filename);
91 void ListKeywords(ostream &);
92
93 bool separateColumnsFromRanges(const vector<string>& list,
94 vector<pair<int, int> >& ranges,
95 vector<string>& listNamesOnly);
96 /// Perform the dumping, based on the current dump list
97 bool Dump(const string &, const vector<string> &list, int);
98 ///Display the selected columns values VS time
99#ifdef PLOTTING_PLEASE
100 int doCurvesDisplay( const vector<string> &list, const string& tableName);
101#endif
102 int doMinMaxPlease(const string& filename, const vector<string>& list, int precision);
103 int doStatsPlease(const string &filename, const vector<string>& list, int precision);
104// void doTBoundary(conf.Get<string>("outfile"), conf.Get<int>("precision"), true);
105 // bool Plot(const vector<string> &list);
106
107public:
108 ///Configures the fitsLoader from the config file and/or command arguments.
109 int ExecConfig(Configuration& conf);
110};
111
112// --------------------------------------------------------------------------
113//
114//! Constructor
115//! @param out
116//! the ostream where to redirect the outputs
117//
118FitsDumper::FitsDumper() : fFile(0), fDotsPlease(false), fNoZeroPlease(false)
119{
120}
121
122// --------------------------------------------------------------------------
123//
124//! Destructor
125//
126FitsDumper::~FitsDumper()
127{
128 if (fFile)
129 delete fFile;
130}
131
132
133string FitsDumper::ValueTypeToStr(char type) const
134{
135 switch (type)
136 {
137 case 'L': return "bool(8)";
138 case 'A': return "char(8)";
139 case 'B': return "byte(8)";
140 case 'I': return "short(16)";
141 case 'J': return "int(32)";
142 case 'K': return "int(64)";
143 case 'E': return "float(32)";
144 case 'D': return "double(64)";
145 default:
146 return "unknown";
147 }
148}
149
150int FitsDumper::ValueTypeToSize(char type) const
151{
152 switch (type)
153 {
154 case 'A': return sizeof(int8_t);
155 case 'L': return sizeof(uint8_t);
156 case 'B': return sizeof(int8_t);
157 case 'I': return sizeof(int16_t);
158 case 'J': return sizeof(int32_t);
159 case 'K': return sizeof(int64_t);
160 case 'E': return sizeof(float);
161 case 'D': return sizeof(double);
162 default:
163 return 0;
164 }
165}
166
167template<class T>
168T FitsDumper::PtrToValue(const unsigned char* &ptr) const
169{
170 T t;
171 reverse_copy(ptr, ptr+sizeof(T), reinterpret_cast<unsigned char*>(&t));
172 ptr += sizeof(T);
173
174 return t;
175}
176// --------------------------------------------------------------------------
177//
178//! Loads the fits file based on the current parameters
179//
180bool FitsDumper::OpenFile(const string &filename, bool force)
181{
182 if (fFile)
183 {
184 fFile->close();
185 delete fFile;
186 }
187
188 try {
189 fFile = new fits(filename, force);
190 }
191 catch (std::runtime_error e)
192 {
193 cout << "Something went wrong while trying to open " << filename;
194 cout << ": " << e.what() << " Aborting dump." << endl;
195 return false;
196 }
197 fFilename = filename;
198
199 const fits::Table::Columns& tCols = fFile->GetColumns();
200
201 for (auto it=tCols.begin(); it != tCols.end(); it++)
202 fColMap.insert(*it);
203
204 const fits::Table::Keys& tkeys = fFile->GetKeys();
205
206 for (auto it=tkeys.begin(); it != tkeys.end(); it++)
207 fKeyMap.insert(*it);
208
209 return true;
210}
211
212bool FitsDumper::OpenTable(const string &)
213{
214 if (!fFile)
215 {
216 cerr << "No file open." << endl;
217 return false;
218 }
219
220
221 return true;
222}
223
224
225void FitsDumper::List()
226{
227 if (!fFile)
228 {
229 cerr << "No file open." << endl;
230 return;
231 }
232
233 cout << "\nFile: " << fFilename << "\n";
234
235 cout << " " << fKeyMap.find("EXTNAME")->second.value << " [";
236 cout << fKeyMap.find("NAXIS2")->second.value << "]\n";
237
238 for (auto it = fColMap.begin(); it != fColMap.end(); it++)
239 {
240 cout << " " << it->first << "[" << it->second.num << "] (" << it->second.unit << ":" << ValueTypeToStr(it->second.type) << ") ";
241 for (auto jt = fKeyMap.begin(); jt != fKeyMap.end(); jt++)
242 if (jt->second.value == it->first)
243 cout << jt->second.comment << endl;
244 }
245
246 cout << endl;
247 cout << flush;
248}
249
250void FitsDumper::ListKeywords(ostream &out)
251{
252 for (auto it=fKeyMap.begin(); it != fKeyMap.end(); it++) {
253 out << "## " << setw(8) << it->first << " = " << setw(10);
254 out << "'" << it->second.value << "'" << " / " << it->second.comment << endl;
255 }
256}
257
258void FitsDumper::ListHeader(const string& filename)
259{
260 if (!fFile)
261 {
262 cerr << "No table open." << endl;
263 return;
264 }
265 ofstream out(filename=="-"?"/dev/stdout":filename);
266 if (!out)
267 {
268 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
269 return;
270 }
271
272 out << "\nTable: " << fKeyMap.find("EXTNAME")->second.value << " (rows=" << fKeyMap.find("NAXIS2")->second.value << ")\n";
273 if (fKeyMap.find("COMMENT") != fKeyMap.end())
274 out << "Comment: \t" << fKeyMap.find("COMMENT")->second.value << "\n";
275
276 ListKeywords(out);
277 out << endl;
278
279}
280
281bool FitsDumper::separateColumnsFromRanges(const vector<string>& list,
282 vector<pair<int, int> >& ranges,
283 vector<string>& listNamesOnly)
284{
285 for (vector<string>::const_iterator it=list.begin(); it!=list.end(); it++)
286 {
287 string columnNameOnly = *it;
288 unsigned long bracketIndex0 = columnNameOnly.find_first_of('[');
289 unsigned long bracketIndex1 = columnNameOnly.find_first_of(']');
290 unsigned long colonIndex = columnNameOnly.find_first_of(':');
291 int columnStart = -1;
292 int columnEnd = -1;
293 if (bracketIndex0 != string::npos)
294 {//there is a range given. Extract the range
295 if (colonIndex != string::npos)
296 {//we have a range here
297 columnStart = atoi(columnNameOnly.substr(bracketIndex0+1, colonIndex-(bracketIndex0+1)).c_str());
298 columnEnd = atoi(columnNameOnly.substr(colonIndex+1, bracketIndex1-(colonIndex+1)).c_str());
299 columnEnd++;
300 }
301 else
302 {//only a single index there
303 columnStart = atoi(columnNameOnly.substr(bracketIndex0+1, bracketIndex1 - (bracketIndex0+1)).c_str());
304 columnEnd = columnStart+1;
305 }
306 columnNameOnly = columnNameOnly.substr(0, bracketIndex0);
307 }
308
309 if (fColMap.find(columnNameOnly) == fColMap.end())
310 {
311 cerr << "ERROR - Column '" << columnNameOnly << "' not found in table." << endl;
312 return false;
313 }
314 fits::Table::Column& cCol = fColMap.find(columnNameOnly)->second;
315 if (bracketIndex0 == string::npos)
316 {//no range given: use the full range
317 ranges.push_back(make_pair(0, cCol.num));
318 columnStart = 0;
319 columnEnd = 1;
320 }
321 else
322 {//use the range extracted earlier
323 if (columnStart < 0)
324 {
325 cerr << "ERROR - Start range for column " << columnNameOnly << " is less than zero (" << columnStart << "). Aborting" << endl;
326 return false;
327 }
328 if (columnEnd>1 && columnEnd > (int)(cCol.num))
329 {
330 cerr << "ERROR - End range for column " << columnNameOnly << " is greater than the last element (" << cCol.num << " vs " << columnEnd << "). Aborting" << endl;
331 return false;
332 }
333 ranges.push_back(make_pair(columnStart, columnEnd));
334 }
335// cout << "Will be exporting from " << columnStart << " to " << columnEnd-1 << " for column " << columnNameOnly << endl;
336 listNamesOnly.push_back(columnNameOnly);
337 }
338 return true;
339}
340// --------------------------------------------------------------------------
341//
342//! Perform the actual dump, based on the current parameters
343//
344bool FitsDumper::Dump(const string &filename, const vector<string> &list, int precision)
345{
346
347 //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
348 vector<pair<int, int> > ranges;
349 vector<string> listNamesOnly;
350
351 if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
352 {
353 cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
354 return false;
355 }
356
357 ofstream out(filename=="-"?"/dev/stdout":filename);
358 if (!out)
359 {
360 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
361 return false;
362 }
363
364 out.precision(precision);
365
366 out << "## --------------------------------------------------------------------------\n";
367 if (filename!="-")
368 out << "## File: \t" << filename << '\n';
369 out << "## Table: \t" << fKeyMap.find("EXTNAME")->second.value << '\n';
370 out << "## Comment: \t" << ((fKeyMap.find("COMMENT") != fKeyMap.end()) ? fKeyMap.find("COMMENT")->second.value : "") << '\n';
371
372 out << "## NumRows: \t" << fFile->GetInt("NAXIS2") << '\n';
373 out << "## --------------------------------------------------------------------------\n";
374 ListKeywords(out);
375 out << "## --------------------------------------------------------------------------\n";
376 out << "#\n";
377
378 auto rangesIt = ranges.begin();
379 for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
380 {
381 const fits::Table::Column& col = fColMap[*it];
382// const MyColumn *col = static_cast<MyColumn*>(fTable->column()[*it]);
383 if (rangesIt->first != 0 || rangesIt->second != (int)(col.num))
384 {
385 out << "#";
386 for (int i=rangesIt->first; i<rangesIt->second; i++)
387 out << " " << *it << "[" << i << "]";
388 out << ": " << col.unit;
389 }
390 else
391 out << "# " << *it << "[" << col.num << "]: " << col.unit;
392
393// FIXME: retrive the column comment
394// if (!col->comment().empty())
395// out << " (" <<col->comment() << ")";
396 out << '\n';
397 }
398 out << "#" << endl;
399
400 while (fFile->GetNextRow())
401 {
402 const size_t row = fFile->GetRow();
403 if (row==fFile->GetNumRows())
404 break;
405
406 rangesIt = ranges.begin();
407
408 for (vector<string>::const_iterator in=listNamesOnly.begin(); in!=listNamesOnly.end(); in++, rangesIt++)
409 {
410 const void *ptr = fFile->SetPtrAddress(*in);
411
412 const fits::Table::Column &cCol = fColMap.find(*in)->second;
413
414 string msg;
415 for (int i=rangesIt->first; i<rangesIt->second; i++)
416 {
417 switch (cCol.type)
418 {
419 case 'A':
420 msg += reinterpret_cast<const char*>(ptr)[i];
421 break;
422 case 'B':
423 out << (unsigned int)reinterpret_cast<const unsigned char*>(ptr)[i] << " ";
424 break;
425 case 'L':
426 out << reinterpret_cast<const bool*>(ptr)[i] << " ";
427 break;
428 case 'I':
429 out << reinterpret_cast<const int16_t*>(ptr)[i] << " ";
430 break;
431 case 'J':
432 out << reinterpret_cast<const int32_t*>(ptr)[i] << " ";
433 break;
434 case 'K':
435 out << reinterpret_cast<const int64_t*>(ptr)[i] << " ";
436 break;
437 case 'E':
438 out << reinterpret_cast<const float*>(ptr)[i] << " ";
439 break;
440 case 'D':
441 out << reinterpret_cast<const double*>(ptr)[i] << " ";
442 break;
443 default:
444 ;
445 }
446 }
447
448 if (cCol.type=='A')
449 out << "'" << msg << "' ";
450 }
451 out << endl;
452 }
453 return true;
454}
455
456// --------------------------------------------------------------------------
457//
458//! Retrieves the configuration parameters
459//! @param conf
460//! the configuration object
461//
462int FitsDumper::ExecConfig(Configuration& conf)
463{
464 if (conf.Has("fitsfile"))
465 {
466 if (!OpenFile(conf.Get<string>("fitsfile"), conf.Get<bool>("force")))
467 return -1;
468 }
469#ifdef PLOTTING_PLEASE
470 if (conf.Get<bool>("stat") && conf.Get<bool>("graph"))
471 {
472 cout << "Invalid conbination of options: cannot graph stats. Aborting" << endl;
473 return -1;
474 }
475#endif
476 if (conf.Get<bool>("minmax") && conf.Get<bool>("stat"))
477 {
478 cout << "Invalid combination of options: cannot do stats and minmax. Aborting" << endl;
479 return -1;
480 }
481#ifdef PLOTTING_PLEASE
482 if (conf.Get<bool>("minmax") && conf.Get<bool>("graph"))
483 {
484 cout << "Invalid combination of options: cannot graph minmax. Aborting" << endl;
485 return -1;
486 }
487#endif
488 if (conf.Get<bool>("stat") && conf.Get<bool>("nozero"))
489 {
490 cout << "Invalid combination of options: nozero only works with minmax. Aborting" << endl;
491 return -1;
492 }
493
494 if (conf.Get<bool>("nozero"))
495 {
496 fNoZeroPlease = true;
497 }
498
499 if (conf.Has("tablename"))
500 {
501 if (!OpenTable(conf.Get<string>("tablename")))
502 return -1;
503 }
504
505 if (conf.Get<bool>("list"))
506 List();
507 /*
508 if (conf.Get<bool>("tstart"))
509 {
510 doTBoundary(conf.Get<string>("outfile"), conf.Get<int>("precision"), true);
511 }
512 if (conf.Get<bool>("tstop"))
513 {
514 doTBoundary(conf.Get<string>("outfile"), conf.Get<int>("precision"), false);
515 }
516 */
517 if (conf.Get<bool>("minmax"))
518 {
519 if (!conf.Has("col"))
520 {
521 cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl;
522 return 0;
523 }
524 doMinMaxPlease(conf.Get<string>("outfile"), conf.Get<vector<string>>("col"), conf.Get<int>("precision"));
525 return 0;
526 }
527
528 if (conf.Get<bool>("stat"))
529 {
530 if (!conf.Has("col"))
531 {
532 cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl;
533 return 0;
534 }
535 doStatsPlease(conf.Get<string>("outfile"), conf.Get<vector<string>>("col"), conf.Get<int>("precision"));
536 return 0;
537 }
538
539#ifdef PLOTTING_PLEASE
540 if (conf.Get<bool>("graph"))
541 {
542 if (!conf.Has("col"))
543 {
544 cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl;
545 return 0;
546 }
547 if (conf.Get<bool>("dots"))
548 fDotsPlease = true;
549 doCurvesDisplay(conf.Get<vector<string>>("col"),
550 conf.Get<string>("tablename"));
551 return 1;
552 }
553#endif
554
555 if (conf.Get<bool>("header"))
556 ListHeader(conf.Get<string>("outfile"));
557
558 if (conf.Get<bool>("header") || conf.Get<bool>("list"))
559 return 1;
560
561 if (conf.Has("outfile"))
562 {
563 if (!conf.Has("col"))
564 {
565 cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl;
566 return 0;
567 }
568 if (!Dump(conf.Get<string>("outfile"),
569 conf.Get<vector<string>>("col"),
570 conf.Get<int>("precision")))
571 return -1;
572 }
573
574
575 return 0;
576}
577
578void PrintUsage()
579{
580 cout <<
581 "fitsdump is a tool to dump data from a FITS table as ascii.\n"
582 "\n"
583 "Usage: fitsdump [OPTIONS] fitsfile col col ... \n"
584 " or: fitsdump [OPTIONS]\n";
585 cout << endl;
586}
587
588void PrintHelp()
589{
590 //
591}
592
593struct minMaxStruct {
594 double min;
595 double max;
596 double average;
597 double squared;
598 long numValues;
599 minMaxStruct() : min(1e10), max(-1e10), average(0), numValues(0) {}
600};
601
602int FitsDumper::doMinMaxPlease(const string& filename, const vector<string>& list, int precision)
603{
604 //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
605 vector<pair<int, int> > ranges;
606 vector<string> names;
607
608 if (!separateColumnsFromRanges(list, ranges, names))
609 {
610 cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
611 return false;
612 }
613
614 ofstream out(filename=="-"?"/dev/stdout":filename);
615 if (!out)
616 {
617 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
618 return false;
619 }
620
621 out.precision(precision);
622
623 vector<minMaxStruct> statData(names.size());
624
625 // Loop over all columns in our list of requested columns
626 while (fFile->GetNextRow())
627 {
628 const size_t row = fFile->GetRow();
629 if (row==fFile->GetNumRows())
630 break;
631
632 auto rangesIt = ranges.begin();
633 auto statsIt = statData.begin();
634
635 for (auto in=names.begin(); in!=names.end(); in++, rangesIt++, statsIt++)
636 {
637 const void *ptr = fFile->SetPtrAddress(*in);
638
639 const fits::Table::Column &cCol = fColMap.find(*in)->second;
640
641 if (*in=="UnixTimeUTC" || *in=="PCTime")
642 {
643 const uint32_t *val = reinterpret_cast<const uint32_t*>(ptr);
644 if (fNoZeroPlease && val[0]==0 && val[1]==0)
645 continue;
646
647 const double cValue = Time(val[0], val[1]).Mjd();
648
649 statsIt->average += cValue;
650 statsIt->squared += cValue*cValue;
651
652 if (cValue < statsIt->min)
653 statsIt->min = cValue;
654
655 if (cValue > statsIt->max)
656 statsIt->max = cValue;
657
658 statsIt->numValues++;
659
660 continue;
661 }
662
663 for (int i=rangesIt->first; i<rangesIt->second; i++)
664 {
665 double cValue = 0;
666 switch (cCol.type)
667 {
668 case 'L':
669 cValue = reinterpret_cast<const bool*>(ptr)[i];
670 break;
671 case 'B':
672 cValue = reinterpret_cast<const int8_t*>(ptr)[i];
673 break;
674 case 'I':
675 cValue = reinterpret_cast<const int16_t*>(ptr)[i];
676 break;
677 case 'J':
678 cValue = reinterpret_cast<const int32_t*>(ptr)[i];
679 break;
680 case 'K':
681 cValue = reinterpret_cast<const int64_t*>(ptr)[i];
682 break;
683 case 'E':
684 cValue = reinterpret_cast<const float*>(ptr)[i];
685 break;
686 case 'D':
687 cValue = reinterpret_cast<const double*>(ptr)[i];
688 break;
689 default:
690 ;
691 }
692
693 if (fNoZeroPlease && cValue == 0)
694 continue;
695
696 statsIt->average += cValue;
697 statsIt->squared += cValue*cValue;
698
699 if (cValue < statsIt->min)
700 statsIt->min = cValue;
701
702 if (cValue > statsIt->max)
703 statsIt->max = cValue;
704
705 statsIt->numValues++;
706 }
707 }
708 }
709
710 // okay. So now I've got ALL the data, loaded.
711 // let's do the summing and averaging in a safe way (i.e. avoid overflow
712 // of variables as much as possible)
713 auto statsIt = statData.begin();
714
715 for (auto it=names.begin(); it!=names.end(); it++, statsIt++)
716 {
717 cout << "[" << *it << "]" << endl;
718 if (statsIt->numValues == 0)
719 {
720 out << "Min: -\nMax: -\nAvg: -\nRms: -" << endl;
721 continue;
722 }
723
724 double &avg = statsIt->average;
725 double &rms = statsIt->squared;
726 long &num = statsIt->numValues;
727
728 avg /= num;
729 rms = sqrt(rms/num - avg*avg);
730
731 out << "Min: " << statsIt->min << '\n';
732 out << "Max: " << statsIt->max << '\n';
733 out << "Avg: " << avg << '\n';
734 out << "Rms: " << rms << endl;
735 }
736
737 /*
738 vector<pair<char, char*> > columnsData;
739 vector<minMaxStruct> statData;
740 int numRows = fFile->GetInt("NAXIS2");
741 auto rangesIt = ranges.begin();
742 for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
743 {
744 fits::Table::Column& cCol = fColMap.find(*it)->second;
745 columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
746// minMaxStuct initData;
747 statData.push_back(minMaxStruct());
748 fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
749 }
750
751 int row = 0;
752 long UTCvalue0=0;
753 long UTCvalue1=0;
754 while (fFile->GetNextRow() && row < numRows)
755 {
756 rangesIt = ranges.begin();
757 auto statsIt = statData.begin();
758 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++)
759 {
760 double cValue = 0;
761 for (int i=rangesIt->first; i<rangesIt->second; i++)
762 {
763 switch (it->first) {
764 case 'L':
765 cValue = reinterpret_cast<unsigned char*>(it->second)[i];
766 break;
767 case 'B':
768 cValue = reinterpret_cast<bool*>(it->second)[i];
769 break;
770 case 'I':
771 cValue = reinterpret_cast<int16_t*>(it->second)[i];
772 break;
773 case 'J':
774 cValue = reinterpret_cast<int32_t*>(it->second)[i];
775 break;
776 case 'K':
777 cValue = reinterpret_cast<int64_t*>(it->second)[i];
778 break;
779 case 'E':
780 cValue = reinterpret_cast<float*>(it->second)[i];
781 break;
782 case 'D':
783 cValue = reinterpret_cast<double*>(it->second)[i];
784 break;
785 default:
786 ;
787 }
788 if (list.size() == 1 && (list[0] == "UnixTimeUTC" || list[0] == "PCTime"))
789 {
790 if (i==0)
791 {
792 UTCvalue0 = cValue;
793 }
794 else
795 {
796 UTCvalue1 = cValue;
797 boost::posix_time::ptime unixTimeT( boost::gregorian::date(1970, boost::gregorian::Jan, 1),
798 boost::posix_time::seconds(UTCvalue0) + boost::posix_time::microsec(UTCvalue1));
799
800 Time mjdTime(unixTimeT);
801 cValue = mjdTime.Mjd();
802 if (!fNoZeroPlease || cValue != 0)
803 {
804 statsIt->average += cValue;
805 if (cValue < statsIt->min)
806 statsIt->min = cValue;
807 if (cValue > statsIt->max)
808 statsIt->max = cValue;
809 statsIt->numValues++;
810 }
811 }
812 }
813 else {
814 if (!fNoZeroPlease || cValue != 0)
815 {
816 statsIt->average += cValue;
817 if (cValue < statsIt->min)
818 statsIt->min = cValue;
819 if (cValue > statsIt->max)
820 statsIt->max = cValue;
821 statsIt->numValues++;
822 }
823 }
824 }
825 }
826 row++;
827 }
828 for (auto it = columnsData.begin(); it != columnsData.end(); it++)
829 delete[] it->second;
830*/
831 return true;
832}
833
834/*
835void FitsDumper::doTBoundary(const string& filename, int precision, bool tStop)
836{
837
838 //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
839 vector<pair<int, int> > ranges;
840 vector<string> listNamesOnly;
841
842 if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
843 {
844 cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
845 return false;
846 }
847
848 ofstream out(filename=="-"?"/dev/stdout":filename);
849 if (!out)
850 {
851 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
852 return false;
853 }
854
855 out.precision(precision);
856
857 // Loop over all columns in our list of requested columns
858 vector<pair<char, char*> > columnsData;
859 vector<minMaxStruct> statData;
860 int numRows = fFile->GetInt("NAXIS2");
861 auto rangesIt = ranges.begin();
862 for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
863 {
864 fits::Table::Column& cCol = fColMap.find(*it)->second;
865 columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
866 // minMaxStuct initData;
867 statData.push_back(minMaxStruct());
868 fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
869 }
870
871 int row = 0;
872 while (fFile->GetNextRow() && row < numRows)
873 {
874 rangesIt = ranges.begin();
875 auto statsIt = statData.begin();
876 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++)
877 {
878 double cValue = 0;
879 for (int i=rangesIt->first; i<rangesIt->second; i++)
880 {
881 switch (it->first) {
882 case 'L':
883 cValue = reinterpret_cast<bool*>(it->second)[i];
884 break;
885 case 'B':
886 cValue = reinterpret_cast<bool*>(it->second)[i];
887 break;
888 case 'I':
889 cValue = reinterpret_cast<int16_t*>(it->second)[i];
890 break;
891 case 'J':
892 cValue = reinterpret_cast<int32_t*>(it->second)[i];
893 break;
894 case 'K':
895 cValue = reinterpret_cast<int64_t*>(it->second)[i];
896 break;
897 case 'E':
898 cValue = reinterpret_cast<float*>(it->second)[i];
899 break;
900 case 'D':
901 cValue = reinterpret_cast<double*>(it->second)[i];
902 break;
903 default:
904 ;
905 }
906 if (!fNoZeroPlease || cValue != 0)
907 {
908 statsIt->average += cValue;
909 if (cValue < statsIt->min)
910 statsIt->min = cValue;
911 if (cValue > statsIt->max)
912 statsIt->max = cValue;
913 statsIt->numValues++;
914 }
915 }
916 }
917 row++;
918 }
919 for (auto it = columnsData.begin(); it != columnsData.end(); it++)
920 delete[] it->second;
921
922 //okay. So now I've got ALL the data, loaded.
923 //let's do the summing and averaging in a safe way (i.e. avoid overflow of variables as much as possible)
924 rangesIt = ranges.begin();
925 auto statsIt = statData.begin();
926
927 auto nameIt = listNamesOnly.begin();
928 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++, nameIt++)
929 {
930 int span = rangesIt->second - rangesIt->first;
931 cout << *nameIt << ": " << endl;
932 if (statsIt->numValues != 0)
933 {
934 statsIt->average /= statsIt->numValues;
935 out << "min: " << statsIt->min << endl;
936 out << "max: " << statsIt->max << endl;
937 out << "mea: " << statsIt->average << endl;
938 }
939 else
940 {
941 out << "min: 0" << endl << "max: 0" << endl << "mea: " << endl;
942 }
943
944 }
945 return true;
946
947}
948*/
949template<typename T>
950void displayStats(char* array, int numElems, ofstream& out)
951{
952 if (numElems == 0)
953 {
954 out << "Min: -\nMax: -\nMed: -\nAvg: -\nRms: -" << endl;
955 return;
956 }
957
958 T *val = reinterpret_cast<T*>(array);
959
960 sort(val, val+numElems);
961
962 out << "Min: " << double(val[0]) << '\n';
963 out << "Max: " << double(val[numElems-1]) << '\n';
964
965 if (numElems>2)
966 {
967 if (numElems%2 == 0)
968 out << "Med: " << double(val[numElems/2] + val[numElems/2+1])/2 << '\n';
969 else
970 out << "Med: " << double(val[numElems/2+1]) << '\n';
971 }
972
973 double avg = 0;
974 double rms = 0;
975 for (int i=0;i<numElems;i++)
976 {
977 avg += val[i];
978 rms += val[i]*val[i];
979 }
980
981 avg /= numElems;
982 rms = sqrt(rms/numElems - avg*avg);
983
984 out << "Avg: " << avg << '\n';
985 out << "Rms: " << rms << endl;
986
987}
988int FitsDumper::doStatsPlease(const string &filename, const vector<string>& list, int precision)
989{
990
991 //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
992 vector<pair<int, int> > ranges;
993 vector<string> names;
994
995 if (!separateColumnsFromRanges(list, ranges, names))
996 {
997 cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
998 return false;
999 }
1000
1001 ofstream out(filename=="-"?"/dev/stdout":filename);
1002 if (!out)
1003 {
1004 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
1005 return false;
1006 }
1007
1008 out.precision(precision);
1009
1010 // Loop over all columns in our list of requested columns
1011 vector<char*> statData;
1012
1013 auto rangesIt = ranges.begin();
1014 for (auto in=names.begin(); in!=names.end(); in++, rangesIt++)
1015 {
1016 cout << *in << endl;
1017 const fits::Table::Column &cCol = fColMap.find(*in)->second;
1018 statData.push_back(new char[(rangesIt->second - rangesIt->first)*cCol.size*fFile->GetNumRows()]);
1019 }
1020
1021 while (fFile->GetNextRow())
1022 {
1023 const size_t row = fFile->GetRow();
1024 if (row==fFile->GetNumRows())
1025 break;
1026
1027 rangesIt = ranges.begin();
1028 auto statsIt = statData.begin();
1029
1030 for (auto in=names.begin(); in!=names.end(); in++, rangesIt++, statsIt++)
1031 {
1032 const void *ptr = fFile->SetPtrAddress(*in);
1033
1034 const fits::Table::Column &cCol = fColMap.find(*in)->second;
1035
1036 const int span = rangesIt->second - rangesIt->first;
1037
1038 for (int i=rangesIt->first; i<rangesIt->second; i++)
1039 {
1040 switch (cCol.type)
1041 {
1042 case 'L':
1043 reinterpret_cast<bool*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const bool*>(ptr)[i];
1044 break;
1045 case 'B':
1046 reinterpret_cast<char*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const char*>(ptr)[i];
1047 break;
1048 case 'I':
1049 reinterpret_cast<int16_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const int16_t*>(ptr)[i];
1050 break;
1051 case 'J':
1052 reinterpret_cast<int32_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const int32_t*>(ptr)[i];
1053 break;
1054 case 'K':
1055 reinterpret_cast<int64_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const int64_t*>(ptr)[i];
1056 break;
1057 case 'E':
1058 reinterpret_cast<float*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const float*>(ptr)[i];
1059 break;
1060 case 'D':
1061 reinterpret_cast<double*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const double*>(ptr)[i];
1062 break;
1063 default:
1064 ;
1065 }
1066 }
1067 }
1068 }
1069
1070 //okay. So now I've got ALL the data, loaded.
1071 //let's do the summing and averaging in a safe way (i.e. avoid overflow of variables as much as possible)
1072 rangesIt = ranges.begin();
1073 auto statsIt = statData.begin();
1074
1075 for (auto it=list.begin(); it!=list.end(); it++, rangesIt++, statsIt++)
1076 {
1077 const fits::Table::Column &cCol = fColMap.find(*it)->second;
1078
1079 const int span = rangesIt->second - rangesIt->first;
1080
1081 out << "[" << *it << "]" << endl;
1082 switch (cCol.type)
1083 {
1084 case 'L':
1085 displayStats<bool>(*statsIt, fFile->GetNumRows()*span, out);
1086 break;
1087 case 'B':
1088 displayStats<char>(*statsIt, fFile->GetNumRows()*span, out);
1089 break;
1090 case 'I':
1091 displayStats<int16_t>(*statsIt, fFile->GetNumRows()*span, out);
1092 break;
1093 case 'J':
1094 displayStats<int32_t>(*statsIt, fFile->GetNumRows()*span, out);
1095 break;
1096 case 'K':
1097 displayStats<int64_t>(*statsIt, fFile->GetNumRows()*span, out);
1098 break;
1099 case 'E':
1100 displayStats<float>(*statsIt, fFile->GetNumRows()*span, out);
1101 break;
1102 case 'D':
1103 displayStats<double>(*statsIt, fFile->GetNumRows()*span, out);
1104 break;
1105 default:
1106 ;
1107 }
1108 }
1109
1110 return true;
1111}
1112#ifdef PLOTTING_PLEASE
1113int FitsDumper::doCurvesDisplay( const vector<string> &list, const string& tableName)
1114{
1115 //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
1116 vector<pair<int, int> > ranges;
1117 vector<string> listNamesOnly;
1118 if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
1119 {
1120 cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
1121 return false;
1122 }
1123 vector<string> curvesNames;
1124 stringstream str;
1125 for (auto it=ranges.begin(), jt=listNamesOnly.begin(); it != ranges.end(); it++, jt++)
1126 {
1127 for (int i=it->first; i<it->second;i++)
1128 {
1129 str.str("");
1130 str << *jt << "[" << i << "]";
1131 curvesNames.push_back(str.str());
1132 }
1133 }
1134 char* handle = new char[17];
1135 sprintf(handle,"FitsDump Display");
1136// Qt::HANDLE h = *handle;//NULL
1137 int argc = 1;
1138 char ** argv = &handle;
1139 QApplication a(argc, argv);
1140
1141
1142
1143 QwtPlot* plot = new QwtPlot();
1144 QwtPlotGrid* grid = new QwtPlotGrid;
1145 grid->enableX(false);
1146 grid->enableY(true);
1147 grid->enableXMin(false);
1148 grid->enableYMin(false);
1149 grid->setMajPen(QPen(Qt::black, 0, Qt::DotLine));
1150 grid->attach(plot);
1151 plot->setAutoReplot(true);
1152 string title = tableName;
1153 plot->setAxisScaleDraw( QwtPlot::xBottom, new TimeScaleDraw());
1154
1155 QWidget window;
1156 QHBoxLayout* layout = new QHBoxLayout(&window);
1157 layout->setContentsMargins(0,0,0,0);
1158 layout->addWidget(plot);
1159
1160 QwtPlotZoomer zoom(plot->canvas());
1161 zoom.setRubberBandPen(QPen(Qt::gray, 2, Qt::DotLine));
1162 zoom.setTrackerPen(QPen(Qt::gray));
1163 int totalSize = 0;
1164 for (unsigned int i=0;i<list.size();i++)
1165 totalSize += ranges[i].second - ranges[i].first;
1166
1167 vector<QwtPlotCurve*> curves(totalSize);
1168 int ii=0;
1169 for (auto it = curves.begin(), jt=curvesNames.begin(); it != curves.end(); it++, jt++)
1170 {
1171 *it = new QwtPlotCurve(jt->c_str());
1172 switch (ii%6)
1173 {
1174 case 0:
1175 (*it)->setPen(QColor(255,0,0));
1176 break;
1177 case 1:
1178 (*it)->setPen(QColor(0,255,0));
1179 break;
1180 case 2:
1181 (*it)->setPen(QColor(0,0,255));
1182 break;
1183 case 3:
1184 (*it)->setPen(QColor(255,255,0));
1185 break;
1186 case 4:
1187 (*it)->setPen(QColor(0,255,255));
1188 break;
1189 case 5:
1190 (*it)->setPen(QColor(255,0,255));
1191 break;
1192 default:
1193 (*it)->setPen(QColor(0,0,0));
1194 };
1195 ii++;
1196 if (fDotsPlease)
1197 (*it)->setStyle(QwtPlotCurve::Dots);
1198 else
1199 (*it)->setStyle(QwtPlotCurve::Lines);
1200 (*it)->attach(plot);
1201 }
1202 plot->insertLegend(new QwtLegend(), QwtPlot::RightLegend);
1203
1204
1205 vector<pair<char, char*> > columnsData;
1206
1207 for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++)
1208 {
1209 fits::Table::Column& cCol = fColMap.find(*it)->second;
1210 columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
1211 fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
1212 }
1213
1214 //add the time column to the given columns
1215 if (fColMap.find("Time") == fColMap.end() &&
1216 fColMap.find("UnixTimeUTC") == fColMap.end())
1217 {
1218 cerr << "Error: time column could not be found in given table. Aborting" << endl;
1219 return false;
1220 }
1221 const fits::Table::Column& timeCol = (fColMap.find("Time") != fColMap.end()) ? fColMap.find("Time")->second : fColMap.find("UnixTimeUTC")->second;
1222 bool unixTime = (fColMap.find("Time") == fColMap.end());
1223 if (unixTime)
1224 ranges.push_back(make_pair(0,2));
1225 else
1226 ranges.push_back(make_pair(0,1));
1227 columnsData.push_back(make_pair(timeCol.type, new char[timeCol.num*timeCol.size]));
1228 fFile->SetPtrAddress(unixTime ? "UnixTimeUTC" : "Time", columnsData[columnsData.size()-1].second);
1229
1230// stringstream str;
1231 str.str("");
1232
1233
1234 vector<double*> xValues(totalSize);
1235 double* yValues;
1236 cout.precision(10);
1237 str.precision(20);
1238 for (auto it=xValues.begin(); it!=xValues.end(); it++)
1239 *it = new double[fFile->GetInt("NAXIS2")];
1240
1241 yValues = new double[fFile->GetInt("NAXIS2")];
1242
1243 cout.precision(3);
1244 int endIndex = 0;
1245 int numRows = fFile->GetInt("NAXIS2");
1246 for (int i=1;i<numRows;i++)
1247 {
1248 fFile->GetNextRow();
1249 cout << "\r" << "Constructing graph " << ((float)(endIndex)/(float)(fFile->GetInt("NAXIS2")))*100.0 << "%";
1250 endIndex++;
1251 auto rangesIt = ranges.begin();
1252 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++)
1253 {
1254 for (int j=rangesIt->first; j<rangesIt->second; j++)
1255 {
1256 switch (it->first) {
1257 case 'L':
1258 str << reinterpret_cast<bool*>(it->second)[j] << " ";
1259 break;
1260 case 'B':
1261 str << reinterpret_cast<char*>(it->second)[j] << " ";
1262 break;
1263 case 'I':
1264 str << reinterpret_cast<int16_t*>(it->second)[j] << " ";
1265 break;
1266 case 'J':
1267 str << reinterpret_cast<int32_t*>(it->second)[j] << " ";
1268 break;
1269 case 'K':
1270 str << reinterpret_cast<int64_t*>(it->second)[j] << " ";
1271 break;
1272 case 'E':
1273 str << reinterpret_cast<float*>(it->second)[j] << " ";
1274 break;
1275 case 'D':
1276 str << reinterpret_cast<double*>(it->second)[j] << " ";
1277 break;
1278 default:
1279 ;
1280 }
1281 }
1282 }
1283
1284 for (auto it=xValues.begin(); it!= xValues.end(); it++)
1285 {
1286 str >> (*it)[i-1];
1287 }
1288 if (unixTime)
1289 {
1290 long u1, u2;
1291 str >> u1 >> u2;
1292
1293 boost::posix_time::ptime unixTimeT( boost::gregorian::date(1970, boost::gregorian::Jan, 1),
1294 boost::posix_time::seconds(u1) + boost::posix_time::microsec(u2));
1295
1296 Time mjdTime(unixTimeT);
1297 yValues[i-1] = mjdTime.Mjd();
1298 if (yValues[i-1] < 40587)
1299 yValues[i-1] += 40587;
1300 }
1301 else
1302 {
1303 str >> yValues[i-1];
1304 if (yValues[i-1] < 40587)
1305 yValues[i-1] += 40587;
1306
1307 Time t(yValues[i-1]);
1308 string time = t.GetAsStr("%H:%M:%S%F");
1309 while (time[time.size()-1] == '0' && time.size() > 2)
1310 {
1311 time = time.substr(0, time.size()-1);
1312 }
1313 }
1314 if (i==1)
1315 {
1316 Time t(yValues[0]);
1317 title += " - " + t.GetAsStr("%Y-%m-%d");
1318 plot->setTitle(title.c_str());
1319 }
1320 }
1321 //set the actual data.
1322 auto jt = xValues.begin();
1323 for (auto it=curves.begin(); it != curves.end(); it++, jt++)
1324 (*it)->setRawData(yValues, *jt, endIndex-1);
1325
1326 QStack<QRectF> stack;
1327 double minX, minY, maxX, maxY;
1328 minX = minY = 1e10;
1329 maxX = maxY = -1e10;
1330 QRectF rect;
1331 QPointF point;
1332 for (auto it=curves.begin(); it!= curves.end(); it++)
1333 {
1334 rect = (*it)->boundingRect();
1335 point = rect.bottomRight();
1336 if (point.x() < minX) minX = point.x();
1337 if (point.y() < minY) minY = point.y();
1338 if (point.x() > maxX) maxX = point.x();
1339 if (point.y() > maxY) maxY = point.y();
1340 point = rect.topLeft();
1341 if (point.x() < minX) minX = point.x();
1342 if (point.y() < minY) minY = point.y();
1343 if (point.x() > maxX) maxX = point.x();
1344 if (point.y() > maxY) maxY = point.y();
1345 }
1346 QPointF bottomRight(maxX, minY);
1347 QPointF topLeft(minX, maxY);
1348 QPointF center((bottomRight+topLeft)/2.f);
1349 stack.push(QRectF(topLeft + (topLeft-center)*(.5f),bottomRight + (bottomRight-center)*(.5f)));
1350 zoom.setZoomStack(stack);
1351
1352// delete[] fitsBuffer;
1353 for (auto it = columnsData.begin(); it != columnsData.end(); it++)
1354 delete[] it->second;
1355 window.resize(600, 400);
1356 window.show();
1357
1358 a.exec();
1359
1360
1361 for (auto it = curves.begin(); it != curves.end(); it++)
1362 {
1363 (*it)->detach();
1364 delete *it;
1365 }
1366 grid->detach();
1367 for (auto it = xValues.begin(); it != xValues.end(); it++)
1368 delete[] *it;
1369
1370 delete[] yValues;
1371
1372 delete[] handle;
1373
1374 return 0;
1375}
1376#endif
1377
1378void SetupConfiguration(Configuration& conf)
1379{
1380 po::options_description configs("Fitsdump options");
1381 configs.add_options()
1382 ("fitsfile,f", var<string>()
1383#if BOOST_VERSION >= 104200
1384 ->required()
1385#endif
1386 , "Name of FITS file")
1387 ("tablename,t", var<string>("DATA")
1388#if BOOST_VERSION >= 104200
1389 ->required()
1390#endif
1391 , "Name of input table")
1392 ("col,c", vars<string>(), "List of columns to dump\narg is a list of columns, separated by a space.\nAdditionnally, a list of sub-columns can be added\ne.g. Data[3] will dump sub-column 3 of column Data\nData[3:4] will dump sub-columns 3 and 4\nOmitting this argument dump the entire column\nnota: all indices start at zero")
1393 ("outfile,o", var<string>("/dev/stdout"), "Name of output file (-:/dev/stdout)")
1394 ("precision,p", var<int>(20), "Precision of ofstream")
1395 ("list,l", po_switch(), "List all tables and columns in file")
1396 ("header,h", po_switch(), "Dump header of given table")
1397 ("stat,s", po_switch(), "Perform statistics instead of dump")
1398 ("minmax,m", po_switch(), "Calculates min and max of data")
1399 ("nozero,z", po_switch(), "skip 0 values for stats")
1400 ("tstart,a", po_switch(), "Give the mjdStart from reading the file data")
1401 ("tstop,b", po_switch(), "Give the mjdStop from reading the file data")
1402 ("force", po_switch(), "Force reading the fits file even if END key is missing")
1403#ifdef PLOTTING_PLEASE
1404 ("graph,g", po_switch(), "Plot the columns instead of dumping them")
1405 ("dots,d", po_switch(), "Plot using dots instead of lines")
1406#endif
1407 ;
1408
1409 po::positional_options_description p;
1410 p.add("fitsfile", 1); // The first positional options
1411 p.add("col", -1); // All others
1412
1413 conf.AddOptions(configs);
1414 conf.SetArgumentPositions(p);
1415}
1416
1417int main(int argc, const char** argv)
1418{
1419 Configuration conf(argv[0]);
1420 conf.SetPrintUsage(PrintUsage);
1421 SetupConfiguration(conf);
1422
1423 if (!conf.DoParse(argc, argv, PrintHelp))
1424 return -1;
1425
1426 FitsDumper loader;
1427 return loader.ExecConfig(conf);
1428}
Note: See TracBrowser for help on using the repository browser.