Changeset 16814 for trunk/Mars/mcore
- Timestamp:
- 06/12/13 15:56:12 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Mars/mcore/zfits.h
r16811 r16814 13 13 #include "fits.h" 14 14 #include "huffman.h" 15 16 17 #define FACT_RAW 0x0 18 #define FACT_SMOOTHING 0x1 19 #define FACT_HUFFMAN16 0x2 20 21 #define FACT_COL_MAJOR 'C' 22 #define FACT_ROW_MAJOR 'R' 15 23 16 24 … … 79 87 80 88 private: 89 90 //Structure helper for reading tiles headers 81 91 typedef struct TileHeader 82 92 { … … 84 94 uint32_t numRows; 85 95 uint64_t size; 86 TileHeader(uint32_t nRows=0, 87 uint64_t s=0) : id({'T', 'I', 'L', 'E'}), 96 97 TileHeader() {} 98 99 TileHeader(uint32_t nRows, 100 uint64_t s) : id({'T', 'I', 'L', 'E'}), 88 101 numRows(nRows), 89 102 size(s) 90 103 { }; 91 friend ostream& operator << (ostream& out, const TileHeader& h)92 {93 out << h.id[0] << h.id[1] << h.id[2] << h.id[3] << " num Rows: " << h.numRows << ", tile size: " << h.size;94 return out;95 }96 104 } __attribute__((__packed__)) TileHeader; 105 106 //Structure helper for reading blocks headers and compresion schemes 107 typedef struct BlockHeader 108 { 109 uint64_t size; 110 char ordering; 111 unsigned char numProcs; 112 uint16_t processings[]; 113 114 BlockHeader(uint64_t s=0, 115 char o=FACT_ROW_MAJOR, 116 unsigned char n=1) : size(s), 117 ordering(o), 118 numProcs(n) 119 {} 120 } __attribute__((__packed__)) BlockHeader; 97 121 98 122 // Do what it takes to initialize the compressed structured … … 103 127 return; 104 128 129 if (fTable.isCompressed) 130 for (auto it=fTable.sortedCols.begin(); it!= fTable.sortedCols.end(); it++) 131 if (it->comp != FACT) 132 { 133 #ifdef __EXCEPTIONS 134 throw runtime_error("ERROR: Only the FACT compression scheme is handled by this reader."); 135 #else 136 gLog << ___err___ << "ERROR: Only the FACT compression scheme is handled by this reader." << endl; 137 return; 138 #endif 139 } 140 141 fColumnOrdering.resize(fTable.sortedCols.size()); 142 for (auto it=fColumnOrdering.begin(); it != fColumnOrdering.end(); it++) 143 (*it) = FACT_ROW_MAJOR; 105 144 //Get compressed specific keywords 106 145 fNumTiles = fTable.isCompressed ? GetInt("NAXIS2") : 0; … … 129 168 vector<char> fTransposedBuffer; ///<intermediate buffer to transpose the rows 130 169 vector<char> fCompressedBuffer; ///<compressed rows 170 vector<char> fColumnOrdering; ///< ordering of the column's rows 131 171 132 172 size_t fNumTiles; ///< Total number of tiles … … 137 177 138 178 vector<vector<pair<int64_t, int64_t> > > fCatalog;///< Catalog, i.e. the main table that points to the compressed data. 139 vector<size_t> fTileSize; ///< size in bytes of each compressed tile140 vector<vector<size_t> > fTileOffsets; ///< offset from start of tile of a given compressed column179 vector<size_t> fTileSize; ///< size in bytes of each compressed tile 180 vector<vector<size_t> > fTileOffsets; ///< offset from start of tile of a given compressed column 141 181 142 182 void AllocateBuffers() … … 148 188 149 189 fTransposedBuffer.resize(fTable.bytes_per_row*fNumRowsPerTile); 150 fCompressedBuffer.resize(fTable.bytes_per_row*fNumRowsPerTile + fTable.num_cols ); //use a bit more memory for compression flags190 fCompressedBuffer.resize(fTable.bytes_per_row*fNumRowsPerTile + fTable.num_cols*(sizeof(BlockHeader)+256)); //use a bit more memory for block headers 151 191 } 152 192 … … 266 306 const char *src = fTransposedBuffer.data(); 267 307 268 for (auto it=fTable.sortedCols.begin(); it!=fTable.sortedCols.end(); it++) 308 uint32_t i=0; 309 for (auto it=fTable.sortedCols.begin(); it!=fTable.sortedCols.end(); it++, i++) 269 310 { 270 311 char *buffer = fBuffer.data() + it->offset; // pointer to column (destination buffer) 271 312 272 switch ( it->comp)313 switch (fColumnOrdering[i]) 273 314 { 274 case UNCOMPRESSED: 275 case SMOOTHMAN: 276 // regular, "semi-transposed" copy 277 for (char *dest=buffer; dest<buffer+thisRoundNumRows*fTable.bytes_per_row; dest+=fTable.bytes_per_row) // row-by-row 278 { 279 memcpy(dest, src, it->bytes); 280 src += it->bytes; // next column 281 } 282 break; 283 284 default: 285 // transposed copy 286 for (char *elem=buffer; elem<buffer+it->bytes; elem+=it->size) // element-by-element (arrays) 287 { 288 for (char *dest=elem; dest<elem+thisRoundNumRows*fTable.bytes_per_row; dest+=fTable.bytes_per_row) // row-by-row 315 case FACT_ROW_MAJOR: 316 // regular, "semi-transposed" copy 317 for (char *dest=buffer; dest<buffer+thisRoundNumRows*fTable.bytes_per_row; dest+=fTable.bytes_per_row) // row-by-row 289 318 { 290 memcpy(dest, src, it-> size);291 src += it-> size; // next element319 memcpy(dest, src, it->bytes); 320 src += it->bytes; // next column 292 321 } 293 } 322 break; 323 324 case FACT_COL_MAJOR: 325 // transposed copy 326 for (char *elem=buffer; elem<buffer+it->bytes; elem+=it->size) // element-by-element (arrays) 327 { 328 for (char *dest=elem; dest<elem+thisRoundNumRows*fTable.bytes_per_row; dest+=fTable.bytes_per_row) // row-by-row 329 { 330 memcpy(dest, src, it->size); 331 src += it->size; // next element 332 } 333 } 334 break; 335 default: 336 #ifdef __EXCEPTIONS 337 throw runtime_error("Unkown column ordering scheme"); 338 #else 339 gLog << ___err___ << "ERROR - unkown column ordering scheme" << endl; 340 return; 341 #endif 294 342 break; 295 343 }; … … 305 353 uint32_t UncompressUNCOMPRESSED(char* dest, 306 354 const char* src, 307 uint32_t numRows, 308 uint32_t sizeOfElems, 309 uint32_t numRowElems) 310 { 311 memcpy(dest, src, numRows*sizeOfElems*numRowElems); 312 return numRows*sizeOfElems*numRowElems; 355 uint32_t numElems, 356 uint32_t sizeOfElems) 357 { 358 memcpy(dest, src, numElems*sizeOfElems); 359 return numElems*sizeOfElems; 313 360 } 314 361 315 362 // Read a bunch of data compressed with the Huffman algorithm 316 uint32_t UncompressHUFFMAN(char* dest, 317 const char* src, 318 uint32_t , 319 uint32_t sizeOfElems, 320 uint32_t numRowElems) 321 { 322 if (sizeOfElems < 2) 323 { 324 cout << "Error, Huffman only works on shorts or longer types. (here: " << sizeOfElems << "). Aborting." << endl; 325 return -1; 326 } 327 363 uint32_t UncompressHUFFMAN16(char* dest, 364 const char* src, 365 uint32_t numChunks) 366 { 328 367 vector<uint16_t> uncompressed; 329 368 330 369 //read compressed sizes (one per row) 331 370 const uint32_t* compressedSizes = reinterpret_cast<const uint32_t*>(src); 332 src += sizeof(uint32_t)*num RowElems;371 src += sizeof(uint32_t)*numChunks; 333 372 334 373 //uncompress the rows, one by one 335 374 uint32_t sizeWritten = 0; 336 for (uint32_t j=0;j<num RowElems;j++)375 for (uint32_t j=0;j<numChunks;j++) 337 376 { 338 377 Huffman::Decode(reinterpret_cast<const unsigned char*>(src), compressedSizes[j], uncompressed); … … 347 386 } 348 387 349 //Read a bunch of data compressed with the smoothman algorithm 350 uint32_t UncompressSMOOTHMAN(int16_t* dest, 351 const char* src, 352 uint32_t numRows, 353 uint32_t sizeOfElems, 354 uint32_t numRowElems) 355 { 356 //call huffman transposed 357 const uint32_t sizeWritten = UncompressHUFFMAN(reinterpret_cast<char*>(dest), src, numRowElems, sizeOfElems, numRows); 358 388 uint32_t UnApplySMOOTHING(int16_t* data, 389 uint32_t numElems) 390 { 359 391 //un-do the integer smoothing 360 for (uint32_t j=2;j<numRowElems*numRows;j++) 361 dest[j] = dest[j] + (dest[j-1]+dest[j-2])/2; 362 363 return sizeWritten; 364 } 365 392 for (uint32_t j=2;j<numElems;j++) 393 data[j] = data[j] + (data[j-1]+data[j-2])/2; 394 395 return numElems*sizeof(uint16_t); 396 } 366 397 // Data has been read from disk. Uncompress it ! 367 398 void UncompressBuffer(const uint32_t &catalogCurrentRow, const uint32_t &thisRoundNumRows) … … 377 408 378 409 //get the compression flag 379 const int64_t compressedOffset = fTileOffsets[catalogCurrentRow][i];//fCatalog[catalogCurrentRow][i].second - fCatalog[catalogCurrentRow][0].second; 380 const char compressedFlag = fCompressedBuffer[compressedOffset]; 381 382 //#define COMPRESSED_FLAG 0x1 383 //#define UNCOMPRESSED_FLAG 0x0 384 385 const char *src = fCompressedBuffer.data()+compressedOffset+1; 386 387 //if this bunch of data is not compressed, modify the compression flag 388 const uint32_t compression = compressedFlag==0 ? UNCOMPRESSED : col.comp; 389 switch (compression) 390 { 391 case UNCOMPRESSED: 392 dest += UncompressUNCOMPRESSED(dest, src, thisRoundNumRows, col.size, col.num); 393 break; 394 395 case SMOOTHMAN: 396 dest += UncompressSMOOTHMAN(reinterpret_cast<int16_t*>(dest), src, thisRoundNumRows, col.size, col.num); 397 break; 398 399 default: 400 ; 410 const int64_t compressedOffset = fTileOffsets[catalogCurrentRow][i]; 411 412 BlockHeader* head = reinterpret_cast<BlockHeader*>(&fCompressedBuffer[compressedOffset]); 413 414 fColumnOrdering[i] = head->ordering; 415 416 uint32_t numRows = (head->ordering==FACT_ROW_MAJOR) ? thisRoundNumRows : col.num; 417 uint32_t numCols = (head->ordering==FACT_COL_MAJOR) ? thisRoundNumRows : col.num; 418 419 const char *src = fCompressedBuffer.data()+compressedOffset+sizeof(BlockHeader)+sizeof(uint16_t)*head->numProcs; 420 421 for (uint32_t j=head->numProcs;j != 0; j--) 422 { 423 uint32_t sizeWritten=0; 424 425 switch (head->processings[j-1]) 426 { 427 case FACT_RAW: 428 if (head->numProcs == 1) 429 sizeWritten = UncompressUNCOMPRESSED(dest, src, numRows*numCols, col.size); 430 break; 431 case FACT_SMOOTHING: 432 sizeWritten = UnApplySMOOTHING(reinterpret_cast<int16_t*>(dest), numRows*numCols); 433 break; 434 case FACT_HUFFMAN16: 435 sizeWritten = UncompressHUFFMAN16(dest, src, numRows); 436 break; 437 default: 438 #ifdef __EXCEPTIONS 439 throw runtime_error("Unknown processing applied to data. Aborting"); 440 #else 441 gLog << ___err___ << "ERROR - Unknown processing applied to data. Aborting" << endl; 442 return; 443 #endif 444 break; 445 } 446 //increment destination counter only when processing done. 447 if (j==1) dest+= sizeWritten; 401 448 } 402 449 }
Note:
See TracChangeset
for help on using the changeset viewer.