Sunday, 23 November 2014

Decompressing the Image

The image decompression is actually pretty routine; it's also fairly specific to this data set.

At the core of the decompressor is a Huffman Coding engine, and decompression has three stages: Extract the compression histogram, initialise the decompressor, and then unpack the lines.

Taking these in turn....

Extract the histogram

This is the compression histogram we talked about last time. It's a set of 511 integers, each of which is a 32 bit frequency count value, corresponding to the frequency of the given offset value.

A simple extraction routine to turn a QByteArray into a histogram would be something like:

QList<QByteArray> data
...
QList<uint32_t> out;
QByteArray p;
...
  p = data.at(i);
  for (int j=0; j < p.size(); j+= 4)
  {
  char b[4];
    b[0] = p.at(j);
    b[1] = p.at(j+1);
    b[2] = p.at(j+2);
    b[3] = p.at(j+3);
    out.append(
              ((unsigned char)b[3] << 24) |
              ((unsigned char)b[2] << 16) |
              ((unsigned char)b[1] << 8) |
              ((unsigned char)b[0]));
  }
i.e. - take the byte arrays that correspond to the records, and for each one then extract the 32 bit value and append it to a new output array.

At the end of this process we have a new QByteArray containing all the records from the histogram. Since we know the offset in the file from the encoding record label, and we know the size is always three records we can do something along the lines of:

QList<uint32_t> _enchistogram;
...
ptr = raw->GetLabel("^ENCODING_HISTOGRAM");
int rn = ptr.toInt() -1; // file count isn't zero based...
...
QList<QByteArray> h;
  h.append(_data[rn++]);
  h.append(_data[rn++]);
  h.append(_data[rn++]);
 _enchistogram = loadhistogram(h);
And our histogram is a QList of integers.

Initialise the decompressor

This comes down to building the binary tree we will use to decode the bitstream in the compressed file.

We simply create a set of Huffman Tree Nodes, then sort them into order, and then go through combining the two lowest frequency nodes into a new node, and then re-sorting the list until we have only one node left.

This is literally the algorithm in this rosetta code article, and also described on this geeks for geeks page.

So we define a simple Huffman node 
class HNode
{
public:
...
  class HNode* right;
  class HNode* left;
  int value;
  int frq;
};

Where the "value" field is the leaf token, and "freq" is a frequency count. The "value" of -1 indicates an internal node.

So we can build the initial node set from the input histogram
for (int i=0; i< list.size(); i++)
  {
  HNode* newNode;
    newNode = new(HNode);
    newNode->frq = list[i];
    newNode->value = 1+i;
    nd.append(newNode);
  }

And then the sorting. We use the Qt qStableSort() function (although this is deprecated, and we should be using std::stable_sort). It's important to use this specific sort, since the way we deal with "equivalent" frequency nodes is important in building the tree, and using an alternative sort will mess that; to quite the Qt docs the property we care about is:
 The item that appeared before the other in the original container will still appear first after the sort.
So we do this:
HNode* _tree;
...
qStableSort(nd.begin(), nd.end(), cmpfreq);
while (nd.size() >1)
{
HNode* src;
  src = new HNode;
  src->right = nd.takeFirst();
  src->left = nd.takeFirst();
  src->frq = src->left->frq + src->right->frq;
  nd.insert(0,src);
  qStableSort(nd.begin(), nd.end(), cmpfreq);
}
_tree = nd.takeFirst();
where the comparison function is simply:

bool cmpfreq(const HNode* p1, const HNode* p2)
{
    return p1->frq < p2->frq;
}
And at the end of this the _tree variable holds the root of the binary tree we can use for the decode.

Unpacking the lines

We can find the location of the line records with the relevant pointer record:
 ptr = raw->GetLabel("^IMAGE");
And then we can take the next 800 records, and unpack them line by line into a new byte array.

Each line unpack is:
  • The first byte gives us the initial pixel value on the line
  • The next bytes are (Huffman coded) offsets from the "previous pixel value"

So we take the first value and then add the coded offset value to each entry as we go along the line.

The unpacking is a simple tree walk, 
QByteArray output;
HNode* tr;
char current;
int pos;

...
  current = in.at(0);
  output[0] = current;
  pos = 1;
  for (int i=1; i< in.size(); i++)
  {
  unsigned char v = in.at(i);
    for (int p=0; p<8; p++)
    {// Walk this byte, updating our tree position
      unsigned char m = 0x80 >> p;
      if (m&v)
        tr = tr->left;
      else
        tr = tr->right;
      if (tr->value != -1)
      {// Hit a leaf - use the value
        current -= tr->value + 256;
        output[pos++] =current;
        tr = _tree;
      }
    }
  }

This is a fairly simple walk, so I won't go into any more detail right now...

Turning it into an image

The simplest way is we can use the QImage() class. This will give us a simple way of writing out a variety of image formats. We can simply create this with
  _image = QImage(width, height, QImage::Format_RGB32);
Then set the image values from the incoming byte array with:
  for (int y=0; y < height; y++)
    for (int x=0; x < height; x++)
    {
    QRgb value;
    unsigned char v;
      v = (unsigned char) im[pos];
      value = qRgb(v,v,v);
      _image.setPixel(x,y, value);
     pos++
    }

(where "im" is the image data from the lines we extracted in the line decompression). We can then just write this out with:
_image.save("out.png");

And the result is:

From Nasa's C4394416.IMQ. See Nasa for copyright/usage information.
This is a Voyager 2 image from the wide angle camera, taken with the orange filter on 24/08/1981 at 02:16:47. 800x800 pixels at 8bpp. And it's pretty clearly Saturn. "If it's not right it'll do until right comes along".
 Interestingly the  original Voyager docs state
on a VAX-750 it takes about one minute of CPU time to decompress an image.
Our code is very inefficient, and has massive room for improvement (all those flexible byte array copies and data duplication). However on my laptop it's about 200-250mS per image. In the "minute of CPU time" that the VAX decompression took, we can get over 250 Images from compressed to saved as a .png on disk. Isn't progress great?

The grid of dots is the reseau grid; it's described on wikipedia, but in summary it's a grid of black dots in known positions placed in the image capture path on the probe. Since we know where these dots "should" be, we can correct for image distortion by post-processing the image to put the dots "back" into the right positions.

There's are also a couple of stuck pixels and distortions from the original capture, some of which persist across images.

Next Steps

Rather than deal with the image undistortion for now we'll leave these Voyager images for a bit, and take a look at the "Vicar" format data records which are used by Cassini and calibrated Voyager images. We can then come back and take a closer look at what we have to do to get undistorted images, and also look at some PDS3 files from ESA's Venus Express.

So stay tuned for something on that next time....