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

Last change on this file since 12733 was 12729, checked in by lyard, 13 years ago
added more options
File size: 42.3 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#include "Time.h"
30#endif
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 while (fFile->GetNextRow() && row < numRows)
633 {
634 rangesIt = ranges.begin();
635 auto statsIt = statData.begin();
636 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++)
637 {
638 double cValue = 0;
639 for (int i=rangesIt->first; i<rangesIt->second; i++)
640 {
641 switch (it->first) {
642 case 'L':
643 cValue = reinterpret_cast<bool*>(it->second)[i];
644 break;
645 case 'B':
646 cValue = reinterpret_cast<bool*>(it->second)[i];
647 break;
648 case 'I':
649 cValue = reinterpret_cast<int16_t*>(it->second)[i];
650 break;
651 case 'J':
652 cValue = reinterpret_cast<int32_t*>(it->second)[i];
653 break;
654 case 'K':
655 cValue = reinterpret_cast<int64_t*>(it->second)[i];
656 break;
657 case 'E':
658 cValue = reinterpret_cast<float*>(it->second)[i];
659 break;
660 case 'D':
661 cValue = reinterpret_cast<double*>(it->second)[i];
662 break;
663 default:
664 ;
665 }
666 if (!fNoZeroPlease || cValue != 0)
667 {
668 statsIt->average += cValue;
669 if (cValue < statsIt->min)
670 statsIt->min = cValue;
671 if (cValue > statsIt->max)
672 statsIt->max = cValue;
673 statsIt->numValues++;
674 }
675 }
676 }
677 row++;
678 }
679 for (auto it = columnsData.begin(); it != columnsData.end(); it++)
680 delete[] it->second;
681
682 //okay. So now I've got ALL the data, loaded.
683 //let's do the summing and averaging in a safe way (i.e. avoid overflow of variables as much as possible)
684 rangesIt = ranges.begin();
685 auto statsIt = statData.begin();
686
687 auto nameIt = listNamesOnly.begin();
688 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++, nameIt++)
689 {
690 int span = rangesIt->second - rangesIt->first;
691 cout << *nameIt << ": " << endl;
692 if (statsIt->numValues != 0)
693 {
694 statsIt->average /= statsIt->numValues;
695 out << "min: " << statsIt->min << endl;
696 out << "max: " << statsIt->max << endl;
697 out << "mea: " << statsIt->average << endl;
698 }
699 else
700 {
701 out << "min: 0" << endl << "max: 0" << endl << "mea: " << endl;
702 }
703
704 }
705 return true;
706}
707
708/*
709void FitsDumper::doTBoundary(const string& filename, int precision, bool tStop)
710{
711
712 //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
713 vector<pair<int, int> > ranges;
714 vector<string> listNamesOnly;
715
716 if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
717 {
718 cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
719 return false;
720 }
721
722 ofstream out(filename=="-"?"/dev/stdout":filename);
723 if (!out)
724 {
725 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
726 return false;
727 }
728
729 out.precision(precision);
730
731 // Loop over all columns in our list of requested columns
732 vector<pair<char, char*> > columnsData;
733 vector<minMaxStruct> statData;
734 int numRows = fFile->GetInt("NAXIS2");
735 auto rangesIt = ranges.begin();
736 for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
737 {
738 fits::Table::Column& cCol = fColMap.find(*it)->second;
739 columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
740 // minMaxStuct initData;
741 statData.push_back(minMaxStruct());
742 fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
743 }
744
745 int row = 0;
746 while (fFile->GetNextRow() && row < numRows)
747 {
748 rangesIt = ranges.begin();
749 auto statsIt = statData.begin();
750 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++)
751 {
752 double cValue = 0;
753 for (int i=rangesIt->first; i<rangesIt->second; i++)
754 {
755 switch (it->first) {
756 case 'L':
757 cValue = reinterpret_cast<bool*>(it->second)[i];
758 break;
759 case 'B':
760 cValue = reinterpret_cast<bool*>(it->second)[i];
761 break;
762 case 'I':
763 cValue = reinterpret_cast<int16_t*>(it->second)[i];
764 break;
765 case 'J':
766 cValue = reinterpret_cast<int32_t*>(it->second)[i];
767 break;
768 case 'K':
769 cValue = reinterpret_cast<int64_t*>(it->second)[i];
770 break;
771 case 'E':
772 cValue = reinterpret_cast<float*>(it->second)[i];
773 break;
774 case 'D':
775 cValue = reinterpret_cast<double*>(it->second)[i];
776 break;
777 default:
778 ;
779 }
780 if (!fNoZeroPlease || cValue != 0)
781 {
782 statsIt->average += cValue;
783 if (cValue < statsIt->min)
784 statsIt->min = cValue;
785 if (cValue > statsIt->max)
786 statsIt->max = cValue;
787 statsIt->numValues++;
788 }
789 }
790 }
791 row++;
792 }
793 for (auto it = columnsData.begin(); it != columnsData.end(); it++)
794 delete[] it->second;
795
796 //okay. So now I've got ALL the data, loaded.
797 //let's do the summing and averaging in a safe way (i.e. avoid overflow of variables as much as possible)
798 rangesIt = ranges.begin();
799 auto statsIt = statData.begin();
800
801 auto nameIt = listNamesOnly.begin();
802 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++, nameIt++)
803 {
804 int span = rangesIt->second - rangesIt->first;
805 cout << *nameIt << ": " << endl;
806 if (statsIt->numValues != 0)
807 {
808 statsIt->average /= statsIt->numValues;
809 out << "min: " << statsIt->min << endl;
810 out << "max: " << statsIt->max << endl;
811 out << "mea: " << statsIt->average << endl;
812 }
813 else
814 {
815 out << "min: 0" << endl << "max: 0" << endl << "mea: " << endl;
816 }
817
818 }
819 return true;
820
821}
822*/
823template<typename T>
824void displayStats(char* array, int numElems, ofstream& out)
825{
826 if (numElems == 0)
827 {
828 out << "Min: 0" << endl;
829 out << "Max: 0" << endl;
830 out << "Med: 0" << endl;
831 out << "Mea: 0" << endl;
832 return;
833 }
834 sort(reinterpret_cast<T*>(array), &reinterpret_cast<T*>(array)[numElems]);
835 out << "Min: " << reinterpret_cast<T*>(array)[0] << endl;
836 out << "Max: " << reinterpret_cast<T*>(array)[numElems-1] << endl;
837 if (numElems%2 == 0)
838 out << "Med: " << (reinterpret_cast<T*>(array)[numElems/2] + reinterpret_cast<int16_t*>(array)[numElems/2+1])/2.f << endl;
839 else
840 out << "Med: " << reinterpret_cast<T*>(array)[numElems/2+1] << endl;
841
842 double average = 0;
843 for (int i=0;i<numElems;i++)
844 average += reinterpret_cast<T*>(array)[i];
845 out << "Mea: " << average/(double)numElems << endl;
846
847}
848int FitsDumper::doStatsPlease(const string &filename, const vector<string>& list, int precision)
849{
850
851 //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
852 vector<pair<int, int> > ranges;
853 vector<string> listNamesOnly;
854
855 if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
856 {
857 cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
858 return false;
859 }
860
861 ofstream out(filename=="-"?"/dev/stdout":filename);
862 if (!out)
863 {
864 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
865 return false;
866 }
867
868 out.precision(precision);
869
870 // Loop over all columns in our list of requested columns
871 vector<pair<char, char*> > columnsData;
872 vector<char*> statData;
873 int numRows = fFile->GetInt("NAXIS2");
874 auto rangesIt = ranges.begin();
875 for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
876 {
877 fits::Table::Column& cCol = fColMap.find(*it)->second;
878 columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
879 statData.push_back(new char[(rangesIt->second - rangesIt->first)*cCol.size*numRows]);
880 fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
881 }
882
883 int row = 0;
884 while (fFile->GetNextRow() && row < numRows)
885 {
886 rangesIt = ranges.begin();
887 auto statsIt = statData.begin();
888 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++)
889 {
890 int span = rangesIt->second - rangesIt->first;
891 for (int i=rangesIt->first; i<rangesIt->second; i++)
892 {
893 switch (it->first) {
894 case 'L':
895 reinterpret_cast<bool*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<bool*>(it->second)[i];
896 break;
897 case 'B':
898 reinterpret_cast<char*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<char*>(it->second)[i];
899 break;
900 case 'I':
901 reinterpret_cast<int16_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<int16_t*>(it->second)[i];
902 break;
903 case 'J':
904 reinterpret_cast<int32_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<int32_t*>(it->second)[i];
905 break;
906 case 'K':
907 reinterpret_cast<int64_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<int64_t*>(it->second)[i];
908 break;
909 case 'E':
910 reinterpret_cast<float*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<float*>(it->second)[i];
911 break;
912 case 'D':
913 reinterpret_cast<double*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<double*>(it->second)[i];
914 break;
915 default:
916 ;
917 }
918 }
919 }
920 row++;
921 }
922 for (auto it = columnsData.begin(); it != columnsData.end(); it++)
923 delete[] it->second;
924
925 //okay. So now I've got ALL the data, loaded.
926 //let's do the summing and averaging in a safe way (i.e. avoid overflow of variables as much as possible)
927 rangesIt = ranges.begin();
928 auto statsIt = statData.begin();
929
930 auto nameIt = list.begin();
931 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++, nameIt++)
932 {
933 int span = rangesIt->second - rangesIt->first;
934 out << *nameIt << ": " << endl;
935 switch (it->first) {
936 case 'L':
937 displayStats<bool>(*statsIt, numRows*span, out);
938 break;
939 case 'B':
940 displayStats<char>(*statsIt, numRows*span, out);
941 break;
942 case 'I':
943 displayStats<int16_t>(*statsIt, numRows*span, out);
944 break;
945 case 'J':
946 displayStats<int32_t>(*statsIt, numRows*span, out);
947 break;
948 case 'K':
949 displayStats<int64_t>(*statsIt, numRows*span, out);
950 break;
951 case 'E':
952 displayStats<float>(*statsIt, numRows*span, out);
953 break;
954 case 'D':
955 displayStats<double>(*statsIt, numRows*span, out);
956 break;
957 default:
958 ;
959 }
960 }
961 return true;
962}
963#ifdef PLOTTING_PLEASE
964int FitsDumper::doCurvesDisplay( const vector<string> &list, const string& tableName)
965{
966 //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
967 vector<pair<int, int> > ranges;
968 vector<string> listNamesOnly;
969 if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
970 {
971 cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
972 return false;
973 }
974 vector<string> curvesNames;
975 stringstream str;
976 for (auto it=ranges.begin(), jt=listNamesOnly.begin(); it != ranges.end(); it++, jt++)
977 {
978 for (int i=it->first; i<it->second;i++)
979 {
980 str.str("");
981 str << *jt << "[" << i << "]";
982 curvesNames.push_back(str.str());
983 }
984 }
985 char* handle = new char[17];
986 sprintf(handle,"FitsDump Display");
987// Qt::HANDLE h = *handle;//NULL
988 int argc = 1;
989 char ** argv = &handle;
990 QApplication a(argc, argv);
991
992
993
994 QwtPlot* plot = new QwtPlot();
995 QwtPlotGrid* grid = new QwtPlotGrid;
996 grid->enableX(false);
997 grid->enableY(true);
998 grid->enableXMin(false);
999 grid->enableYMin(false);
1000 grid->setMajPen(QPen(Qt::black, 0, Qt::DotLine));
1001 grid->attach(plot);
1002 plot->setAutoReplot(true);
1003 string title = tableName;
1004 plot->setAxisScaleDraw( QwtPlot::xBottom, new TimeScaleDraw());
1005
1006 QWidget window;
1007 QHBoxLayout* layout = new QHBoxLayout(&window);
1008 layout->setContentsMargins(0,0,0,0);
1009 layout->addWidget(plot);
1010
1011 QwtPlotZoomer zoom(plot->canvas());
1012 zoom.setRubberBandPen(QPen(Qt::gray, 2, Qt::DotLine));
1013 zoom.setTrackerPen(QPen(Qt::gray));
1014 int totalSize = 0;
1015 for (unsigned int i=0;i<list.size();i++)
1016 totalSize += ranges[i].second - ranges[i].first;
1017
1018 vector<QwtPlotCurve*> curves(totalSize);
1019 int ii=0;
1020 for (auto it = curves.begin(), jt=curvesNames.begin(); it != curves.end(); it++, jt++)
1021 {
1022 *it = new QwtPlotCurve(jt->c_str());
1023 switch (ii%6)
1024 {
1025 case 0:
1026 (*it)->setPen(QColor(255,0,0));
1027 break;
1028 case 1:
1029 (*it)->setPen(QColor(0,255,0));
1030 break;
1031 case 2:
1032 (*it)->setPen(QColor(0,0,255));
1033 break;
1034 case 3:
1035 (*it)->setPen(QColor(255,255,0));
1036 break;
1037 case 4:
1038 (*it)->setPen(QColor(0,255,255));
1039 break;
1040 case 5:
1041 (*it)->setPen(QColor(255,0,255));
1042 break;
1043 default:
1044 (*it)->setPen(QColor(0,0,0));
1045 };
1046 ii++;
1047 if (fDotsPlease)
1048 (*it)->setStyle(QwtPlotCurve::Dots);
1049 else
1050 (*it)->setStyle(QwtPlotCurve::Lines);
1051 (*it)->attach(plot);
1052 }
1053 plot->insertLegend(new QwtLegend(), QwtPlot::RightLegend);
1054
1055
1056 vector<pair<char, char*> > columnsData;
1057
1058 for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++)
1059 {
1060 fits::Table::Column& cCol = fColMap.find(*it)->second;
1061 columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
1062 fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
1063 }
1064
1065 //add the time column to the given columns
1066 if (fColMap.find("Time") == fColMap.end() && fColMap.find("UnixTimeUTC") == fColMap.end())
1067 {
1068 cerr << "Error: time column could not be found in given table. Aborting" << endl;
1069 return false;
1070 }
1071 const fits::Table::Column& timeCol = (fColMap.find("Time") != fColMap.end()) ? fColMap.find("Time")->second : fColMap.find("UnixTimeUTC")->second;
1072 bool unixTime = (fColMap.find("Time") == fColMap.end());
1073 if (unixTime)
1074 ranges.push_back(make_pair(0,2));
1075 else
1076 ranges.push_back(make_pair(0,1));
1077 columnsData.push_back(make_pair(timeCol.type, new char[timeCol.num*timeCol.size]));
1078 fFile->SetPtrAddress(unixTime ? "UnixTimeUTC" : "Time", columnsData[columnsData.size()-1].second);
1079
1080// stringstream str;
1081 str.str("");
1082
1083
1084 vector<double*> xValues(totalSize);
1085 double* yValues;
1086 cout.precision(10);
1087 str.precision(20);
1088 for (auto it=xValues.begin(); it!=xValues.end(); it++)
1089 *it = new double[fFile->GetInt("NAXIS2")];
1090
1091 yValues = new double[fFile->GetInt("NAXIS2")];
1092
1093 cout.precision(3);
1094 int endIndex = 0;
1095 int numRows = fFile->GetInt("NAXIS2");
1096 for (int i=1;i<numRows;i++)
1097 {
1098 fFile->GetNextRow();
1099 cout << "\r" << "Constructing graph " << ((float)(endIndex)/(float)(fFile->GetInt("NAXIS2")))*100.0 << "%";
1100 endIndex++;
1101 auto rangesIt = ranges.begin();
1102 for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++)
1103 {
1104 for (int j=rangesIt->first; j<rangesIt->second; j++)
1105 {
1106 switch (it->first) {
1107 case 'L':
1108 str << reinterpret_cast<bool*>(it->second)[j] << " ";
1109 break;
1110 case 'B':
1111 str << reinterpret_cast<char*>(it->second)[j] << " ";
1112 break;
1113 case 'I':
1114 str << reinterpret_cast<int16_t*>(it->second)[j] << " ";
1115 break;
1116 case 'J':
1117 str << reinterpret_cast<int32_t*>(it->second)[j] << " ";
1118 break;
1119 case 'K':
1120 str << reinterpret_cast<int64_t*>(it->second)[j] << " ";
1121 break;
1122 case 'E':
1123 str << reinterpret_cast<float*>(it->second)[j] << " ";
1124 break;
1125 case 'D':
1126 str << reinterpret_cast<double*>(it->second)[j] << " ";
1127 break;
1128 default:
1129 ;
1130 }
1131 }
1132 }
1133
1134 for (auto it=xValues.begin(); it!= xValues.end(); it++)
1135 {
1136 str >> (*it)[i-1];
1137 }
1138 if (unixTime)
1139 {
1140 long u1, u2;
1141 str >> u1 >> u2;
1142
1143 boost::posix_time::ptime unixTimeT( boost::gregorian::date(1970, boost::gregorian::Jan, 1),
1144 boost::posix_time::seconds(u1) + boost::posix_time::microsec(u2));
1145
1146 Time mjdTime(unixTimeT);
1147 yValues[i-1] = mjdTime.Mjd();
1148 }
1149 else
1150 {
1151 str >> yValues[i-1];
1152 if (yValues[i-1] < 40587)
1153 yValues[i-1] += 40587;
1154
1155 Time t(yValues[i-1]);
1156 string time = t.GetAsStr("%H:%M:%S%F");
1157 while (time[time.size()-1] == '0' && time.size() > 2)
1158 {
1159 time = time.substr(0, time.size()-1);
1160 }
1161 }
1162 if (i==1)
1163 {
1164 Time t(yValues[0]);
1165 title += " - " + t.GetAsStr("%Y-%m-%d");
1166 plot->setTitle(title.c_str());
1167 }
1168 }
1169 //set the actual data.
1170 auto jt = xValues.begin();
1171 for (auto it=curves.begin(); it != curves.end(); it++, jt++)
1172 (*it)->setRawData(yValues, *jt, endIndex-1);
1173
1174 QStack<QRectF> stack;
1175 double minX, minY, maxX, maxY;
1176 minX = minY = 1e10;
1177 maxX = maxY = -1e10;
1178 QRectF rect;
1179 QPointF point;
1180 for (auto it=curves.begin(); it!= curves.end(); it++)
1181 {
1182 rect = (*it)->boundingRect();
1183 point = rect.bottomRight();
1184 if (point.x() < minX) minX = point.x();
1185 if (point.y() < minY) minY = point.y();
1186 if (point.x() > maxX) maxX = point.x();
1187 if (point.y() > maxY) maxY = point.y();
1188 point = rect.topLeft();
1189 if (point.x() < minX) minX = point.x();
1190 if (point.y() < minY) minY = point.y();
1191 if (point.x() > maxX) maxX = point.x();
1192 if (point.y() > maxY) maxY = point.y();
1193 }
1194 QPointF bottomRight(maxX, minY);
1195 QPointF topLeft(minX, maxY);
1196 QPointF center((bottomRight+topLeft)/2.f);
1197 stack.push(QRectF(topLeft + (topLeft-center)*(.5f),bottomRight + (bottomRight-center)*(.5f)));
1198 zoom.setZoomStack(stack);
1199
1200// delete[] fitsBuffer;
1201 for (auto it = columnsData.begin(); it != columnsData.end(); it++)
1202 delete[] it->second;
1203 window.resize(600, 400);
1204 window.show();
1205
1206 a.exec();
1207
1208
1209 for (auto it = curves.begin(); it != curves.end(); it++)
1210 {
1211 (*it)->detach();
1212 delete *it;
1213 }
1214 grid->detach();
1215 for (auto it = xValues.begin(); it != xValues.end(); it++)
1216 delete[] *it;
1217
1218 delete[] yValues;
1219
1220 delete[] handle;
1221
1222 return 0;
1223}
1224#endif
1225
1226void SetupConfiguration(Configuration& conf)
1227{
1228 po::options_description configs("Fitsdump options");
1229 configs.add_options()
1230 ("fitsfile,f", var<string>()
1231#if BOOST_VERSION >= 104200
1232 ->required()
1233#endif
1234 , "Name of FITS file")
1235 ("tablename,t", var<string>("DATA")
1236#if BOOST_VERSION >= 104200
1237 ->required()
1238#endif
1239 , "Name of input table")
1240 ("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")
1241 ("outfile,o", var<string>("/dev/stdout"), "Name of output file (-:/dev/stdout)")
1242 ("precision,p", var<int>(20), "Precision of ofstream")
1243 ("list,l", po_switch(), "List all tables and columns in file")
1244 ("header,h", po_switch(), "Dump header of given table")
1245 ("stat,s", po_switch(), "Perform statistics instead of dump")
1246 ("minmax,m", po_switch(), "Calculates min and max of data")
1247 ("nozero,z", po_switch(), "skip 0 values for stats")
1248 ("tstart,a", po_switch(), "Give the mjdStart from reading the file data")
1249 ("tstop,b", po_switch(), "Give the mjdStop from reading the file data")
1250#ifdef PLOTTING_PLEASE
1251 ("graph,g", po_switch(), "Plot the columns instead of dumping them")
1252 ("dots,d", po_switch(), "Plot using dots instead of lines")
1253#endif
1254 ;
1255
1256 po::positional_options_description p;
1257 p.add("fitsfile", 1); // The first positional options
1258 p.add("col", -1); // All others
1259
1260 conf.AddOptions(configs);
1261 conf.SetArgumentPositions(p);
1262}
1263
1264int main(int argc, const char** argv)
1265{
1266 Configuration conf(argv[0]);
1267 conf.SetPrintUsage(PrintUsage);
1268 SetupConfiguration(conf);
1269
1270 if (!conf.DoParse(argc, argv, PrintHelp))
1271 return -1;
1272
1273 FitsDumper loader;
1274 return loader.ExecConfig(conf);
1275}
Note: See TracBrowser for help on using the repository browser.