Changeset 20036 for trunk/FACT++
- Timestamp:
- 01/13/21 19:49:18 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/FACT++/src/fixfits.cc
r20035 r20036 36 36 { 37 37 ostringstream str; 38 str << setfill(' ');39 38 str << value; 40 39 … … 42 41 str.str(""); 43 42 44 str << setw(8-keyword.size()) << keyword; 43 str << setfill(' ') << left; 44 str << setw(8) << keyword; 45 45 46 46 str << "= "; 47 47 48 str << setw(20 -num_length) << value << " / " << comment;48 str << setw(20) << value << " / " << comment; 49 49 50 50 str << setw(80-str.str().size()) << ' '; 51 51 52 cout << '\n'<< str.str() << '\n';52 cout << keyword.size() << ":" << num_length << ":" << str.str().size() << ":" << str.str() << '\n'; 53 53 54 54 if (str.str().size() != 80) … … 68 68 control.add_options() 69 69 ("file", var<string>()->required(), "Input file name. Fits file to be repaired.") 70 ("out,o", var<string>(), "Output file name (+'.repaired' if empty)") 70 ("out,o", var<string>(), "Output file name (+'.repaired' if empty), dry-run if not provided") 71 ("no-progress", po_switch(), "Suppress progress output") 72 ("dry-run", po_switch(), "Do not write any output") 71 73 ; 72 74 … … 105 107 const string input_name = conf.Get<string>("file"); 106 108 const string output_name = conf.Has("out") ? conf.Get<string>("out") : 107 (fs::path(input_name)/".repaired").string(); 108 109 cout << "Reading: " << input_name << endl; 109 (input_name+".repaired"); 110 111 const bool no_progress = conf.Get<bool>("no-progress"); 112 const bool dry_run = conf.Get<bool>("dry-run"); 113 114 cout << "\nReading: " << input_name << endl; 110 115 111 116 ifstream input(input_name.c_str(), ios_base::binary); 112 113 cout << "Looking for beginning of data table." << endl; 117 if (!input) 118 { 119 cerr << "Error opening file: " << strerror(errno) << endl; 120 return 6; 121 } 122 123 cout << "Seeking beginning of data table." << endl; 114 124 115 125 // ------------------------------------------------------------ … … 137 147 num_bytes+=80; 138 148 139 cout << "\rByte " << num_bytes << " "; 140 } 149 if (!no_progress) 150 cout << "\rByte " << num_bytes << " "; 151 } 152 if (no_progress) 153 cout << "Reached " << num_bytes << '\n'; 141 154 cout << "\nFound beginning of data table." << endl; 142 155 … … 155 168 num_bytes+=80; 156 169 157 cout << "\rByte " << num_bytes << " "; 158 } 170 if (!no_progress) 171 cout << "\rByte " << num_bytes << " "; 172 } 173 if (no_progress) 174 cout << "Reached " << num_bytes << '\n'; 159 175 160 176 cout << "\nReached end of data table header" << endl; … … 168 184 num_bytes+=80; 169 185 170 cout << "\rByte " << num_bytes << " "; 171 } 186 if (!no_progress) 187 cout << "\rByte " << num_bytes << " "; 188 } 189 if (no_progress) 190 cout << "Reached " << num_bytes << '\n'; 172 191 173 192 … … 180 199 input.read(buffer, 4); 181 200 num_bytes += 4; 182 cout << "\rByte " << num_bytes << " "; 183 } 201 202 if (!no_progress) 203 cout << "\rByte " << num_bytes << " "; 204 } 205 if (no_progress) 206 cout << "Reached " << num_bytes << '\n'; 184 207 185 208 input.seekg(num_bytes-4); 186 cout << "\nFound start of very first compressed tile" << endl;209 cout << "\nFound start of first compressed tile," << endl; 187 210 188 211 // ------------------------------------------------------------ … … 196 219 BlockHeader block_head; 197 220 198 input.read( (char*)(&tile_head), sizeof(TileHeader));199 200 cout << "First tile is " << tile_head.size << " bytes long " << endl;221 input.read(reinterpret_cast<char*>(&tile_head), sizeof(TileHeader)); 222 223 cout << "First tile is " << tile_head.size << " bytes long." << endl; 201 224 202 225 // ------------------------------------------------------------ … … 205 228 size_t num_blocks = 0; 206 229 207 while ( input.tellg() < byte_where_tile_start+tile_head.size)208 { 209 input.read( (char*)(&block_head), sizeof(BlockHeader));230 while (size_t(input.tellg()) < byte_where_tile_start+tile_head.size) 231 { 232 input.read(reinterpret_cast<char*>(&block_head), sizeof(BlockHeader)); 210 233 block_bytes += block_head.size; 211 234 … … 220 243 } 221 244 222 if ( input.tellg() != byte_where_tile_start+tile_head.size)245 if (size_t(input.tellg()) != byte_where_tile_start+tile_head.size) 223 246 { 224 247 cout << "Discrepancy between tile size and sum of blocks size." << endl; … … 226 249 } 227 250 228 cout << "There are " << num_blocks << " blocks in each tile. Now looping over all tiles" << endl;251 cout << "There are " << num_blocks << " blocks in each tile.\nLooping over all tiles" << endl; 229 252 230 253 // ------------------------------------------------------------ … … 253 276 } 254 277 255 input.read( (char*)(&tile_head), sizeof(TileHeader));278 input.read(reinterpret_cast<char*>(&tile_head), sizeof(TileHeader)); 256 279 if (memcmp(tile_head.id, "TILE", 4) || input.eof()) 257 280 { … … 263 286 catalog.emplace_back(); 264 287 265 for ( int i=0;i<num_blocks;i++)288 for (size_t i=0; i<num_blocks; i++) 266 289 { 267 290 input.read((char*)(&block_head), sizeof(BlockHeader)); … … 285 308 catalog.back().emplace_back(0,0); 286 309 287 cout << "\rByte " << input.tellg() << " "; 310 if (!no_progress) 311 cout << "\rByte " << input.tellg() << " "; 288 312 num_tiles++; 289 313 } … … 304 328 305 329 //Now we have everything that we need. Go again through the input file and write the output, fixed file 306 cout << " Data recovered.\n Writing output file: " << output_name << endl;330 cout << "Recovery finished!\n\nWriting: " << output_name << endl; 307 331 308 332 // ============================================================ 309 333 310 ofstream output(output_name.c_str(), ios_base::binary); 311 input.seekg(0); 312 313 //copy everything until the catalog start and update header values on our way 314 counter = 0; 315 for (uint64_t i=0; i<output_catalog_start; i+=80) 316 { 317 input.read(buffer, 80); 318 319 if (!strcmp(buffer, "XTENSION= 'BINTABLE' / binary table extension ")) 320 counter++; 321 322 if (counter == 2) 323 { 324 bool rc = true; 325 if (!memcmp(buffer, "NAXIS2", 6)) 326 rc &= update_header_key(buffer, "NAXIS2", num_tiles, "number of rows in table"); 327 if (!memcmp(buffer, "ZNAXIS2", 7)) 328 rc &= update_header_key(buffer, "ZNAXIS2", num_tiles, "Number of uncompressed rows"); 329 if (!memcmp(buffer, "PCOUNT", 6)) 330 rc &= update_header_key(buffer, "PCOUNT", output_heap_end - output_catalog_end, "size of special data area"); 331 if (!memcmp(buffer, "THEAP", 5)) 332 rc &= update_header_key(buffer, "THEAP", theap, "Pointer to heap"); 333 if (!memcmp(buffer, "ZHEAPPTR", 8)) 334 rc &= update_header_key(buffer, "ZHEAPPTR", zheapptr, "Pointer to data"); 335 336 if (!rc) 337 return 5; 338 } 339 340 output.write(buffer, 80); 341 342 cout << "\rByte " << input.tellg() << " "; 343 } 344 345 // ------------------------------------------------------------ 346 347 // now write the updated catalog 348 const uint64_t reserved_num_tiles = zheapptr / (num_cols*16); 349 const uint32_t one_catalog_row_size = num_cols*16; 350 const uint32_t total_catalog_size = reserved_num_tiles*one_catalog_row_size; 351 352 vector<char> swapped_catalog(total_catalog_size); 353 354 uint32_t shift = 0; 355 for (auto it=catalog.cbegin(); it!=catalog.cend(); it++) 356 { 357 revcpy<sizeof(int64_t)>(swapped_catalog.data() + shift, (char*)(it->data()), num_cols*2); 358 shift += one_catalog_row_size; 359 } 360 361 if (num_tiles < reserved_num_tiles) 362 memset(swapped_catalog.data()+shift, 0, total_catalog_size - shift); 363 364 output.write(swapped_catalog.data(), total_catalog_size); 365 366 // ------------------------------------------------------------ 367 368 // we are now at the beginning of the heap in the output file. Move there in the input file 369 input.seekg(output_heap_start); 370 371 uint64_t num_heap_bytes_written = 0; 372 uint64_t total_heap_bytes_to_write = output_heap_end - output_heap_start; 373 374 cout << endl << "Catalog updated. Copying " << total_heap_bytes_to_write << " of data now." << endl; 375 376 while (num_heap_bytes_written != total_heap_bytes_to_write) 377 { 378 uint64_t bytes_to_write = 2048; 379 if (total_heap_bytes_to_write - num_heap_bytes_written < 2048) 380 bytes_to_write = total_heap_bytes_to_write - num_heap_bytes_written; 381 382 input.read(buffer, bytes_to_write); 383 output.write(buffer, bytes_to_write); 384 num_heap_bytes_written += bytes_to_write; 385 cout << "\rByte " << input.tellg() << " "; 386 } 387 388 cout << "\nData written. Adding FITS padding. " << endl; 389 390 // ------------------------------------------------------------ 391 392 //now add the FITS padding at the end 393 while (output.tellp() % 2880 != 0) 394 { 395 output.write("\0", 1); 396 } 334 if (!dry_run) 335 { 336 ofstream output(output_name.c_str(), ios_base::binary); 337 input.seekg(0); 338 339 //copy everything until the catalog start and update header values on our way 340 counter = 0; 341 for (uint64_t i=0; i<output_catalog_start; i+=80) 342 { 343 input.read(buffer, 80); 344 345 if (!strcmp(buffer, "XTENSION= 'BINTABLE' / binary table extension ")) 346 counter++; 347 348 if (counter == 2) 349 { 350 bool rc = true; 351 352 if (!memcmp(buffer, "NAXIS2", 6)) 353 rc &= update_header_key(buffer, "NAXIS2", num_tiles, "number of rows in table"); 354 if (!memcmp(buffer, "ZNAXIS2", 7)) 355 rc &= update_header_key(buffer, "ZNAXIS2", num_tiles, "Number of uncompressed rows"); 356 if (!memcmp(buffer, "PCOUNT", 6)) 357 rc &= update_header_key(buffer, "PCOUNT", output_heap_end - output_catalog_end, "size of special data area"); 358 if (!memcmp(buffer, "THEAP", 5)) 359 rc &= update_header_key(buffer, "THEAP", theap, "Pointer to heap"); 360 if (!memcmp(buffer, "ZHEAPPTR", 8)) 361 rc &= update_header_key(buffer, "ZHEAPPTR", zheapptr, "Pointer to data"); 362 363 if (!rc) 364 return 5; 365 } 366 367 output.write(buffer, 80); 368 369 if (!no_progress) 370 cout << "\rByte " << input.tellg() << " "; 371 } 372 373 if (no_progress) 374 cout << "Reached " << input.tellg() << '\n'; 375 376 // ------------------------------------------------------------ 377 378 // now write the updated catalog 379 const uint64_t reserved_num_tiles = zheapptr / (num_cols*16); 380 const uint32_t one_catalog_row_size = num_cols*16; 381 const uint32_t total_catalog_size = reserved_num_tiles*one_catalog_row_size; 382 383 vector<char> swapped_catalog(total_catalog_size); 384 385 uint32_t shift = 0; 386 for (auto it=catalog.cbegin(); it!=catalog.cend(); it++) 387 { 388 revcpy<sizeof(int64_t)>(swapped_catalog.data() + shift, (char*)(it->data()), num_cols*2); 389 shift += one_catalog_row_size; 390 } 391 392 if (num_tiles < reserved_num_tiles) 393 memset(swapped_catalog.data()+shift, 0, total_catalog_size - shift); 394 395 output.write(swapped_catalog.data(), total_catalog_size); 396 397 // ------------------------------------------------------------ 398 399 // we are now at the beginning of the heap in the output file. Move there in the input file 400 input.seekg(output_heap_start); 401 402 uint64_t num_heap_bytes_written = 0; 403 uint64_t total_heap_bytes_to_write = output_heap_end - output_heap_start; 404 405 cout << "\nCatalog updated.\nCopying " << total_heap_bytes_to_write << " of data now." << endl; 406 407 while (num_heap_bytes_written != total_heap_bytes_to_write) 408 { 409 uint64_t bytes_to_write = 2048; 410 if (total_heap_bytes_to_write - num_heap_bytes_written < 2048) 411 bytes_to_write = total_heap_bytes_to_write - num_heap_bytes_written; 412 413 input.read(buffer, bytes_to_write); 414 output.write(buffer, bytes_to_write); 415 num_heap_bytes_written += bytes_to_write; 416 417 if (!no_progress) 418 cout << "\rByte " << input.tellg() << " "; 419 } 420 if (no_progress) 421 cout << "Reached " << input.tellg() << '\n'; 422 423 cout << "\nData written. Adding FITS padding. " << endl; 424 425 // ------------------------------------------------------------ 426 427 //now add the FITS padding at the end 428 while (output.tellp() % 2880 != 0) 429 output.write("\0", 1); 430 } 431 else 432 cout << " ... skipped ...\n"; 433 397 434 cout << "File recovered!\n"; 398 435 399 input.close(); 400 output.close(); 401 402 // ------------------------------------------------------------ 403 404 cout << "\nValid num tiles: " << num_tiles; 405 cout << "\nCatalog start : " << output_catalog_start; 406 cout << "\nCatalog end : " << output_catalog_end; 407 cout << "\nHeap start : " << output_heap_start; 408 cout << "\nHeap end : " << output_heap_end; 409 cout << "\nZHEAPPTR : " << zheapptr; 410 cout << "\nTHEAP : " << theap; 411 cout << "\nPCOUNT : " << output_heap_end - output_catalog_end; 412 cout << "\nNAXIS2 : " << num_tiles; 413 cout << "\nZNAXIS2 : " << num_tiles; 414 436 // ------------------------------------------------------------ 437 438 cout << "\nValid num tiles : " << setw(12) << num_tiles; 439 cout << "\nNAXIS2 : " << setw(12) << num_tiles; 440 cout << "\nZNAXIS2 : " << setw(12) << num_tiles; 441 cout << "\nCatalog start : " << setw(12) << output_catalog_start; 442 cout << "\nCatalog end : " << setw(12) << output_catalog_end; 443 cout << "\nZHEAPPTR : " << setw(12) << zheapptr; 444 cout << "\nTHEAP : " << setw(12) << theap; 445 cout << "\nHeap start : " << setw(12) << output_heap_start; 446 cout << "\nHeap end : " << setw(12) << output_heap_end; 447 cout << "\nPCOUNT : " << setw(12) << output_heap_end - output_catalog_end; 448 cout << '\n'; 449 cout << "\nInput file size : " << setw(12) << fs::file_size(input_name); 450 if (!dry_run) 451 { 452 cout << "\nOutput file size : " << setw(12) << fs::file_size(output_name); 453 cout << "\nLost bytes : " << setw(12) << fs::file_size(input_name)-fs::file_size(output_name); 454 cout << "\nLost fraction : " << setw(12) << 1.-fs::file_size(output_name)/fs::file_size(input_name); 455 } 415 456 cout << "\n\nYou should now:\n"; 416 457 cout << "- update TSTOPI and TSTOPF using \"fitsdump -c UnixTimeUTC --minmax --nozero\" and then fv\n";
Note:
See TracChangeset
for help on using the changeset viewer.