/* * fitsCompressor.cc * * Created on: May 7, 2013 * Author: lyard */ #include "Configuration.h" #include "../externals/factfits.h" #include "../externals/ofits.h" #include "../externals/checksum.h" #include #include #include #include using namespace std; typedef struct TileHeader { char id[4]; uint32_t numRows; uint64_t size; TileHeader(uint32_t nRows=0, uint64_t s=0) : id({'T', 'I', 'L', 'E'}), numRows(nRows), size(s) { }; } __attribute__((__packed__)) TileHeader; typedef struct BlockHeader { uint64_t size; char ordering; unsigned char numProcs; BlockHeader(uint64_t s=0, char o=FACT_ROW_MAJOR, unsigned char n=1) : size(s), ordering(o), numProcs(n) {} } __attribute__((__packed__)) BlockHeader; class CompressedFitsFile { public: class HeaderEntry { public: /** * Default constructor */ HeaderEntry(): _key(""), _value(""), _comment(""), _fitsString("") { } /** * Regular constructor. * @param key the name of the keyword entry * @param value its value * @param comment an optionnal comment to be placed after the value */ template HeaderEntry(const string& k, const T& val, const string& comm) : _key(k), _value(""), _comment(comm), _fitsString("") { setValue(val); } /** * From fits.h */ string Trim(const string &str, char c=' ') { // Trim Both leading and trailing spaces const size_t pstart = str.find_first_not_of(c); // Find the first character position after excluding leading blank spaces const size_t pend = str.find_last_not_of(c); // Find the first character position from reverse af // if all spaces or empty return an empty string if (string::npos==pstart || string::npos==pend) return string(); return str.substr(pstart, pend-pstart+1); } /** * Constructor from the original fits entry */ HeaderEntry(const string& line) { _fitsString = line.data(); //parse the line _key = Trim(line.substr(0,8)); //COMMENT and/or HISTORY values if (line.substr(8,2)!= "= ") { _value = ""; _comment = Trim(line.substr(10)); return; } string next=line.substr(10); const size_t slash = next.find_first_of('/'); _value = Trim(Trim(Trim(next.substr(0, slash)), '\'')); _comment = Trim(next.substr(slash+1)); } /** * Alternative constroctor from the fits entry */ HeaderEntry(const vector& lineVec) { HeaderEntry(string(lineVec.data())); } /** * Destructor */ virtual ~HeaderEntry(){} const string& value() const {return _value;} const string& key() const {return _key;} const string& comment() const {return _comment;} const string& fitsString() const {return _fitsString;} /** * Set a keyword value. * @param value The value to be set * @param update whether the value already exist or not. To be modified soon. */ template void setValue(const T& val) { ostringstream str; str << val; _value = str.str(); buildFitsString(); }; /** * Set the comment for a given entry */ void setComment(const string& comm) { _comment = comm; buildFitsString(); } private: /** * Construct the FITS header string from the key, value and comment */ void buildFitsString() { ostringstream str; unsigned int totSize = 0; // Tuncate the key if required if (_key.length() > 8) { str << _key.substr(0, 8); totSize += 8; } else { str << _key; totSize += _key.length(); } // Append space if key is less than 8 chars long for (int i=totSize; i<8;i++) { str << " "; totSize++; } // Add separator str << "= "; totSize += 2; // Format value if (_value.length() < 20) for (;totSize<30-_value.length();totSize++) str << " "; if (_value.length() > 70) { str << _value.substr(0,70); totSize += 70; } else { str << _value; totSize += _value.size(); } // If there is space remaining, add comment area if (totSize < 77) { str << " / "; totSize += 3; if (totSize < 80) { unsigned int commentSize = 80 - totSize; if (_comment.length() > commentSize) { str << _comment.substr(0,commentSize); totSize += commentSize; } else { str << _comment; totSize += _comment.length(); } } } // If there is yet some free space, fill up the entry with spaces for (int i=totSize; i<80;i++) str << " "; _fitsString = str.str(); // Check for correct completion if (_fitsString.length() != 80) cout << "Error |" << _fitsString << "| is not of length 80" << endl; } string _key; ///< the key (name) of the header entry string _value; ///< the value of the header entry string _comment; ///< the comment associated to the header entry string _fitsString; ///< the string that will be written to the fits file }; /** * Supported compressions */ typedef enum { UNCOMPRESSED, SMOOTHMAN } FitsCompression; /** * Columns class */ class ColumnEntry { public: /** * Default constructor */ ColumnEntry(); /** * Default destructor */ virtual ~ColumnEntry(){} /** * Constructor from values * @param n the column name * @param t the column type * @param numof the number of entries in the column * @param comp the compression type for this column */ ColumnEntry(const string& n, char t, int numOf, BlockHeader& head, vector& seq) : _name(n), _num(numOf), _typeSize(0), _offset(0), _type(t), _description(""), _header(head), _compSequence(seq) { switch (t) { case 'L': case 'A': case 'B': _typeSize = 1; break; case 'I': _typeSize = 2; break; case 'J': case 'E': _typeSize = 4; break; case 'K': case 'D': _typeSize = 8; break; default: cout << "Error: typename " << t << " missing in the current implementation" << endl; }; ostringstream str; str << "data format of field: "; switch (t) { case 'L': str << "1-byte BOOL"; break; case 'A': str << "1-byte CHAR"; break; case 'B': str << "BYTE"; break; case 'I': str << "2-byte INTEGER"; break; case 'J': str << "4-byte INTEGER"; break; case 'K': str << "8-byte INTEGER"; break; case 'E': str << "4-byte FLOAT"; break; case 'D': str << "8-byte FLOAT"; break; } _description = str.str(); } const string& name() const { return _name;}; int width() const { return _num*_typeSize;}; int offset() const { return _offset; }; int numElems() const { return _num; }; int sizeOfElems() const { return _typeSize;}; void setOffset(int off) { _offset = off;}; char type() const { return _type;}; string getDescription() const { return _description;} BlockHeader& getBlockHeader() { return _header;} const vector& getCompressionSequence() const { return _compSequence;} const char& getColumnOrdering() const { return _header.ordering;} string getCompressionString() const { return "FACT"; /* ostringstream str; for (uint32_t i=0;i<_compSequence.size();i++) switch (_compSequence[i]) { case FACT_RAW: if (str.str().size() == 0) str << "RAW"; break; case FACT_SMOOTHING: str << "SMOOTHING "; break; case FACT_HUFFMAN16: str << "HUFFMAN16 "; break; }; return str.str();*/ } private: string _name; ///< name of the column int _num; ///< number of elements contained in one row of this column int _typeSize; ///< the number of bytes taken by one element int _offset; ///< the offset of the column, in bytes, from the beginning of one row char _type; ///< the type of the column, as specified by the fits documentation string _description; ///< a description for the column. It will be placed in the header BlockHeader _header; vector _compSequence; }; public: ///@brief default constructor. Assigns a default number of rows and tiles CompressedFitsFile(uint32_t numTiles=100, uint32_t numRowsPerTile=100); ///@brief default destructor virtual ~CompressedFitsFile(); ///@brief get the header of the file vector& getHeaderEntries() { return _header;} protected: ///@brief protected function to allocate the intermediate buffers bool reallocateBuffers(); //FITS related stuff vector _header; ///< Header keys vector _columns; ///< Columns in the file uint32_t _numTiles; ///< Number of tiles (i.e. groups of rows) uint32_t _numRowsPerTile; ///< Number of rows per tile uint32_t _totalNumRows; ///< Total number of raws uint32_t _rowWidth; ///< Total number of bytes in one row bool _headerFlushed; ///< Flag telling whether the header record is synchronized with the data on disk char* _buffer; ///< Memory buffer to store rows while they are not compressed Checksum _checksum; ///< Checksum for asserting the consistency of the data fstream _file; ///< The actual file streamer for accessing disk data //compression related stuff typedef pair CatalogEntry; typedef vector CatalogRow; typedef vector CatalogType; CatalogType _catalog; ///< Catalog, i.e. the main table that points to the compressed data. uint64_t _heapPtr; ///< the address in the file of the heap area vector _transposedBuffer; ///< Memory buffer to store rows while they are transposed vector _compressedBuffer; ///< Memory buffer to store rows while they are compressed //thread related stuff uint32_t _numThreads; ///< The number of threads that will be used to compress uint32_t _threadIndex; ///< A variable to assign threads indices vector _thread; ///< The thread handler of the compressor vector _threadNumRows; ///< Total number of rows for thread to compress vector _threadStatus; ///< Flag telling whether the buffer to be transposed (and compressed) is full or empty //thread states. Not all used, but they do not hurt static const uint32_t _THREAD_WAIT_; ///< Thread doing nothing static const uint32_t _THREAD_COMPRESS_; ///< Thread working, compressing static const uint32_t _THREAD_DECOMPRESS_; ///< Thread working, decompressing static const uint32_t _THREAD_WRITE_; ///< Thread writing data to disk static const uint32_t _THREAD_READ_; ///< Thread reading data from disk static const uint32_t _THREAD_EXIT_; ///< Thread exiting static HeaderEntry _dummyHeaderEntry; ///< Dummy entry for returning if requested on is not found }; class CompressedFitsWriter : public CompressedFitsFile { public: ///@brief Default constructor. 100 tiles of 100 rows each are assigned by default CompressedFitsWriter(uint32_t numTiles=100, uint32_t numRowsPerTile=100); ///@brief default destructor virtual ~CompressedFitsWriter(); ///@brief add one column to the file bool addColumn(const ColumnEntry& column); ///@brief sets a given header key bool setHeaderKey(const HeaderEntry&); ///@brief open a new fits file bool open(const string& fileName, const string& tableName="Data"); ///@brief close the opened file bool close(); ///@brief write one row of data, already placed in bufferToWrite. Does the byte-swapping bool writeBinaryRow(const char* bufferToWrite); ///@brief assign a given (already loaded) drs calibration void setDrsCalib(int16_t* data); ///@brief set the number of worker threads compressing the data. bool setNumWorkingThreads(uint32_t num); private: ///@brief compresses one buffer of data, the one given by threadIndex uint64_t compressBuffer(uint32_t threadIndex); ///@brief writes an already compressed buffer to disk bool writeCompressedDataToDisk(uint32_t threadID, uint32_t sizeToWrite); ///@brief add the header checksum to the datasum void addHeaderChecksum(Checksum& checksum); ///@brief write the header. If closingFile is set to true, checksum is calculated void writeHeader(bool closingFile = false); ///@brief write the compressed data catalog. If closingFile is set to true, checksum is calculated void writeCatalog(bool closingFile=false); /// FIXME this was a bad idea. Move everything to the regular header vector _defaultHeader; /// the main function compressing the data static void* threadFunction(void* context); /// Write the drs calibration to disk, if any void writeDrsCalib(); /// Copy and transpose (or semi-transpose) one tile of data void copyTransposeTile(uint32_t index); /// Specific compression functions uint32_t compressUNCOMPRESSED(char* dest, const char* src, uint32_t numRows, uint32_t sizeOfElems, uint32_t numRowElems); uint32_t compressHUFFMAN(char* dest, const char* src, uint32_t numRows, uint32_t sizeOfElems, uint32_t numRowElems); uint32_t compressSMOOTHMAN(char* dest, char* src, uint32_t numRows, uint32_t sizeOfElems, uint32_t numRowElems); uint32_t applySMOOTHING(char* dest, char* src, uint32_t numRows, uint32_t sizeOfElems, uint32_t numRowElems); int32_t _checkOffset; ///< offset to the data pointer to calculate the checksum int16_t* _drsCalibData; ///< array of the Drs baseline mean int32_t _threadLooper; ///< Which thread will deal with the upcoming bunch of data ? pthread_mutex_t _mutex; ///< mutex for compressing threads static string _emptyBlock; ///< an empty block to be apened at the end of a file so that its length is a multiple of 2880 static string _fitsHeader; ///< the default header to be written in every fits file }; const uint32_t CompressedFitsFile::_THREAD_WAIT_ = 0; const uint32_t CompressedFitsFile::_THREAD_COMPRESS_ = 1; const uint32_t CompressedFitsFile::_THREAD_DECOMPRESS_= 2; const uint32_t CompressedFitsFile::_THREAD_WRITE_ = 3; const uint32_t CompressedFitsFile::_THREAD_READ_ = 4; const uint32_t CompressedFitsFile::_THREAD_EXIT_ = 5; template<> void CompressedFitsFile::HeaderEntry::setValue(const string& v) { string val = v; if (val.size() > 2 && val[0] == '\'') { size_t pos = val.find_last_of("'"); if (pos != string::npos && pos != 0) val = val.substr(1, pos-1); } ostringstream str; str << "'" << val << "'"; for (int i=str.str().length(); i<20;i++) str << " "; _value = str.str(); buildFitsString(); } /** * Default header to be written in all fits files */ string CompressedFitsWriter::_fitsHeader = "SIMPLE = T / file does conform to FITS standard " "BITPIX = 8 / number of bits per data pixel " "NAXIS = 0 / number of data axes " "EXTEND = T / FITS dataset may contain extensions " "CHECKSUM= '4AcB48bA4AbA45bA' / Checksum for the whole HDU " "DATASUM = ' 0' / Checksum for the data block " "COMMENT FITS (Flexible Image Transport System) format is defined in 'Astronomy" "COMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H " "END " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " "; /** * Empty block to be appenned at the end of files, so that the length matches multiple of 2880 bytes * */ string CompressedFitsWriter::_emptyBlock = " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " "; CompressedFitsFile::HeaderEntry CompressedFitsFile::_dummyHeaderEntry; /**************************************************************** * SUPER CLASS DEFAULT CONSTRUCTOR ****************************************************************/ CompressedFitsFile::CompressedFitsFile(uint32_t numTiles, uint32_t numRowsPerTile) : _header(), _columns(), _numTiles(numTiles), _numRowsPerTile(numRowsPerTile), _totalNumRows(0), _rowWidth(0), _headerFlushed(false), _buffer(NULL), _checksum(0), _heapPtr(0), _transposedBuffer(1), _compressedBuffer(1), _numThreads(1), _threadIndex(0), _thread(1), _threadNumRows(1), _threadStatus(1) { _catalog.resize(_numTiles); _transposedBuffer[0] = NULL; _compressedBuffer[0] = NULL; _threadStatus[0] = _THREAD_WAIT_; _threadNumRows[0] = 0; } /**************************************************************** * SUPER CLASS DEFAULT DESTRUCTOR ****************************************************************/ CompressedFitsFile::~CompressedFitsFile() { if (_buffer != NULL) { delete[] _buffer; _buffer = NULL; for (uint32_t i=0;i<_numThreads;i++) { _compressedBuffer[i] = _compressedBuffer[i]-4; delete[] _transposedBuffer[i]; delete[] _compressedBuffer[i]; _transposedBuffer[i] = NULL; _compressedBuffer[i] = NULL; } } if (_file.is_open()) _file.close(); } /**************************************************************** * REALLOCATE BUFFER ****************************************************************/ bool CompressedFitsFile::reallocateBuffers() { if (_buffer) { delete[] _buffer; for (uint32_t i=0;i<_compressedBuffer.size();i++) { _compressedBuffer[i] = _compressedBuffer[i]-4; delete[] _transposedBuffer[i]; delete[] _compressedBuffer[i]; } } _buffer = new char[_rowWidth*_numRowsPerTile]; if (_buffer == NULL) return false; if (_compressedBuffer.size() != _numThreads) { _transposedBuffer.resize(_numThreads); _compressedBuffer.resize(_numThreads); } for (uint32_t i=0;i<_numThreads;i++) { _transposedBuffer[i] = new char[_rowWidth*_numRowsPerTile]; _compressedBuffer[i] = new char[_rowWidth*_numRowsPerTile + _columns.size() + sizeof(TileHeader)]; //use a bit more memory for compression flags if (_transposedBuffer[i] == NULL || _compressedBuffer[i] == NULL) return false; //shift the compressed buffer by 4 bytes, for checksum calculation memset(_compressedBuffer[i], 0, 4); _compressedBuffer[i] = _compressedBuffer[i]+4; //initialize the tile header TileHeader tileHeader; memcpy(_compressedBuffer[i], &tileHeader, sizeof(TileHeader)); } return true; } /**************************************************************** * DEFAULT WRITER CONSTRUCTOR ****************************************************************/ CompressedFitsWriter::CompressedFitsWriter(uint32_t numTiles, uint32_t numRowsPerTile) : CompressedFitsFile(numTiles, numRowsPerTile), _checkOffset(0), _drsCalibData(NULL), _threadLooper(0) { _defaultHeader.push_back(HeaderEntry("XTENSION", "'BINTABLE' ", "binary table extension")); _defaultHeader.push_back(HeaderEntry("BITPIX", 8, "8-bit bytes")); _defaultHeader.push_back(HeaderEntry("NAXIS", 2, "2-dimensional binary table")); _defaultHeader.push_back(HeaderEntry("NAXIS1", _rowWidth, "width of table in bytes")); _defaultHeader.push_back(HeaderEntry("NAXIS2", numTiles, "num of rows in table")); _defaultHeader.push_back(HeaderEntry("PCOUNT", 0, "size of special data area")); _defaultHeader.push_back(HeaderEntry("GCOUNT", 1, "one data group (required keyword)")); _defaultHeader.push_back(HeaderEntry("TFIELDS", _columns.size(), "number of fields in each row")); _defaultHeader.push_back(HeaderEntry("CHECKSUM", "'0000000000000000' ", "Checksum for the whole HDU")); _defaultHeader.push_back(HeaderEntry("DATASUM", " 0", "Checksum for the data block")); //compression stuff _defaultHeader.push_back(HeaderEntry("ZTABLE", "T", "Table is compressed")); _defaultHeader.push_back(HeaderEntry("ZNAXIS1", 0, "Width of uncompressed rows")); _defaultHeader.push_back(HeaderEntry("ZNAXIS2", 0, "Number of uncompressed rows")); _defaultHeader.push_back(HeaderEntry("ZPCOUNT", 0, "")); _defaultHeader.push_back(HeaderEntry("ZHEAPPTR", 0, "")); _defaultHeader.push_back(HeaderEntry("ZTILELEN", numRowsPerTile, "Number of rows per tile")); _defaultHeader.push_back(HeaderEntry("THEAP", 0, "")); pthread_mutex_init(&_mutex, NULL); } /**************************************************************** * DEFAULT DESTRUCTOR ****************************************************************/ CompressedFitsWriter::~CompressedFitsWriter() { pthread_mutex_destroy(&_mutex); } /**************************************************************** * SET THE POINTER TO THE DRS CALIBRATION ****************************************************************/ void CompressedFitsWriter::setDrsCalib(int16_t* data) { _drsCalibData = data; } /**************************************************************** * SET NUM WORKING THREADS ****************************************************************/ bool CompressedFitsWriter::setNumWorkingThreads(uint32_t num) { if (_file.is_open()) return false; if (num < 1 || num > 64) { cout << "ERROR: num threads must be between 1 and 64. Ignoring" << endl; return false; } _numThreads = num; _transposedBuffer[0] = NULL; _compressedBuffer[0] = NULL; _threadStatus.resize(num); _thread.resize(num); _threadNumRows.resize(num); for (uint32_t i=0;i header; header.push_back(HeaderEntry("XTENSION", "'BINTABLE' ", "binary table extension")); header.push_back(HeaderEntry("BITPIX" , 8 , "8-bit bytes")); header.push_back(HeaderEntry("NAXIS" , 2 , "2-dimensional binary table")); header.push_back(HeaderEntry("NAXIS1" , 1024*1440*2 , "width of table in bytes")); header.push_back(HeaderEntry("NAXIS2" , 1 , "number of rows in table")); header.push_back(HeaderEntry("PCOUNT" , 0 , "size of special data area")); header.push_back(HeaderEntry("GCOUNT" , 1 , "one data group (required keyword)")); header.push_back(HeaderEntry("TFIELDS" , 1 , "number of fields in each row")); header.push_back(HeaderEntry("CHECKSUM", "'0000000000000000' ", "Checksum for the whole HDU")); header.push_back(HeaderEntry("DATASUM" , " 0" , "Checksum for the data block")); header.push_back(HeaderEntry("EXTNAME" , "'ZDrsCellOffsets' ", "name of this binary table extension")); header.push_back(HeaderEntry("TTYPE1" , "'OffsetCalibration' ", "label for field 1")); header.push_back(HeaderEntry("TFORM1" , "'1474560I' ", "data format of field: 2-byte INTEGER")); for (uint32_t i=0;i(swappedBytes)[j]; reinterpret_cast(swappedBytes)[j] = reinterpret_cast(swappedBytes)[j+1]; reinterpret_cast(swappedBytes)[j+1] = inter; } _file.write(reinterpret_cast(swappedBytes), 2048); checksum.add(reinterpret_cast(swappedBytes), 2048); } uint64_t whereDidIStop = _file.tellp(); delete[] swappedBytes; //No need to pad the data, as (1440*1024*2)%2880==0 //calculate the checksum from the header ostringstream str; str << checksum.val(); header[9] = HeaderEntry("DATASUM", str.str(), "Checksum for the data block"); for (vector::iterator it=header.begin();it!=header.end(); it++) checksum.add(it->fitsString().c_str(), 80); string end("END "); string space(" "); checksum.add(end.c_str(), 80); int headerRowsLeft = 36 - (header.size() + 1)%36; for (int i=0;i::iterator it=_columns.begin(); it != _columns.end(); it++) { if (it->name() == column.name()) { cout << "Warning: column already exist (" << column.name() << "). Ignoring" << endl; return false; } } _columns.push_back(column); _columns.back().setOffset(_rowWidth); _rowWidth += column.width(); reallocateBuffers(); ostringstream str, str2, str3; str << "TTYPE" << _columns.size(); str2 << column.name(); str3 << "label for field "; if (_columns.size() < 10) str3 << " "; if (_columns.size() < 100) str3 << " "; str3 << _columns.size(); setHeaderKey(HeaderEntry(str.str(), str2.str(), str3.str())); str.str(""); str2.str(""); str3.str(""); str << "TFORM" << _columns.size(); str2 << "1QB"; str3 << "data format of field " << _columns.size(); setHeaderKey(HeaderEntry(str.str(), str2.str(), str3.str())); str.str(""); str2.str(""); str3.str(""); str << "ZFORM" << _columns.size(); str2 << column.numElems() << column.type(); str3 << "Original format of field " << _columns.size(); setHeaderKey(HeaderEntry(str.str(), str2.str(), str3.str())); str.str(""); str2.str(""); str3.str(""); str << "ZCTYP" << _columns.size(); str2 << column.getCompressionString(); str3 << "Comp. Scheme of field " << _columns.size(); setHeaderKey(HeaderEntry(str.str(), str2.str(), str3.str())); //resize the catalog vector accordingly for (uint32_t i=0;i<_numTiles;i++) { _catalog[i].resize(_columns.size()); for (uint32_t j=0;j<_catalog[i].size();j++) _catalog[i][j] = make_pair(0,0); } return true; } /**************************************************************** * SET HEADER KEY ****************************************************************/ bool CompressedFitsWriter::setHeaderKey(const HeaderEntry& entry) { HeaderEntry ent = entry; for (vector::iterator it=_header.begin(); it != _header.end(); it++) { if (it->key() == entry.key()) { if (entry.comment() == "") ent.setComment(it->comment()); (*it) = ent; _headerFlushed = false; return true; } } for (vector::iterator it=_defaultHeader.begin(); it != _defaultHeader.end(); it++) { if (it->key() == entry.key()) { if (entry.comment() == "") ent.setComment(it->comment()); (*it) = ent; _headerFlushed = false; return true; } } if (_totalNumRows != 0) { cout << "Error: new header keys (" << entry.key() << ") must be set before the first row is written. Ignoring." << endl; return false; } _header.push_back(entry); _headerFlushed = false; return true; } /**************************************************************** * OPEN ****************************************************************/ bool CompressedFitsWriter::open(const string& fileName, const string& tableName) { _file.open(fileName.c_str(), ios_base::out); if (!_file.is_open()) { cout << "Error: Could not open the file (" << fileName << ")." << endl; return false; } _defaultHeader.push_back(HeaderEntry("EXTNAME", tableName, "name of this binary table extension")); _headerFlushed = false; _threadIndex = 0; //create requested number of threads for (uint32_t i=0;i<_numThreads;i++) pthread_create(&(_thread[i]), NULL, threadFunction, this); //wait for the threads to start while (_numThreads != _threadIndex) usleep(1000); //set the writing fence to the last thread _threadIndex = _numThreads-1; return (_file.good()); } /**************************************************************** * WRITE HEADER ****************************************************************/ void CompressedFitsWriter::writeHeader(bool closingFile) { if (_headerFlushed) return; if (!_file.is_open()) return; long cPos = _file.tellp(); _file.seekp(0); _file.write(_fitsHeader.c_str(), 2880); //Write the DRS calib table here ! writeDrsCalib(); //we are now at the beginning of the main table. Write its header for (vector::iterator it=_defaultHeader.begin(); it != _defaultHeader.end(); it++) _file.write(it->fitsString().c_str(), 80); for (vector::iterator it=_header.begin(); it != _header.end(); it++) _file.write(it->fitsString().c_str(), 80); _file.write("END ", 80); long here = _file.tellp(); if (here%2880) _file.write(_emptyBlock.c_str(), 2880 - here%2880); _headerFlushed = true; here = _file.tellp(); if (here%2880) cout << "Error: seems that header did not finish at the end of a block." << endl; if (here > cPos && cPos != 0) { cout << "Error, entries were added after the first row was written. This is not supposed to happen." << endl; return; } here = _file.tellp(); writeCatalog(closingFile); here = _file.tellp() - here; _heapPtr = here; if (cPos != 0) _file.seekp(cPos); } /**************************************************************** * WRITE CATALOG * WARNING: writeCatalog is only meant to be used by writeHeader. * external usage will most likely corrupt the file ****************************************************************/ void CompressedFitsWriter::writeCatalog(bool closingFile) { uint32_t sizeWritten = 0; for (uint32_t i=0;i<_catalog.size();i++) { for (uint32_t j=0;j<_catalog[i].size();j++) { //swap the bytes int8_t swappedEntry[16]; swappedEntry[0] = reinterpret_cast(&_catalog[i][j].first)[7]; swappedEntry[1] = reinterpret_cast(&_catalog[i][j].first)[6]; swappedEntry[2] = reinterpret_cast(&_catalog[i][j].first)[5]; swappedEntry[3] = reinterpret_cast(&_catalog[i][j].first)[4]; swappedEntry[4] = reinterpret_cast(&_catalog[i][j].first)[3]; swappedEntry[5] = reinterpret_cast(&_catalog[i][j].first)[2]; swappedEntry[6] = reinterpret_cast(&_catalog[i][j].first)[1]; swappedEntry[7] = reinterpret_cast(&_catalog[i][j].first)[0]; swappedEntry[8] = reinterpret_cast(&_catalog[i][j].second)[7]; swappedEntry[9] = reinterpret_cast(&_catalog[i][j].second)[6]; swappedEntry[10] = reinterpret_cast(&_catalog[i][j].second)[5]; swappedEntry[11] = reinterpret_cast(&_catalog[i][j].second)[4]; swappedEntry[12] = reinterpret_cast(&_catalog[i][j].second)[3]; swappedEntry[13] = reinterpret_cast(&_catalog[i][j].second)[2]; swappedEntry[14] = reinterpret_cast(&_catalog[i][j].second)[1]; swappedEntry[15] = reinterpret_cast(&_catalog[i][j].second)[0]; if (closingFile) { _checksum.add(reinterpret_cast(swappedEntry), 16); } _file.write(reinterpret_cast(&swappedEntry[0]), 2*sizeof(int64_t)); sizeWritten += 2*sizeof(int64_t); } } //we do not reserve space for now because fverify does not like that. //TODO bug should be fixed in the new version. Install it on the cluster and restor space reservation return ; //write the padding so that the HEAP section starts at a 2880 bytes boundary if (sizeWritten % 2880 != 0) { vector nullVec(2880 - sizeWritten%2880, 0); _file.write(nullVec.data(), 2880 - sizeWritten%2880); } } /**************************************************************** * ADD HEADER CHECKSUM ****************************************************************/ void CompressedFitsWriter::addHeaderChecksum(Checksum& checksum) { for (vector::iterator it=_defaultHeader.begin();it!=_defaultHeader.end(); it++) _checksum.add(it->fitsString().c_str(), 80); for (vector::iterator it=_header.begin(); it != _header.end(); it++) _checksum.add(it->fitsString().c_str(), 80); string end("END "); string space(" "); checksum.add(end.c_str(), 80); int headerRowsLeft = 36 - (_defaultHeader.size() + _header.size() + 1)%36; for (int i=0;i 0) { setHeaderKey(HeaderEntry("NAXIS1", 2*sizeof(int64_t)*_catalog[0].size(), "width of table in bytes")); setHeaderKey(HeaderEntry("NAXIS2", _numTiles, "")); setHeaderKey(HeaderEntry("TFIELDS", _columns.size(), "number of fields in each row")); int64_t heapSize = 0; int64_t compressedOffset = 0; for (uint32_t i=0;i<_catalog.size();i++) { compressedOffset += sizeof(TileHeader); heapSize += sizeof(TileHeader); for (uint32_t j=0;j<_catalog[i].size();j++) { heapSize += _catalog[i][j].first; // cout << "heapSize: " << heapSize << endl; //set the catalog offsets to their actual values _catalog[i][j].second = compressedOffset; compressedOffset += _catalog[i][j].first; //special case if entry has zero length if (_catalog[i][j].first == 0) _catalog[i][j].second = 0; } } setHeaderKey(HeaderEntry("PCOUNT", heapSize, "size of special data area")); } else { setHeaderKey(HeaderEntry("NAXIS1", _rowWidth, "width of table in bytes")); setHeaderKey(HeaderEntry("NAXIS2", 0, "")); setHeaderKey(HeaderEntry("TFIELDS", _columns.size(), "number of fields in each row")); setHeaderKey(HeaderEntry("PCOUNT", 0, "size of special data area")); } writeHeader(true); ostringstream str; str << _checksum.val(); setHeaderKey(HeaderEntry("DATASUM", str.str(), "")); addHeaderChecksum(_checksum); setHeaderKey(HeaderEntry("CHECKSUM", _checksum.str(), "")); //update header value writeHeader(); //update file length long here = _file.tellp(); if (here%2880) { vector nullVec(2880 - here%2880, 0); _file.write(nullVec.data(), 2880 - here%2880); } _file.close(); return true; } /**************************************************************** * COPY TRANSPOSE TILE ****************************************************************/ void CompressedFitsWriter::copyTransposeTile(uint32_t index) { uint32_t thisRoundNumRows = (_totalNumRows%_numRowsPerTile) ? _totalNumRows%_numRowsPerTile : _numRowsPerTile; //copy the tile and transpose it uint32_t offset = 0; for (uint32_t i=0;i<_columns.size();i++) { switch (_columns[i].getColumnOrdering())//getCompression()) { case FACT_ROW_MAJOR: for (uint32_t k=0;k(&src[j*sizeOfElems*numRows]), numRows*(sizeOfElems/2)); reinterpret_cast(&dest[huffmanOffset])[0] = huffmanOutput.size() - previousHuffmanSize; huffmanOffset += sizeof(uint32_t); previousHuffmanSize = huffmanOutput.size(); } const size_t totalSize = huffmanOutput.size() + huffmanOffset; //only copy if not larger than not-compressed size if (totalSize < numRows*sizeOfElems*numRowElems) memcpy(&dest[huffmanOffset], huffmanOutput.data(), huffmanOutput.size()); return totalSize; } /**************************************************************** * COMPRESS BUFFER ****************************************************************/ uint32_t CompressedFitsWriter::compressSMOOTHMAN(char* dest, char* src, uint32_t numRows, uint32_t sizeOfElems, uint32_t numRowElems) { uint32_t colWidth = numRowElems; for (int j=colWidth*numRows-1;j>1;j--) reinterpret_cast(src)[j] = reinterpret_cast(src)[j] - (reinterpret_cast(src)[j-1]+reinterpret_cast(src)[j-2])/2; //call the huffman transposed return compressHUFFMAN(dest, src, numRowElems, sizeOfElems, numRows); } uint32_t CompressedFitsWriter::applySMOOTHING(char* , char* src, uint32_t numRows, uint32_t sizeOfElems, uint32_t numRowElems) { uint32_t colWidth = numRowElems; for (int j=colWidth*numRows-1;j>1;j--) reinterpret_cast(src)[j] = reinterpret_cast(src)[j] - (reinterpret_cast(src)[j-1]+reinterpret_cast(src)[j-2])/2; return numRows*sizeOfElems*numRowElems; } /**************************************************************** * COMPRESS BUFFER ****************************************************************/ uint64_t CompressedFitsWriter::compressBuffer(uint32_t threadIndex) { uint32_t thisRoundNumRows = (_threadNumRows[threadIndex]%_numRowsPerTile) ? _threadNumRows[threadIndex]%_numRowsPerTile : _numRowsPerTile; uint32_t offset=0; uint32_t currentCatalogRow = (_threadNumRows[threadIndex]-1)/_numRowsPerTile; uint64_t compressedOffset = sizeof(TileHeader); //skip the 'TILE' marker and tile size entry //now compress each column one by one by calling compression on arrays for (uint32_t i=0;i<_columns.size();i++) { _catalog[currentCatalogRow][i].second = compressedOffset; if (_columns[i].numElems() == 0) continue; BlockHeader& head = _columns[i].getBlockHeader(); const vector& sequence = _columns[i].getCompressionSequence(); //set the default byte telling if uncompressed the compressed Flag uint64_t previousOffset = compressedOffset; //skip header data compressedOffset += sizeof(BlockHeader) + sizeof(uint16_t)*sequence.size(); for (uint32_t j=0;j _columns[i].sizeOfElems()*_columns[i].numElems()*thisRoundNumRows+sizeof(BlockHeader)+sizeof(uint16_t)*sequence.size()) {//if so set flag and redo it uncompressed cout << "REDOING UNCOMPRESSED" << endl; compressedOffset = previousOffset + sizeof(BlockHeader) + 1; compressedOffset += compressUNCOMPRESSED(&(_compressedBuffer[threadIndex][compressedOffset]), &(_transposedBuffer[threadIndex][offset]), thisRoundNumRows, _columns[i].sizeOfElems(), _columns[i].numElems()); BlockHeader he; he.size = compressedOffset - previousOffset; he.numProcs = 1; he.ordering = FACT_ROW_MAJOR; memcpy(&(_compressedBuffer[threadIndex][previousOffset]), (char*)(&he), sizeof(BlockHeader)); _compressedBuffer[threadIndex][previousOffset+sizeof(BlockHeader)] = FACT_RAW; offset += thisRoundNumRows*_columns[i].sizeOfElems()*_columns[i].numElems(); _catalog[currentCatalogRow][i].first = compressedOffset - _catalog[currentCatalogRow][i].second; continue; } head.size = compressedOffset - previousOffset; memcpy(&(_compressedBuffer[threadIndex][previousOffset]), (char*)(&head), sizeof(BlockHeader)); memcpy(&(_compressedBuffer[threadIndex][previousOffset+sizeof(BlockHeader)]), sequence.data(), sizeof(uint16_t)*sequence.size()); offset += thisRoundNumRows*_columns[i].sizeOfElems()*_columns[i].numElems(); _catalog[currentCatalogRow][i].first = compressedOffset - _catalog[currentCatalogRow][i].second; } TileHeader tHead(thisRoundNumRows, compressedOffset); memcpy(_compressedBuffer[threadIndex], &tHead, sizeof(TileHeader)); return compressedOffset; } /**************************************************************** * WRITE COMPRESS DATA TO DISK ****************************************************************/ bool CompressedFitsWriter::writeCompressedDataToDisk(uint32_t threadID, uint32_t sizeToWrite) { char* checkSumPointer = _compressedBuffer[threadID]; int32_t extraBytes = 0; uint32_t sizeToChecksum = sizeToWrite; if (_checkOffset != 0) {//should we extend the array to the left ? sizeToChecksum += _checkOffset; checkSumPointer -= _checkOffset; memset(checkSumPointer, 0, _checkOffset); } if (sizeToChecksum%4 != 0) {//should we extend the array to the right ? extraBytes = 4 - (sizeToChecksum%4); memset(checkSumPointer+sizeToChecksum, 0,extraBytes); sizeToChecksum += extraBytes; } //do the checksum _checksum.add(checkSumPointer, sizeToChecksum); _checkOffset = (4 - extraBytes)%4; //write data to disk _file.write(_compressedBuffer[threadID], sizeToWrite); return _file.good(); } /**************************************************************** * WRITER THREAD LOOP ****************************************************************/ void* CompressedFitsWriter::threadFunction(void* context) { CompressedFitsWriter* myself =static_cast(context); uint32_t myID = 0; pthread_mutex_lock(&(myself->_mutex)); myID = myself->_threadIndex++; pthread_mutex_unlock(&(myself->_mutex)); uint32_t threadToWaitForBeforeWriting = (myID == 0) ? myself->_numThreads-1 : myID-1; while (myself->_threadStatus[myID] != _THREAD_EXIT_) { while (myself->_threadStatus[myID] == _THREAD_WAIT_) usleep(100000); if (myself->_threadStatus[myID] != _THREAD_COMPRESS_) continue; uint32_t numBytes = myself->compressBuffer(myID); myself->_threadStatus[myID] = _THREAD_WRITE_; //wait for the previous data to be written while (myself->_threadIndex != threadToWaitForBeforeWriting) usleep(1000); //do the actual writing to disk pthread_mutex_lock(&(myself->_mutex)); myself->writeCompressedDataToDisk(myID, numBytes); myself->_threadIndex = myID; pthread_mutex_unlock(&(myself->_mutex)); myself->_threadStatus[myID] = _THREAD_WAIT_; } return NULL; } /**************************************************************** * PRINT USAGE ****************************************************************/ void printUsage() { cout << endl; cout << "The FACT-Fits compressor reads an input Fits file from FACT" " and compresses it.\n It can use a drs calibration in order to" " improve the compression ratio. If so, the input calibration" " is embedded into the compressed file.\n" " By default, the Data column will be compressed using SMOOTHMAN (Thomas' algorithm)" " while other columns will be compressed with the AMPLITUDE coding (Veritas)" "Usage: Compressed_Fits_Test "; cout << endl; } /**************************************************************** * PRINT HELP ****************************************************************/ void printHelp() { cout << endl; cout << "The inputFile is required. It must have fits in its filename and the compressed file will be written in the same folder. " "The fz extension will be added, replacing the .gz one if required \n" "If output is specified, then it will replace the automatically generated output filename\n" "If --drs, followed by a drs calib then it will be applied to the data before compressing\n" "rowPerTile can be used to tune how many rows are in each tile. Default is 100\n" "threads gives the number of threads to use. Cannot be less than the default (1)\n" "compression explicitely gives the compression scheme to use for a given column. The syntax is:\n" "= with one of the following:\n" "UNCOMPRESSED\n" "AMPLITUDE\n" "HUFFMAN\n" "SMOOTHMAN\n" "INT_WAVELET\n" "\n" "--quiet removes any textual output, except error messages\n" "--verify makes the compressor check the compressed data. It will read it back, and compare the reconstructed CHECKSUM and DATASUM with the original file values." ; cout << endl << endl; } /**************************************************************** * SETUP CONFIGURATION ****************************************************************/ void setupConfiguration(Configuration& conf) { po::options_description configs("FitsCompressor options"); configs.add_options() ("inputFile,i", vars(), "Input file") ("drs,d", var(), "Input drs calibration file") ("rowPerTile,r", var(), "Number of rows per tile. Default is 100") ("output,o", var(), "Output file. If empty, .fz is appened to the original name") ("threads,t", var(), "Number of threads to use for compression") ("compression,c", vars(), "which compression to use for which column. Syntax =") ("quiet,q", po_switch(), "Should the program display any text at all ?") ("verify,v", po_switch(), "Should we verify the data that has been compressed ?") ; po::positional_options_description positional; positional.add("inputFile", -1); conf.AddOptions(configs); conf.SetArgumentPositions(positional); } /**************************************************************** * MAIN ****************************************************************/ int main(int argc, const char** argv) { Configuration conf(argv[0]); conf.SetPrintUsage(printUsage); setupConfiguration(conf); if (!conf.DoParse(argc, argv, printHelp)) return -1; //initialize the file names to nothing. string fileNameIn = ""; string fileNameOut = ""; string drsFileName = ""; uint32_t numRowsPerTile = 100; bool displayText=true; //parse configuration if (conf.Get("quiet")) displayText = false; const vector inputFileNameVec = conf.Vec("inputFile"); if (inputFileNameVec.size() != 1) { cout << "Error: "; if (inputFileNameVec.size() == 0) cout << "no"; else cout << inputFileNameVec.size(); cout << " input file(s) given. Expected one. Aborting. Input:" << endl;; for (unsigned int i=0;i("drs"); //Should we verify the data ? bool verifyDataPlease = false; if (conf.Has("verify")) verifyDataPlease = conf.Get("verify"); //should we use a specific output filename ? if (conf.Has("output")) fileNameOut = conf.Get("output"); else { size_t pos = fileNameIn.find(".fits"); if (pos == string::npos) { cout << "ERROR: input file does not seems ot be fits. Aborting." << endl; return -1; } fileNameOut = fileNameIn.substr(0, pos) + ".fz"; } //should we use specific compression on some columns ? const vector columnsCompression = conf.Vec("compression"); //split up values between column names and compression scheme vector> compressions; for (unsigned int i=0;i("rowPerTile"); /************************************************************************************ * Done reading configuration. Open relevant files ************************************************************************************/ //Open input's fits file factfits inFile(fileNameIn); if (inFile.IsCompressedFITS()) { cout << "ERROR: input file is already a compressed fits. Cannot be compressed again: Aborting." << endl; return -1; } //decide how many tiles should be put in the compressed file uint32_t originalNumRows = inFile.GetNumRows(); uint32_t numTiles = (originalNumRows%numRowsPerTile) ? (originalNumRows/numRowsPerTile)+1 : originalNumRows/numRowsPerTile; CompressedFitsWriter outFile(numTiles, numRowsPerTile); //should we use a specific number of threads for compressing ? unsigned int numThreads = 1; if (conf.Has("threads")) { numThreads = conf.Get("threads"); outFile.setNumWorkingThreads(numThreads); } if (!outFile.open(fileNameOut)) { cout << "Error: could not open " << fileNameOut << " for writing" << endl; return -1; } //Because the file to open MUST be given by the constructor, I must use a pointer instead factfits* drsFile = NULL; //try to open the Drs file. If any. if (drsFileName != "") { try { drsFile = new factfits(drsFileName); } catch (...) { cout << "Error: could not open " << drsFileName << " for calibration" << endl; return -1; } } if (displayText) { cout << endl; cout << "**********************" << endl; cout << "Will compress from : " << fileNameIn << endl; cout << "to : " << fileNameOut << endl; if (drsFileName != "") cout << "while calibrating with: " << drsFileName << endl; cout << "Compression will use : " << numThreads << " worker threads" << endl; cout << "Data will be verified : "; if (verifyDataPlease) cout << "yes" << endl; else cout << "no (WARNING !)" << endl; cout << "**********************" << endl; cout << endl; } /************************************************************************************ * Done opening input files. Allocate memory and configure output file ************************************************************************************/ //allocate the buffer for temporary storage of each read/written row uint32_t rowWidth = inFile.GetUInt("NAXIS1"); char* buffer = new char[rowWidth]; //get the source columns const fits::Table::Columns& columns = inFile.GetColumns(); const fits::Table::SortedColumns& sortedColumns = inFile.GetSortedColumns(); if (displayText) cout << "Input file has " << columns.size() << " columns and " << inFile.GetNumRows() << " rows" << endl; //Add columns. for (uint32_t i=0;i "; } //get header structures BlockHeader rawHeader; BlockHeader smoothmanHeader(0, FACT_ROW_MAJOR, 2); vector rawProcessings(1); rawProcessings[0] = FACT_RAW; vector smoothmanProcessings(2); smoothmanProcessings[0] = FACT_SMOOTHING; smoothmanProcessings[1] = FACT_HUFFMAN16; // smoothmanProcessings[2] = FACT_RAW; //first lets see if we do have an explicit request bool explicitRequest = false; for (unsigned int j=0;jfirst;//header[i].key(); if (k == "XTENSION" || k == "BITPIX" || k == "PCOUNT" || k == "GCOUNT" || k == "TFIELDS") continue; if (k == "CHECKSUM") { outFile.setHeaderKey(CompressedFitsFile::HeaderEntry("ZCHKSUM", i->second.value, i->second.comment)); continue; } if (k == "DATASUM") { outFile.setHeaderKey(CompressedFitsFile::HeaderEntry("ZDTASUM", i->second.value, i->second.comment)); continue; } k = k.substr(0,5); if (k == "TTYPE" || k == "TFORM") { string tmpKey = i->second.fitsString; tmpKey[0] = 'Z'; outFile.setHeaderKey(tmpKey); continue; } if (k == "NAXIS") continue; outFile.setHeaderKey(i->second.fitsString); // cout << i->first << endl; } //deal with the DRS calib int16_t* drsCalib16 = NULL; //load the drs calib. data int32_t startCellOffset = -1; if (drsFileName != "") { drsCalib16 = new int16_t[1440*1024]; float* drsCalibFloat = NULL; try { drsCalibFloat = reinterpret_cast(drsFile->SetPtrAddress("BaselineMean")); } catch (...) { cout << "Could not find column BaselineMean in drs calibration file " << drsFileName << ". Aborting" << endl; return -1; } //read the calibration and calculate its integer value drsFile->GetNextRow(); for (uint32_t i=0;i<1440*1024;i++) drsCalib16[i] = (int16_t)(drsCalibFloat[i]*4096.f/2000.f); //get the start cells offsets for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end(); it++) if (it->first == "StartCellData") { startCellOffset = it->second.offset; if (it->second.type != 'I') { cout << "Wrong type for the StartCellData Column: " << it->second.type << " instead of I expected"<< endl; return -1; } } if (startCellOffset == -1) { cout << "WARNING: Could not find StartCellData in input file " << fileNameIn << ". Doing it uncalibrated"<< endl; } else { //assign it to the ouput file outFile.setDrsCalib(drsCalib16); } } /************************************************************************************ * Done configuring compression. Do the real job now ! ************************************************************************************/ if (displayText) cout << "Converting file..." << endl; int numSlices = -1; int32_t dataOffset = -1; //Get the pointer to the column that must be drs-calibrated for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end(); it++) if (it->first == "Data") { numSlices = it->second.num; dataOffset = it->second.offset; } if (numSlices % 1440 != 0) { cout << "seems like the number of samples is not a multiple of 1440. Aborting." << endl; return -1; } if (dataOffset == -1) { cout << "Could not find the column Data in the input file. Aborting." << endl; return -1; } numSlices /= 1440; //set pointers to the readout data to later be able to gather it to "buffer". vector readPointers; vector readOffsets; vector readElemSize; vector readNumElems; for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end(); it++) { readPointers.push_back(inFile.SetPtrAddress(it->first)); readOffsets.push_back(it->second.offset); readElemSize.push_back(it->second.size); readNumElems.push_back(it->second.num); } //Convert each row one after the other for (uint32_t i=0;i(&buffer[startCellOffset])[j]; if (thisStartCell < 0) continue; for (int k=0;k(&buffer[dataOffset])[numSlices*j + k] -= drsCalib16[1024*j + (thisStartCell+k)%1024]; } } outFile.writeBinaryRow(buffer); }; //Get table name for later use in case the compressed file is to be verified string tableName = inFile.GetStr("EXTNAME"); if (displayText) cout << endl << "Done. Flushing output file..." << endl; inFile.close(); if (!outFile.close()) { cout << "Something went wrong while writing the catalog: negative index" << endl; return false; } if (drsCalib16 != NULL) delete[] drsCalib16; if (displayText) cout << "Done." << endl; /************************************************************************************ * Actual job done. Should we verify what we did ? ************************************************************************************/ if (verifyDataPlease) { if (displayText) cout << "Now verify data..." << endl; } else return 0; //get a compressed reader //TEMP try to copy the file too // string copyName("/gpfs/scratch/fact/etienne/copyFile.fz"); // string copyName(fileNameOut+".copy"); string copyName(""); factfits verifyFile(fileNameOut, copyName, tableName, false); //and the header of the compressed file const fits::Table::Keys& header2 = verifyFile.GetKeys(); //get a non-compressed writer ofits reconstructedFile; //figure out its name: /dev/null unless otherwise specified string reconstructedName = fileNameOut+".recons"; // reconstructedName = "/dev/null"; reconstructedFile.open(reconstructedName.c_str(), false); //reconstruct the original columns from the compressed file. string origChecksumStr; string origDatasum; //reset tablename value so that it is re-read from compressed table's header tableName = ""; /************************************************************************************ * Reconstruction setup done. Rebuild original header ************************************************************************************/ //re-tranlate the keys for (fits::Table::Keys::const_iterator it=header2.begin(); it!= header2.end(); it++) { string k = it->first; if (k == "XTENSION" || k == "BITPIX" || k == "PCOUNT" || k == "GCOUNT" || k == "TFIELDS" || k == "ZTABLE" || k == "ZNAXIS1" || k == "ZNAXIS2" || k == "ZHEAPPTR" || k == "ZPCOUNT" || k == "ZTILELEN" || k == "THEAP" || k == "CHECKSUM" || k == "DATASUM" || k == "FCTCPVER") { continue; } if (k == "ZCHKSUM") { reconstructedFile.SetKeyComment("CHECKSUM", it->second.comment); origChecksumStr = it->second.value; continue; } if (k == "ZDTASUM") { reconstructedFile.SetKeyComment("DATASUM", it->second.comment); origDatasum = it->second.value; continue; } if (k == "EXTNAME") { tableName = it->second.value; } k = k.substr(0,5); if (k == "TTYPE") {//we have an original column name here. //manually deal with these in order to preserve the ordering (easier than re-constructing yet another list on the fly) continue; } if (k == "TFORM" || k == "NAXIS" || k == "ZCTYP" ) { continue; } if (k == "ZFORM" || k == "ZTYPE") { string tmpKey = it->second.fitsString; tmpKey[0] = 'T'; reconstructedFile.SetKeyFromFitsString(tmpKey); continue; } reconstructedFile.SetKeyFromFitsString(it->second.fitsString); } if (tableName == "") { cout << "Error: table name from file " << fileNameOut << " could not be found. Aborting" << endl; return -1; } //Restore the original columns for (uint32_t numCol=1; numCol<10000; numCol++) { ostringstream str; str << numCol; if (!verifyFile.HasKey("TTYPE"+str.str())) break; string ttype = verifyFile.GetStr("TTYPE"+str.str()); string tform = verifyFile.GetStr("ZFORM"+str.str()); char type = tform[tform.size()-1]; string number = tform.substr(0, tform.size()-1); int numElems = atoi(number.c_str()); if (number == "") numElems=1; reconstructedFile.AddColumn(numElems, type, ttype, "", "", false); } reconstructedFile.WriteTableHeader(tableName.c_str()); /************************************************************************************ * Original header restored. Do the data ************************************************************************************/ //set pointers to the readout data to later be able to gather it to "buffer". readPointers.clear(); readOffsets.clear(); readElemSize.clear(); readNumElems.clear(); for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end(); it++) { readPointers.push_back(verifyFile.SetPtrAddress(it->first)); readOffsets.push_back(it->second.offset); readElemSize.push_back(it->second.size); readNumElems.push_back(it->second.num); } //do the actual reconstruction work uint32_t i=1; while (i<=verifyFile.GetNumRows() && verifyFile.GetNextRow()) { int count=0; for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end();it++) { memcpy(&buffer[readOffsets[count]], readPointers[count], readElemSize[count]*readNumElems[count]); count++; } if (displayText) cout << "\r Row " << i << flush; reconstructedFile.WriteRow(buffer, rowWidth); i++; } if (displayText) cout << endl; //close reconstruction input and output // Do NOT close the verify file, otherwise data cannot be flushed to copy file // verifyFile.close(); reconstructedFile.close(); //get original and reconstructed checksum and datasum std::pair origChecksum = make_pair(origChecksumStr, atoi(origDatasum.c_str())); std::pair newChecksum = reconstructedFile.GetChecksumData(); //verify that no mistake was made if (origChecksum.second != newChecksum.second) { cout << "ERROR: datasums are NOT identical: " << (uint32_t)(origChecksum.second) << " vs " << (uint32_t)(newChecksum.second) << endl; return -1; } if (origChecksum.first != newChecksum.first) { cout << "WARNING: checksums are NOT Identical: " << origChecksum.first << " vs " << newChecksum.first << endl; } else { if (true) cout << "Ok" << endl; } delete[] buffer; return 0; }