Index: trunk/FACT++/src/fitsdump.cc =================================================================== --- trunk/FACT++/src/fitsdump.cc (revision 12716) +++ trunk/FACT++/src/fitsdump.cc (revision 12723) @@ -11,7 +11,7 @@ #include -#include - -//#define PLOTTING_PLEASE +#include "externals/fits.h" + +#define PLOTTING_PLEASE #ifdef PLOTTING_PLEASE @@ -32,15 +32,4 @@ using namespace std; -class MyColumn : public CCfits::Column -{ -public: - const string &comment() const { return CCfits::Column::comment(); } - long width() const - {//the width() method seems to be buggy (or empty ?) as it always return 1... redo it ourselves. - string inter = format(); - inter = inter.substr(0, inter.size()-1); - return atoi(inter.c_str()); - } -}; #ifdef PLOTTING_PLEASE @@ -68,17 +57,19 @@ private: - - CCfits::FITS *fFile; /// FITS pointer - CCfits::Table *fTable; /// Table pointer - map fColMap; /// map between the column names and their CCfits objects + bool fDotsPlease; + fits* fFile; + string fFilename; + + fits::Table::Columns fColMap; + fits::Table::Keys fKeyMap; // Convert CCfits::ValueType into a human readable string - string ValueTypeToStr(CCfits::ValueType type) const; + string ValueTypeToStr(char type) const; // Convert CCfits::ValueType into a number of associated bytes - int ValueTypeToSize(CCfits::ValueType type) const; + int ValueTypeToSize(char type) const; /// Calculate the buffer size required to read a row of the fits table, as well as the offsets to each column - vector CalculateOffsets() const; +// vector CalculateOffsets() const; template @@ -89,5 +80,5 @@ /// Write a single row of the selected data - int WriteRow(ostream &, const vector &, const vector &, unsigned char *, const vector >&) const; + // int WriteRow(ostream &, const vector &, const vector &, unsigned char *, const vector >&) const; bool OpenFile(const string &); /// Open a file @@ -108,4 +99,5 @@ int doCurvesDisplay( const vector &list, const string& tableName); #endif + int doStatsPlease(const string &filename, const vector& list, int precision); // bool Plot(const vector &list); @@ -121,5 +113,5 @@ //! the ostream where to redirect the outputs // -FitsDumper::FitsDumper() : fFile(0), fTable(0) +FitsDumper::FitsDumper() : fFile(0), fDotsPlease(false) { } @@ -136,19 +128,15 @@ -string FitsDumper::ValueTypeToStr(CCfits::ValueType type) const +string FitsDumper::ValueTypeToStr(char type) const { switch (type) { - case CCfits::Tbyte: return "uint8_t"; - case CCfits::Tushort: return "uint16_t"; - case CCfits::Tshort: return "int16_t"; - case CCfits::Tuint: return "uint32_t"; - case CCfits::Tint: return "int32_t"; - case CCfits::Tulong: return "uint32_t"; - case CCfits::Tlong: return "int32_t"; - case CCfits::Tlonglong: return "int64_t"; - case CCfits::Tfloat: return "float"; - case CCfits::Tdouble: return "double"; - + case 'L': return "bool(8)"; + case 'B': return "byte(8)"; + case 'I': return "short(16)"; + case 'J': return "int(32)"; + case 'K': return "int(64)"; + case 'E': return "float(32)"; + case 'D': return "double(64)"; default: return "unknwown"; @@ -156,19 +144,15 @@ } -int FitsDumper::ValueTypeToSize(CCfits::ValueType type) const +int FitsDumper::ValueTypeToSize(char type) const { switch (type) { - case CCfits::Tbyte: return sizeof(uint8_t); - case CCfits::Tushort: - case CCfits::Tshort: return sizeof(uint16_t); - case CCfits::Tuint: - case CCfits::Tint: - case CCfits::Tulong: - case CCfits::Tlong: return sizeof(uint32_t); - case CCfits::Tlonglong: return sizeof(uint64_t); - case CCfits::Tfloat: return sizeof(float); - case CCfits::Tdouble: return sizeof(double); - + case 'L': return sizeof(uint8_t); + case 'B': return sizeof(int8_t); + case 'I': return sizeof(int16_t); + case 'J': return sizeof(int32_t); + case 'K': return sizeof(int64_t); + case 'E': return sizeof(float); + case 'D': return sizeof(double); default: return 0; @@ -185,175 +169,4 @@ return t; } -/* -template -double FitsDumper::PtrToDouble(const unsigned char *ptr) const -{ - T t; - reverse_copy(ptr, ptr+sizeof(T), reinterpret_cast(&t)); - return t; -} - -double FitsDumper::PtrToDouble(const unsigned char *ptr, CCfits::ValueType type) const -{ - switch (type) - { - case CCfits::Tbyte: return PtrToDouble (ptr); - case CCfits::Tushort: return PtrToDouble(ptr); - case CCfits::Tuint: - case CCfits::Tulong: return PtrToDouble(ptr); - case CCfits::Tshort: return PtrToDouble (ptr); - case CCfits::Tint: - case CCfits::Tlong: return PtrToDouble (ptr); - case CCfits::Tlonglong: return PtrToDouble (ptr); - case CCfits::Tfloat: return PtrToDouble (ptr); - case CCfits::Tdouble: return PtrToDouble (ptr); - - default: - throw runtime_error("Data type not implemented yet."); - } -} -*/ - -// -------------------------------------------------------------------------- -// -//! Writes a single row of the selected FITS data to the output file. -//! -//! Data type \b not yet implemented: -//! CCfits::Tnull, CCfits::Tbit, CCfits::Tlogical, CCfits::Tstring, -//! CCfits::Tcomplex, CCfits::Tdblcomplex, CCfits::VTbit, -//! CCfits::VTbyte, CCfits::VTlogical, CCfits::VTushort, -//! CCfits::VTshort, CCfits::VTuint, CCfits::VTint, CCfits::VTulong, -//! CCfits::VTlong, CCfits::VTlonglong, CCfits::VTfloat, -//! CCfits::VTdouble, CCfits::VTcomplex, CCfits::VTdblcomplex -//! -//! @param offsets -//! a vector containing the offsets to the columns (in bytes) -//! @param targetFile -//! the ofstream where to write to -//! @param fitsBuffer -//! the memory were the row has been loaded by cfitsio -// -int FitsDumper::WriteRow(ostream &out, const vector &list, const vector &offsets, unsigned char* fitsBuffer, const vector >& ranges) const -{ - int cnt = 0; - vector >::const_iterator jt = ranges.begin(); - for (vector::const_iterator it=list.begin(); it!=list.end(); it++, jt++) - { - if (jt == ranges.end()) - { - cout << "ERROR: END OF RANGE POINTER REACHED" << endl; - return false; - } - const MyColumn *col = *it; - - // CCfits starts counting at 1 not 0 - const int offset = offsets[col->index()-1]; - - // Get the pointer to the array which we are supposed to print - const unsigned char *ptr = fitsBuffer + offset; - - // Loop over all array entries - int sizeToSkip = 0; - switch (col->type()) - { - case CCfits::Tbyte: sizeToSkip = sizeof(uint8_t); break; - case CCfits::Tushort: sizeToSkip = sizeof(uint16_t); break; - case CCfits::Tuint: - case CCfits::Tulong: sizeToSkip = sizeof(uint32_t); break; - case CCfits::Tshort: sizeToSkip = sizeof(int16_t); break; - case CCfits::Tint: - case CCfits::Tlong: sizeToSkip = sizeof(int32_t); break; - case CCfits::Tlonglong: sizeToSkip = sizeof(int64_t); break; - case CCfits::Tfloat: sizeToSkip = sizeof(float); break; - case CCfits::Tdouble: sizeToSkip = sizeof(double); break; - default: - cerr << "Data type not implemented yet." << endl; - return 0; - } - ptr += sizeToSkip*jt->first; - for (int width=jt->first; widthsecond; width++) - { - switch (col->type()) - { - case CCfits::Tbyte: out << PtrToValue (ptr); break; - case CCfits::Tushort: out << PtrToValue(ptr); break; - case CCfits::Tuint: - case CCfits::Tulong: out << PtrToValue(ptr); break; - case CCfits::Tshort: out << PtrToValue (ptr); break; - case CCfits::Tint: - case CCfits::Tlong: out << PtrToValue (ptr); break; - case CCfits::Tlonglong: out << PtrToValue (ptr); break; - case CCfits::Tfloat: out << PtrToValue (ptr); break; - case CCfits::Tdouble: out << PtrToValue (ptr); break; - - default: - cerr << "Data type not implemented yet." << endl; - return 0; - } - - out << " "; - cnt++; - } - } - - if (cnt>0) - out << endl; - - if (out.fail()) - { - cerr << "ERROR - writing output: " << strerror(errno) << endl; - return -1; - } - - return cnt; -} - -// -------------------------------------------------------------------------- -// -//! Calculates the required buffer size for reading one row of the current -//! table. Also calculates the offsets to all the columns -// -vector FitsDumper::CalculateOffsets() const -{ - map sizes; - - for (map::const_iterator it=fColMap.begin(); - it!=fColMap.end(); it++) - { - const int &width = it->second->width(); - const int &idx = it->second->index(); - - const int size = ValueTypeToSize(it->second->type()); - if (size==0) - { - cerr << "Data type " << (int)it->second->type() << " not implemented yet." << endl; - return vector(); - } -//cout << "data size: " << size << " width: " << width << " index: " << idx << endl; - int realwidth = (width == 0)? 1 : width; - sizes[idx] = size*realwidth; - } - - //calculate the offsets in the vector. - vector result(1, 0); - - int size = 0; - int idx = 0; - - for (map::const_iterator it=sizes.begin(); it!=sizes.end(); it++) - { - size += it->second; - result.push_back(size); - - if (it->first == ++idx) - continue; - - cerr << "Expected index " << idx << ", but found " << it->first << endl; - return vector(); - } - - return result; -} - // -------------------------------------------------------------------------- // @@ -363,16 +176,29 @@ { if (fFile) + { + fFile->close(); delete fFile; - - ostringstream str; - try - { - fFile = new CCfits::FITS(filename); - } - catch (CCfits::FitsException e) - { - cerr << "Could not open FITS file " << filename << " reason: " << e.message() << endl; + } + + try { + fFile = new fits(filename); + } + catch (std::runtime_error e) + { + cout << "Something went wrong while trying to open " << filename; + cout << ": " << e.what() << " Aborting dump." << endl; return false; } + fFilename = filename; + + const fits::Table::Columns& tCols = fFile->getColumns(); + + for (auto it=tCols.begin(); it != tCols.end(); it++) + fColMap.insert(*it); + + const fits::Table::Keys& tkeys = fFile->getKeys(); + + for (auto it=tkeys.begin(); it != tkeys.end(); it++) + fKeyMap.insert(*it); return true; @@ -387,30 +213,4 @@ } - const multimap< string, CCfits::ExtHDU * > extMap = fFile->extension(); - if (extMap.find(tablename) == extMap.end()) - { - cerr << "Table '" << tablename << "' not found." << endl; - return false; - } - - fTable = dynamic_cast(extMap.find(tablename)->second); - if (!fTable) - { - cerr << "Object '" << tablename << "' returned not a CCfits::Table." << endl; - return false; - } - - map& tCols = fTable->column(); - for (map::const_iterator it = tCols.begin(); it != tCols.end(); it++) - { - fColMap.insert(make_pair(it->first, static_cast(it->second))); - } -// fColMap = fTable->column(); - - fTable->makeThisCurrent(); - - fTable->getComments(); - fTable->getHistory(); - fTable->readAllKeys(); return true; @@ -426,85 +226,26 @@ } - cout << "\nFile: " << fFile->name() << "\n"; - - const multimap< string, CCfits::ExtHDU * > extMap = fFile->extension(); - for (std::multimap::const_iterator it=extMap.begin(); it != extMap.end(); it++) - { - - CCfits::Table *table = dynamic_cast(extMap.find(it->first)->second); - - table->makeThisCurrent(); - table->readData(); - - cout << " " << it->first << " [" << table->rows() << "]\n"; - - const map &cols = table->column(); - -// for (map::const_iterator id = cols.begin(); ib != cols.end(); ib++) -// { -// TFORM -// } - - for (map::const_iterator ic=cols.begin(); - ic != cols.end(); ic++) - { - const MyColumn *col = static_cast(ic->second); - - cout << " " << col->name() << "[" ; - if (col->width() == 0) - cout << 1; - else - cout << col->width(); - cout << "] * " << col->scale() << " (" << col->unit() << ":" << ValueTypeToStr(col->type())<< ") " << col->comment() << "\n"; - /* - inline size_t Column::repeat () const - inline bool Column::varLength () const - inline double Column::zero () const - inline const String& Column::display () const - inline const String& Column::dimen () const - inline const String& Column::TBCOL () - inline const String& Column::TTYPE () - inline const String& Column::TFORM () - inline const String& Column::TDISP () - inline const String& Column::TZERO () - inline const String& Column::TDIM () - inline const String& Column::TNULL () - inline const String& Column::TLMIN () - inline const String& Column::TLMAX () - inline const String& Column::TDMAX () - inline const String& Column::TDMIN () - */ - } - cout << '\n'; - } + cout << "\nFile: " << fFilename << "\n"; + + cout << " " << fKeyMap.find("EXTNAME")->second.value << " ["; + cout << fKeyMap.find("NAXIS2")->second.value << "]\n"; + + for (auto it = fColMap.begin(); it != fColMap.end(); it++) + { + cout << " " << it->first << "[" << it->second.num << "] (" << it->second.unit << ":" << ValueTypeToStr(it->second.type) << ") "; + for (auto jt = fKeyMap.begin(); jt != fKeyMap.end(); jt++) + if (jt->second.value == it->first) + cout << jt->second.comment << endl; + } + + cout << endl; cout << flush; } -class MyKeyword : public CCfits::Keyword -{ -public: - CCfits::ValueType keytype() const { return CCfits::Keyword::keytype(); } -}; - void FitsDumper::ListKeywords(ostream &out) { - map keys = fTable->keyWord(); - for (map::const_iterator it=keys.begin(); - it!=keys.end(); it++) - { - string str; - double d; - int l; - - const MyKeyword *kw = static_cast(it->second); - kw->keytype(); - out << "## " << setw(8) << kw->name() << " = " << setw(10); - if (kw->keytype()==16) - out << ("'"+kw->value(str)+"'"); - if (kw->keytype()==31) - out << kw->value(l); - if (kw->keytype()==82) - out << kw->value(d); - out << " / " << kw->comment() << endl; + for (auto it=fKeyMap.begin(); it != fKeyMap.end(); it++) { + out << "## " << setw(8) << it->first << " = " << setw(10); + out << "'" << it->second.value << "'" << " / " << it->second.comment << endl; } } @@ -512,5 +253,5 @@ void FitsDumper::ListHeader() { - if (!fTable) + if (!fFile) { cerr << "No table open." << endl; @@ -518,12 +259,11 @@ } - cout << "\nTable: " << fTable->name() << " (rows=" << fTable->rows() << ")\n"; - if (!fTable->comment().empty()) - cout << "Comment: \t" << fTable->comment() << '\n'; - if (!fTable->history().empty()) - cout << "History: \t" << fTable->history() << '\n'; + cout << "\nTable: " << fKeyMap.find("EXTNAME")->second.value << " (rows=" << fKeyMap.find("NAXIS2")->second.value << ")\n"; + if (fKeyMap.find("COMMENT") != fKeyMap.end()) + cout << "Comment: \t" << fKeyMap.find("COMMENT")->second.value << "\n"; ListKeywords(cout); cout << endl; + } @@ -538,5 +278,4 @@ unsigned long bracketIndex1 = columnNameOnly.find_first_of(']'); unsigned long colonIndex = columnNameOnly.find_first_of(':'); -// cout << bracketIndex0 << " " << bracketIndex1 << " " << colonIndex << endl; int columnStart = -1; int columnEnd = -1; @@ -553,5 +292,4 @@ columnStart = atoi(columnNameOnly.substr(bracketIndex0+1, bracketIndex1 - (bracketIndex0+1)).c_str()); columnEnd = columnStart+1; -// cout << "Cstart " << columnStart << " end: " << columnEnd << endl; } columnNameOnly = columnNameOnly.substr(0, bracketIndex0); @@ -563,9 +301,8 @@ return false; } -// cout << "The column name is: " << columnNameOnly << endl; - MyColumn *cCol = static_cast(fColMap.find(columnNameOnly)->second); + fits::Table::Column& cCol = fColMap.find(columnNameOnly)->second; if (bracketIndex0 == string::npos) {//no range given: use the full range - ranges.push_back(make_pair(0, cCol->width())); + ranges.push_back(make_pair(0, cCol.num)); columnStart = 0; columnEnd = 1; @@ -578,7 +315,7 @@ return false; } - if (columnEnd>1 && columnEnd > cCol->width()) + if (columnEnd>1 && columnEnd > (int)(cCol.num)) { - cerr << "ERROR - End range for column " << columnNameOnly << " is greater than the last element (" << cCol->width() << " vs " << columnEnd << "). Aborting" << endl; + cerr << "ERROR - End range for column " << columnNameOnly << " is greater than the last element (" << cCol.num << " vs " << columnEnd << "). Aborting" << endl; return false; } @@ -607,12 +344,4 @@ } - // FIXME: Maybe do this when opening a table? - const vector offsets = CalculateOffsets(); - // for (int i=0;iname() << '\n'; - if (!fTable->comment().empty()) - out << "## Comment: \t" << fTable->comment() << '\n'; - if (!fTable->history().empty()) - out << "## History: \t" << fTable->history() << '\n'; - out << "## NumRows: \t" << fTable->rows() << '\n'; + out << "## Table: \t" << fKeyMap.find("EXTNAME")->second.value << '\n'; + out << "## Comment: \t" << ((fKeyMap.find("COMMENT") != fKeyMap.end()) ? fKeyMap.find("COMMENT")->second.value : "") << '\n'; + + out << "## NumRows: \t" << fFile->GetInt("NAXIS2") << '\n'; out << "## --------------------------------------------------------------------------\n"; ListKeywords(out); out << "## --------------------------------------------------------------------------\n"; out << "#\n"; - //vector >::const_iterator rangesIt = ranges.begin(); - //the above is soooo yesterday ;) let's try the new auto feature of c++0x + auto rangesIt = ranges.begin(); for (vector::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++) { - const MyColumn *col = static_cast(fTable->column()[*it]); - if (rangesIt->first != 0 || rangesIt->second != col->width()) + const fits::Table::Column& col = fColMap[*it]; +// const MyColumn *col = static_cast(fTable->column()[*it]); + if (rangesIt->first != 0 || rangesIt->second != col.num) { out << "#"; for (int i=rangesIt->first; isecond; i++) - out << " " << col->name() << "[" << i << "]"; - out << ": " << col->unit(); + out << " " << *it << "[" << i << "]"; + out << ": " << col.unit; } else - out << "# " << col->name() << "[" << col->width() << "]: " << col->unit(); - - if (!col->comment().empty()) - out << " (" <comment() << ")"; + out << "# " << *it << "[" << col.num << "]: " << col.unit; + +// FIXME: retrive the column comment +// if (!col->comment().empty()) +// out << " (" <comment() << ")"; out << '\n'; } @@ -661,108 +389,54 @@ // Loop over all columns in our list of requested columns - vector columns; + vector > columnsData; + for (vector::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++) { - MyColumn *cCol = static_cast(fColMap.find(*it)->second); - columns.push_back(cCol); - } - const int size = offsets[offsets.size()-1]; - unsigned char* fitsBuffer = new unsigned char[size]; - - int status = 0; - int endIndex = (fColMap.find("Time") == fColMap.end()) ? fTable->rows()-1 : fTable->rows(); - for (int i=1; i<=endIndex; i++) - { - fits_read_tblbytes(fFile->fitsPointer(), i, 1, size, fitsBuffer, &status); - if (status) + fits::Table::Column& cCol = fColMap.find(*it)->second; + columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size])); + fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second); + } + int numRows = fFile->GetInt("NAXIS2"); + int row = 0; + while (fFile->GetNextRow() && row < numRows) + { + row++; + rangesIt = ranges.begin(); + for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++) { - char text[30];//max length of cfitsio error strings (from doc) - fits_get_errstatus(status, text); - - cerr << "Reading row " << i << " failed: " << text << " (fits_insert_rows,rc=" << status << ")" << endl; - break; - } - if (WriteRow(out, columns, offsets, fitsBuffer, ranges)<0) - { - status=1; - break; - } - } - delete[] fitsBuffer; - - return status==0; -} - -// -------------------------------------------------------------------------- -// -//! Perform the actual dump, based on the current parameters -// -/* -bool FitsDumper::Plot(const vector &list) -{ - for (vector::const_iterator it=list.begin(); it!=list.end(); it++) - if (fColMap.find(*it) == fColMap.end()) - { - cerr << "WARNING - Column '" << *it << "' not found in table." << endl; - return false; - } - - // FIXME: Maybe do this when opening a table? - const vector offsets = CalculateOffsets(); - if (offsets.size()==0) - return false; - - // Loop over all columns in our list of requested columns - const CCfits::Column *col[3] = - { - list.size()>0 ? fColMap.find(list[0])->second : 0, - list.size()>1 ? fColMap.find(list[1])->second : 0, - list.size()>2 ? fColMap.find(list[2])->second : 0 - }; - - const int size = offsets[offsets.size()-1]; - unsigned char* fitsBuffer = new unsigned char[size]; - - const int idx = 0; - - // CCfits starts counting at 1 not 0 - const size_t pos[3] = - { - col[0] ? offsets[col[0]->index()-1] + idx*col[0]->width()*ValueTypeToSize(col[0]->type()) : 0, - col[1] ? offsets[col[1]->index()-1] + idx*col[1]->width()*ValueTypeToSize(col[1]->type()) : 0, - col[2] ? offsets[col[2]->index()-1] + idx*col[2]->width()*ValueTypeToSize(col[2]->type()) : 0 - }; - - int status = 0; - for (int i=1; i<=fTable->rows(); i++) - { - fits_read_tblbytes(fFile->fitsPointer(), i, 1, size, fitsBuffer, &status); - if (status) - { - cerr << "An error occurred while reading fits row #" << i << " error code: " << status << endl; - break; - } - - try - { - const double x[3] = + for (int i=rangesIt->first; isecond; i++) { - col[0] ? PtrToDouble(fitsBuffer+pos[0], col[0]->type()) : 0, - col[1] ? PtrToDouble(fitsBuffer+pos[1], col[1]->type()) : 0, - col[2] ? PtrToDouble(fitsBuffer+pos[2], col[2]->type()) : 0 - }; - } - catch (const runtime_error &e) - { - cerr << e.what() << endl; - status=1; - break; - } - } - delete[] fitsBuffer; - - return status==0; -} -*/ + switch (it->first) { + case 'L': + out << reinterpret_cast(it->second)[i] << " "; + break; + case 'B': + out << reinterpret_cast(it->second)[i] << " "; + break; + case 'I': + out << reinterpret_cast(it->second)[i] << " "; + break; + case 'J': + out << reinterpret_cast(it->second)[i] << " "; + break; + case 'K': + out << reinterpret_cast(it->second)[i] << " "; + break; + case 'E': + out << reinterpret_cast(it->second)[i] << " "; + break; + case 'D': + out << reinterpret_cast(it->second)[i] << " "; + break; + default: + ; + } + } + } + } + for (auto it = columnsData.begin(); it != columnsData.end(); it++) + delete[] it->second; + return true; +} // -------------------------------------------------------------------------- @@ -780,4 +454,10 @@ } + if (conf.Get("stat") && conf.Get("graph")) + { + cout << "Invalid conbination of options: cannot graph stats. Aborting" << endl; + return -1; + } + if (conf.Get("list")) List(); @@ -789,6 +469,5 @@ } -#ifdef PLOTTING_PLEASE - if (conf.Get("graph")) + if (conf.Get("stat")) { if (!conf.Has("col")) @@ -797,4 +476,18 @@ return 0; } + doStatsPlease(conf.Get("outfile"), conf.Get>("col"), conf.Get("precision")); + return 0; + } + +#ifdef PLOTTING_PLEASE + if (conf.Get("graph")) + { + if (!conf.Has("col")) + { + cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl; + return 0; + } + if (conf.Get("dots")) + fDotsPlease = true; doCurvesDisplay(conf.Get>("col"), conf.Get("tablename")); @@ -839,4 +532,137 @@ { // +} + +template +void displayStats(char* array, int numElems) +{ + sort(reinterpret_cast(array), &reinterpret_cast(array)[numElems]); + cout << "Min: " << reinterpret_cast(array)[0] << endl; + cout << "Max: " << reinterpret_cast(array)[numElems-1] << endl; + if (numElems%2 == 0) + cout << "Med: " << (reinterpret_cast(array)[numElems/2] + reinterpret_cast(array)[numElems/2+1])/2.f << endl; + else + cout << "Med: " << reinterpret_cast(array)[numElems/2+1] << endl; + + double average = 0; + for (int i=0;i(array)[i]; + cout << "Mea: " << average/(double)numElems << endl; + +} +int FitsDumper::doStatsPlease(const string &filename, const vector& list, int precision) +{ + + //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file + vector > ranges; + vector listNamesOnly; + + if (!separateColumnsFromRanges(list, ranges, listNamesOnly)) + { + cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl; + return false; + } + + ofstream out(filename=="-"?"/dev/stdout":filename); + if (!out) + { + cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl; + return false; + } + + out.precision(precision); + + // Loop over all columns in our list of requested columns + vector > columnsData; + vector statData; + int numRows = fFile->GetInt("NAXIS2"); + auto rangesIt = ranges.begin(); + for (vector::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++) + { + fits::Table::Column& cCol = fColMap.find(*it)->second; + columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size])); + statData.push_back(new char[(rangesIt->second - rangesIt->first)*cCol.size*numRows]); + fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second); + } + + int row = 0; + while (fFile->GetNextRow() && row < numRows) + { + rangesIt = ranges.begin(); + auto statsIt = statData.begin(); + for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++) + { + int span = rangesIt->second - rangesIt->first; + for (int i=rangesIt->first; isecond; i++) + { + switch (it->first) { + case 'L': + reinterpret_cast(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast(it->second)[i]; + break; + case 'B': + reinterpret_cast(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast(it->second)[i]; + break; + case 'I': + reinterpret_cast(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast(it->second)[i]; + break; + case 'J': + reinterpret_cast(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast(it->second)[i]; + break; + case 'K': + reinterpret_cast(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast(it->second)[i]; + break; + case 'E': + reinterpret_cast(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast(it->second)[i]; + break; + case 'D': + reinterpret_cast(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast(it->second)[i]; + break; + default: + ; + } + } + } + row++; + } + for (auto it = columnsData.begin(); it != columnsData.end(); it++) + delete[] it->second; + + //okay. So now I've got ALL the data, loaded. + //let's do the summing and averaging in a safe way (i.e. avoid overflow of variables as much as possible) + rangesIt = ranges.begin(); + auto statsIt = statData.begin(); + int round; + auto nameIt = listNamesOnly.begin(); + for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++, nameIt++) + { + int span = rangesIt->second - rangesIt->first; + cout << *nameIt << ": " << endl; + switch (it->first) { + case 'L': + displayStats(*statsIt, numRows*span); + break; + case 'B': + displayStats(*statsIt, numRows*span); + break; + case 'I': + displayStats(*statsIt, numRows*span); + break; + case 'J': + displayStats(*statsIt, numRows*span); + break; + case 'K': + displayStats(*statsIt, numRows*span); + break; + case 'E': + displayStats(*statsIt, numRows*span); + break; + case 'D': + displayStats(*statsIt, numRows*span); + break; + default: + ; + } + } + return true; } #ifdef PLOTTING_PLEASE @@ -894,5 +720,5 @@ for (unsigned int i=0;i curves(totalSize); int ii=0; @@ -924,21 +750,22 @@ }; ii++; - (*it)->setStyle(QwtPlotCurve::Lines); + if (fDotsPlease) + (*it)->setStyle(QwtPlotCurve::Dots); + else + (*it)->setStyle(QwtPlotCurve::Lines); (*it)->attach(plot); } plot->insertLegend(new QwtLegend(), QwtPlot::RightLegend); - const vector offsets = CalculateOffsets(); - if (offsets.size()==0) - return false; - - - // Loop over all columns in our list of requested columns - vector columns; + + vector > columnsData; + for (vector::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++) { - MyColumn *cCol = static_cast(fColMap.find(*it)->second); - columns.push_back(cCol); + fits::Table::Column& cCol = fColMap.find(*it)->second; + columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size])); + fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second); } + //add the time column to the given columns if (fColMap.find("Time") == fColMap.end() && fColMap.find("UnixTimeUTC") == fColMap.end()) @@ -947,24 +774,12 @@ return false; } - MyColumn* timeCol; - bool unixTime = false; - int endIndex = 0; - if (fColMap.find("Time") != fColMap.end()) - { - timeCol = static_cast(fColMap.find("Time")->second); + const fits::Table::Column& timeCol = (fColMap.find("Time") != fColMap.end()) ? fColMap.find("Time")->second : fColMap.find("UnixTimeUTC")->second; + bool unixTime = (fColMap.find("Time") == fColMap.end()); + if (unixTime) + ranges.push_back(make_pair(0,2)); + else ranges.push_back(make_pair(0,1)); - endIndex = fTable->rows(); - } - if (fColMap.find("UnixTimeUTC") != fColMap.end()) - { - timeCol = static_cast(fColMap.find("UnixTimeUTC")->second); - ranges.push_back(make_pair(0,2)); - endIndex = fTable->rows()-1; - unixTime = true; - } - columns.push_back(timeCol); - ///// - const int size = offsets[offsets.size()-1]; - unsigned char* fitsBuffer = new unsigned char[size]; + columnsData.push_back(make_pair(timeCol.type, new char[timeCol.num*timeCol.size])); + fFile->SetPtrAddress(unixTime ? "UnixTimeUTC" : "Time", columnsData[columnsData.size()-1].second); // stringstream str; @@ -977,29 +792,52 @@ str.precision(20); for (auto it=xValues.begin(); it!=xValues.end(); it++) - *it = new double[fTable->rows()]; - - yValues = new double[fTable->rows()]; + *it = new double[fFile->GetInt("NAXIS2")]; + + yValues = new double[fFile->GetInt("NAXIS2")]; cout.precision(3); - for (int i=1; i<=endIndex; i++) + int endIndex = 0; + int numRows = fFile->GetInt("NAXIS2"); + for (int i=1;ifitsPointer(), i, 1, size, fitsBuffer, &status); - if (status) + fFile->GetNextRow(); + cout << "\r" << "Constructing graph " << ((float)(endIndex)/(float)(fFile->GetInt("NAXIS2")))*100.0 << "%"; + endIndex++; + auto rangesIt = ranges.begin(); + for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++) { - cerr << "An error occurred while reading fits row #" << i << " error code: " << status << endl; - break; + for (int i=rangesIt->first; isecond; i++) + { + switch (it->first) { + case 'L': + str << reinterpret_cast(it->second)[i] << " "; + break; + case 'B': + str << reinterpret_cast(it->second)[i] << " "; + break; + case 'I': + str << reinterpret_cast(it->second)[i] << " "; + break; + case 'J': + str << reinterpret_cast(it->second)[i] << " "; + break; + case 'K': + str << reinterpret_cast(it->second)[i] << " "; + break; + case 'E': + str << reinterpret_cast(it->second)[i] << " "; + break; + case 'D': + str << reinterpret_cast(it->second)[i] << " "; + break; + default: + ; + } + } } - if (WriteRow(str, columns, offsets, fitsBuffer, ranges)<0) - { - status=1; - cerr << "An Error occured while reading the fits row " << i << endl; - return -1; - } -// yValues[i-1] = i; + for (auto it=xValues.begin(); it!= xValues.end(); it++) { str >> (*it)[i-1]; -// cout << (*it)[i-1] << " "; } if (unixTime) @@ -1007,5 +845,5 @@ long u1, u2; str >> u1 >> u2; - // cout << u1 << " " << u2; + boost::posix_time::ptime unixTimeT( boost::gregorian::date(1970, boost::gregorian::Jan, 1), boost::posix_time::seconds(u1) + boost::posix_time::microsec(u2)); @@ -1013,5 +851,4 @@ Time mjdTime(unixTimeT); yValues[i-1] = mjdTime.Mjd(); - // cout << " " << mjdTime.Mjd() << endl; } else @@ -1020,4 +857,11 @@ if (yValues[i-1] < 40587) yValues[i-1] += 40587; + + Time t(yValues[i-1]); + string time = t.GetAsStr("%H:%M:%S%F"); + while (time[time.size()-1] == '0' && time.size() > 2) + { + time = time.substr(0, time.size()-1); + } } if (i==1) @@ -1027,11 +871,9 @@ plot->setTitle(title.c_str()); } -// cout << yValues[i-1] << " "; -// cout << endl; } //set the actual data. auto jt = xValues.begin(); for (auto it=curves.begin(); it != curves.end(); it++, jt++) - (*it)->setRawData(yValues, *jt, endIndex); + (*it)->setRawData(yValues, *jt, endIndex-1); QStack stack; @@ -1061,5 +903,7 @@ zoom.setZoomStack(stack); - delete[] fitsBuffer; +// delete[] fitsBuffer; + for (auto it = columnsData.begin(); it != columnsData.end(); it++) + delete[] it->second; window.resize(600, 400); window.show(); @@ -1101,6 +945,8 @@ ("list,l", po_switch(), "List all tables and columns in file") ("header,h", po_switch(), "Dump header of given table") + ("stat,s", po_switch(), "Perform statistics instead of dump") #ifdef PLOTTING_PLEASE ("graph,g", po_switch(), "Plot the columns instead of dumping them") + ("dots,d", po_switch(), "Plot using dots instead of lines") #endif ;