int Decompress	(	const string &	ifile,
		const string &	ofile
	)

Definition at line 246 of file zfits.cc.

References Huffman::Decode(), fits::GetColumns(), fits::GetUInt(), size, start(), and Time::UnixTime().

Referenced by main().

{
    // open a fits file
    ifstream fin(ifile);

    // open output file
    ofstream fout(ofile);

    // get and check magic number
    unsigned char m[2];
    fin.read((char*)m, 2);
    if (m[0]!='z'+128 || m[1]!='f'+128)
        throw runtime_error("File not a compressed fits file.");

    // get length of additional header information
    size_t hlen = 0;
    fin.read((char*)&hlen, sizeof(size_t));
    if (hlen>0)
        throw runtime_error("Only Version-zero files supported.");

    // get size of FITS header
    size_t hs = 0;
    fin.read((char*)&hs, sizeof(size_t));
    if (!fin)
        throw runtime_error("Could not access header size.");

    // copy the header from the input to the output file
    // and prefix the output file as a compressed fits file
    string header;
    header.resize(hs);

    fin.read((char*)header.c_str(), header.size());
    fout.write((char*)header.c_str(), header.size());
    if (!fin)
        throw runtime_error("Could not read full header");

    string templ("tmpXXXXXX");
    int fd = mkstemp((char*)templ.c_str());
    const ssize_t rc = write(fd, header.c_str(), header.size());
    close(fd);

    if (rc<0)
        throw runtime_error("Could not write to temporary file: "+string(strerror(errno)));

    // open the output file to get the header parsed
    fits info(templ);

    remove(templ.c_str());

    // get the maximum size of one row and a list
    // of all columns ordered by their offset
    size_t row_tot = 0;
    const fits::Table::Columns &cols = info.GetColumns();
    map<size_t, fits::Table::Column> columns;
    for (auto it=cols.begin(); it!=cols.end(); it++)
    {
        columns[it->second.offset] = it->second;
        row_tot += it->second.num*it->second.size;
    }

    // very simple timer
    double sec = 0;
    double frac = 0;

    size_t com = 2+hs+2*sizeof(size_t);
    size_t tot = hs;

    const size_t masklen = cols.size()/8+1;

    // loop over all rows
    vector<char> buf(row_tot+masklen);
    vector<char> swap(row_tot);
    uint32_t offset = 0;

    const uint64_t nrows = info.GetUInt("NAXIS2");

    Time start;

    cout << fixed;
    for (uint32_t irow=0; irow<nrows; irow++)
    {
        fin.read(buf.data()+offset, buf.size()-offset);

        const uint8_t *mask = reinterpret_cast<uint8_t*>(buf.data());
        offset = masklen;

        char *ptr = swap.data();

        uint32_t icol = 0;
        for (auto it=columns.begin(); it!= columns.end(); it++, icol++)
        {
            const size_t &num  = it->second.num;
            const size_t &size = it->second.size;

            if (mask[icol/8]&(1<<(icol%8)))
            {
                Time now;

                vector<uint16_t> out(num*size/2);
                int len = Huffman::Decode((uint8_t*)buf.data()+offset, buf.size()-offset, out);
                if (len<0)
                    throw runtime_error("Decoding failed.");

                sec += Time().UnixTime()-now.UnixTime();

                offset += len;

                revcpy<2>(ptr, (char*)out.data(), num);
            }
            else
            {
                switch (size)
                {
                case 1: memcpy   (ptr, buf.data()+offset, num*size); break;
                case 2: revcpy<2>(ptr, buf.data()+offset, num);      break;
                case 4: revcpy<4>(ptr, buf.data()+offset, num);      break;
                case 8: revcpy<8>(ptr, buf.data()+offset, num);      break;
                }

                offset += num*size;
            }

            ptr += num*size;
        }

        com += offset+masklen;
        tot += row_tot;

        fout.write((char*)swap.data(), swap.size());

        memmove(buf.data(), buf.data()+offset, buf.size()-offset);
        offset = buf.size()-offset;

        if (!fout)
            throw runtime_error("Error writing to output file");

        const float proc = float(irow)/nrows;
        if (proc>frac)
        {
            const double elep = Time().UnixTime()-start.UnixTime();
            cout << "\r" << setprecision(0) << setw(3) << 100*proc << "% [" << setprecision(1) << setw(5) << 100.*com/tot << "%] cpu:" << sec << "s out:" << tot/1000000/elep << "MB/s" << flush;
            frac += 0.01;
        }
    }

    const double elep = Time().UnixTime()-start.UnixTime();
    cout << setprecision(0) << "\r100% [" << setprecision(1) << setw(5) << 100.*com/tot << "%] cpu:" << sec << "s out:" << tot/1000000/elep << "MB/s" << endl;

    return 0;
}

Here is the call graph for this function:

Here is the caller graph for this function: