source: trunk/FACT++/src/fitsCompressor.cc@ 16773

Last change on this file since 16773 was 16506, checked in by lyard, 13 years ago
changed output text of fitsCompressor
File size: 80.0 KB
Line 
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
20using namespace std;
21
22class CompressedFitsFile
23{
24public:
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
375class 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
448const uint32_t CompressedFitsFile::_THREAD_WAIT_ = 0;
449const uint32_t CompressedFitsFile::_THREAD_COMPRESS_ = 1;
450const uint32_t CompressedFitsFile::_THREAD_DECOMPRESS_= 2;
451const uint32_t CompressedFitsFile::_THREAD_WRITE_ = 3;
452const uint32_t CompressedFitsFile::_THREAD_READ_ = 4;
453const uint32_t CompressedFitsFile::_THREAD_EXIT_ = 5;
454
455#define COMPRESSED_FLAG 0x1
456#define UNCOMPRESSED_FLAG 0x0
457
458template<>
459void 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 */
480string 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 */
521string 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
558CompressedFitsFile::HeaderEntry CompressedFitsFile::_dummyHeaderEntry;
559/****************************************************************
560 * SUPER CLASS DEFAULT CONSTRUCTOR
561 ****************************************************************/
562CompressedFitsFile::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 ****************************************************************/
591CompressedFitsFile::~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 ****************************************************************/
613bool 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 + _columns.size()]; //use a bit more memory for compression flags
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 ****************************************************************/
648CompressedFitsWriter::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 ****************************************************************/
679CompressedFitsWriter::~CompressedFitsWriter()
680{
681 pthread_mutex_destroy(&_mutex);
682}
683
684/****************************************************************
685 * SET THE POINTER TO THE DRS CALIBRATION
686 ****************************************************************/
687void CompressedFitsWriter::setDrsCalib(int16_t* data)
688{
689 _drsCalibData = data;
690}
691
692/****************************************************************
693 * SET NUM WORKING THREADS
694 ****************************************************************/
695bool 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 ****************************************************************/
721void 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" , "'ZDrsCellOffsets' ", "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 ****************************************************************/
793bool 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 ****************************************************************/
855bool 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 ****************************************************************/
893bool 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 ****************************************************************/
918void 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 ****************************************************************/
974void 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 ****************************************************************/
1024void 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 ****************************************************************/
1041bool 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 int64_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 if (compressedOffset < 0) return false;
1080 _catalog[i][j].second = compressedOffset;
1081 compressedOffset += _catalog[i][j].first;
1082 //special case if entry has zero length
1083 if (_catalog[i][j].first == 0) _catalog[i][j].second = 0;
1084 }
1085 setHeaderKey(HeaderEntry("PCOUNT", heapSize, "size of special data area"));
1086 }
1087 writeHeader(true);
1088 ostringstream str;
1089 str << _checksum.val();
1090 setHeaderKey(HeaderEntry("DATASUM", str.str(), ""));
1091 addHeaderChecksum(_checksum);
1092 setHeaderKey(HeaderEntry("CHECKSUM", _checksum.str(), ""));
1093 //update header value
1094 writeHeader();
1095 //update file length
1096 long here = _file.tellp();
1097 if (here%2880)
1098 {
1099 vector<char> nullVec(2880 - here%2880, 0);
1100 _file.write(nullVec.data(), 2880 - here%2880);
1101 }
1102 _file.close();
1103 return true;
1104}
1105
1106/****************************************************************
1107 * COPY TRANSPOSE TILE
1108 ****************************************************************/
1109void CompressedFitsWriter::copyTransposeTile(uint32_t index)
1110{
1111 uint32_t thisRoundNumRows = (_totalNumRows%_numRowsPerTile) ? _totalNumRows%_numRowsPerTile : _numRowsPerTile;
1112 //copy the tile and transpose it
1113 uint32_t offset = 0;
1114 for (uint32_t i=0;i<_columns.size();i++)
1115 {
1116 switch (_columns[i].getCompression())
1117 {
1118 case UNCOMPRESSED:
1119 case SMOOTHMAN:
1120 for (uint32_t k=0;k<thisRoundNumRows;k++)
1121 {//regular, "semi-transposed" copy
1122 memcpy(&(_transposedBuffer[index][offset]), &_buffer[k*_rowWidth + _columns[i].offset()], _columns[i].sizeOfElems()*_columns[i].numElems());
1123 offset += _columns[i].sizeOfElems()*_columns[i].numElems();
1124 }
1125 break;
1126 default :
1127 for (int j=0;j<_columns[i].numElems();j++)
1128 for (uint32_t k=0;k<thisRoundNumRows;k++)
1129 {//transposed copy
1130 memcpy(&(_transposedBuffer[index][offset]), &_buffer[k*_rowWidth + _columns[i].offset() + _columns[i].sizeOfElems()*j], _columns[i].sizeOfElems());
1131 offset += _columns[i].sizeOfElems();
1132 }
1133 break;
1134
1135 };
1136 }
1137}
1138
1139/****************************************************************
1140 * WRITE BINARY ROW
1141 ****************************************************************/
1142bool CompressedFitsWriter::writeBinaryRow(const char* bufferToWrite)
1143{
1144 if (_totalNumRows == 0)
1145 writeHeader();
1146
1147 memcpy(&_buffer[_rowWidth*(_totalNumRows%_numRowsPerTile)], bufferToWrite, _rowWidth);
1148 _totalNumRows++;
1149 if (_totalNumRows%_numRowsPerTile == 0)
1150 {
1151 //which is the next thread that we should use ?
1152 while (_threadStatus[_threadLooper] == _THREAD_COMPRESS_)
1153 usleep(100000);
1154
1155 copyTransposeTile(_threadLooper);
1156
1157 while (_threadStatus[_threadLooper] != _THREAD_WAIT_)
1158 usleep(100000);
1159 _threadNumRows[_threadLooper] = _totalNumRows;
1160 _threadStatus[_threadLooper] = _THREAD_COMPRESS_;
1161 _threadLooper = (_threadLooper+1)%_numThreads;
1162 }
1163 return _file.good();
1164}
1165
1166/****************************************************************
1167 * COMPRESS BUFFER
1168 ****************************************************************/
1169uint32_t CompressedFitsWriter::compressUNCOMPRESSED(char* dest, const char* src, uint32_t numRows, uint32_t sizeOfElems, uint32_t numRowElems)
1170{
1171 memcpy(dest, src, numRows*sizeOfElems*numRowElems);
1172 return numRows*sizeOfElems*numRowElems;
1173}
1174
1175/****************************************************************
1176 * COMPRESS BUFFER
1177 ****************************************************************/
1178uint32_t CompressedFitsWriter::compressHUFFMAN(char* dest, const char* src, uint32_t numRows, uint32_t sizeOfElems, uint32_t numRowElems)
1179{
1180 string huffmanOutput;
1181 uint32_t previousHuffmanSize = 0;
1182 if (numRows < 2)
1183 {//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.
1184 return numRows*sizeOfElems*numRowElems + 1000;
1185 }
1186 if (sizeOfElems < 2 )
1187 {
1188 cout << "Fatal ERROR: HUFMANN can only encode short or longer types" << endl;
1189 return 0;
1190 }
1191 uint32_t huffmanOffset = 0;
1192 for (uint32_t j=0;j<numRowElems;j++)
1193 {
1194 Huffman::Encode(huffmanOutput,
1195 reinterpret_cast<const uint16_t*>(&src[j*sizeOfElems*numRows]),
1196 numRows*(sizeOfElems/2));
1197 reinterpret_cast<uint32_t*>(&dest[huffmanOffset])[0] = huffmanOutput.size() - previousHuffmanSize;
1198 huffmanOffset += sizeof(uint32_t);
1199 previousHuffmanSize = huffmanOutput.size();
1200 }
1201 const size_t totalSize = huffmanOutput.size() + huffmanOffset;
1202
1203 //only copy if not larger than not-compressed size
1204 if (totalSize < numRows*sizeOfElems*numRowElems)
1205 memcpy(&dest[huffmanOffset], huffmanOutput.data(), huffmanOutput.size());
1206
1207 return totalSize;
1208}
1209
1210/****************************************************************
1211 * COMPRESS BUFFER
1212 ****************************************************************/
1213uint32_t CompressedFitsWriter::compressSMOOTHMAN(char* dest, char* src, uint32_t numRows, uint32_t sizeOfElems, uint32_t numRowElems)
1214{
1215 uint32_t colWidth = numRowElems;
1216 for (int j=colWidth*numRows-1;j>1;j--)
1217 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;
1218 //call the huffman transposed
1219 return compressHUFFMAN(dest, src, numRowElems, sizeOfElems, numRows);
1220}
1221
1222/****************************************************************
1223 * COMPRESS BUFFER
1224 ****************************************************************/
1225uint32_t CompressedFitsWriter::compressBuffer(uint32_t threadIndex)
1226{
1227 uint32_t thisRoundNumRows = (_threadNumRows[threadIndex]%_numRowsPerTile) ? _threadNumRows[threadIndex]%_numRowsPerTile : _numRowsPerTile;
1228 uint32_t offset=0;
1229 uint32_t currentCatalogRow = (_threadNumRows[threadIndex]-1)/_numRowsPerTile;
1230 int64_t compressedOffset = 0;
1231
1232 //now compress each column one by one by calling compression on arrays
1233 for (uint32_t i=0;i<_columns.size();i++)
1234 {
1235 _catalog[currentCatalogRow][i].second = compressedOffset;
1236
1237 uint32_t compression = _columns[i].getCompression();
1238
1239 if (_columns[i].numElems() == 0) continue;
1240 //set the default byte telling if uncompressed the compressed Flag
1241 int64_t previousOffset = compressedOffset;
1242 _compressedBuffer[threadIndex][compressedOffset++] = COMPRESSED_FLAG;
1243 switch (compression)
1244 {
1245 case UNCOMPRESSED:
1246 compressedOffset += compressUNCOMPRESSED(&(_compressedBuffer[threadIndex][compressedOffset]), &(_transposedBuffer[threadIndex][offset]), thisRoundNumRows, _columns[i].sizeOfElems(), _columns[i].numElems());
1247 break;
1248 case SMOOTHMAN:
1249 compressedOffset += compressSMOOTHMAN(&(_compressedBuffer[threadIndex][compressedOffset]), &(_transposedBuffer[threadIndex][offset]), thisRoundNumRows, _columns[i].sizeOfElems(), _columns[i].numElems());
1250 break;
1251
1252 default:
1253 ;
1254 };
1255 //check if compressed size is larger than uncompressed
1256 if (compression != UNCOMPRESSED &&
1257 compressedOffset - previousOffset > _columns[i].sizeOfElems()*_columns[i].numElems()*thisRoundNumRows+1)
1258 {//if so set flag and redo it uncompressed
1259 compressedOffset = previousOffset;
1260 _compressedBuffer[threadIndex][compressedOffset++] = UNCOMPRESSED_FLAG;
1261 compressedOffset += compressUNCOMPRESSED(&(_compressedBuffer[threadIndex][compressedOffset]), &(_transposedBuffer[threadIndex][offset]), thisRoundNumRows, _columns[i].sizeOfElems(), _columns[i].numElems());
1262 }
1263 offset += thisRoundNumRows*_columns[i].sizeOfElems()*_columns[i].numElems();
1264 _catalog[currentCatalogRow][i].first = compressedOffset - _catalog[currentCatalogRow][i].second;
1265 }
1266 return compressedOffset;
1267}
1268
1269/****************************************************************
1270 * WRITE COMPRESS DATA TO DISK
1271 ****************************************************************/
1272bool CompressedFitsWriter::writeCompressedDataToDisk(uint32_t threadID, uint32_t sizeToWrite)
1273{
1274 char* checkSumPointer = _compressedBuffer[threadID];
1275 int32_t extraBytes = 0;
1276 uint32_t sizeToChecksum = sizeToWrite;
1277 if (_checkOffset != 0)
1278 {//should we extend the array to the left ?
1279 sizeToChecksum += _checkOffset;
1280 checkSumPointer -= _checkOffset;
1281 memset(checkSumPointer, 0, _checkOffset);
1282 }
1283 if (sizeToChecksum%4 != 0)
1284 {//should we extend the array to the right ?
1285 extraBytes = 4 - (sizeToChecksum%4);
1286 memset(checkSumPointer+sizeToChecksum, 0,extraBytes);
1287 sizeToChecksum += extraBytes;
1288 }
1289 //do the checksum
1290 _checksum.add(checkSumPointer, sizeToChecksum);
1291 _checkOffset = (4 - extraBytes)%4;
1292 //write data to disk
1293 _file.write(_compressedBuffer[threadID], sizeToWrite);
1294 return _file.good();
1295}
1296
1297/****************************************************************
1298 * WRITER THREAD LOOP
1299 ****************************************************************/
1300void* CompressedFitsWriter::threadFunction(void* context)
1301{
1302 CompressedFitsWriter* myself =static_cast<CompressedFitsWriter*>(context);
1303
1304 uint32_t myID = 0;
1305 pthread_mutex_lock(&(myself->_mutex));
1306 myID = myself->_threadIndex++;
1307 pthread_mutex_unlock(&(myself->_mutex));
1308 uint32_t threadToWaitForBeforeWriting = (myID == 0) ? myself->_numThreads-1 : myID-1;
1309 while (myself->_threadStatus[myID] != _THREAD_EXIT_)
1310 {
1311 while (myself->_threadStatus[myID] == _THREAD_WAIT_)
1312 usleep(100000);
1313 if (myself->_threadStatus[myID] != _THREAD_COMPRESS_)
1314 continue;
1315 uint32_t numBytes = myself->compressBuffer(myID);
1316 myself->_threadStatus[myID] = _THREAD_WRITE_;
1317
1318 //wait for the previous data to be written
1319 while (myself->_threadIndex != threadToWaitForBeforeWriting)
1320 usleep(1000);
1321 //do the actual writing to disk
1322 pthread_mutex_lock(&(myself->_mutex));
1323 myself->writeCompressedDataToDisk(myID, numBytes);
1324 myself->_threadIndex = myID;
1325 pthread_mutex_unlock(&(myself->_mutex));
1326 myself->_threadStatus[myID] = _THREAD_WAIT_;
1327 }
1328 return NULL;
1329}
1330
1331/****************************************************************
1332 * PRINT USAGE
1333 ****************************************************************/
1334void printUsage()
1335{
1336 cout << endl;
1337 cout << "The FACT-Fits compressor reads an input Fits file from FACT"
1338 " and compresses it.\n It can use a drs calibration in order to"
1339 " improve the compression ratio. If so, the input calibration"
1340 " is embedded into the compressed file.\n"
1341 " By default, the Data column will be compressed using SMOOTHMAN (Thomas' algorithm)"
1342 " while other columns will be compressed with the AMPLITUDE coding (Veritas)"
1343 "Usage: Compressed_Fits_Test <inputFile>";
1344 cout << endl;
1345}
1346
1347/****************************************************************
1348 * PRINT HELP
1349 ****************************************************************/
1350void printHelp()
1351{
1352 cout << endl;
1353 cout << "The inputFile is required. It must have fits in its filename and the compressed file will be written in the same folder. "
1354 "The fz extension will be added, replacing the .gz one if required \n"
1355 "If output is specified, then it will replace the automatically generated output filename\n"
1356 "If --drs, followed by a drs calib then it will be applied to the data before compressing\n"
1357 "rowPerTile can be used to tune how many rows are in each tile. Default is 100\n"
1358 "threads gives the number of threads to use. Cannot be less than the default (1)\n"
1359 "compression explicitely gives the compression scheme to use for a given column. The syntax is:\n"
1360 "<ColumnName>=<CompressionScheme> with <CompressionScheme> one of the following:\n"
1361 "UNCOMPRESSED\n"
1362 "AMPLITUDE\n"
1363 "HUFFMAN\n"
1364 "SMOOTHMAN\n"
1365 "INT_WAVELET\n"
1366 "\n"
1367 "--quiet removes any textual output, except error messages\n"
1368 "--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."
1369 ;
1370 cout << endl << endl;
1371}
1372
1373/****************************************************************
1374 * SETUP CONFIGURATION
1375 ****************************************************************/
1376void setupConfiguration(Configuration& conf)
1377{
1378 po::options_description configs("FitsCompressor options");
1379 configs.add_options()
1380 ("inputFile,i", vars<string>(), "Input file")
1381 ("drs,d", var<string>(), "Input drs calibration file")
1382 ("rowPerTile,r", var<unsigned int>(), "Number of rows per tile. Default is 100")
1383 ("output,o", var<string>(), "Output file. If empty, .fz is appened to the original name")
1384 ("threads,t", var<unsigned int>(), "Number of threads to use for compression")
1385 ("compression,c", vars<string>(), "which compression to use for which column. Syntax <colName>=<compressionScheme>")
1386 ("quiet,q", po_switch(), "Should the program display any text at all ?")
1387 ("verify,v", po_switch(), "Should we verify the data that has been compressed ?")
1388 ;
1389 po::positional_options_description positional;
1390 positional.add("inputFile", -1);
1391 conf.AddOptions(configs);
1392 conf.SetArgumentPositions(positional);
1393}
1394
1395/****************************************************************
1396 * MAIN
1397 ****************************************************************/
1398int main(int argc, const char** argv)
1399{
1400 Configuration conf(argv[0]);
1401 conf.SetPrintUsage(printUsage);
1402 setupConfiguration(conf);
1403
1404 if (!conf.DoParse(argc, argv, printHelp))
1405 return -1;
1406
1407 //initialize the file names to nothing.
1408 string fileNameIn = "";
1409 string fileNameOut = "";
1410 string drsFileName = "";
1411 uint32_t numRowsPerTile = 100;
1412 bool displayText=true;
1413
1414 //parse configuration
1415 if (conf.Get<bool>("quiet")) displayText = false;
1416 const vector<string> inputFileNameVec = conf.Vec<string>("inputFile");
1417 if (inputFileNameVec.size() != 1)
1418 {
1419 cout << "Error: ";
1420 if (inputFileNameVec.size() == 0) cout << "no";
1421 else cout << inputFileNameVec.size();
1422 cout << " input file(s) given. Expected one. Aborting. Input:" << endl;;
1423 for (unsigned int i=0;i<inputFileNameVec.size(); i++)
1424 cout << inputFileNameVec[i] << endl;
1425 return -1;
1426 }
1427
1428 //Assign the input filename
1429 fileNameIn = inputFileNameVec[0];
1430
1431 //Check if we have a drs calib too
1432 if (conf.Has("drs")) drsFileName = conf.Get<string>("drs");
1433
1434 //Should we verify the data ?
1435 bool verifyDataPlease = false;
1436 if (conf.Has("verify")) verifyDataPlease = conf.Get<bool>("verify");
1437
1438
1439 //should we use a specific output filename ?
1440 if (conf.Has("output"))
1441 fileNameOut = conf.Get<string>("output");
1442 else
1443 {
1444 size_t pos = fileNameIn.find(".fits");
1445 if (pos == string::npos)
1446 {
1447 cout << "ERROR: input file does not seems ot be fits. Aborting." << endl;
1448 return -1;
1449 }
1450 fileNameOut = fileNameIn.substr(0, pos) + ".fz";
1451 }
1452
1453 //should we use specific compression on some columns ?
1454 const vector<string> columnsCompression = conf.Vec<string>("compression");
1455
1456 //split up values between column names and compression scheme
1457 vector<std::pair<string, string>> compressions;
1458 for (unsigned int i=0;i<columnsCompression.size();i++)
1459 {
1460 size_t pos = columnsCompression[i].find_first_of("=");
1461 if (pos == string::npos)
1462 {
1463 cout << "ERROR: Something wrong occured while parsing " << columnsCompression[i] << ". Aborting." << endl;
1464 return -1;
1465 }
1466 string comp = columnsCompression[i].substr(pos+1);
1467 if (comp != "UNCOMPRESSED" && comp != "AMPLITUDE" && comp != "HUFFMAN" &&
1468 comp != "SMOOTHMAN" && comp != "INT_WAVELET")
1469 {
1470 cout << "Unkown compression scheme requested (" << comp << "). Aborting." << endl;
1471 return -1;
1472 }
1473 compressions.push_back(make_pair(columnsCompression[i].substr(0, pos), comp));
1474 }
1475
1476 //How many rows per tile should we use ?
1477 if (conf.Has("rowPerTile")) numRowsPerTile = conf.Get<unsigned int>("rowPerTile");
1478
1479 /************************************************************************************
1480 * Done reading configuration. Open relevant files
1481 ************************************************************************************/
1482
1483 //Open input's fits file
1484 factfits inFile(fileNameIn);
1485
1486 if (inFile.IsCompressedFITS())
1487 {
1488 cout << "ERROR: input file is already a compressed fits. Cannot be compressed again: Aborting." << endl;
1489 return -1;
1490 }
1491
1492 //decide how many tiles should be put in the compressed file
1493 uint32_t originalNumRows = inFile.GetNumRows();
1494 uint32_t numTiles = (originalNumRows%numRowsPerTile) ? (originalNumRows/numRowsPerTile)+1 : originalNumRows/numRowsPerTile;
1495 CompressedFitsWriter outFile(numTiles, numRowsPerTile);
1496
1497 //should we use a specific number of threads for compressing ?
1498 unsigned int numThreads = 1;
1499 if (conf.Has("threads"))
1500 {
1501 numThreads = conf.Get<unsigned int>("threads");
1502 outFile.setNumWorkingThreads(numThreads);
1503 }
1504
1505 if (!outFile.open(fileNameOut))
1506 {
1507 cout << "Error: could not open " << fileNameOut << " for writing" << endl;
1508 return -1;
1509 }
1510
1511 //Because the file to open MUST be given by the constructor, I must use a pointer instead
1512 factfits* drsFile = NULL;
1513 //try to open the Drs file. If any.
1514 if (drsFileName != "")
1515 {
1516 try
1517 {
1518 drsFile = new factfits(drsFileName);
1519 }
1520 catch (...)
1521 {
1522 cout << "Error: could not open " << drsFileName << " for calibration" << endl;
1523 return -1;
1524 }
1525 }
1526
1527 if (displayText)
1528 {
1529 cout << endl;
1530 cout << "**********************" << endl;
1531 cout << "Will compress from : " << fileNameIn << endl;
1532 cout << "to : " << fileNameOut << endl;
1533 if (drsFileName != "")
1534 cout << "while calibrating with: " << drsFileName << endl;
1535 cout << "Compression will use : " << numThreads << " worker threads" << endl;
1536 cout << "Data will be verified : ";
1537 if (verifyDataPlease)
1538 cout << "yes" << endl;
1539 else
1540 cout << "no (WARNING !)" << endl;
1541 cout << "**********************" << endl;
1542 cout << endl;
1543 }
1544
1545 /************************************************************************************
1546 * Done opening input files. Allocate memory and configure output file
1547 ************************************************************************************/
1548
1549 //allocate the buffer for temporary storage of each read/written row
1550 uint32_t rowWidth = inFile.GetUInt("NAXIS1");
1551 char* buffer = new char[rowWidth];
1552
1553 //get the source columns
1554 const fits::Table::Columns& columns = inFile.GetColumns();
1555 const fits::Table::SortedColumns& sortedColumns = inFile.GetSortedColumns();
1556 if (displayText)
1557 cout << "Input file has " << columns.size() << " columns and " << inFile.GetNumRows() << " rows" << endl;
1558
1559 //Add columns.
1560 for (uint32_t i=0;i<sortedColumns.size(); i++)
1561 {
1562 //get column name
1563 ostringstream str;
1564 str << "TTYPE" << i+1;
1565 string colName = inFile.GetStr(str.str());
1566 if (displayText)
1567 {
1568 cout << "Column " << colName;
1569 for (uint32_t j=colName.size();j<21;j++)
1570 cout << " ";
1571 cout << " -> ";
1572 }
1573
1574 //first lets see if we do have an explicit request
1575 bool explicitRequest = false;
1576 for (unsigned int j=0;j<compressions.size();j++)
1577 {
1578 if (compressions[j].first == colName)
1579 {
1580 explicitRequest = true;
1581 if (displayText) cout << compressions[j].second << endl;
1582 if (compressions[j].second == "UNCOMPRESSED")
1583 outFile.addColumn(CompressedFitsFile::ColumnEntry(colName, sortedColumns[i].type, sortedColumns[i].num, CompressedFitsFile::UNCOMPRESSED));
1584 if (compressions[j].second == "SMOOTHMAN")
1585 outFile.addColumn(CompressedFitsFile::ColumnEntry(colName, sortedColumns[i].type, sortedColumns[i].num, CompressedFitsFile::SMOOTHMAN));
1586 break;
1587 }
1588 }
1589
1590 if (explicitRequest) continue;
1591
1592 if (colName != "Data")
1593 {
1594 if (displayText) cout << "UNCOMPRESSED" << endl;
1595 outFile.addColumn(CompressedFitsFile::ColumnEntry(colName, sortedColumns[i].type, sortedColumns[i].num, CompressedFitsFile::UNCOMPRESSED));
1596 }
1597 else
1598 {
1599 if (displayText) cout << "SMOOTHMAN" << endl;
1600 outFile.addColumn(CompressedFitsFile::ColumnEntry(colName, sortedColumns[i].type, sortedColumns[i].num, CompressedFitsFile::SMOOTHMAN));
1601 }
1602 }
1603
1604 //translate original header entries to their Z-version
1605 const fits::Table::Keys& header = inFile.GetKeys();
1606 for (fits::Table::Keys::const_iterator i=header.begin(); i!= header.end(); i++)
1607 {
1608 string k = i->first;//header[i].key();
1609 if (k == "XTENSION" || k == "BITPIX" || k == "PCOUNT" || k == "GCOUNT" || k == "TFIELDS")
1610 continue;
1611 if (k == "CHECKSUM")
1612 {
1613 outFile.setHeaderKey(CompressedFitsFile::HeaderEntry("ZCHKSUM", i->second.value, i->second.comment));
1614 continue;
1615 }
1616 if (k == "DATASUM")
1617 {
1618 outFile.setHeaderKey(CompressedFitsFile::HeaderEntry("ZDTASUM", i->second.value, i->second.comment));
1619 continue;
1620 }
1621 k = k.substr(0,5);
1622 if (k == "TTYPE" || k == "TFORM")
1623 {
1624 string tmpKey = i->second.fitsString;
1625 tmpKey[0] = 'Z';
1626 outFile.setHeaderKey(tmpKey);
1627 continue;
1628 }
1629 if (k == "NAXIS")
1630 continue;
1631 outFile.setHeaderKey(i->second.fitsString);
1632 }
1633
1634 //deal with the DRS calib
1635 int16_t* drsCalib16 = NULL;
1636
1637 //load the drs calib. data
1638 int32_t startCellOffset = -1;
1639 if (drsFileName != "")
1640 {
1641 drsCalib16 = new int16_t[1440*1024];
1642 float* drsCalibFloat = NULL;
1643 try
1644 {
1645 drsCalibFloat = reinterpret_cast<float*>(drsFile->SetPtrAddress("BaselineMean"));
1646 }
1647 catch (...)
1648 {
1649 cout << "Could not find column BaselineMean in drs calibration file " << drsFileName << ". Aborting" << endl;
1650 return -1;
1651 }
1652
1653 //read the calibration and calculate its integer value
1654 drsFile->GetNextRow();
1655 for (uint32_t i=0;i<1440*1024;i++)
1656 drsCalib16[i] = (int16_t)(drsCalibFloat[i]*4096.f/2000.f);
1657
1658 //assign it to the ouput file
1659 outFile.setDrsCalib(drsCalib16);
1660
1661 //get the start cells offsets
1662 for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end(); it++)
1663 if (it->first == "StartCellData")
1664 {
1665 startCellOffset = it->second.offset;
1666 if (it->second.type != 'I')
1667 {
1668 cout << "Wrong type for the StartCellData Column: " << it->second.type << " instead of I expected"<< endl;
1669 return -1;
1670 }
1671 }
1672
1673 if (startCellOffset == -1)
1674 {
1675 cout << "Could not find StartCellData in input file " << fileNameIn << ". Aborting."<< endl;
1676 return -1;
1677 }
1678 }
1679
1680 /************************************************************************************
1681 * Done configuring compression. Do the real job now !
1682 ************************************************************************************/
1683
1684 if (displayText) cout << "Converting file..." << endl;
1685
1686 int numSlices = -1;
1687 int32_t dataOffset = -1;
1688
1689 //Get the pointer to the column that must be drs-calibrated
1690 for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end(); it++)
1691 if (it->first == "Data")
1692 {
1693 numSlices = it->second.num;
1694 dataOffset = it->second.offset;
1695 }
1696 if (numSlices % 1440 != 0)
1697 {
1698 cout << "seems like the number of samples is not a multiple of 1440. Aborting." << endl;
1699 return -1;
1700 }
1701 if (dataOffset == -1)
1702 {
1703 cout << "Could not find the column Data in the input file. Aborting." << endl;
1704 return -1;
1705 }
1706
1707 numSlices /= 1440;
1708
1709 //set pointers to the readout data to later be able to gather it to "buffer".
1710 vector<void*> readPointers;
1711 vector<int32_t> readOffsets;
1712 vector<int32_t> readElemSize;
1713 vector<int32_t> readNumElems;
1714 for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end(); it++)
1715 {
1716 readPointers.push_back(inFile.SetPtrAddress(it->first));
1717 readOffsets.push_back(it->second.offset);
1718 readElemSize.push_back(it->second.size);
1719 readNumElems.push_back(it->second.num);
1720 }
1721
1722 //Convert each row one after the other
1723 for (uint32_t i=0;i<inFile.GetNumRows();i++)
1724 {
1725 if (displayText) cout << "\r Row " << i+1 << flush;
1726 inFile.GetNextRow();
1727 //copy from inFile internal structures to buffer
1728 int32_t count=0;
1729 for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end();it++)
1730 {
1731 memcpy(&buffer[readOffsets[count]], readPointers[count], readElemSize[count]*readNumElems[count]);
1732 count++;
1733 }
1734
1735 if (startCellOffset != -1)
1736 {//drs calibrate this data
1737 for (int j=0;j<1440;j++)
1738 {
1739 const int thisStartCell = reinterpret_cast<int16_t*>(&buffer[startCellOffset])[j];
1740 if (thisStartCell < 0) continue;
1741 for (int k=0;k<numSlices;k++)
1742 reinterpret_cast<int16_t*>(&buffer[dataOffset])[numSlices*j + k] -= drsCalib16[1024*j + (thisStartCell+k)%1024];
1743 }
1744 }
1745 outFile.writeBinaryRow(buffer);
1746 };
1747
1748 //Get table name for later use in case the compressed file is to be verified
1749 string tableName = inFile.GetStr("EXTNAME");
1750
1751 if (displayText) cout << endl << "Done. Flushing output file..." << endl;
1752
1753 inFile.close();
1754 if (!outFile.close())
1755 {
1756 cout << "Something went wrong while writing the catalog: negative index" << endl;
1757 return false;
1758 }
1759
1760 delete[] drsCalib16;
1761
1762 if (displayText) cout << "Done." << endl;
1763
1764 /************************************************************************************
1765 * Actual job done. Should we verify what we did ?
1766 ************************************************************************************/
1767
1768 if (verifyDataPlease)
1769 {
1770 if (displayText) cout << "Now verify data..." << endl;
1771 }
1772 else
1773 return 0;
1774
1775 //get a compressed reader
1776//TEMP try to copy the file too
1777// string copyName("/scratch/copyFile.fz");
1778 string copyName("");
1779 factfits verifyFile(fileNameOut, copyName, tableName, false);
1780
1781 //and the header of the compressed file
1782 const fits::Table::Keys& header2 = verifyFile.GetKeys();
1783
1784 //get a non-compressed writer
1785 ofits reconstructedFile;
1786
1787 //figure out its name: /dev/null unless otherwise specified
1788 string reconstructedName = fileNameOut+".recons";
1789 reconstructedName = "/dev/null";
1790 reconstructedFile.open(reconstructedName.c_str(), false);
1791
1792 //reconstruct the original columns from the compressed file.
1793 string origChecksumStr;
1794 string origDatasum;
1795
1796 //reset tablename value so that it is re-read from compressed table's header
1797 tableName = "";
1798
1799 /************************************************************************************
1800 * Reconstruction setup done. Rebuild original header
1801 ************************************************************************************/
1802
1803 //re-tranlate the keys
1804 for (fits::Table::Keys::const_iterator it=header2.begin(); it!= header2.end(); it++)
1805 {
1806 string k = it->first;
1807 if (k == "XTENSION" || k == "BITPIX" || k == "PCOUNT" || k == "GCOUNT" ||
1808 k == "TFIELDS" || k == "ZTABLE" || k == "ZNAXIS1" || k == "ZNAXIS2" ||
1809 k == "ZHEAPPTR" || k == "ZPCOUNT" || k == "ZTILELEN" || k == "THEAP" ||
1810 k == "CHECKSUM" || k == "DATASUM" || k == "FCTCPVER")
1811 {
1812 continue;
1813 }
1814
1815 if (k == "ZCHKSUM")
1816 {
1817 reconstructedFile.SetKeyComment("CHECKSUM", it->second.comment);
1818 origChecksumStr = it->second.value;
1819 continue;
1820 }
1821
1822 if (k == "ZDTASUM")
1823 {
1824 reconstructedFile.SetKeyComment("DATASUM", it->second.comment);
1825 origDatasum = it->second.value;
1826 continue;
1827 }
1828
1829 if (k == "EXTNAME")
1830 {
1831 tableName = it->second.value;
1832 }
1833
1834 k = k.substr(0,5);
1835
1836 if (k == "TTYPE")
1837 {//we have an original column name here.
1838 //manually deal with these in order to preserve the ordering (easier than re-constructing yet another list on the fly)
1839 continue;
1840 }
1841
1842 if (k == "TFORM" || k == "NAXIS" || k == "ZCTYP" )
1843 {
1844 continue;
1845 }
1846
1847 if (k == "ZFORM" || k == "ZTYPE")
1848 {
1849 string tmpKey = it->second.fitsString;
1850 tmpKey[0] = 'T';
1851 reconstructedFile.SetKeyFromFitsString(tmpKey);
1852 continue;
1853 }
1854
1855 reconstructedFile.SetKeyFromFitsString(it->second.fitsString);
1856 }
1857
1858 if (tableName == "")
1859 {
1860 cout << "Error: table name from file " << fileNameOut << " could not be found. Aborting" << endl;
1861 return -1;
1862 }
1863
1864 //Restore the original columns
1865 for (uint32_t numCol=1; numCol<10000; numCol++)
1866 {
1867 ostringstream str;
1868 str << numCol;
1869 if (!verifyFile.HasKey("TTYPE"+str.str())) break;
1870
1871 string ttype = verifyFile.GetStr("TTYPE"+str.str());
1872 string tform = verifyFile.GetStr("ZFORM"+str.str());
1873 char type = tform[tform.size()-1];
1874 string number = tform.substr(0, tform.size()-1);
1875 int numElems = atoi(number.c_str());
1876
1877 if (number == "") numElems=1;
1878
1879 reconstructedFile.AddColumn(numElems, type, ttype, "", "", false);
1880 }
1881
1882 reconstructedFile.WriteTableHeader(tableName.c_str());
1883
1884 /************************************************************************************
1885 * Original header restored. Do the data
1886 ************************************************************************************/
1887
1888 //set pointers to the readout data to later be able to gather it to "buffer".
1889 readPointers.clear();
1890 readOffsets.clear();
1891 readElemSize.clear();
1892 readNumElems.clear();
1893 for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end(); it++)
1894 {
1895 readPointers.push_back(verifyFile.SetPtrAddress(it->first));
1896 readOffsets.push_back(it->second.offset);
1897 readElemSize.push_back(it->second.size);
1898 readNumElems.push_back(it->second.num);
1899 }
1900
1901 //do the actual reconstruction work
1902 uint32_t i=1;
1903 while (i<=verifyFile.GetNumRows() && verifyFile.GetNextRow())
1904 {
1905 int count=0;
1906 for (fits::Table::Columns::const_iterator it=columns.begin(); it!= columns.end();it++)
1907 {
1908 memcpy(&buffer[readOffsets[count]], readPointers[count], readElemSize[count]*readNumElems[count]);
1909 count++;
1910 }
1911 if (displayText) cout << "\r Row " << i << flush;
1912 reconstructedFile.WriteRow(buffer, rowWidth);
1913 i++;
1914 }
1915
1916 if (displayText) cout << endl;
1917
1918 //close reconstruction input and output
1919 verifyFile.close();
1920 reconstructedFile.close();
1921
1922 //get original and reconstructed checksum and datasum
1923 std::pair<string, int> origChecksum = make_pair(origChecksumStr, atoi(origDatasum.c_str()));
1924 std::pair<string, int> newChecksum = reconstructedFile.GetChecksumData();
1925
1926 //verify that no mistake was made
1927 if (origChecksum.second != newChecksum.second)
1928 {
1929 cout << "ERROR: datasums are NOT identical: " << (uint32_t)(origChecksum.second) << " vs " << (uint32_t)(newChecksum.second) << endl;
1930 return -1;
1931 }
1932 if (origChecksum.first != newChecksum.first)
1933 {
1934 cout << "WARNING: checksums are NOT Identical: " << origChecksum.first << " vs " << newChecksum.first << endl;
1935 }
1936 else
1937 {
1938 if (displayText) cout << "Ok" << endl;
1939 }
1940
1941 delete[] buffer;
1942 return 0;
1943
1944}
1945
Note: See TracBrowser for help on using the repository browser.