Index: /trunk/FACT++/src/fixfits.cc
===================================================================
--- /trunk/FACT++/src/fixfits.cc	(revision 20035)
+++ /trunk/FACT++/src/fixfits.cc	(revision 20035)
@@ -0,0 +1,422 @@
+//Compile with g++ -o FixFits fixFits.cpp -O2 -std=c++11
+//This is a "FITS file fixer" program for FACT.
+//It goes through the file, relying only on basic header data and compressed blocks headers
+//It will reconstruct the catalog data too in the output file a strip all corrupt data
+//The actual header keywords must be edited with the FTOOLS afterwards
+
+#include <string>
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <cstring>
+#include <vector>
+#include <algorithm>
+#include <sstream>
+
+#include <boost/filesystem.hpp>
+
+#include "FITS.h"
+#include "Time.h"
+#include "Configuration.h"
+
+using namespace std;
+using namespace FITS;
+namespace fs = boost::filesystem;
+
+
+template<size_t N>
+    void revcpy(char *dest, const char *src, int num)
+{
+    const char *pend = src + num*N;
+    for (const char *ptr = src; ptr<pend; ptr+=N, dest+=N)
+        std::reverse_copy(ptr, ptr+N, dest);
+}
+
+bool update_header_key(char buffer[], const string& keyword, uint64_t value, const string& comment)
+{
+    ostringstream str;
+    str << setfill(' ');
+    str << value;
+
+    const auto num_length = str.str().size();
+    str.str("");
+
+    str << setw(8-keyword.size()) << keyword;
+
+    str << "= ";
+
+    str << setw(20-num_length) << value << " / " << comment;
+
+    str << setw(80-str.str().size()) << ' ';
+
+    cout << '\n' << str.str() << '\n';
+
+    if (str.str().size() != 80)
+    {
+        cerr << "Wrong " << keyword << " string size." << endl;
+        return false;
+    }
+
+    memcpy(buffer, str.str().c_str(), 80);
+
+    return true;
+}
+
+void SetupConfiguration(Configuration &conf)
+{
+    po::options_description control("Root to SQL");
+    control.add_options()
+        ("file",  var<string>()->required(),  "Input file name. Fits file to be repaired.")
+        ("out,o", var<string>(),              "Output file name (+'.repaired' if empty)")
+        ;
+
+    po::positional_options_description p;
+    p.add("file", 1); // All positional options
+    p.add("out",  1); // All positional options
+
+    conf.AddOptions(control);
+    conf.SetArgumentPositions(p);
+}
+
+void PrintUsage()
+{
+    cout <<
+        "fixfits - Fix the header of a none closed raw-data file\n"
+        "\n"
+        "Usage: fixfits input-file-name [output-file-name]\n"
+        "\n"
+        ;
+    cout << endl;
+}
+
+
+int main(int argc, const char* argv[])
+{
+    Time start;
+
+    Configuration conf(argv[0]);
+    conf.SetPrintUsage(PrintUsage);
+    SetupConfiguration(conf);
+
+    if (!conf.DoParse(argc, argv))
+        return 127;
+
+    // ----------------------------- Evaluate options --------------------------
+    const string input_name  = conf.Get<string>("file");
+    const string output_name = conf.Has("out") ? conf.Get<string>("out") :
+        (fs::path(input_name)/".repaired").string();
+
+    cout << "Reading: " << input_name << endl;
+
+    ifstream input(input_name.c_str(), ios_base::binary);
+
+    cout << "Looking for beginning of data table." << endl;
+
+    // ------------------------------------------------------------
+
+    //first off, find the beginning of the data table
+    char buffer[2048];
+    buffer[80] = '\0';
+
+    size_t counter = 0;
+    size_t num_bytes = 80;
+    input.read(buffer, 80);
+
+    while (counter != 2)
+    {
+        if (!strcmp(buffer, "XTENSION= 'BINTABLE'           / binary table extension                         "))
+            counter++;
+
+        if (input.eof())
+        {
+            cout << "Reached end-of-file but could not find start of data table." << endl;
+            return 2;
+        }
+
+        input.read(buffer, 80);
+        num_bytes+=80;
+
+        cout << "\rByte " << num_bytes << "      ";
+    }
+    cout << "\nFound beginning of data table." << endl;
+
+    // ------------------------------------------------------------
+
+    //we just read the start of the data table. continue until the end of the table header
+    while (strcmp(buffer, "END                                                                             "))
+    {
+        if (input.eof())
+        {
+            cout << "\nData table header was not closed properly. Aborting." << endl;
+            return 3;
+        }
+
+        input.read(buffer, 80);
+        num_bytes+=80;
+
+        cout << "\rByte " << num_bytes << "       ";
+    }
+
+    cout << "\nReached end of data table header" << endl;
+
+    // ------------------------------------------------------------
+
+    //now skip the fits padding
+    while (num_bytes % 2880 != 0)
+    {
+        input.read(buffer, 80);
+        num_bytes+=80;
+
+        cout << "\rByte " << num_bytes << "       ";
+    }
+
+
+    const size_t output_catalog_start = input.tellg();
+
+    //now we are in the catalog. skip it until we find the string 'TILE'
+    buffer[4] = '\0';
+    while (strcmp(buffer, "TILE"))
+    {
+        input.read(buffer, 4);
+        num_bytes += 4;
+        cout << "\rByte " << num_bytes << "      ";
+    }
+
+    input.seekg(num_bytes-4);
+    cout << "\nFound start of very first compressed tile" << endl;
+
+    // ------------------------------------------------------------
+
+    const size_t output_heap_start = input.tellg();
+
+    size_t byte_where_tile_start = input.tellg();
+
+    //now figure out how many blocks are in each tile
+    TileHeader tile_head;
+    BlockHeader block_head;
+
+    input.read((char*)(&tile_head), sizeof(TileHeader));
+
+    cout << "First tile is " << tile_head.size << " bytes long" << endl;
+
+    // ------------------------------------------------------------
+
+    size_t block_bytes = 0;
+    size_t num_blocks = 0;
+
+    while (input.tellg() < byte_where_tile_start+tile_head.size)
+    {
+        input.read((char*)(&block_head), sizeof(BlockHeader));
+        block_bytes += block_head.size;
+
+        input.seekg(byte_where_tile_start + sizeof(TileHeader) + block_bytes);
+        num_blocks ++;
+
+        if (input.eof())
+        {
+            cout << "First tile seems broken unable to recover anything." << endl;
+            return 4;
+        }
+    }
+
+    if (input.tellg() != byte_where_tile_start+tile_head.size)
+    {
+        cout << "Discrepancy between tile size and sum of blocks size." << endl;
+        return 1;
+    }
+
+    cout << "There are " << num_blocks << " blocks in each tile. Now looping over all tiles" << endl;
+
+    // ------------------------------------------------------------
+
+    size_t num_tiles = 0;
+    input.seekg(output_heap_start);
+
+    size_t output_heap_end = input.tellg();
+    size_t previous_heap_end = input.tellg();
+
+    vector<vector<pair<int64_t, int64_t> > > catalog;
+
+    uint64_t offset_in_heap = 0;
+    while (!input.eof())
+    {
+        previous_heap_end = output_heap_end;
+        output_heap_end = input.tellg();
+
+        byte_where_tile_start = input.tellg();
+        block_bytes = 0;
+
+        if (input.peek() == EOF)
+        {
+            cout << "Reached end-of-file." << endl;
+            break;
+        }
+
+        input.read((char*)(&tile_head), sizeof(TileHeader));
+        if (memcmp(tile_head.id, "TILE", 4) || input.eof())
+        {
+            cout << "Reached end-of-table padding or end-of-file." << endl;
+            break;
+        }
+
+        offset_in_heap += sizeof(TileHeader);
+        catalog.emplace_back();
+
+        for (int i=0;i<num_blocks;i++)
+        {
+            input.read((char*)(&block_head), sizeof(BlockHeader));
+            if (input.eof())
+            {
+                cout << "Reached end-of-file." << endl;
+                break;
+            }
+            block_bytes += block_head.size;
+            input.seekg(byte_where_tile_start + sizeof(TileHeader) + block_bytes);
+            if (input.eof())
+            {
+                cout << "Reached end-of-file." << endl;
+                break;
+            }
+            catalog.back().emplace_back((int64_t)(block_head.size),offset_in_heap);
+            offset_in_heap += block_head.size;
+        }
+
+        //one last catalog entry for the column that has no data
+        catalog.back().emplace_back(0,0);
+
+        cout << "\rByte " << input.tellg() << "      ";
+        num_tiles++;
+    }
+
+    input.clear();
+
+    // ------------------------------------------------------------
+
+    //ignore last read tile as we have no idea if it looked ok or not
+    num_tiles--;
+    output_heap_end = previous_heap_end;
+
+    const uint64_t num_cols = num_blocks+1;//+1 because of not-used drs-time-marker column
+    const uint64_t output_catalog_end = output_catalog_start + num_tiles*num_cols*16;
+
+    const uint64_t zheapptr = output_heap_start - output_catalog_start;
+    const uint64_t theap    = output_catalog_end - output_catalog_start;
+
+    //Now we have everything that we need. Go again through the input file and write the output, fixed file
+    cout << "Data recovered.\n Writing output file: " << output_name << endl;
+
+    // ============================================================
+
+    ofstream output(output_name.c_str(), ios_base::binary);
+    input.seekg(0);
+
+    //copy everything until the catalog start and update header values on our way
+    counter = 0;
+    for (uint64_t i=0; i<output_catalog_start; i+=80)
+    {
+        input.read(buffer, 80);
+
+        if (!strcmp(buffer, "XTENSION= 'BINTABLE'           / binary table extension                         "))
+            counter++;
+
+        if (counter == 2)
+        {
+            bool rc = true;
+            if (!memcmp(buffer, "NAXIS2", 6))
+                rc &= update_header_key(buffer, "NAXIS2", num_tiles, "number of rows in table");
+            if (!memcmp(buffer, "ZNAXIS2", 7))
+                rc &= update_header_key(buffer, "ZNAXIS2", num_tiles, "Number of uncompressed rows");
+            if (!memcmp(buffer, "PCOUNT", 6))
+                rc &= update_header_key(buffer, "PCOUNT", output_heap_end - output_catalog_end, "size of special data area");
+            if (!memcmp(buffer, "THEAP", 5))
+                rc &= update_header_key(buffer, "THEAP", theap, "Pointer to heap");
+            if (!memcmp(buffer, "ZHEAPPTR", 8))
+                rc &= update_header_key(buffer, "ZHEAPPTR", zheapptr, "Pointer to data");
+
+            if (!rc)
+                return 5;
+        }
+
+        output.write(buffer, 80);
+
+        cout << "\rByte " << input.tellg() << "       ";
+    }
+
+    // ------------------------------------------------------------
+
+    // now write the updated catalog
+    const uint64_t reserved_num_tiles   = zheapptr / (num_cols*16);
+    const uint32_t one_catalog_row_size = num_cols*16;
+    const uint32_t total_catalog_size   = reserved_num_tiles*one_catalog_row_size;
+
+    vector<char> swapped_catalog(total_catalog_size);
+
+    uint32_t shift = 0;
+    for (auto it=catalog.cbegin(); it!=catalog.cend(); it++)
+    {
+        revcpy<sizeof(int64_t)>(swapped_catalog.data() + shift, (char*)(it->data()), num_cols*2);
+        shift += one_catalog_row_size;
+    }
+
+    if (num_tiles < reserved_num_tiles)
+        memset(swapped_catalog.data()+shift, 0, total_catalog_size - shift);
+
+    output.write(swapped_catalog.data(), total_catalog_size);
+
+    // ------------------------------------------------------------
+
+    // we are now at the beginning of the heap in the output file. Move there in the input file
+    input.seekg(output_heap_start);
+
+    uint64_t num_heap_bytes_written = 0;
+    uint64_t total_heap_bytes_to_write = output_heap_end - output_heap_start;
+
+    cout << endl << "Catalog updated. Copying " << total_heap_bytes_to_write << " of data now." << endl;
+
+    while (num_heap_bytes_written != total_heap_bytes_to_write)
+    {
+        uint64_t bytes_to_write = 2048;
+        if (total_heap_bytes_to_write - num_heap_bytes_written < 2048)
+            bytes_to_write = total_heap_bytes_to_write - num_heap_bytes_written;
+
+        input.read(buffer, bytes_to_write);
+        output.write(buffer, bytes_to_write);
+        num_heap_bytes_written += bytes_to_write;
+        cout << "\rByte " << input.tellg() << "           ";
+    }
+
+    cout << "\nData written. Adding FITS padding. " << endl;
+
+    // ------------------------------------------------------------
+
+    //now add the FITS padding at the end
+    while (output.tellp() % 2880 != 0)
+    {
+        output.write("\0", 1);
+    }
+    cout << "File recovered!\n";
+
+    input.close();
+    output.close();
+
+    // ------------------------------------------------------------
+
+    cout << "\nValid num tiles: " << num_tiles;
+    cout << "\nCatalog start  : " << output_catalog_start;
+    cout << "\nCatalog end    : " << output_catalog_end;
+    cout << "\nHeap start     : " << output_heap_start;
+    cout << "\nHeap end       : " << output_heap_end;
+    cout << "\nZHEAPPTR       : " << zheapptr;
+    cout << "\nTHEAP          : " << theap;
+    cout << "\nPCOUNT         : " << output_heap_end - output_catalog_end;
+    cout << "\nNAXIS2         : " << num_tiles;
+    cout << "\nZNAXIS2        : " << num_tiles;
+
+    cout << "\n\nYou should now:\n";
+    cout << "- update TSTOPI and TSTOPF using \"fitsdump -c UnixTimeUTC --minmax --nozero\" and then fv\n";
+    cout << "- update DATE-END using TSTOP given to the MjDtoISO program in ~fact_arc/ingestScripts and then fv again\n";
+    cout << "- update the checksum of the file with \"fchecksum update+ <filename>\"\n" << endl;
+
+    return 0;
+
+}