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

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