1 | /*
|
---|
2 | * fitsCompressor.cc
|
---|
3 | *
|
---|
4 | * Created on: May 7, 2013
|
---|
5 | * Author: lyard
|
---|
6 | */
|
---|
7 |
|
---|
8 |
|
---|
9 | #include "Configuration.h"
|
---|
10 | #include "../externals/factfits.h"
|
---|
11 | #include "../externals/ofits.h"
|
---|
12 | #include "../externals/checksum.h"
|
---|
13 |
|
---|
14 | #include <map>
|
---|
15 | #include <fstream>
|
---|
16 | #include <sstream>
|
---|
17 | #include <iostream>
|
---|
18 |
|
---|
19 |
|
---|
20 | using namespace std;
|
---|
21 |
|
---|
22 | class CompressedFitsFile
|
---|
23 | {
|
---|
24 | public:
|
---|
25 | class HeaderEntry
|
---|
26 | {
|
---|
27 | public:
|
---|
28 | /**
|
---|
29 | * Default constructor
|
---|
30 | */
|
---|
31 | HeaderEntry(): _key(""),
|
---|
32 | _value(""),
|
---|
33 | _comment(""),
|
---|
34 | _fitsString("")
|
---|
35 | {
|
---|
36 | }
|
---|
37 |
|
---|
38 | /**
|
---|
39 | * Regular constructor.
|
---|
40 | * @param key the name of the keyword entry
|
---|
41 | * @param value its value
|
---|
42 | * @param comment an optionnal comment to be placed after the value
|
---|
43 | */
|
---|
44 | template<typename T>
|
---|
45 | HeaderEntry(const string& k,
|
---|
46 | const T& val,
|
---|
47 | const string& comm) : _key(k),
|
---|
48 | _value(""),
|
---|
49 | _comment(comm),
|
---|
50 | _fitsString("")
|
---|
51 | {
|
---|
52 | setValue(val);
|
---|
53 | }
|
---|
54 |
|
---|
55 | /**
|
---|
56 | * From fits.h
|
---|
57 | */
|
---|
58 | string Trim(const string &str, char c=' ')
|
---|
59 | {
|
---|
60 | // Trim Both leading and trailing spaces
|
---|
61 | const size_t pstart = str.find_first_not_of(c); // Find the first character position after excluding leading blank spaces
|
---|
62 | const size_t pend = str.find_last_not_of(c); // Find the first character position from reverse af
|
---|
63 |
|
---|
64 | // if all spaces or empty return an empty string
|
---|
65 | if (string::npos==pstart || string::npos==pend)
|
---|
66 | return string();
|
---|
67 |
|
---|
68 | return str.substr(pstart, pend-pstart+1);
|
---|
69 | }
|
---|
70 | /**
|
---|
71 | * Constructor from the original fits entry
|
---|
72 | */
|
---|
73 | HeaderEntry(const string& line)
|
---|
74 | {
|
---|
75 | _fitsString = line.data();
|
---|
76 |
|
---|
77 | //parse the line
|
---|
78 | _key = Trim(line.substr(0,8));
|
---|
79 | //COMMENT and/or HISTORY values
|
---|
80 | if (line.substr(8,2)!= "= ")
|
---|
81 | {
|
---|
82 | _value = "";
|
---|
83 | _comment = Trim(line.substr(10));
|
---|
84 | return;
|
---|
85 | }
|
---|
86 | string next=line.substr(10);
|
---|
87 | const size_t slash = next.find_first_of('/');
|
---|
88 | _value = Trim(Trim(Trim(next.substr(0, slash)), '\''));
|
---|
89 | _comment = Trim(next.substr(slash+1));
|
---|
90 | }
|
---|
91 | /**
|
---|
92 | * Alternative constroctor from the fits entry
|
---|
93 | */
|
---|
94 | HeaderEntry(const vector<char>& lineVec)
|
---|
95 | {
|
---|
96 | HeaderEntry(string(lineVec.data()));
|
---|
97 | }
|
---|
98 | /**
|
---|
99 | * Destructor
|
---|
100 | */
|
---|
101 | virtual ~HeaderEntry(){}
|
---|
102 |
|
---|
103 | const string& value() const {return _value;}
|
---|
104 | const string& key() const {return _key;}
|
---|
105 | const string& comment() const {return _comment;}
|
---|
106 | const string& fitsString() const {return _fitsString;}
|
---|
107 |
|
---|
108 | /**
|
---|
109 | * Set a keyword value.
|
---|
110 | * @param value The value to be set
|
---|
111 | * @param update whether the value already exist or not. To be modified soon.
|
---|
112 | */
|
---|
113 | template<typename T>
|
---|
114 | void setValue(const T& val)
|
---|
115 | {
|
---|
116 | ostringstream str;
|
---|
117 | str << val;
|
---|
118 | _value = str.str();
|
---|
119 | buildFitsString();
|
---|
120 | };
|
---|
121 | /**
|
---|
122 | * Set the comment for a given entry
|
---|
123 | */
|
---|
124 | void setComment(const string& comm)
|
---|
125 | {
|
---|
126 | _comment = comm;
|
---|
127 | buildFitsString();
|
---|
128 | }
|
---|
129 |
|
---|
130 | private:
|
---|
131 | /**
|
---|
132 | * Construct the FITS header string from the key, value and comment
|
---|
133 | */
|
---|
134 | void buildFitsString()
|
---|
135 | {
|
---|
136 | ostringstream str;
|
---|
137 | unsigned int totSize = 0;
|
---|
138 |
|
---|
139 | // Tuncate the key if required
|
---|
140 | if (_key.length() > 8)
|
---|
141 | {
|
---|
142 | str << _key.substr(0, 8);
|
---|
143 | totSize += 8;
|
---|
144 | }
|
---|
145 | else
|
---|
146 | {
|
---|
147 | str << _key;
|
---|
148 | totSize += _key.length();
|
---|
149 | }
|
---|
150 |
|
---|
151 | // Append space if key is less than 8 chars long
|
---|
152 | for (int i=totSize; i<8;i++)
|
---|
153 | {
|
---|
154 | str << " ";
|
---|
155 | totSize++;
|
---|
156 | }
|
---|
157 |
|
---|
158 | // Add separator
|
---|
159 | str << "= ";
|
---|
160 | totSize += 2;
|
---|
161 |
|
---|
162 | // Format value
|
---|
163 | if (_value.length() < 20)
|
---|
164 | for (;totSize<30-_value.length();totSize++)
|
---|
165 | str << " ";
|
---|
166 |
|
---|
167 | if (_value.length() > 70)
|
---|
168 | {
|
---|
169 | str << _value.substr(0,70);
|
---|
170 | totSize += 70;
|
---|
171 | }
|
---|
172 | else
|
---|
173 | {
|
---|
174 | str << _value;
|
---|
175 | totSize += _value.size();
|
---|
176 | }
|
---|
177 |
|
---|
178 | // If there is space remaining, add comment area
|
---|
179 | if (totSize < 77)
|
---|
180 | {
|
---|
181 | str << " / ";
|
---|
182 | totSize += 3;
|
---|
183 | if (totSize < 80)
|
---|
184 | {
|
---|
185 | unsigned int commentSize = 80 - totSize;
|
---|
186 | if (_comment.length() > commentSize)
|
---|
187 | {
|
---|
188 | str << _comment.substr(0,commentSize);
|
---|
189 | totSize += commentSize;
|
---|
190 | }
|
---|
191 | else
|
---|
192 | {
|
---|
193 | str << _comment;
|
---|
194 | totSize += _comment.length();
|
---|
195 | }
|
---|
196 | }
|
---|
197 | }
|
---|
198 |
|
---|
199 | // If there is yet some free space, fill up the entry with spaces
|
---|
200 | for (int i=totSize; i<80;i++)
|
---|
201 | str << " ";
|
---|
202 |
|
---|
203 | _fitsString = str.str();
|
---|
204 |
|
---|
205 | // Check for correct completion
|
---|
206 | if (_fitsString.length() != 80)
|
---|
207 | cout << "Error |" << _fitsString << "| is not of length 80" << endl;
|
---|
208 | }
|
---|
209 |
|
---|
210 | string _key; ///< the key (name) of the header entry
|
---|
211 | string _value; ///< the value of the header entry
|
---|
212 | string _comment; ///< the comment associated to the header entry
|
---|
213 | string _fitsString; ///< the string that will be written to the fits file
|
---|
214 | };
|
---|
215 |
|
---|
216 | /**
|
---|
217 | * Supported compressions
|
---|
218 | */
|
---|
219 | typedef enum
|
---|
220 | {
|
---|
221 | UNCOMPRESSED,
|
---|
222 | SMOOTHMAN
|
---|
223 | } FitsCompression;
|
---|
224 |
|
---|
225 | /**
|
---|
226 | * Columns class
|
---|
227 | */
|
---|
228 | class ColumnEntry
|
---|
229 | {
|
---|
230 | public:
|
---|
231 | /**
|
---|
232 | * Default constructor
|
---|
233 | */
|
---|
234 | ColumnEntry();
|
---|
235 |
|
---|
236 | /**
|
---|
237 | * Default destructor
|
---|
238 | */
|
---|
239 | virtual ~ColumnEntry(){}
|
---|
240 |
|
---|
241 | /**
|
---|
242 | * Constructor from values
|
---|
243 | * @param n the column name
|
---|
244 | * @param t the column type
|
---|
245 | * @param numof the number of entries in the column
|
---|
246 | * @param comp the compression type for this column
|
---|
247 | */
|
---|
248 | ColumnEntry(const string& n,
|
---|
249 | char t,
|
---|
250 | int numOf,
|
---|
251 | uint32_t comp) : _name(n),
|
---|
252 | _num(numOf),
|
---|
253 | _typeSize(0),
|
---|
254 | _offset(0),
|
---|
255 | _type(t),
|
---|
256 | _description(""),
|
---|
257 | _compression(comp)
|
---|
258 | {
|
---|
259 | switch (t)
|
---|
260 | {
|
---|
261 | case 'L':
|
---|
262 | case 'A':
|
---|
263 | case 'B': _typeSize = 1; break;
|
---|
264 | case 'I': _typeSize = 2; break;
|
---|
265 | case 'J':
|
---|
266 | case 'E': _typeSize = 4; break;
|
---|
267 | case 'K':
|
---|
268 | case 'D': _typeSize = 8; break;
|
---|
269 | default:
|
---|
270 | cout << "Error: typename " << t << " missing in the current implementation" << endl;
|
---|
271 | };
|
---|
272 |
|
---|
273 | ostringstream str;
|
---|
274 | str << "data format of field: ";
|
---|
275 |
|
---|
276 | switch (t)
|
---|
277 | {
|
---|
278 | case 'L': str << "1-byte BOOL"; break;
|
---|
279 | case 'A': str << "1-byte CHAR"; break;
|
---|
280 | case 'B': str << "BYTE"; break;
|
---|
281 | case 'I': str << "2-byte INTEGER"; break;
|
---|
282 | case 'J': str << "4-byte INTEGER"; break;
|
---|
283 | case 'K': str << "8-byte INTEGER"; break;
|
---|
284 | case 'E': str << "4-byte FLOAT"; break;
|
---|
285 | case 'D': str << "8-byte FLOAT"; break;
|
---|
286 | }
|
---|
287 |
|
---|
288 | _description = str.str();
|
---|
289 | }
|
---|
290 |
|
---|
291 | const string& name() const { return _name;};
|
---|
292 | int width() const { return _num*_typeSize;};
|
---|
293 | int offset() const { return _offset; };
|
---|
294 | int numElems() const { return _num; };
|
---|
295 | int sizeOfElems() const { return _typeSize;};
|
---|
296 | void setOffset(int off) { _offset = off;};
|
---|
297 | char type() const { return _type;};
|
---|
298 | string getDescription() const { return _description;}
|
---|
299 | uint32_t getCompression() const { return _compression;}
|
---|
300 |
|
---|
301 | string getCompressionString() const
|
---|
302 | {
|
---|
303 | switch (_compression)
|
---|
304 | {
|
---|
305 | case UNCOMPRESSED: return "RAW";
|
---|
306 | case SMOOTHMAN: return "SMOO";
|
---|
307 | };
|
---|
308 | return "UNKNOWN";
|
---|
309 | }
|
---|
310 |
|
---|
311 | private:
|
---|
312 |
|
---|
313 | string _name; ///< name of the column
|
---|
314 | int _num; ///< number of elements contained in one row of this column
|
---|
315 | int _typeSize; ///< the number of bytes taken by one element
|
---|
316 | int _offset; ///< the offset of the column, in bytes, from the beginning of one row
|
---|
317 | char _type; ///< the type of the column, as specified by the fits documentation
|
---|
318 | string _description; ///< a description for the column. It will be placed in the header
|
---|
319 | uint32_t _compression; ///< the compression of the column. Only looked for if used by the compressed writer
|
---|
320 | };
|
---|
321 |
|
---|
322 | public:
|
---|
323 | ///@brief default constructor. Assigns a default number of rows and tiles
|
---|
324 | CompressedFitsFile(uint32_t numTiles=100, uint32_t numRowsPerTile=100);
|
---|
325 |
|
---|
326 | ///@brief default destructor
|
---|
327 | virtual ~CompressedFitsFile();
|
---|
328 |
|
---|
329 | ///@brief get the header of the file
|
---|
330 | vector<HeaderEntry>& getHeaderEntries() { return _header;}
|
---|
331 |
|
---|
332 | protected:
|
---|
333 | ///@brief protected function to allocate the intermediate buffers
|
---|
334 | bool reallocateBuffers();
|
---|
335 |
|
---|
336 | //FITS related stuff
|
---|
337 | vector<HeaderEntry> _header; ///< Header keys
|
---|
338 | vector<ColumnEntry> _columns; ///< Columns in the file
|
---|
339 | uint32_t _numTiles; ///< Number of tiles (i.e. groups of rows)
|
---|
340 | uint32_t _numRowsPerTile; ///< Number of rows per tile
|
---|
341 | uint32_t _totalNumRows; ///< Total number of raws
|
---|
342 | uint32_t _rowWidth; ///< Total number of bytes in one row
|
---|
343 | bool _headerFlushed; ///< Flag telling whether the header record is synchronized with the data on disk
|
---|
344 | char* _buffer; ///< Memory buffer to store rows while they are not compressed
|
---|
345 | Checksum _checksum; ///< Checksum for asserting the consistency of the data
|
---|
346 | fstream _file; ///< The actual file streamer for accessing disk data
|
---|
347 |
|
---|
348 | //compression related stuff
|
---|
349 | typedef pair<int64_t, int64_t> CatalogEntry;
|
---|
350 | typedef vector<CatalogEntry> CatalogRow;
|
---|
351 | typedef vector<CatalogRow> CatalogType;
|
---|
352 | CatalogType _catalog; ///< Catalog, i.e. the main table that points to the compressed data.
|
---|
353 | uint64_t _heapPtr; ///< the address in the file of the heap area
|
---|
354 | vector<char*> _transposedBuffer; ///< Memory buffer to store rows while they are transposed
|
---|
355 | vector<char*> _compressedBuffer; ///< Memory buffer to store rows while they are compressed
|
---|
356 |
|
---|
357 | //thread related stuff
|
---|
358 | uint32_t _numThreads; ///< The number of threads that will be used to compress
|
---|
359 | uint32_t _threadIndex; ///< A variable to assign threads indices
|
---|
360 | vector<pthread_t> _thread; ///< The thread handler of the compressor
|
---|
361 | vector<uint32_t> _threadNumRows; ///< Total number of rows for thread to compress
|
---|
362 | vector<uint32_t> _threadStatus; ///< Flag telling whether the buffer to be transposed (and compressed) is full or empty
|
---|
363 |
|
---|
364 | //thread states. Not all used, but they do not hurt
|
---|
365 | static const uint32_t _THREAD_WAIT_; ///< Thread doing nothing
|
---|
366 | static const uint32_t _THREAD_COMPRESS_; ///< Thread working, compressing
|
---|
367 | static const uint32_t _THREAD_DECOMPRESS_; ///< Thread working, decompressing
|
---|
368 | static const uint32_t _THREAD_WRITE_; ///< Thread writing data to disk
|
---|
369 | static const uint32_t _THREAD_READ_; ///< Thread reading data from disk
|
---|
370 | static const uint32_t _THREAD_EXIT_; ///< Thread exiting
|
---|
371 |
|
---|
372 | static HeaderEntry _dummyHeaderEntry; ///< Dummy entry for returning if requested on is not found
|
---|
373 | };
|
---|
374 |
|
---|
375 | class CompressedFitsWriter : public CompressedFitsFile
|
---|
376 | {
|
---|
377 | public:
|
---|
378 | ///@brief Default constructor. 100 tiles of 100 rows each are assigned by default
|
---|
379 | CompressedFitsWriter(uint32_t numTiles=100, uint32_t numRowsPerTile=100);
|
---|
380 |
|
---|
381 | ///@brief default destructor
|
---|
382 | virtual ~CompressedFitsWriter();
|
---|
383 |
|
---|
384 | ///@brief add one column to the file
|
---|
385 | bool addColumn(const ColumnEntry& column);
|
---|
386 |
|
---|
387 | ///@brief sets a given header key
|
---|
388 | bool setHeaderKey(const HeaderEntry&);
|
---|
389 |
|
---|
390 | ///@brief open a new fits file
|
---|
391 | bool open(const string& fileName, const string& tableName="Data");
|
---|
392 |
|
---|
393 | ///@brief close the opened file
|
---|
394 | bool close();
|
---|
395 |
|
---|
396 | ///@brief write one row of data, already placed in bufferToWrite. Does the byte-swapping
|
---|
397 | bool writeBinaryRow(const char* bufferToWrite);
|
---|
398 |
|
---|
399 | ///@brief assign a given (already loaded) drs calibration
|
---|
400 | void setDrsCalib(int16_t* data);
|
---|
401 |
|
---|
402 | ///@brief set the number of worker threads compressing the data.
|
---|
403 | bool setNumWorkingThreads(uint32_t num);
|
---|
404 |
|
---|
405 | private:
|
---|
406 | ///@brief compresses one buffer of data, the one given by threadIndex
|
---|
407 | uint32_t compressBuffer(uint32_t threadIndex);
|
---|
408 |
|
---|
409 | ///@brief writes an already compressed buffer to disk
|
---|
410 | bool writeCompressedDataToDisk(uint32_t threadID, uint32_t sizeToWrite);
|
---|
411 |
|
---|
412 | ///@brief add the header checksum to the datasum
|
---|
413 | void addHeaderChecksum(Checksum& checksum);
|
---|
414 |
|
---|
415 | ///@brief write the header. If closingFile is set to true, checksum is calculated
|
---|
416 | void writeHeader(bool closingFile = false);
|
---|
417 |
|
---|
418 | ///@brief write the compressed data catalog. If closingFile is set to true, checksum is calculated
|
---|
419 | void writeCatalog(bool closingFile=false);
|
---|
420 |
|
---|
421 | /// FIXME this was a bad idea. Move everything to the regular header
|
---|
422 | vector<HeaderEntry> _defaultHeader;
|
---|
423 |
|
---|
424 | /// the main function compressing the data
|
---|
425 | static void* threadFunction(void* context);
|
---|
426 |
|
---|
427 | /// Write the drs calibration to disk, if any
|
---|
428 | void writeDrsCalib();
|
---|
429 |
|
---|
430 | /// Copy and transpose (or semi-transpose) one tile of data
|
---|
431 | void copyTransposeTile(uint32_t index);
|
---|
432 |
|
---|
433 | /// Specific compression functions
|
---|
434 | uint32_t compressUNCOMPRESSED(char* dest, const char* src, uint32_t numRows, uint32_t sizeOfElems, uint32_t numRowElems);
|
---|
435 | uint32_t compressHUFFMAN(char* dest, const char* src, uint32_t numRows, uint32_t sizeOfElems, uint32_t numRowElems);
|
---|
436 | uint32_t compressSMOOTHMAN(char* dest, char* src, uint32_t numRows, uint32_t sizeOfElems, uint32_t numRowElems);
|
---|
437 |
|
---|
438 | int32_t _checkOffset; ///< offset to the data pointer to calculate the checksum
|
---|
439 | int16_t* _drsCalibData; ///< array of the Drs baseline mean
|
---|
440 | int32_t _threadLooper; ///< Which thread will deal with the upcoming bunch of data ?
|
---|
441 | pthread_mutex_t _mutex; ///< mutex for compressing threads
|
---|
442 |
|
---|
443 |
|
---|
444 | static string _emptyBlock; ///< an empty block to be apened at the end of a file so that its length is a multiple of 2880
|
---|
445 | static string _fitsHeader; ///< the default header to be written in every fits file
|
---|
446 | };
|
---|
447 |
|
---|
448 | const uint32_t CompressedFitsFile::_THREAD_WAIT_ = 0;
|
---|
449 | const uint32_t CompressedFitsFile::_THREAD_COMPRESS_ = 1;
|
---|
450 | const uint32_t CompressedFitsFile::_THREAD_DECOMPRESS_= 2;
|
---|
451 | const uint32_t CompressedFitsFile::_THREAD_WRITE_ = 3;
|
---|
452 | const uint32_t CompressedFitsFile::_THREAD_READ_ = 4;
|
---|
453 | const uint32_t CompressedFitsFile::_THREAD_EXIT_ = 5;
|
---|
454 |
|
---|
455 | #define COMPRESSED_FLAG 0x1
|
---|
456 | #define UNCOMPRESSED_FLAG 0x0
|
---|
457 |
|
---|
458 | template<>
|
---|
459 | void CompressedFitsFile::HeaderEntry::setValue(const string& v)
|
---|
460 | {
|
---|
461 | string val = v;
|
---|
462 | if (val.size() > 2 && val[0] == '\'')
|
---|
463 | {
|
---|
464 | size_t pos = val.find_last_of("'");
|
---|
465 | if (pos != string::npos && pos != 0)
|
---|
466 | val = val.substr(1, pos-1);
|
---|
467 | }
|
---|
468 | ostringstream str;
|
---|
469 |
|
---|
470 | str << "'" << val << "'";
|
---|
471 | for (int i=str.str().length(); i<20;i++)
|
---|
472 | str << " ";
|
---|
473 | _value = str.str();
|
---|
474 | buildFitsString();
|
---|
475 | }
|
---|
476 |
|
---|
477 | /**
|
---|
478 | * Default header to be written in all fits files
|
---|
479 | */
|
---|
480 | string CompressedFitsWriter::_fitsHeader = "SIMPLE = T / file does conform to FITS standard "
|
---|
481 | "BITPIX = 8 / number of bits per data pixel "
|
---|
482 | "NAXIS = 0 / number of data axes "
|
---|
483 | "EXTEND = T / FITS dataset may contain extensions "
|
---|
484 | "CHECKSUM= '4AcB48bA4AbA45bA' / Checksum for the whole HDU "
|
---|
485 | "DATASUM = ' 0' / Checksum for the data block "
|
---|
486 | "COMMENT FITS (Flexible Image Transport System) format is defined in 'Astronomy"
|
---|
487 | "COMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H "
|
---|
488 | "END "
|
---|
489 | " "
|
---|
490 | " "
|
---|
491 | " "
|
---|
492 | " "
|
---|
493 | " "
|
---|
494 | " "
|
---|
495 | " "
|
---|
496 | " "
|
---|
497 | " "
|
---|
498 | " "
|
---|
499 | " "
|
---|
500 | " "
|
---|
501 | " "
|
---|
502 | " "
|
---|
503 | " "
|
---|
504 | " "
|
---|
505 | " "
|
---|
506 | " "
|
---|
507 | " "
|
---|
508 | " "
|
---|
509 | " "
|
---|
510 | " "
|
---|
511 | " "
|
---|
512 | " "
|
---|
513 | " "
|
---|
514 | " "
|
---|
515 | " ";
|
---|
516 |
|
---|
517 | /**
|
---|
518 | * Empty block to be appenned at the end of files, so that the length matches multiple of 2880 bytes
|
---|
519 | *
|
---|
520 | */
|
---|
521 | string CompressedFitsWriter::_emptyBlock = " "
|
---|
522 | " "
|
---|
523 | " "
|
---|
524 | " "
|
---|
525 | " "
|
---|
526 | " "
|
---|
527 | " "
|
---|
528 | " "
|
---|
529 | " "
|
---|
530 | " "
|
---|
531 | " "
|
---|
532 | " "
|
---|
533 | " "
|
---|
534 | " "
|
---|
535 | " "
|
---|
536 | " "
|
---|
537 | " "
|
---|
538 | " "
|
---|
539 | " "
|
---|
540 | " "
|
---|
541 | " "
|
---|
542 | " "
|
---|
543 | " "
|
---|
544 | " "
|
---|
545 | " "
|
---|
546 | " "
|
---|
547 | " "
|
---|
548 | " "
|
---|
549 | " "
|
---|
550 | " "
|
---|
551 | " "
|
---|
552 | " "
|
---|
553 | " "
|
---|
554 | " "
|
---|
555 | " "
|
---|
556 | " ";
|
---|
557 |
|
---|
558 | CompressedFitsFile::HeaderEntry CompressedFitsFile::_dummyHeaderEntry;
|
---|
559 | /****************************************************************
|
---|
560 | * SUPER CLASS DEFAULT CONSTRUCTOR
|
---|
561 | ****************************************************************/
|
---|
562 | CompressedFitsFile::CompressedFitsFile(uint32_t numTiles,
|
---|
563 | uint32_t numRowsPerTile) : _header(),
|
---|
564 | _columns(),
|
---|
565 | _numTiles(numTiles),
|
---|
566 | _numRowsPerTile(numRowsPerTile),
|
---|
567 | _totalNumRows(0),
|
---|
568 | _rowWidth(0),
|
---|
569 | _headerFlushed(false),
|
---|
570 | _buffer(NULL),
|
---|
571 | _checksum(0),
|
---|
572 | _heapPtr(0),
|
---|
573 | _transposedBuffer(1),
|
---|
574 | _compressedBuffer(1),
|
---|
575 | _numThreads(1),
|
---|
576 | _threadIndex(0),
|
---|
577 | _thread(1),
|
---|
578 | _threadNumRows(1),
|
---|
579 | _threadStatus(1)
|
---|
580 | {
|
---|
581 | _catalog.resize(_numTiles);
|
---|
582 | _transposedBuffer[0] = NULL;
|
---|
583 | _compressedBuffer[0] = NULL;
|
---|
584 | _threadStatus[0] = _THREAD_WAIT_;
|
---|
585 | _threadNumRows[0] = 0;
|
---|
586 | }
|
---|
587 |
|
---|
588 | /****************************************************************
|
---|
589 | * SUPER CLASS DEFAULT DESTRUCTOR
|
---|
590 | ****************************************************************/
|
---|
591 | CompressedFitsFile::~CompressedFitsFile()
|
---|
592 | {
|
---|
593 | if (_buffer != NULL)
|
---|
594 | {
|
---|
595 | delete[] _buffer;
|
---|
596 | _buffer = NULL;
|
---|
597 | for (uint32_t i=0;i<_numThreads;i++)
|
---|
598 | {
|
---|
599 | _compressedBuffer[i] = _compressedBuffer[i]-4;
|
---|
600 | delete[] _transposedBuffer[i];
|
---|
601 | delete[] _compressedBuffer[i];
|
---|
602 | _transposedBuffer[i] = NULL;
|
---|
603 | _compressedBuffer[i] = NULL;
|
---|
604 | }
|
---|
605 | }
|
---|
606 | if (_file.is_open())
|
---|
607 | _file.close();
|
---|
608 | }
|
---|
609 |
|
---|
610 | /****************************************************************
|
---|
611 | * REALLOCATE BUFFER
|
---|
612 | ****************************************************************/
|
---|
613 | bool CompressedFitsFile::reallocateBuffers()
|
---|
614 | {
|
---|
615 | if (_buffer)
|
---|
616 | {
|
---|
617 | delete[] _buffer;
|
---|
618 | for (uint32_t i=0;i<_compressedBuffer.size();i++)
|
---|
619 | {
|
---|
620 | _compressedBuffer[i] = _compressedBuffer[i]-4;
|
---|
621 | delete[] _transposedBuffer[i];
|
---|
622 | delete[] _compressedBuffer[i];
|
---|
623 | }
|
---|
624 | }
|
---|
625 | _buffer = new char[_rowWidth*_numRowsPerTile];
|
---|
626 | if (_buffer == NULL) return false;
|
---|
627 | if (_compressedBuffer.size() != _numThreads)
|
---|
628 | {
|
---|
629 | _transposedBuffer.resize(_numThreads);
|
---|
630 | _compressedBuffer.resize(_numThreads);
|
---|
631 | }
|
---|
632 | for (uint32_t i=0;i<_numThreads;i++)
|
---|
633 | {
|
---|
634 | _transposedBuffer[i] = new char[_rowWidth*_numRowsPerTile];
|
---|
635 | _compressedBuffer[i] = new char[_rowWidth*_numRowsPerTile + 1024*1024]; //use a bit more memory, in case the compression algorithms uses more
|
---|
636 | if (_transposedBuffer[i] == NULL || _compressedBuffer[i] == NULL)
|
---|
637 | return false;
|
---|
638 | //shift the compressed buffer by 4 bytes, for checksum calculation
|
---|
639 | memset(_compressedBuffer[i], 0, 4);
|
---|
640 | _compressedBuffer[i] = _compressedBuffer[i]+4;
|
---|
641 | }
|
---|
642 | return true;
|
---|
643 | }
|
---|
644 |
|
---|
645 | /****************************************************************
|
---|
646 | * DEFAULT WRITER CONSTRUCTOR
|
---|
647 | ****************************************************************/
|
---|
648 | CompressedFitsWriter::CompressedFitsWriter(uint32_t numTiles,
|
---|
649 | uint32_t numRowsPerTile) : CompressedFitsFile(numTiles, numRowsPerTile),
|
---|
650 | _checkOffset(0),
|
---|
651 | _drsCalibData(NULL),
|
---|
652 | _threadLooper(0)
|
---|
653 | {
|
---|
654 | _defaultHeader.push_back(HeaderEntry("XTENSION", "'BINTABLE' ", "binary table extension"));
|
---|
655 | _defaultHeader.push_back(HeaderEntry("BITPIX", 8, "8-bit bytes"));
|
---|
656 | _defaultHeader.push_back(HeaderEntry("NAXIS", 2, "2-dimensional binary table"));
|
---|
657 | _defaultHeader.push_back(HeaderEntry("NAXIS1", _rowWidth, "width of table in bytes"));
|
---|
658 | _defaultHeader.push_back(HeaderEntry("NAXIS2", numTiles, "num of rows in table"));
|
---|
659 | _defaultHeader.push_back(HeaderEntry("PCOUNT", 0, "size of special data area"));
|
---|
660 | _defaultHeader.push_back(HeaderEntry("GCOUNT", 1, "one data group (required keyword)"));
|
---|
661 | _defaultHeader.push_back(HeaderEntry("TFIELDS", _columns.size(), "number of fields in each row"));
|
---|
662 | _defaultHeader.push_back(HeaderEntry("CHECKSUM", "'0000000000000000' ", "Checksum for the whole HDU"));
|
---|
663 | _defaultHeader.push_back(HeaderEntry("DATASUM", " 0", "Checksum for the data block"));
|
---|
664 | //compression stuff
|
---|
665 | _defaultHeader.push_back(HeaderEntry("ZTABLE", "T", "Table is compressed"));
|
---|
666 | _defaultHeader.push_back(HeaderEntry("ZNAXIS1", 0, "Width of uncompressed rows"));
|
---|
667 | _defaultHeader.push_back(HeaderEntry("ZNAXIS2", 0, "Number of uncompressed rows"));
|
---|
668 | _defaultHeader.push_back(HeaderEntry("ZPCOUNT", 0, ""));
|
---|
669 | _defaultHeader.push_back(HeaderEntry("ZHEAPPTR", 0, ""));
|
---|
670 | _defaultHeader.push_back(HeaderEntry("ZTILELEN", numRowsPerTile, "Number of rows per tile"));
|
---|
671 | _defaultHeader.push_back(HeaderEntry("THEAP", 0, ""));
|
---|
672 |
|
---|
673 | pthread_mutex_init(&_mutex, NULL);
|
---|
674 | }
|
---|
675 |
|
---|
676 | /****************************************************************
|
---|
677 | * DEFAULT DESTRUCTOR
|
---|
678 | ****************************************************************/
|
---|
679 | CompressedFitsWriter::~CompressedFitsWriter()
|
---|
680 | {
|
---|
681 | pthread_mutex_destroy(&_mutex);
|
---|
682 | }
|
---|
683 |
|
---|
684 | /****************************************************************
|
---|
685 | * SET THE POINTER TO THE DRS CALIBRATION
|
---|
686 | ****************************************************************/
|
---|
687 | void CompressedFitsWriter::setDrsCalib(int16_t* data)
|
---|
688 | {
|
---|
689 | _drsCalibData = data;
|
---|
690 | }
|
---|
691 |
|
---|
692 | /****************************************************************
|
---|
693 | * SET NUM WORKING THREADS
|
---|
694 | ****************************************************************/
|
---|
695 | bool CompressedFitsWriter::setNumWorkingThreads(uint32_t num)
|
---|
696 | {
|
---|
697 | if (_file.is_open())
|
---|
698 | return false;
|
---|
699 | if (num < 1 || num > 64)
|
---|
700 | {
|
---|
701 | cout << "ERROR: num threads must be between 1 and 64. Ignoring" << endl;
|
---|
702 | return false;
|
---|
703 | }
|
---|
704 | _numThreads = num;
|
---|
705 | _transposedBuffer[0] = NULL;
|
---|
706 | _compressedBuffer[0] = NULL;
|
---|
707 | _threadStatus.resize(num);
|
---|
708 | _thread.resize(num);
|
---|
709 | _threadNumRows.resize(num);
|
---|
710 | for (uint32_t i=0;i<num;i++)
|
---|
711 | {
|
---|
712 | _threadNumRows[i] = 0;
|
---|
713 | _threadStatus[i] = _THREAD_WAIT_;
|
---|
714 | }
|
---|
715 | return reallocateBuffers();
|
---|
716 | }
|
---|
717 |
|
---|
718 | /****************************************************************
|
---|
719 | * WRITE DRS CALIBRATION TO FILE
|
---|
720 | ****************************************************************/
|
---|
721 | void CompressedFitsWriter::writeDrsCalib()
|
---|
722 | {
|
---|
723 | //if file was not loaded, ignore
|
---|
724 | if (_drsCalibData == NULL)
|
---|
725 | return;
|
---|
726 | uint64_t whereDidIStart = _file.tellp();
|
---|
727 | vector<HeaderEntry> header;
|
---|
728 | header.push_back(HeaderEntry("XTENSION", "'BINTABLE' ", "binary table extension"));
|
---|
729 | header.push_back(HeaderEntry("BITPIX" , 8 , "8-bit bytes"));
|
---|
730 | header.push_back(HeaderEntry("NAXIS" , 2 , "2-dimensional binary table"));
|
---|
731 | header.push_back(HeaderEntry("NAXIS1" , 1024*1440*2 , "width of table in bytes"));
|
---|
732 | header.push_back(HeaderEntry("NAXIS2" , 1 , "number of rows in table"));
|
---|
733 | header.push_back(HeaderEntry("PCOUNT" , 0 , "size of special data area"));
|
---|
734 | header.push_back(HeaderEntry("GCOUNT" , 1 , "one data group (required keyword)"));
|
---|
735 | header.push_back(HeaderEntry("TFIELDS" , 1 , "number of fields in each row"));
|
---|
736 | header.push_back(HeaderEntry("CHECKSUM", "'0000000000000000' ", "Checksum for the whole HDU"));
|
---|
737 | header.push_back(HeaderEntry("DATASUM" , " 0" , "Checksum for the data block"));
|
---|
738 | header.push_back(HeaderEntry("EXTNAME" , "'DrsCalib' ", "name of this binary table extension"));
|
---|
739 | header.push_back(HeaderEntry("TTYPE1" , "'OffsetCalibration' ", "label for field 1"));
|
---|
740 | header.push_back(HeaderEntry("TFORM1" , "'1474560I' ", "data format of field: 2-byte INTEGER"));
|
---|
741 |
|
---|
742 | for (uint32_t i=0;i<header.size();i++)
|
---|
743 | _file.write(header[i].fitsString().c_str(), 80);
|
---|
744 | //End the header
|
---|
745 | _file.write("END ", 80);
|
---|
746 | long here = _file.tellp();
|
---|
747 | if (here%2880)
|
---|
748 | _file.write(_emptyBlock.c_str(), 2880 - here%2880);
|
---|
749 | //now write the data itself
|
---|
750 | int16_t* swappedBytes = new int16_t[1024];
|
---|
751 | Checksum checksum;
|
---|
752 | for (int32_t i=0;i<1440;i++)
|
---|
753 | {
|
---|
754 | memcpy(swappedBytes, &(_drsCalibData[i*1024]), 2048);
|
---|
755 | for (int32_t j=0;j<2048;j+=2)
|
---|
756 | {
|
---|
757 | int8_t inter;
|
---|
758 | inter = reinterpret_cast<int8_t*>(swappedBytes)[j];
|
---|
759 | reinterpret_cast<int8_t*>(swappedBytes)[j] = reinterpret_cast<int8_t*>(swappedBytes)[j+1];
|
---|
760 | reinterpret_cast<int8_t*>(swappedBytes)[j+1] = inter;
|
---|
761 | }
|
---|
762 | _file.write(reinterpret_cast<char*>(swappedBytes), 2048);
|
---|
763 | checksum.add(reinterpret_cast<char*>(swappedBytes), 2048);
|
---|
764 | }
|
---|
765 | uint64_t whereDidIStop = _file.tellp();
|
---|
766 | delete[] swappedBytes;
|
---|
767 | //No need to pad the data, as (1440*1024*2)%2880==0
|
---|
768 |
|
---|
769 | //calculate the checksum from the header
|
---|
770 | ostringstream str;
|
---|
771 | str << checksum.val();
|
---|
772 | header[9] = HeaderEntry("DATASUM", str.str(), "Checksum for the data block");
|
---|
773 | for (vector<HeaderEntry>::iterator it=header.begin();it!=header.end(); it++)
|
---|
774 | checksum.add(it->fitsString().c_str(), 80);
|
---|
775 | string end("END ");
|
---|
776 | string space(" ");
|
---|
777 | checksum.add(end.c_str(), 80);
|
---|
778 | int headerRowsLeft = 36 - (header.size() + 1)%36;
|
---|
779 | for (int i=0;i<headerRowsLeft;i++)
|
---|
780 | checksum.add(space.c_str(), 80);
|
---|
781 | //udpate the checksum keyword
|
---|
782 | header[8] = HeaderEntry("CHECKSUM", checksum.str(), "Checksum for the whole HDU");
|
---|
783 | //and eventually re-write the header data
|
---|
784 | _file.seekp(whereDidIStart);
|
---|
785 | for (uint32_t i=0;i<header.size();i++)
|
---|
786 | _file.write(header[i].fitsString().c_str(), 80);
|
---|
787 | _file.seekp(whereDidIStop);
|
---|
788 | }
|
---|
789 |
|
---|
790 | /****************************************************************
|
---|
791 | * ADD COLUMN
|
---|
792 | ****************************************************************/
|
---|
793 | bool CompressedFitsWriter::addColumn(const ColumnEntry& column)
|
---|
794 | {
|
---|
795 | if (_totalNumRows != 0)
|
---|
796 | {
|
---|
797 | cout << "Error: cannot add new columns once first row has been written" << endl;
|
---|
798 | return false;
|
---|
799 | }
|
---|
800 | for (vector<ColumnEntry>::iterator it=_columns.begin(); it != _columns.end(); it++)
|
---|
801 | {
|
---|
802 | if (it->name() == column.name())
|
---|
803 | {
|
---|
804 | cout << "Warning: column already exist (" << column.name() << "). Ignoring" << endl;
|
---|
805 | return false;
|
---|
806 | }
|
---|
807 | }
|
---|
808 | _columns.push_back(column);
|
---|
809 | _columns.back().setOffset(_rowWidth);
|
---|
810 | _rowWidth += column.width();
|
---|
811 | reallocateBuffers();
|
---|
812 |
|
---|
813 | ostringstream str, str2, str3;
|
---|
814 | str << "TTYPE" << _columns.size();
|
---|
815 | str2 << column.name();
|
---|
816 | str3 << "label for field ";
|
---|
817 | if (_columns.size() < 10) str3 << " ";
|
---|
818 | if (_columns.size() < 100) str3 << " ";
|
---|
819 | str3 << _columns.size();
|
---|
820 | setHeaderKey(HeaderEntry(str.str(), str2.str(), str3.str()));
|
---|
821 | str.str("");
|
---|
822 | str2.str("");
|
---|
823 | str3.str("");
|
---|
824 | str << "TFORM" << _columns.size();
|
---|
825 | str2 << "1QB";
|
---|
826 | str3 << "data format of field " << _columns.size();
|
---|
827 | setHeaderKey(HeaderEntry(str.str(), str2.str(), str3.str()));
|
---|
828 | str.str("");
|
---|
829 | str2.str("");
|
---|
830 | str3.str("");
|
---|
831 | str << "ZFORM" << _columns.size();
|
---|
832 | str2 << column.numElems() << column.type();
|
---|
833 | str3 << "Original format of field " << _columns.size();
|
---|
834 | setHeaderKey(HeaderEntry(str.str(), str2.str(), str3.str()));
|
---|
835 | str.str("");
|
---|
836 | str2.str("");
|
---|
837 | str3.str("");
|
---|
838 | str << "ZCTYP" << _columns.size();
|
---|
839 | str2 << column.getCompressionString();
|
---|
840 | str3 << "Comp. Scheme of field " << _columns.size();
|
---|
841 | setHeaderKey(HeaderEntry(str.str(), str2.str(), str3.str()));
|
---|
842 | //resize the catalog vector accordingly
|
---|
843 | for (uint32_t i=0;i<_numTiles;i++)
|
---|
844 | {
|
---|
845 | _catalog[i].resize(_columns.size());
|
---|
846 | for (uint32_t j=0;j<_catalog[i].size();j++)
|
---|
847 | _catalog[i][j] = make_pair(0,0);
|
---|
848 | }
|
---|
849 | return true;
|
---|
850 | }
|
---|
851 |
|
---|
852 | /****************************************************************
|
---|
853 | * SET HEADER KEY
|
---|
854 | ****************************************************************/
|
---|
855 | bool CompressedFitsWriter::setHeaderKey(const HeaderEntry& entry)
|
---|
856 | {
|
---|
857 | HeaderEntry ent = entry;
|
---|
858 | for (vector<HeaderEntry>::iterator it=_header.begin(); it != _header.end(); it++)
|
---|
859 | {
|
---|
860 | if (it->key() == entry.key())
|
---|
861 | {
|
---|
862 | if (entry.comment() == "")
|
---|
863 | ent.setComment(it->comment());
|
---|
864 | (*it) = ent;
|
---|
865 | _headerFlushed = false;
|
---|
866 | return true;
|
---|
867 | }
|
---|
868 | }
|
---|
869 | for (vector<HeaderEntry>::iterator it=_defaultHeader.begin(); it != _defaultHeader.end(); it++)
|
---|
870 | {
|
---|
871 | if (it->key() == entry.key())
|
---|
872 | {
|
---|
873 | if (entry.comment() == "")
|
---|
874 | ent.setComment(it->comment());
|
---|
875 | (*it) = ent;
|
---|
876 | _headerFlushed = false;
|
---|
877 | return true;
|
---|
878 | }
|
---|
879 | }
|
---|
880 | if (_totalNumRows != 0)
|
---|
881 | {
|
---|
882 | cout << "Error: new header keys (" << entry.key() << ") must be set before the first row is written. Ignoring." << endl;
|
---|
883 | return false;
|
---|
884 | }
|
---|
885 | _header.push_back(entry);
|
---|
886 | _headerFlushed = false;
|
---|
887 | return true;
|
---|
888 | }
|
---|
889 |
|
---|
890 | /****************************************************************
|
---|
891 | * OPEN
|
---|
892 | ****************************************************************/
|
---|
893 | bool CompressedFitsWriter::open(const string& fileName, const string& tableName)
|
---|
894 | {
|
---|
895 | _file.open(fileName.c_str(), ios_base::out);
|
---|
896 | if (!_file.is_open())
|
---|
897 | {
|
---|
898 | cout << "Error: Could not open the file (" << fileName << ")." << endl;
|
---|
899 | return false;
|
---|
900 | }
|
---|
901 | _defaultHeader.push_back(HeaderEntry("EXTNAME", tableName, "name of this binary table extension"));
|
---|
902 | _headerFlushed = false;
|
---|
903 | _threadIndex = 0;
|
---|
904 | //create requested number of threads
|
---|
905 | for (uint32_t i=0;i<_numThreads;i++)
|
---|
906 | pthread_create(&(_thread[i]), NULL, threadFunction, this);
|
---|
907 | //wait for the threads to start
|
---|
908 | while (_numThreads != _threadIndex)
|
---|
909 | usleep(1000);
|
---|
910 | //set the writing fence to the last thread
|
---|
911 | _threadIndex = _numThreads-1;
|
---|
912 | return (_file.good());
|
---|
913 | }
|
---|
914 |
|
---|
915 | /****************************************************************
|
---|
916 | * WRITE HEADER
|
---|
917 | ****************************************************************/
|
---|
918 | void CompressedFitsWriter::writeHeader(bool closingFile)
|
---|
919 | {
|
---|
920 | if (_headerFlushed)
|
---|
921 | return;
|
---|
922 | if (!_file.is_open())
|
---|
923 | return;
|
---|
924 |
|
---|
925 | long cPos = _file.tellp();
|
---|
926 |
|
---|
927 | _file.seekp(0);
|
---|
928 |
|
---|
929 | _file.write(_fitsHeader.c_str(), 2880);
|
---|
930 |
|
---|
931 | //Write the DRS calib table here !
|
---|
932 | writeDrsCalib();
|
---|
933 |
|
---|
934 | //we are now at the beginning of the main table. Write its header
|
---|
935 | for (vector<HeaderEntry>::iterator it=_defaultHeader.begin(); it != _defaultHeader.end(); it++)
|
---|
936 | _file.write(it->fitsString().c_str(), 80);
|
---|
937 |
|
---|
938 | for (vector<HeaderEntry>::iterator it=_header.begin(); it != _header.end(); it++)
|
---|
939 | _file.write(it->fitsString().c_str(), 80);
|
---|
940 |
|
---|
941 | _file.write("END ", 80);
|
---|
942 | long here = _file.tellp();
|
---|
943 | if (here%2880)
|
---|
944 | _file.write(_emptyBlock.c_str(), 2880 - here%2880);
|
---|
945 |
|
---|
946 | _headerFlushed = true;
|
---|
947 |
|
---|
948 | here = _file.tellp();
|
---|
949 |
|
---|
950 | if (here%2880)
|
---|
951 | cout << "Error: seems that header did not finish at the end of a block." << endl;
|
---|
952 |
|
---|
953 | if (here > cPos && cPos != 0)
|
---|
954 | {
|
---|
955 | cout << "Error, entries were added after the first row was written. This is not supposed to happen." << endl;
|
---|
956 | return;
|
---|
957 | }
|
---|
958 |
|
---|
959 | here = _file.tellp();
|
---|
960 | writeCatalog(closingFile);
|
---|
961 |
|
---|
962 | here = _file.tellp() - here;
|
---|
963 | _heapPtr = here;
|
---|
964 |
|
---|
965 | if (cPos != 0)
|
---|
966 | _file.seekp(cPos);
|
---|
967 | }
|
---|
968 |
|
---|
969 | /****************************************************************
|
---|
970 | * WRITE CATALOG
|
---|
971 | * WARNING: writeCatalog is only meant to be used by writeHeader.
|
---|
972 | * external usage will most likely corrupt the file
|
---|
973 | ****************************************************************/
|
---|
974 | void CompressedFitsWriter::writeCatalog(bool closingFile)
|
---|
975 | {
|
---|
976 | uint32_t sizeWritten = 0;
|
---|
977 | for (uint32_t i=0;i<_catalog.size();i++)
|
---|
978 | {
|
---|
979 | for (uint32_t j=0;j<_catalog[i].size();j++)
|
---|
980 | {
|
---|
981 | //swap the bytes
|
---|
982 | int8_t swappedEntry[16];
|
---|
983 | swappedEntry[0] = reinterpret_cast<int8_t*>(&_catalog[i][j].first)[7];
|
---|
984 | swappedEntry[1] = reinterpret_cast<int8_t*>(&_catalog[i][j].first)[6];
|
---|
985 | swappedEntry[2] = reinterpret_cast<int8_t*>(&_catalog[i][j].first)[5];
|
---|
986 | swappedEntry[3] = reinterpret_cast<int8_t*>(&_catalog[i][j].first)[4];
|
---|
987 | swappedEntry[4] = reinterpret_cast<int8_t*>(&_catalog[i][j].first)[3];
|
---|
988 | swappedEntry[5] = reinterpret_cast<int8_t*>(&_catalog[i][j].first)[2];
|
---|
989 | swappedEntry[6] = reinterpret_cast<int8_t*>(&_catalog[i][j].first)[1];
|
---|
990 | swappedEntry[7] = reinterpret_cast<int8_t*>(&_catalog[i][j].first)[0];
|
---|
991 |
|
---|
992 | swappedEntry[8] = reinterpret_cast<int8_t*>(&_catalog[i][j].second)[7];
|
---|
993 | swappedEntry[9] = reinterpret_cast<int8_t*>(&_catalog[i][j].second)[6];
|
---|
994 | swappedEntry[10] = reinterpret_cast<int8_t*>(&_catalog[i][j].second)[5];
|
---|
995 | swappedEntry[11] = reinterpret_cast<int8_t*>(&_catalog[i][j].second)[4];
|
---|
996 | swappedEntry[12] = reinterpret_cast<int8_t*>(&_catalog[i][j].second)[3];
|
---|
997 | swappedEntry[13] = reinterpret_cast<int8_t*>(&_catalog[i][j].second)[2];
|
---|
998 | swappedEntry[14] = reinterpret_cast<int8_t*>(&_catalog[i][j].second)[1];
|
---|
999 | swappedEntry[15] = reinterpret_cast<int8_t*>(&_catalog[i][j].second)[0];
|
---|
1000 | if (closingFile)
|
---|
1001 | {
|
---|
1002 | _checksum.add(reinterpret_cast<char*>(swappedEntry), 16);
|
---|
1003 | }
|
---|
1004 | _file.write(reinterpret_cast<char*>(&swappedEntry[0]), 2*sizeof(int64_t));
|
---|
1005 | sizeWritten += 2*sizeof(int64_t);
|
---|
1006 | }
|
---|
1007 | }
|
---|
1008 |
|
---|
1009 | //we do not reserve space for now because fverify does not like that.
|
---|
1010 | //TODO bug should be fixed in the new version. Install it on the cluster and restor space reservation
|
---|
1011 | return ;
|
---|
1012 |
|
---|
1013 | //write the padding so that the HEAP section starts at a 2880 bytes boundary
|
---|
1014 | if (sizeWritten % 2880 != 0)
|
---|
1015 | {
|
---|
1016 | vector<char> nullVec(2880 - sizeWritten%2880, 0);
|
---|
1017 | _file.write(nullVec.data(), 2880 - sizeWritten%2880);
|
---|
1018 | }
|
---|
1019 | }
|
---|
1020 |
|
---|
1021 | /****************************************************************
|
---|
1022 | * ADD HEADER CHECKSUM
|
---|
1023 | ****************************************************************/
|
---|
1024 | void CompressedFitsWriter::addHeaderChecksum(Checksum& checksum)
|
---|
1025 | {
|
---|
1026 | for (vector<HeaderEntry>::iterator it=_defaultHeader.begin();it!=_defaultHeader.end(); it++)
|
---|
1027 | _checksum.add(it->fitsString().c_str(), 80);
|
---|
1028 | for (vector<HeaderEntry>::iterator it=_header.begin(); it != _header.end(); it++)
|
---|
1029 | _checksum.add(it->fitsString().c_str(), 80);
|
---|
1030 | string end("END ");
|
---|
1031 | string space(" ");
|
---|
1032 | checksum.add(end.c_str(), 80);
|
---|
1033 | int headerRowsLeft = 36 - (_defaultHeader.size() + _header.size() + 1)%36;
|
---|
1034 | for (int i=0;i<headerRowsLeft;i++)
|
---|
1035 | checksum.add(space.c_str(), 80);
|
---|
1036 | }
|
---|
1037 |
|
---|
1038 | /****************************************************************
|
---|
1039 | * CLOSE
|
---|
1040 | ****************************************************************/
|
---|
1041 | bool CompressedFitsWriter::close()
|
---|
1042 | {
|
---|
1043 | for (uint32_t i=0;i<_numThreads;i++)
|
---|
1044 | while (_threadStatus[i] != _THREAD_WAIT_)
|
---|
1045 | usleep(100000);
|
---|
1046 | for (uint32_t i=0;i<_numThreads;i++)
|
---|
1047 | _threadStatus[i] = _THREAD_EXIT_;
|
---|
1048 | for (uint32_t i=0;i<_numThreads;i++)
|
---|
1049 | pthread_join(_thread[i], NULL);
|
---|
1050 | //flush the rows that were not written yet
|
---|
1051 | if (_totalNumRows%_numRowsPerTile != 0)
|
---|
1052 | {
|
---|
1053 | copyTransposeTile(0);
|
---|
1054 |
|
---|
1055 | _threadNumRows[0] = _totalNumRows;
|
---|
1056 | uint32_t numBytes = compressBuffer(0);
|
---|
1057 | writeCompressedDataToDisk(0, numBytes);
|
---|
1058 | }
|
---|
1059 | //compression stuff
|
---|
1060 | setHeaderKey(HeaderEntry("ZNAXIS1", _rowWidth, "Width of uncompressed rows"));
|
---|
1061 | setHeaderKey(HeaderEntry("ZNAXIS2", _totalNumRows, "Number of uncompressed rows"));
|
---|
1062 | //TODO calculate the real offset from the main table to the start of the HEAP data area
|
---|
1063 | setHeaderKey(HeaderEntry("ZHEAPPTR", _heapPtr, ""));
|
---|
1064 | setHeaderKey(HeaderEntry("THEAP", _heapPtr, ""));
|
---|
1065 |
|
---|
1066 | //regular stuff
|
---|
1067 | if (_catalog.size() > 0)
|
---|
1068 | {
|
---|
1069 | setHeaderKey(HeaderEntry("NAXIS1", 2*sizeof(int64_t)*_catalog[0].size(), "width of table in bytes"));
|
---|
1070 | setHeaderKey(HeaderEntry("NAXIS2", _numTiles, ""));
|
---|
1071 | setHeaderKey(HeaderEntry("TFIELDS", _columns.size(), "number of fields in each row"));
|
---|
1072 | int64_t heapSize = 0;
|
---|
1073 | uint32_t compressedOffset = 0;
|
---|
1074 | for (uint32_t i=0;i<_catalog.size();i++)
|
---|
1075 | for (uint32_t j=0;j<_catalog[i].size();j++)
|
---|
1076 | {
|
---|
1077 | heapSize += _catalog[i][j].first;
|
---|
1078 | //set the catalog offsets to their actual values
|
---|
1079 | _catalog[i][j].second = compressedOffset;
|
---|
1080 | compressedOffset += _catalog[i][j].first;
|
---|
1081 | //special case if entry has zero length
|
---|
1082 | if (_catalog[i][j].first == 0) _catalog[i][j].second = 0;
|
---|
1083 | }
|
---|
1084 | setHeaderKey(HeaderEntry("PCOUNT", heapSize, "size of special data area"));
|
---|
1085 | }
|
---|
1086 | writeHeader(true);
|
---|
1087 | ostringstream str;
|
---|
1088 | str << _checksum.val();
|
---|
1089 | setHeaderKey(HeaderEntry("DATASUM", str.str(), ""));
|
---|
1090 | addHeaderChecksum(_checksum);
|
---|
1091 | setHeaderKey(HeaderEntry("CHECKSUM", _checksum.str(), ""));
|
---|
1092 | //update header value
|
---|
1093 | writeHeader();
|
---|
1094 | //update file length
|
---|
1095 | long here = _file.tellp();
|
---|
1096 | if (here%2880)
|
---|
1097 | {
|
---|
1098 | vector<char> nullVec(2880 - here%2880, 0);
|
---|
1099 | _file.write(nullVec.data(), 2880 - here%2880);
|
---|
1100 | }
|
---|
1101 | _file.close();
|
---|
1102 | return true;
|
---|
1103 | }
|
---|
1104 |
|
---|
1105 | /****************************************************************
|
---|
1106 | * COPY TRANSPOSE TILE
|
---|
1107 | ****************************************************************/
|
---|
1108 | void CompressedFitsWriter::copyTransposeTile(uint32_t index)
|
---|
1109 | {
|
---|
1110 | uint32_t thisRoundNumRows = (_totalNumRows%_numRowsPerTile) ? _totalNumRows%_numRowsPerTile : _numRowsPerTile;
|
---|
1111 | //copy the tile and transpose it
|
---|
1112 | uint32_t offset = 0;
|
---|
1113 | for (uint32_t i=0;i<_columns.size();i++)
|
---|
1114 | {
|
---|
1115 | switch (_columns[i].getCompression())
|
---|
1116 | {
|
---|
1117 | case UNCOMPRESSED:
|
---|
1118 | case SMOOTHMAN:
|
---|
1119 | for (uint32_t k=0;k<thisRoundNumRows;k++)
|
---|
1120 | {//regular, "semi-transposed" copy
|
---|
1121 | memcpy(&(_transposedBuffer[index][offset]), &_buffer[k*_rowWidth + _columns[i].offset()], _columns[i].sizeOfElems()*_columns[i].numElems());
|
---|
1122 | offset += _columns[i].sizeOfElems()*_columns[i].numElems();
|
---|
1123 | }
|
---|
1124 | break;
|
---|
1125 | default :
|
---|
1126 | for (int j=0;j<_columns[i].numElems();j++)
|
---|
1127 | for (uint32_t k=0;k<thisRoundNumRows;k++)
|
---|
1128 | {//transposed copy
|
---|
1129 | memcpy(&(_transposedBuffer[index][offset]), &_buffer[k*_rowWidth + _columns[i].offset() + _columns[i].sizeOfElems()*j], _columns[i].sizeOfElems());
|
---|
1130 | offset += _columns[i].sizeOfElems();
|
---|
1131 | }
|
---|
1132 | break;
|
---|
1133 |
|
---|
1134 | };
|
---|
1135 | }
|
---|
1136 | }
|
---|
1137 |
|
---|
1138 | /****************************************************************
|
---|
1139 | * WRITE BINARY ROW
|
---|
1140 | ****************************************************************/
|
---|
1141 | bool CompressedFitsWriter::writeBinaryRow(const char* bufferToWrite)
|
---|
1142 | {
|
---|
1143 | if (_totalNumRows == 0)
|
---|
1144 | writeHeader();
|
---|
1145 |
|
---|
1146 | memcpy(&_buffer[_rowWidth*(_totalNumRows%_numRowsPerTile)], bufferToWrite, _rowWidth);
|
---|
1147 | _totalNumRows++;
|
---|
1148 | if (_totalNumRows%_numRowsPerTile == 0)
|
---|
1149 | {
|
---|
1150 | //which is the next thread that we should use ?
|
---|
1151 | while (_threadStatus[_threadLooper] == _THREAD_COMPRESS_)
|
---|
1152 | usleep(100000);
|
---|
1153 |
|
---|
1154 | copyTransposeTile(_threadLooper);
|
---|
1155 |
|
---|
1156 | while (_threadStatus[_threadLooper] != _THREAD_WAIT_)
|
---|
1157 | usleep(100000);
|
---|
1158 | _threadNumRows[_threadLooper] = _totalNumRows;
|
---|
1159 | _threadStatus[_threadLooper] = _THREAD_COMPRESS_;
|
---|
1160 | _threadLooper = (_threadLooper+1)%_numThreads;
|
---|
1161 | }
|
---|
1162 | return _file.good();
|
---|
1163 | }
|
---|
1164 |
|
---|
1165 | /****************************************************************
|
---|
1166 | * COMPRESS BUFFER
|
---|
1167 | ****************************************************************/
|
---|
1168 | uint32_t CompressedFitsWriter::compressUNCOMPRESSED(char* dest, const char* src, uint32_t numRows, uint32_t sizeOfElems, uint32_t numRowElems)
|
---|
1169 | {
|
---|
1170 | memcpy(dest, src, numRows*sizeOfElems*numRowElems);
|
---|
1171 | return numRows*sizeOfElems*numRowElems;
|
---|
1172 | }
|
---|
1173 |
|
---|
1174 | /****************************************************************
|
---|
1175 | * COMPRESS BUFFER
|
---|
1176 | ****************************************************************/
|
---|
1177 | uint32_t CompressedFitsWriter::compressHUFFMAN(char* dest, const char* src, uint32_t numRows, uint32_t sizeOfElems, uint32_t numRowElems)
|
---|
1178 | {
|
---|
1179 | string huffmanOutput;
|
---|
1180 | uint32_t previousHuffmanSize = 0;
|
---|
1181 | if (numRows < 2)
|
---|
1182 | {//if we have less than 2 elems to compress, Huffman encoder does not work (and has no point). Just return larger size than uncompressed to trigger the raw storage.
|
---|
1183 | return numRows*sizeOfElems*numRowElems + 1000;
|
---|
1184 | }
|
---|
1185 | if (sizeOfElems < 2 )
|
---|
1186 | {
|
---|
1187 | cout << "Fatal ERROR: HUFMANN can only encode short or longer types" << endl;
|
---|
1188 | return 0;
|
---|
1189 | }
|
---|
1190 | uint32_t huffmanOffset = 0;
|
---|
1191 | for (uint32_t j=0;j<numRowElems;j++)
|
---|
1192 | {
|
---|
1193 | Huffman::Encode(huffmanOutput,
|
---|
1194 | reinterpret_cast<const uint16_t*>(&src[j*sizeOfElems*numRows]),
|
---|
1195 | numRows*(sizeOfElems/2));
|
---|
1196 | reinterpret_cast<uint32_t*>(&dest[huffmanOffset])[0] = huffmanOutput.size() - previousHuffmanSize;
|
---|
1197 | huffmanOffset += sizeof(uint32_t);
|
---|
1198 | previousHuffmanSize = huffmanOutput.size();
|
---|
1199 | }
|
---|
1200 | memcpy(&dest[huffmanOffset], huffmanOutput.data(), huffmanOutput.size());
|
---|
1201 | return huffmanOutput.size() + huffmanOffset;
|
---|
1202 | }
|
---|
1203 |
|
---|
1204 | /****************************************************************
|
---|
1205 | * COMPRESS BUFFER
|
---|
1206 | ****************************************************************/
|
---|
1207 | uint32_t CompressedFitsWriter::compressSMOOTHMAN(char* dest, char* src, uint32_t numRows, uint32_t sizeOfElems, uint32_t numRowElems)
|
---|
1208 | {
|
---|
1209 | uint32_t colWidth = numRowElems;
|
---|
1210 | for (int j=colWidth*numRows-1;j>1;j--)
|
---|
1211 | reinterpret_cast<int16_t*>(src)[j] = reinterpret_cast<int16_t*>(src)[j] - (reinterpret_cast<int16_t*>(src)[j-1]+reinterpret_cast<int16_t*>(src)[j-2])/2;
|
---|
1212 | //call the huffman transposed
|
---|
1213 | return compressHUFFMAN(dest, src, numRowElems, sizeOfElems, numRows);
|
---|
1214 | }
|
---|
1215 |
|
---|
1216 | /****************************************************************
|
---|
1217 | * COMPRESS BUFFER
|
---|
1218 | ****************************************************************/
|
---|
1219 | uint32_t CompressedFitsWriter::compressBuffer(uint32_t threadIndex)
|
---|
1220 | {
|
---|
1221 | uint32_t thisRoundNumRows = (_threadNumRows[threadIndex]%_numRowsPerTile) ? _threadNumRows[threadIndex]%_numRowsPerTile : _numRowsPerTile;
|
---|
1222 | uint32_t offset=0;
|
---|
1223 | uint32_t currentCatalogRow = (_threadNumRows[threadIndex]-1)/_numRowsPerTile;
|
---|
1224 | int64_t compressedOffset = 0;
|
---|
1225 |
|
---|
1226 | //now compress each column one by one by calling compression on arrays
|
---|
1227 | for (uint32_t i=0;i<_columns.size();i++)
|
---|
1228 | {
|
---|
1229 | _catalog[currentCatalogRow][i].second = compressedOffset;
|
---|
1230 |
|
---|
1231 | uint32_t compression = _columns[i].getCompression();
|
---|
1232 |
|
---|
1233 | if (_columns[i].numElems() == 0) continue;
|
---|
1234 | //set the default byte telling if uncompressed the compressed Flag
|
---|
1235 | int64_t previousOffset = compressedOffset;
|
---|
1236 | _compressedBuffer[threadIndex][compressedOffset++] = COMPRESSED_FLAG;
|
---|
1237 | switch (compression)
|
---|
1238 | {
|
---|
1239 | case UNCOMPRESSED:
|
---|
1240 | compressedOffset += compressUNCOMPRESSED(&(_compressedBuffer[threadIndex][compressedOffset]), &(_transposedBuffer[threadIndex][offset]), thisRoundNumRows, _columns[i].sizeOfElems(), _columns[i].numElems());
|
---|
1241 | break;
|
---|
1242 | case SMOOTHMAN:
|
---|
1243 | compressedOffset += compressSMOOTHMAN(&(_compressedBuffer[threadIndex][compressedOffset]), &(_transposedBuffer[threadIndex][offset]), thisRoundNumRows, _columns[i].sizeOfElems(), _columns[i].numElems());
|
---|
1244 | break;
|
---|
1245 |
|
---|
1246 | default:
|
---|
1247 | ;
|
---|
1248 | };
|
---|
1249 | //check if compressed size is larger than uncompressed
|
---|
1250 | if (compression != UNCOMPRESSED &&
|
---|
1251 | compressedOffset - previousOffset > _columns[i].sizeOfElems()*_columns[i].numElems()*thisRoundNumRows+1)
|
---|
1252 | {//if so set flag and redo it uncompressed
|
---|
1253 | compressedOffset = previousOffset;
|
---|
1254 | _compressedBuffer[threadIndex][compressedOffset++] = UNCOMPRESSED_FLAG;
|
---|
1255 | compressedOffset += compressUNCOMPRESSED(&(_compressedBuffer[threadIndex][compressedOffset]), &(_transposedBuffer[threadIndex][offset]), thisRoundNumRows, _columns[i].sizeOfElems(), _columns[i].numElems());
|
---|
1256 | }
|
---|
1257 | offset += thisRoundNumRows*_columns[i].sizeOfElems()*_columns[i].numElems();
|
---|
1258 | _catalog[currentCatalogRow][i].first = compressedOffset - _catalog[currentCatalogRow][i].second;
|
---|
1259 | }
|
---|
1260 | return compressedOffset;
|
---|
1261 | }
|
---|
1262 |
|
---|
1263 | /****************************************************************
|
---|
1264 | * WRITE COMPRESS DATA TO DISK
|
---|
1265 | ****************************************************************/
|
---|
1266 | bool CompressedFitsWriter::writeCompressedDataToDisk(uint32_t threadID, uint32_t sizeToWrite)
|
---|
1267 | {
|
---|
1268 | char* checkSumPointer = _compressedBuffer[threadID];
|
---|
1269 | int32_t extraBytes = 0;
|
---|
1270 | uint32_t sizeToChecksum = sizeToWrite;
|
---|
1271 | if (_checkOffset != 0)
|
---|
1272 | {//should we extend the array to the left ?
|
---|
1273 | sizeToChecksum += _checkOffset;
|
---|
1274 | checkSumPointer -= _checkOffset;
|
---|
1275 | memset(checkSumPointer, 0, _checkOffset);
|
---|
1276 | }
|
---|
1277 | if (sizeToChecksum%4 != 0)
|
---|
1278 | {//should we extend the array to the right ?
|
---|
1279 | extraBytes = 4 - (sizeToChecksum%4);
|
---|
1280 | memset(checkSumPointer+sizeToChecksum, 0,extraBytes);
|
---|
1281 | sizeToChecksum += extraBytes;
|
---|
1282 | }
|
---|
1283 | //do the checksum
|
---|
1284 | _checksum.add(checkSumPointer, sizeToChecksum);
|
---|
1285 | _checkOffset = (4 - extraBytes)%4;
|
---|
1286 | //write data to disk
|
---|
1287 | _file.write(_compressedBuffer[threadID], sizeToWrite);
|
---|
1288 | return _file.good();
|
---|
1289 | }
|
---|
1290 |
|
---|
1291 | /****************************************************************
|
---|
1292 | * WRITER THREAD LOOP
|
---|
1293 | ****************************************************************/
|
---|
1294 | void* CompressedFitsWriter::threadFunction(void* context)
|
---|
1295 | {
|
---|
1296 | CompressedFitsWriter* myself =static_cast<CompressedFitsWriter*>(context);
|
---|
1297 |
|
---|
1298 | uint32_t myID = 0;
|
---|
1299 | pthread_mutex_lock(&(myself->_mutex));
|
---|
1300 | myID = myself->_threadIndex++;
|
---|
1301 | pthread_mutex_unlock(&(myself->_mutex));
|
---|
1302 | uint32_t threadToWaitForBeforeWriting = (myID == 0) ? myself->_numThreads-1 : myID-1;
|
---|
1303 | while (myself->_threadStatus[myID] != _THREAD_EXIT_)
|
---|
1304 | {
|
---|
1305 | while (myself->_threadStatus[myID] == _THREAD_WAIT_)
|
---|
1306 | usleep(100000);
|
---|
1307 | if (myself->_threadStatus[myID] != _THREAD_COMPRESS_)
|
---|
1308 | continue;
|
---|
1309 | uint32_t numBytes = myself->compressBuffer(myID);
|
---|
1310 | myself->_threadStatus[myID] = _THREAD_WRITE_;
|
---|
1311 |
|
---|
1312 | //wait for the previous data to be written
|
---|
1313 | while (myself->_threadIndex != threadToWaitForBeforeWriting)
|
---|
1314 | usleep(1000);
|
---|
1315 | //do the actual writing to disk
|
---|
1316 | pthread_mutex_lock(&(myself->_mutex));
|
---|
1317 | myself->writeCompressedDataToDisk(myID, numBytes);
|
---|
1318 | myself->_threadIndex = myID;
|
---|
1319 | pthread_mutex_unlock(&(myself->_mutex));
|
---|
1320 | myself->_threadStatus[myID] = _THREAD_WAIT_;
|
---|
1321 | }
|
---|
1322 | return NULL;
|
---|
1323 | }
|
---|
1324 |
|
---|
1325 | /****************************************************************
|
---|
1326 | * PRINT USAGE
|
---|
1327 | ****************************************************************/
|
---|
1328 | void printUsage()
|
---|
1329 | {
|
---|
1330 | cout << endl;
|
---|
1331 | cout << "The FACT-Fits compressor reads an input Fits file from FACT"
|
---|
1332 | " and compresses it.\n It can use a drs calibration in order to"
|
---|
1333 | " improve the compression ratio. If so, the input calibration"
|
---|
1334 | " is embedded into the compressed file.\n"
|
---|
1335 | " By default, the Data column will be compressed using SMOOTHMAN (Thomas' algorithm)"
|
---|
1336 | " while other columns will be compressed with the AMPLITUDE coding (Veritas)"
|
---|
1337 | "Usage: Compressed_Fits_Test <inputFile>";
|
---|
1338 | cout << endl;
|
---|
1339 | }
|
---|
1340 |
|
---|
1341 | /****************************************************************
|
---|
1342 | * PRINT HELP
|
---|
1343 | ****************************************************************/
|
---|
1344 | void printHelp()
|
---|
1345 | {
|
---|
1346 | cout << endl;
|
---|
1347 | cout << "The inputFile is required. It must have fits in its filename and the compressed file will be written in the same folder. "
|
---|
1348 | "The fz extension will be added, replacing the .gz one if required \n"
|
---|
1349 | "If output is specified, then it will replace the automatically generated output filename\n"
|
---|
1350 | "If --drs, followed by a drs calib then it will be applied to the data before compressing\n"
|
---|
1351 | "rowPerTile can be used to tune how many rows are in each tile. Default is 100\n"
|
---|
1352 | "threads gives the number of threads to use. Cannot be less than the default (1)\n"
|
---|
1353 | "compression explicitely gives the compression scheme to use for a given column. The syntax is:\n"
|
---|
1354 | "<ColumnName>=<CompressionScheme> with <CompressionScheme> one of the following:\n"
|
---|
1355 | "UNCOMPRESSED\n"
|
---|
1356 | "AMPLITUDE\n"
|
---|
1357 | "HUFFMAN\n"
|
---|
1358 | "SMOOTHMAN\n"
|
---|
1359 | "INT_WAVELET\n"
|
---|
1360 | "\n"
|
---|
1361 | "--quiet removes any textual output, except error messages\n"
|
---|
1362 | "--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."
|
---|
1363 | ;
|
---|
1364 | cout << endl << endl;
|
---|
1365 | }
|
---|
1366 |
|
---|
1367 | /****************************************************************
|
---|
1368 | * SETUP CONFIGURATION
|
---|
1369 | ****************************************************************/
|
---|
1370 | void setupConfiguration(Configuration& conf)
|
---|
1371 | {
|
---|
1372 | po::options_description configs("FitsCompressor options");
|
---|
1373 | configs.add_options()
|
---|
1374 | ("inputFile,i", vars<string>(), "Input file")
|
---|
1375 | ("drs,d", var<string>(), "Input drs calibration file")
|
---|
1376 | ("rowPerTile,r", var<unsigned int>(), "Number of rows per tile. Default is 100")
|
---|
1377 | ("output,o", var<string>(), "Output file. If empty, .fz is appened to the original name")
|
---|
1378 | ("threads,t", var<unsigned int>(), "Number of threads to use for compression")
|
---|
1379 | ("compression,c", vars<string>(), "which compression to use for which column. Syntax <colName>=<compressionScheme>")
|
---|
1380 | ("quiet,q", po_switch(), "Should the program display any text at all ?")
|
---|
1381 | ("verify,v", po_switch(), "Should we verify the data that has been compressed ?")
|
---|
1382 | ;
|
---|
1383 | po::positional_options_description positional;
|
---|
1384 | positional.add("inputFile", -1);
|
---|
1385 | conf.AddOptions(configs);
|
---|
1386 | conf.SetArgumentPositions(positional);
|
---|
1387 | }
|
---|
1388 |
|
---|
1389 | /****************************************************************
|
---|
1390 | * MAIN
|
---|
1391 | ****************************************************************/
|
---|
1392 | int main(int argc, const char** argv)
|
---|
1393 | {
|
---|
1394 | Configuration conf(argv[0]);
|
---|
1395 | conf.SetPrintUsage(printUsage);
|
---|
1396 | setupConfiguration(conf);
|
---|
1397 |
|
---|
1398 | if (!conf.DoParse(argc, argv, printHelp))
|
---|
1399 | return -1;
|
---|
1400 |
|
---|
1401 | //initialize the file names to nothing.
|
---|
1402 | string fileNameIn = "";
|
---|
1403 | string fileNameOut = "";
|
---|
1404 | string drsFileName = "";
|
---|
1405 | uint32_t numRowsPerTile = 100;
|
---|
1406 | bool displayText=true;
|
---|
1407 |
|
---|
1408 | //parse configuration
|
---|
1409 | if (conf.Get<bool>("quiet")) displayText = false;
|
---|
1410 | const vector<string> inputFileNameVec = conf.Vec<string>("inputFile");
|
---|
1411 | if (inputFileNameVec.size() != 1)
|
---|
1412 | {
|
---|
1413 | cout << "Error: ";
|
---|
1414 | if (inputFileNameVec.size() == 0) cout << "no";
|
---|
1415 | else cout << inputFileNameVec.size();
|
---|
1416 | cout << " input file(s) given. Expected one. Aborting. Input:" << endl;;
|
---|
1417 | for (unsigned int i=0;i<inputFileNameVec.size(); i++)
|
---|
1418 | cout << inputFileNameVec[i] << endl;
|
---|
1419 | return -1;
|
---|
1420 | }
|
---|
1421 |
|
---|
1422 | //Assign the input filename
|
---|
1423 | fileNameIn = inputFileNameVec[0];
|
---|
1424 |
|
---|
1425 | //Check if we have a drs calib too
|
---|
1426 | if (conf.Has("drs")) drsFileName = conf.Get<string>("drs");
|
---|
1427 |
|
---|
1428 | //Should we verify the data ?
|
---|
1429 | bool verifyDataPlease = false;
|
---|
1430 | if (conf.Has("verify")) verifyDataPlease = conf.Get<bool>("verify");
|
---|
1431 |
|
---|
1432 |
|
---|
1433 | //should we use a specific output filename ?
|
---|
1434 | if (conf.Has("output"))
|
---|
1435 | fileNameOut = conf.Get<string>("output");
|
---|
1436 | else
|
---|
1437 | {
|
---|
1438 | size_t pos = fileNameIn.find(".fits");
|
---|
1439 | if (pos == string::npos)
|
---|
1440 | {
|
---|
1441 | cout << "ERROR: input file does not seems ot be fits. Aborting." << endl;
|
---|
1442 | return -1;
|
---|
1443 | }
|
---|
1444 | fileNameOut = fileNameIn.substr(0, pos) + ".fz";
|
---|
1445 | }
|
---|
1446 |
|
---|
1447 | //should we use specific compression on some columns ?
|
---|
1448 | const vector<string> columnsCompression = conf.Vec<string>("compression");
|
---|
1449 |
|
---|
1450 | //split up values between column names and compression scheme
|
---|
1451 | vector<std::pair<string, string>> compressions;
|
---|
1452 | for (unsigned int i=0;i<columnsCompression.size();i++)
|
---|
1453 | {
|
---|
1454 | size_t pos = columnsCompression[i].find_first_of("=");
|
---|
1455 | if (pos == string::npos)
|
---|
1456 | {
|
---|
1457 | cout << "ERROR: Something wrong occured while parsing " << columnsCompression[i] << ". Aborting." << endl;
|
---|
1458 | return -1;
|
---|
1459 | }
|
---|
1460 | string comp = columnsCompression[i].substr(pos+1);
|
---|
1461 | if (comp != "UNCOMPRESSED" && comp != "AMPLITUDE" && comp != "HUFFMAN" &&
|
---|
1462 | comp != "SMOOTHMAN" && comp != "INT_WAVELET")
|
---|
1463 | {
|
---|
1464 | cout << "Unkown compression scheme requested (" << comp << "). Aborting." << endl;
|
---|
1465 | return -1;
|
---|
1466 | }
|
---|
1467 | compressions.push_back(make_pair(columnsCompression[i].substr(0, pos), comp));
|
---|
1468 | }
|
---|
1469 |
|
---|
1470 | //How many rows per tile should we use ?
|
---|
1471 | if (conf.Has("rowPerTile")) numRowsPerTile = conf.Get<unsigned int>("rowPerTile");
|
---|
1472 |
|
---|
1473 | /************************************************************************************
|
---|
1474 | * Done reading configuration. Open relevant files
|
---|
1475 | ************************************************************************************/
|
---|
1476 |
|
---|
1477 | //Open input's fits file
|
---|
1478 | factFits inFile(fileNameIn);
|
---|
1479 |
|
---|
1480 | if (inFile.IsCompressedFITS())
|
---|
1481 | {
|
---|
1482 | cout << "ERROR: input file is already a compressed fits. Cannot be compressed again: Aborting." << endl;
|
---|
1483 | return -1;
|
---|
1484 | }
|
---|
1485 |
|
---|
1486 | //decide how many tiles should be put in the compressed file
|
---|
1487 | uint32_t originalNumRows = inFile.GetNumRows();
|
---|
1488 | uint32_t numTiles = (originalNumRows%numRowsPerTile) ? (originalNumRows/numRowsPerTile)+1 : originalNumRows/numRowsPerTile;
|
---|
1489 | CompressedFitsWriter outFile(numTiles, numRowsPerTile);
|
---|
1490 |
|
---|
1491 | //should we use a specific number of threads for compressing ?
|
---|
1492 | unsigned int numThreads = 1;
|
---|
1493 | if (conf.Has("threads"))
|
---|
1494 | {
|
---|
1495 | numThreads = conf.Get<unsigned int>("threads");
|
---|
1496 | outFile.setNumWorkingThreads(numThreads);
|
---|
1497 | }
|
---|
1498 |
|
---|
1499 | if (!outFile.open(fileNameOut))
|
---|
1500 | {
|
---|
1501 | cout << "Error: could not open " << fileNameOut << " for writing" << endl;
|
---|
1502 | return -1;
|
---|
1503 | }
|
---|
1504 |
|
---|
1505 | //Because the file to open MUST be given by the constructor, I must use a pointer instead
|
---|
1506 | factFits* drsFile = NULL;
|
---|
1507 | //try to open the Drs file. If any.
|
---|
1508 | if (drsFileName != "")
|
---|
1509 | {
|
---|
1510 | try
|
---|
1511 | {
|
---|
1512 | drsFile = new factFits(drsFileName);
|
---|
1513 | }
|
---|
1514 | catch (...)
|
---|
1515 | {
|
---|
1516 | cout << "Error: could not open " << drsFileName << " for calibration" << endl;
|
---|
1517 | return -1;
|
---|
1518 | }
|
---|
1519 | }
|
---|
1520 |
|
---|
1521 | if (displayText)
|
---|
1522 | {
|
---|
1523 | cout << endl;
|
---|
1524 | cout << "**********************" << endl;
|
---|
1525 | cout << "Will compress from : " << fileNameIn << endl;
|
---|
1526 | cout << "to : " << fileNameOut << endl;
|
---|
1527 | if (drsFileName != "")
|
---|
1528 | cout << "while calibrating with: " << drsFileName << endl;
|
---|
1529 | cout << "Compression will use : " << numThreads << " worker threads" << endl;
|
---|
1530 | cout << "Data will be verified : ";
|
---|
1531 | if (verifyDataPlease)
|
---|
1532 | cout << "yes" << endl;
|
---|
1533 | else
|
---|
1534 | cout << "no (WARNING !)" << endl;
|
---|
1535 | cout << "**********************" << endl;
|
---|
1536 | cout << endl;
|
---|
1537 | }
|
---|
1538 |
|
---|
1539 | /************************************************************************************
|
---|
1540 | * Done opening input files. Allocate memory and configure output file
|
---|
1541 | ************************************************************************************/
|
---|
1542 |
|
---|
1543 | //allocate the buffer for temporary storage of each read/written row
|
---|
1544 | uint32_t rowWidth = inFile.GetUInt("NAXIS1");
|
---|
1545 | char* buffer = new char[rowWidth];
|
---|
1546 |
|
---|
1547 | //get the source columns
|
---|
1548 | const fits::Table::Columns& columns = inFile.GetColumns();
|
---|
1549 | const fits::Table::SortedColumns& sortedColumns = inFile.GetSortedColumns();
|
---|
1550 | if (displayText)
|
---|
1551 | cout << "Input file has " << columns.size() << " columns and " << inFile.GetNumRows() << " rows" << endl;
|
---|
1552 |
|
---|
1553 | //Add columns.
|
---|
1554 | for (uint32_t i=0;i<sortedColumns.size(); i++)
|
---|
1555 | {
|
---|
1556 | //get column name
|
---|
1557 | ostringstream str;
|
---|
1558 | str << "TTYPE" << i+1;
|
---|
1559 | string colName = inFile.GetStr(str.str());
|
---|
1560 | if (displayText)
|
---|
1561 | {
|
---|
1562 | cout << "Column " << colName;
|
---|
1563 | for (uint32_t j=colName.size();j<21;j++)
|
---|
1564 | cout << " ";
|
---|
1565 | cout << " -> ";
|
---|
1566 | }
|
---|
1567 |
|
---|
1568 | //first lets see if we do have an explicit request
|
---|
1569 | bool explicitRequest = false;
|
---|
1570 | for (unsigned int j=0;j<compressions.size();j++)
|
---|
1571 | {
|
---|
1572 | if (compressions[j].first == colName)
|
---|
1573 | {
|
---|
1574 | explicitRequest = true;
|
---|
1575 | if (displayText) cout << compressions[j].second << endl;
|
---|
1576 | if (compressions[j].second == "UNCOMPRESSED")
|
---|
1577 | outFile.addColumn(CompressedFitsFile::ColumnEntry(colName, sortedColumns[i].type, sortedColumns[i].num, CompressedFitsFile::UNCOMPRESSED));
|
---|
1578 | if (compressions[j].second == "SMOOTHMAN")
|
---|
1579 | outFile.addColumn(CompressedFitsFile::ColumnEntry(colName, sortedColumns[i].type, sortedColumns[i].num, CompressedFitsFile::SMOOTHMAN));
|
---|
1580 | break;
|
---|
1581 | }
|
---|
1582 | }
|
---|
1583 |
|
---|
1584 | if (explicitRequest) continue;
|
---|
1585 |
|
---|
1586 | if (colName != "Data")
|
---|
1587 | {
|
---|
1588 | if (displayText) cout << "UNCOMPRESSED" << endl;
|
---|
1589 | outFile.addColumn(CompressedFitsFile::ColumnEntry(colName, sortedColumns[i].type, sortedColumns[i].num, CompressedFitsFile::UNCOMPRESSED));
|
---|
1590 | }
|
---|
1591 | else
|
---|
1592 | {
|
---|
1593 | if (displayText) cout << "SMOOTHMAN" << endl;
|
---|
1594 | outFile.addColumn(CompressedFitsFile::ColumnEntry(colName, sortedColumns[i].type, sortedColumns[i].num, CompressedFitsFile::SMOOTHMAN));
|
---|
1595 | }
|
---|
1596 | }
|
---|
1597 |
|
---|
1598 | //translate original header entries to their Z-version
|
---|
1599 | const fits::Table::Keys& header = inFile.GetKeys();
|
---|
1600 | for (fits::Table::Keys::const_iterator i=header.begin(); i!= header.end(); i++)
|
---|
1601 | {
|
---|
1602 | string k = i->first;//header[i].key();
|
---|
1603 | if (k == "XTENSION" || k == "BITPIX" || k == "PCOUNT" || k == "GCOUNT" || k == "TFIELDS")
|
---|
1604 | continue;
|
---|
1605 | if (k == "CHECKSUM")
|
---|
1606 | {
|
---|
1607 | outFile.setHeaderKey(CompressedFitsFile::HeaderEntry("ZCHKSUM", i->second.value, i->second.comment));
|
---|
1608 | continue;
|
---|
1609 | }
|
---|
1610 | if (k == "DATASUM")
|
---|
1611 | {
|
---|
1612 | outFile.setHeaderKey(CompressedFitsFile::HeaderEntry("ZDTASUM", i->second.value, i->second.comment));
|
---|
1613 | continue;
|
---|
1614 | }
|
---|
1615 | k = k.substr(0,5);
|
---|
1616 | if (k == "TTYPE" || k == "TFORM")
|
---|
1617 | {
|
---|
1618 | string tmpKey = i->second.fitsString;
|
---|
1619 | tmpKey[0] = 'Z';
|
---|
1620 | outFile.setHeaderKey(tmpKey);
|
---|
1621 | continue;
|
---|
1622 | }
|
---|
1623 | if (k == "NAXIS")
|
---|
1624 | continue;
|
---|
1625 | outFile.setHeaderKey(i->second.fitsString);
|
---|
1626 | }
|
---|
1627 |
|
---|
1628 | //deal with the DRS calib
|
---|
1629 | int16_t* drsCalib16 = NULL;
|
---|
1630 |
|
---|
1631 | //load the drs calib. data
|
---|
1632 | int32_t startCellOffset = -1;
|
---|
1633 | if (drsFileName != "")
|
---|
1634 | {
|
---|
1635 | drsCalib16 = new int16_t[1440*1024];
|
---|
1636 | float* drsCalibFloat = NULL;
|
---|
1637 | try
|
---|
1638 | {
|
---|
1639 | drsCalibFloat = reinterpret_cast<float*>(drsFile->SetPtrAddress("BaselineMean"));
|
---|
1640 | }
|
---|
1641 | catch (...)
|
---|
1642 | {
|
---|
1643 | cout << "Could not find column BaselineMean in drs calibration file " << drsFileName << ". Aborting" << endl;
|
---|
1644 | return -1;
|
---|
1645 | }
|
---|
1646 |
|
---|
1647 | //read the calibration and calculate its integer value
|
---|
1648 | drsFile->GetNextRow();
|
---|
1649 | for (uint32_t i=0;i<1440*1024;i++)
|
---|
1650 | drsCalib16[i] = (int16_t)(drsCalibFloat[i]*4096.f/2000.f);
|
---|
1651 |
|
---|
1652 | //assign it to the ouput file
|
---|
1653 | outFile.setDrsCalib(drsCalib16);
|
---|
1654 |
|
---|
1655 | //get the start cells offsets
|
---|
1656 | for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end(); it++)
|
---|
1657 | if (it->first == "StartCellData")
|
---|
1658 | {
|
---|
1659 | startCellOffset = it->second.offset;
|
---|
1660 | if (it->second.type != 'I')
|
---|
1661 | {
|
---|
1662 | cout << "Wrong type for the StartCellData Column: " << it->second.type << " instead of I expected"<< endl;
|
---|
1663 | return -1;
|
---|
1664 | }
|
---|
1665 | }
|
---|
1666 |
|
---|
1667 | if (startCellOffset == -1)
|
---|
1668 | {
|
---|
1669 | cout << "Could not find StartCellData in input file " << fileNameIn << ". Aborting."<< endl;
|
---|
1670 | return -1;
|
---|
1671 | }
|
---|
1672 | }
|
---|
1673 |
|
---|
1674 | /************************************************************************************
|
---|
1675 | * Done configuring compression. Do the real job now !
|
---|
1676 | ************************************************************************************/
|
---|
1677 |
|
---|
1678 | if (displayText) cout << "Converting file..." << endl;
|
---|
1679 |
|
---|
1680 | int numSlices = -1;
|
---|
1681 | int32_t dataOffset = -1;
|
---|
1682 |
|
---|
1683 | //Get the pointer to the column that must be drs-calibrated
|
---|
1684 | for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end(); it++)
|
---|
1685 | if (it->first == "Data")
|
---|
1686 | {
|
---|
1687 | numSlices = it->second.num;
|
---|
1688 | dataOffset = it->second.offset;
|
---|
1689 | }
|
---|
1690 | if (numSlices % 1440 != 0)
|
---|
1691 | {
|
---|
1692 | cout << "seems like the number of samples is not a multiple of 1440. Aborting." << endl;
|
---|
1693 | return -1;
|
---|
1694 | }
|
---|
1695 | if (dataOffset == -1)
|
---|
1696 | {
|
---|
1697 | cout << "Could not find the column Data in the input file. Aborting." << endl;
|
---|
1698 | return -1;
|
---|
1699 | }
|
---|
1700 |
|
---|
1701 | numSlices /= 1440;
|
---|
1702 |
|
---|
1703 | //set pointers to the readout data to later be able to gather it to "buffer".
|
---|
1704 | vector<void*> readPointers;
|
---|
1705 | vector<int32_t> readOffsets;
|
---|
1706 | vector<int32_t> readElemSize;
|
---|
1707 | vector<int32_t> readNumElems;
|
---|
1708 | for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end(); it++)
|
---|
1709 | {
|
---|
1710 | readPointers.push_back(inFile.SetPtrAddress(it->first));
|
---|
1711 | readOffsets.push_back(it->second.offset);
|
---|
1712 | readElemSize.push_back(it->second.size);
|
---|
1713 | readNumElems.push_back(it->second.num);
|
---|
1714 | }
|
---|
1715 |
|
---|
1716 | //Convert each row one after the other
|
---|
1717 | for (uint32_t i=0;i<inFile.GetNumRows();i++)
|
---|
1718 | {
|
---|
1719 | if (displayText) cout << "\r Row " << i+1 << flush;
|
---|
1720 | inFile.GetNextRow();
|
---|
1721 | //copy from inFile internal structures to buffer
|
---|
1722 | int32_t count=0;
|
---|
1723 | for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end();it++)
|
---|
1724 | {
|
---|
1725 | memcpy(&buffer[readOffsets[count]], readPointers[count], readElemSize[count]*readNumElems[count]);
|
---|
1726 | count++;
|
---|
1727 | }
|
---|
1728 |
|
---|
1729 | if (startCellOffset != -1)
|
---|
1730 | {//drs calibrate this data
|
---|
1731 | for (int j=0;j<1440;j++)
|
---|
1732 | {
|
---|
1733 | int thisStartCell = reinterpret_cast<int16_t*>(&buffer[startCellOffset])[j];
|
---|
1734 | for (int k=0;k<numSlices;k++)
|
---|
1735 | reinterpret_cast<int16_t*>(&buffer[dataOffset])[numSlices*j + k] -= drsCalib16[1024*j + (thisStartCell+k)%1024];
|
---|
1736 | }
|
---|
1737 | }
|
---|
1738 | outFile.writeBinaryRow(buffer);
|
---|
1739 | };
|
---|
1740 |
|
---|
1741 | //Get table name for later use in case the compressed file is to be verified
|
---|
1742 | string tableName = inFile.GetStr("EXTNAME");
|
---|
1743 |
|
---|
1744 | if (displayText) cout << endl << "Done. Flushing output file..." << endl;
|
---|
1745 |
|
---|
1746 | inFile.close();
|
---|
1747 | outFile.close();
|
---|
1748 | delete[] drsCalib16;
|
---|
1749 |
|
---|
1750 | if (displayText) cout << "Done." << endl;
|
---|
1751 |
|
---|
1752 | /************************************************************************************
|
---|
1753 | * Actual job done. Should we verify what we did ?
|
---|
1754 | ************************************************************************************/
|
---|
1755 |
|
---|
1756 | if (verifyDataPlease)
|
---|
1757 | {
|
---|
1758 | if (displayText) cout << "Now verify data..." << endl;
|
---|
1759 | }
|
---|
1760 | else
|
---|
1761 | return 0;
|
---|
1762 |
|
---|
1763 | //get a compressed reader
|
---|
1764 | factFits verifyFile(fileNameOut, tableName, false);
|
---|
1765 |
|
---|
1766 | //and the header of the compressed file
|
---|
1767 | const fits::Table::Keys& header2 = verifyFile.GetKeys();
|
---|
1768 |
|
---|
1769 | //get a non-compressed writer
|
---|
1770 | ofits reconstructedFile;
|
---|
1771 |
|
---|
1772 | //figure out its name: /dev/null unless otherwise specified
|
---|
1773 | string reconstructedName = fileNameOut+".recons";
|
---|
1774 | reconstructedName = "/dev/null";
|
---|
1775 | reconstructedFile.open(reconstructedName.c_str(), false);
|
---|
1776 |
|
---|
1777 | //reconstruct the original columns from the compressed file.
|
---|
1778 | string origChecksumStr;
|
---|
1779 | string origDatasum;
|
---|
1780 |
|
---|
1781 | //reset tablename value so that it is re-read from compressed table's header
|
---|
1782 | tableName = "";
|
---|
1783 |
|
---|
1784 | /************************************************************************************
|
---|
1785 | * Reconstruction setup done. Rebuild original header
|
---|
1786 | ************************************************************************************/
|
---|
1787 |
|
---|
1788 | //re-tranlate the keys
|
---|
1789 | for (fits::Table::Keys::const_iterator it=header2.begin(); it!= header2.end(); it++)
|
---|
1790 | {
|
---|
1791 | string k = it->first;
|
---|
1792 | if (k == "XTENSION" || k == "BITPIX" || k == "PCOUNT" || k == "GCOUNT" ||
|
---|
1793 | k == "TFIELDS" || k == "ZTABLE" || k == "ZNAXIS1" || k == "ZNAXIS2" ||
|
---|
1794 | k == "ZHEAPPTR" || k == "ZPCOUNT" || k == "ZTILELEN" || k == "THEAP" ||
|
---|
1795 | k == "CHECKSUM" || k == "DATASUM" || k == "FCTCPVER")
|
---|
1796 | {
|
---|
1797 | continue;
|
---|
1798 | }
|
---|
1799 |
|
---|
1800 | if (k == "ZCHKSUM")
|
---|
1801 | {
|
---|
1802 | reconstructedFile.SetKeyComment("CHECKSUM", it->second.comment);
|
---|
1803 | origChecksumStr = it->second.value;
|
---|
1804 | continue;
|
---|
1805 | }
|
---|
1806 |
|
---|
1807 | if (k == "ZDTASUM")
|
---|
1808 | {
|
---|
1809 | reconstructedFile.SetKeyComment("DATASUM", it->second.comment);
|
---|
1810 | origDatasum = it->second.value;
|
---|
1811 | continue;
|
---|
1812 | }
|
---|
1813 |
|
---|
1814 | if (k == "EXTNAME")
|
---|
1815 | {
|
---|
1816 | tableName = it->second.value;
|
---|
1817 | }
|
---|
1818 |
|
---|
1819 | k = k.substr(0,5);
|
---|
1820 |
|
---|
1821 | if (k == "TTYPE")
|
---|
1822 | {//we have an original column name here.
|
---|
1823 | //manually deal with these in order to preserve the ordering (easier than re-constructing yet another list on the fly)
|
---|
1824 | continue;
|
---|
1825 | }
|
---|
1826 |
|
---|
1827 | if (k == "TFORM" || k == "NAXIS" || k == "ZCTYP" )
|
---|
1828 | {
|
---|
1829 | continue;
|
---|
1830 | }
|
---|
1831 |
|
---|
1832 | if (k == "ZFORM" || k == "ZTYPE")
|
---|
1833 | {
|
---|
1834 | string tmpKey = it->second.fitsString;
|
---|
1835 | tmpKey[0] = 'T';
|
---|
1836 | reconstructedFile.SetKeyFromFitsString(tmpKey);
|
---|
1837 | continue;
|
---|
1838 | }
|
---|
1839 |
|
---|
1840 | reconstructedFile.SetKeyFromFitsString(it->second.fitsString);
|
---|
1841 | }
|
---|
1842 |
|
---|
1843 | if (tableName == "")
|
---|
1844 | {
|
---|
1845 | cout << "Error: table name could not be found. Aborting" << endl;
|
---|
1846 | return -1;
|
---|
1847 | }
|
---|
1848 |
|
---|
1849 | //Restore the original columns
|
---|
1850 | for (uint32_t numCol=1; numCol<10000; numCol++)
|
---|
1851 | {
|
---|
1852 | ostringstream str;
|
---|
1853 | str << numCol;
|
---|
1854 | if (!verifyFile.HasKey("TTYPE"+str.str())) break;
|
---|
1855 |
|
---|
1856 | string ttype = verifyFile.GetStr("TTYPE"+str.str());
|
---|
1857 | string tform = verifyFile.GetStr("ZFORM"+str.str());
|
---|
1858 | char type = tform[tform.size()-1];
|
---|
1859 | string number = tform.substr(0, tform.size()-1);
|
---|
1860 | int numElems = atoi(number.c_str());
|
---|
1861 |
|
---|
1862 | if (number == "") numElems=1;
|
---|
1863 |
|
---|
1864 | reconstructedFile.AddColumn(numElems, type, ttype, "", "", false);
|
---|
1865 | }
|
---|
1866 |
|
---|
1867 | reconstructedFile.WriteTableHeader(tableName.c_str());
|
---|
1868 |
|
---|
1869 | /************************************************************************************
|
---|
1870 | * Original header restored. Do the data
|
---|
1871 | ************************************************************************************/
|
---|
1872 |
|
---|
1873 | //set pointers to the readout data to later be able to gather it to "buffer".
|
---|
1874 | readPointers.clear();
|
---|
1875 | readOffsets.clear();
|
---|
1876 | readElemSize.clear();
|
---|
1877 | readNumElems.clear();
|
---|
1878 | for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end(); it++)
|
---|
1879 | {
|
---|
1880 | readPointers.push_back(verifyFile.SetPtrAddress(it->first));
|
---|
1881 | readOffsets.push_back(it->second.offset);
|
---|
1882 | readElemSize.push_back(it->second.size);
|
---|
1883 | readNumElems.push_back(it->second.num);
|
---|
1884 | }
|
---|
1885 |
|
---|
1886 | //do the actual reconstruction work
|
---|
1887 | uint32_t i=1;
|
---|
1888 | while (i<=verifyFile.GetNumRows() && verifyFile.GetNextRow())
|
---|
1889 | {
|
---|
1890 | int count=0;
|
---|
1891 | for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end();it++)
|
---|
1892 | {
|
---|
1893 | memcpy(&buffer[readOffsets[count]], readPointers[count], readElemSize[count]*readNumElems[count]);
|
---|
1894 | count++;
|
---|
1895 | }
|
---|
1896 | if (displayText) cout << "\r Row " << i << flush;
|
---|
1897 | reconstructedFile.WriteRow(buffer, rowWidth);
|
---|
1898 | i++;
|
---|
1899 | }
|
---|
1900 |
|
---|
1901 | if (displayText) cout << endl;
|
---|
1902 |
|
---|
1903 | //close reconstruction input and output
|
---|
1904 | verifyFile.close();
|
---|
1905 | reconstructedFile.close();
|
---|
1906 |
|
---|
1907 | //get original and reconstructed checksum and datasum
|
---|
1908 | std::pair<string, int> origChecksum = make_pair(origChecksumStr, atoi(origDatasum.c_str()));
|
---|
1909 | std::pair<string, int> newChecksum = reconstructedFile.GetChecksumData();
|
---|
1910 |
|
---|
1911 | //verify that no mistake was made
|
---|
1912 | if (origChecksum.second != newChecksum.second)
|
---|
1913 | {
|
---|
1914 | cout << "ERROR: datasums are NOT identical: " << (uint32_t)(origChecksum.second) << " vs " << (uint32_t)(newChecksum.second) << endl;
|
---|
1915 | return -1;
|
---|
1916 | }
|
---|
1917 | if (origChecksum.first != newChecksum.first)
|
---|
1918 | {
|
---|
1919 | cout << "WARNING: checksums are NOT Identical: " << origChecksum.first << " vs " << newChecksum.first << endl;
|
---|
1920 | }
|
---|
1921 | else
|
---|
1922 | {
|
---|
1923 | if (displayText) cout << "Data is identical" << endl;
|
---|
1924 | }
|
---|
1925 |
|
---|
1926 | delete[] buffer;
|
---|
1927 | return 0;
|
---|
1928 |
|
---|
1929 | }
|
---|
1930 |
|
---|