Thursday, 18 December 2014

Venus Express

Since this is the week that the Venus Express mission ended then let's look at the format of data from the Venus Monitoring Camera, the VMC.


Getting the data

The Venus Express data sets are collected at the ESA Project page, and the VMC sets are at the associated FTP directory.

Looking under the per-mission catalogues then the DATA/ directories are split into orbits, and image files are in there as .IMG files.


The Sensor


The VMC actually has a single 1024*1024 CCD sensor which has four independent optical channels over the sensor - see the Calibration Report (PDF) for details, but the end result is: we have four fixed channels available - Visible, UV and two IR (centred at .935 & 1.01 uM respectively).

The data file contains a tag "DETECTOR_ID" which is either VEX_VMC_VIS,  VEX_VMC_UV, VEX_VMC_NIR-1 or VEX_VMC_NIR-2 depending on the source. This is also used for the filename so;
  • V*_UV1.IMG is UV
  • V*_VI2.IMG is Visible
  • V*_UV21.IMG is UV #1
  • V*_UV22.IMG is UV #2

The Data

The data is actually the same as the Vicar record from Cassini, but with the addition of a PDS data header . To decode this with our existing code is fairly simple;
  • Check that the file starts with a "PDS_VERSION_ID" value of "PDS3"
  • Parse through a set of (line delimited) LABEL=VALUE records
  • When we hit a single line "END" then we're done
  • At this point we look for RECORD_BYTES and LABEL_RECORDS from the data we ingested
  • Skip "RECORD_BYTES * LABEL_RECORDS" bytes from the start of the file to move over the header
  • Load a Vicar/Cassini style image from that point

Image Level Correction

One minor gotcha is that we have a visible border from the imaging on these pictures - these leave a peak white "Hexagon" frame around the image. We can just ignore this during processing, however when we level correct then these values will swamp our actual image data. We could look at a more complex framing filter on the cv calls, but it's simpler for this case to implement a custom level min/max detect rather than use the stock cv::minMaxLoc() call, and have our version simply ignore these peak 16 bit values (65535) on the input; i.e.
static void getMinMax(cv::Mat img, double* min, double* max)
{
int minVal, maxVal;

  minVal = 65536;
  maxVal = 0;
  for(int i=0; i<img.rows; i++)
    for(int j=0; j<img.cols; j++)
    {
    uint16_t v;
        v = img.at<uint16_t>(i,j);
        if (v < minVal)
          minVal = v;
        else if ((v > maxVal) && (v < 65535))
          maxVal = v;
    }
  *min = (double)minVal;
  *max = (double)maxVal;
}


These are fairly simple tweaks and putting them together: From VMC Orbit 50, Image52, UV2



And we're (more or less) done.


Thursday, 11 December 2014

Automatic Feature Detection

We can use the feature detection algorithms in OpenCV to automatically align some of the image data; this isn't particularly effective on many pictures, but it can work in a few cases.

For this example we'll take three images of the same target taken through different filters and try to align them to produce a single full colour version:
  • W1484593934_1.IMG = Blue
  • W1484593810_1.IMG = Green
  • W1484593662_1.IMG = Red

These show the moon Rhea, and are from the coiss_2009 data set (and are actually the true colour filters).

We follow a standard decode and range optimisation, as for the previous example. However in this case we have to convert the images to 8 bit for the image alignment algorithms to operate on them. This is simply
im16d[i].convertTo(ch8[i], CV_8U, 0.0039);

Where im16d[i] is the source image and ch8[i] the destination. The 0.0039 factor is scaling the range back down from 16 bit to 8 bit values i.e. 256/65536. There's an argument we could have skipped this and just ranged to 8 bits, but this keeps our test code consistent if we always deal with a 16 bit backend.

Following this we align the images. OpenCV offers a standard workflow for this and a number of default algorithms to work with in this case. The basic set of steps is:
  1. Use a feature detector to get the "interesting" points in the image - i.e. features we can match between images. This gets us the KeyPoints List.
  2. Use a descriptor extractor to get a descriptor for each of the KeyPoints in the image, which gives us a Matrix of descriptors, one per keypoint.
  3.  Use a matcher to calculate the matches between our two sets of descriptors. This gives us a set of DMatch matches, containing points in one set and the nearest match in the other. These matches aren't necessarily all correct, just the best available.
  4. Cull the matches, to just use the "strongest" available matches. This gives us a better chance of getting a consistently correct relation between the two images.
  5. Use the findHomography() call to work out the transform between the points of the match.
  6. Apply this transform to the image to be aligned using warpPerspective()
Then we will have a set of aligned images which we can merge together into a single colour version.

For this first attempt we can use the stock feature detector, extractor and matcher components of OpenCV on a pair of images. Breaking the first 5 steps of the matching procedure into a single function to take two images and return the matrix from findHomography() we have:
static cv::Mat GetAlignment(cv::Mat im0, cv::Mat im1)
{
Then we have the three stage processors:
cv::OrbFeatureDetector detector;
cv::OrbDescriptorExtractor extractor;

cv::BFMatcher matcher(cv::NORM_HAMMING, true);

The detector outputs KeyPoint vectors;
std::vector<cv::KeyPoint> keypoints_0, keypoints_1;
and we run this with
  detector.detect(im0, keypoints_0);
  detector.detect(im1, keypoints_1);
The extractor pulls out those KeyPoints into a matrix:
cv::Mat desc_0, desc_1;
  extractor.compute(im0, keypoints_0, desc_0);
  extractor.compute(im1, keypoints_1, desc_1);
And then the matcher is passed the two descriptor vectors, and produces a set of matches;
std::vector<cv::DMatch> matches;
  matcher.match(desc_0, desc_1, matches);
Following this we need to filter the matches to pick out the "best" ones. In this case we look at the "distance" field of each match where a lower value is "better" (i.e. a closer match).
We check the complete set for the min (and  max, if required), and then only allow through matches that are within a given distance range. In this case we tweak the distance we accept manually (i.e. the magic numbers in the check "matches[i].distance" against "(3*min_dist)" || "< 10".
double min_dist = 9999;
  for (uint i=0; i< matches.size(); i++)
  {
    double dist = matches[i].distance;
    if(dist < min_dist)
      min_dist = dist;
  }

  std::vector<cv::DMatch> good_match;
  for (uint i=0; i< matches.size(); i++)
  {
    if ((matches[i].distance < (3*min_dist)) || (matches[i].distance < 10))
      good_match.push_back(matches[i]);
  }
At this point we can visually review the matches by passing the vector, along with the keypoints and source images to drawMatches(), using this piece of boilerplate code:
 cv::Mat imgMatch;
  cv::drawMatches(im0, keypoints_0, im1, keypoints_1, good_match, imgMatch,
    cv::Scalar::all(-1), cv::Scalar::all(-1),
    std::vector<char>(), cv::DrawMatchesFlags::NOT_DRAW_SINGLE_POINTS );
  cv::namedWindow("matches", CV_WINDOW_KEEPRATIO);
  cv::imshow("matches", imgMatch);

  cv::waitKey(0);

For example between Red and Blue we see
This shows that we have a few stray points. Annoyingly these are "strong" matches as far as OpenCV is concerned, and we could be smarter and try to eliminate them, but for the simple example we've actually got a good enough set to work with here given the number of genuine matches, so we can get a translation matrix from findHomography(), operating on the points in this match: i.e.
std::vector<cv::Point2f> pts1;
std::vector<cv::Point2f> pts2;
for(uint i = 0; i < good_match.size(); i++)  {   pts1.push_back(keypoints_0[good_match[i].queryIdx].pt);    pts2.push_back(keypoints_1[good_match[i].trainIdx].pt);
}
return cv::findHomography(pts2, pts1, CV_RANSAC);
So using the above function we can calculate the transform required to align the Red image to the green and blue image respectively. As a result if we have ch8[0] as Blue, ch8[1] as green and ch8[2] as red we can obtain the transform matrices as:
cv::Mat HgB, HgG;
  HgB = GetAlignment(ch8[2], ch8[0]);
  HgG = GetAlignment(ch8[2], ch8[1]);
And then build an image list which contains the Red image, and the Blue and Green images warped through  warpPerspective() to align them with the Red channel image:
std::vector<cv::Mat> imlist(3);
  imlist[2] = im16d[2];
  warpPerspective(im16d[0], imlist[0], HgB, cv::Size(im16d[2].cols, im16d[2].rows));
  warpPerspective(im16d[1], imlist[1], HgG, cv::Size(im16d[2].cols, im16d[2].rows));
And combine this into a single output image as before:
cv::Mat out;
cv::merge(imlist, out);
cv::imshow("combined", out);
cv::waitKey(0);


And... it's Grey. Well, at least we're done.

This isn't great - we have some obvious colour fringes in the output, although this is also partly due to the time between images and relative movement of the moon and the camera. And we ideally want to cull the outlying points. And, well, it's just grey.
But for an automatic alignment using the "stock" calls it's not that bad.... Next up we can look at being a little smarter.

Wednesday, 10 December 2014

Colours - Simple scaling and combining

So far the images have all been retrieved as greyscale, the basic sensor return format.

To support colour capture then the probes have image filters, which are placed in front of the sensor and can then be used to filter for colour, or for various chemical or non-visible wavelengths.

In the case of Voyager the two cameras each have a filter wheel, and for Cassini  each has a pair of image filters. For details on the particular Cassini filter combinations NASA have an FAQ on this topic.

For our purposes then to render colour or false colour images we use the filter information from the tag list, and then combine multiple images of a single target to generate the final image. Usually some minor position shifts are needed to align multiple images though.

There's a few ways to try and stitch the images together. Initially we'll do this manually, since most of OpenCV's stock feature detection algorithms will cause us problems on “plain circles taken with different filters” - we won't get good feature matches from some of this data. Following this we can look at piecing together a colour image of a moon, using feature detection to automate the alignment process slightly.

The first, and simplest thing we can do is to grab three images, one red, one green and one blue. Ideally these images should be taken in sequence and should all be showing the same subject in (roughly) the same position.

For this one we'll take the three files
  • "N1461106394_1.IMG"
  • "N1461106328_1.IMG"
  • "N1461106361_1.IMG
These are blue, green and red images of Saturn respectively.

Since we will have to align the images we want to handle X/Y offsets and scaling of the images relative to each other as a minimum: We actually want to undistort the image frame, but simple offsets will do for now.

To do this under open CV we use an Affine Transform.

Code-wise we provide two triangles – one from the image before the modification and the second from the image afterwards. We can then tell OpenCV to calculate the change in the points required to map from one to the other, and then apply this change to the image.


The simplest example is to think of three corners of the image. If we have the translation for any three corners then it defines the scale and repositioning of the image:

So we can do this in OpenCV with:
  • Create a new (empty) image
  • Get three of the corners of the old image
  • Get the three corners when the image has been scaled and moved
  • Calculate the transform
  • Apply it to the source image, output in the destination image.
In this case I calculated the initially transform manually, to align the rings of Saturn, and taking the (minor) offset for this case into the Red layer information.

So we create a new empty matrix with:
dst = cv::Mat::zeros(_im16.rows, _im16.cols, _im16.type());
And the three original corner points are simply: cv::Point2f st[3];
st[0] = cv::Point2f(0,0);
st[1] = cv::Point2f(_im16.cols-1, 0 );
st[2] = cv::Point2f(0,_im16.rows-1);
Following the scale the new corner point positions will be given by:

double scale;
double offsetX;
double offsetY;
cv::Point2f dt[3];
...
dt[0] = cv::Point2f(offsetX, offsetY);
dt[1] = cv::Point2f((_im16.cols*scale) +offsetX, offsetY);
dt[2] = cv::Point2f(offsetX, (_im16.rows*scale)+offsetY);

And then we can calculate the transform using getAffineTransform(), and apply it with warpAffine(), i.e.:
scalemat = getAffineTransform(st, dt);
cv::warpAffine(_im16, dst, scalemat, dst.size());
And now dst contains our scaled and shifted output. So for our three images we can define a simple structure

#define IM_BLUE (0)
#define IM_GREEN (1)
#define IM_RED (2)
typedef struct
{
   QString nm;
   int color;
   double scale;
   double xshift;
   double yshift;
}CFILES;
CFILES fileinfo[] = {
{"N1461106394_1.IMG", IM_BLUE, 1, 0, 0},
{"N1461106328_1.IMG", IM_GREEN, 1, 0, 0},
{"N1461106361_1.IMG", IM_RED, 1, -3, -1},
};
The BGR ordering is important - it's OpenCV's standard colour ordering for RGB data.
At this point we can simply walk through and decode the three images in turn, placing an appropriately scaled and shifted image in an output Mat.
i.e.: assuming we have:
RawImage* raw[3];
std::vector<cv::Mat> ch(3);
And the basic decode is:

col = fileinfo[i].color;
ch[col] = DecodeAndScale(fileinfo[i].nm);
We can now combine the three colour channels held in the channels vector using the OpenCV merge() function:


cv::Mat out;
cv::merge(ch, out);
And now the Matrix “out” contains our output image; Adding some simple scaling up the colour to use the entire dynamic range (as per the original greyscale scaling code) gets us:


This isn't perfectly aligned, but it's a proof of concept that our simple channel merge for colours is at least sane. We'd need to handle undistortion and colour correction of data to get the colours "right", but it'll do for now.

Sunday, 30 November 2014

Using OpenCV to process the images

We've hit the limits of our simple QImage based approach, so the next step is to use OpenCV, which will allow us to process and store the images without data loss.

Adding OpenCV to the project

(At least on Linux) This is simply a modification to the .pro file, using pkgconfig to provide Qt with the linking information necessary:
CONFIG += link_pkgconfig
PKGCONFIG += opencv

Getting the data into OpenCV

The basic structure for storing image data in OpenCV is cv::Mat, so we have
cv::Mat _im16;
and create our image store with
_im16 = cv::Mat(_height, _width, CV_16UC1);
The CV_16UC1 type tells OpenCV the data in the matrix is "16 bit, unsigned, single channel". At this point we can transfer data from our list of ints, and set a single point in the output with:
 _im16.at<uint16_t>(Point(x,y)) = raw->GetPixel(x,y);
Where  raw->GetPixel(x,y) is accessing data from the integer array in the existing code.
So putting this together, and adding a simple 8 bit handler (which populates an  8UC1 cv::Mat then converts it to 16UC1) we have a fragment like:
  if (_samplesz == 1)
  {
  cv::Mat im8;
    im8 = cv::Mat(_height, _width, CV_8UC1);
    for (x=0; x < _width; x++)
      for (y=0; y < _height; y++)
        im8.at<uint8_t>(Point(x,y)) = raw->GetPixel(x,y);

    im8.convertTo(_im16, CV_16UC1);
  }
  else if (_samplesz == 2)
  {
    _im16 = cv::Mat(_height, _width, CV_16UC1);
    for (x=0; x < _width; x++)
      for (y=0; y < _height; y++)
        _im16.at<uint16_t>(Point(x,y)) = raw->GetPixel(x,y);
  }
  else
  {
   // error case...
  }
At this point we scale the values we have in the matrix to use the full 16 bit range range, so we can use the OpenCV function minMaxLoc() to get the minimum and maximum values in the array, and we can scale the values with convertTo(), i.e.
  double minVal;
 double maxVal;
 double scale;
  minMaxLoc(_im16, &minVal, &maxVal,NULL, NULL);
  scale = 65535.0/(maxVal - minVal);
  _im16 -= minVal;
  _im16.convertTo(_im16, -1, scale);

(Note the scalar subtraction of minVal which reduces the lowest pixel values to 0).
Following this we can write the OpenCV image with
QString name
...
const char* nm;
  nm = name.toLocal8Bit();
  imwrite(nm, _im16);
And we can display it on the screen with the following boilerplate:
  namedWindow( "Display image", WINDOW_AUTOSIZE);
 imshow( "Display image", _im16);
 waitKey(0)
;

More Transforms

Since we now have the image in OpenCV data structures we can apply some basic transforms to the image, such as this crude despeckle to remove some of the pixel glitches, by doing a Morphological Transform pair. Although this is lossy it effectively removes the speckling; if we didn't want data loss we could apply the scaling we derive from this image to the original:
Mat out;
Mat el;

QString nm;
const char* nms;

  el = getStructuringElement(MORPH_ELLIPSE, Size(3,3));

  out = _im16.clone();
// Despeckle
  morphologyEx(out, out, MORPH_OPEN, el);
  morphologyEx(out, out, MORPH_CLOSE, el);
// Rescale image pixels
  minMaxLoc(out, &minVal, &maxVal,NULL, NULL);
  scale = 65535.0/(maxVal - minVal);
  out -= minVal;
  out.convertTo(
out, -1, scale);
  nm = "despeckle_" + QString::number(sz) +".tiff";
  nms = nm.toLocal8Bit();
  imwrite(nms, out);
We use the .clone() in the above code to avoid modifying the base image, and the result of the despeckle and rescale in the output file is:

Which is an improvement on the darker image we had previously, regardless of the loss.... Oddly this process highlights a bright line at the top of the image, which implies we may have a parser issue in the basic data to resolve. However that'll do for today...

Thursday, 27 November 2014

Cassini - Getting the pixel data (and writing an image)

Getting the image data

So, given the information we have from the last step we can extract pixel information. The process is essentially:
  • Get the file data
  • Skip the label, using the LBLSIZE
  • Skip the binary header, using NBB and RECSIZE
  • Extract image records - count given by NL * NB (or N2 * N3)
Where each extracted image record is:
  • Extract RECSIZE data
  • Skip NBB Bytes
  • Append what's left to the image data output
This gives us an array of raw pixel data values in the image data output.

So for the case where we have the data in
QByteArray ba;

And we can access the tags we extracted using the method:
raw->GetLabel("foo")

And we have the output image data:
QByteArray im;

Then we can just do:
  int datalineskip = raw->GetLabel("NBB").toInt();
  int datalinesz = raw->GetLabel("RECSIZE").toInt();
  ba = ba.remove(0, raw->GetLabel("LBLSIZE").toInt());

  for (int i=0 ; i< raw->GetLabel("NLB").toInt() ; i++)
  {
  QByteArray imd;
    imd = ba.left(datalinesz);
    imd.remove(0, datalineskip);
    //    _binaryheader.append(imd);
    ba.remove(0, datalinesz);
  }

int imsz = raw->GetLabel("N2").toInt();
  imsz = imsz * raw->GetLabel("N3").toInt();
  for (int i=0 ; i< imsz; i++)
  {
  QByteArray imd;
    imd = ba.left(datalinesz);
    imd.remove(0, datalineskip);
    im.append(imd);
    ba.remove(0, datalinesz);
  }
And we're done. There's a few assumptions about the way data is organised here, and we're just skipping over the binary header data so we could do that in a simpler manner, but this processing seems to hold up for the records I've used it on.

Turning it into Pixel Data

Now we have the image data we can turn these into pixel values. In this case, since we're dealing with larger sample sizes than a byte, we'll store them as alist of integers:
QList<int> _data;
This is overkill, since the samples are 16bpp, and practically are actually 12bit values from the DAC, but we're still only dealing with 1024*1024 images at most so this isn't enough for me to worry about. We could use an int16_t or similar if we had to though.

So  our line and sample counts are given by:
int lineCount = _labels["NL"].toInt();
int sampleCount = _labels["NS"].toInt();

And we can get the actual pixel data depth from
 if (_labels["FORMAT"] == "'BYTE'")
 {
   _samplesz = 1;
 }
 else if (_labels["FORMAT"] == "'HALF'")
 {
   _samplesz = 2;
 }

... else error...

Then we can just form up the data into the output list either directly, or as byte pairs concatenated to 16bpp values:
int cursor =0;
  if (_samplesz == 1)
  {
    for (int y=0; y < lineCount; y++)
      for (int x=0; x < sampleCount; x++)
      {
        if (cursor >= im.size())
          return;
        _data.append(im.at(cursor++));
      }
  }
  else if (_samplesz == 2)
  {
    for (int y=0; y < lineCount; y++)
      for (int x=0; x < sampleCount; x++)
      {
        int v;
        int vh,vl;

          vh = (unsigned char) im.at(cursor);
          vl = (unsigned char) im.at(cursor+1);
          v = (vh << 8) | vl;
          _data.append(v);
        cursor +=_samplesz;
        if (cursor >= im.size())
        {
        return;
        }
    }

 }

Then making a bitmap

Since the basic bitmap types only hold 8bpp we shift down the 16bpp values by 4 to reduce the range when we generate our preview image. This is a bit rubbish, and we really ought to do a proper min/max calculation and level correction  but it'll do for this example. Otherwise this is basically the same as the QImage process we followed previously.
  _image = QImage(_width, _height, QImage::Format_RGB32);

  cursor =0;
  for (int y=0; y < _height; y++)
  {
    for (int x=0; x < _width; x++)
    {
    QRgb value;
    int v;
      v = _data.at(cursor++);
      if (_samplesz == 2)
        v = v >>4;

      value = qRgb(v,v,v);
      _image.setPixel(x,y, value);
    }
  }

  sname = _labels["MISSION_NAME"];
  sname += " " + _labels["INSTRUMENT_NAME"];
  sname += " " + _labels["IMAGE_TIME"];
  sname += " " + _labels["TARGET_NAME"];
  sname += " " + _labels["FILTER_NAME"];
  sname += ".png";

  _image.save(sname);

And the output is

For this case we have the output

From Nasa's N1767153499_1.IMG. See Nasa for copyright/usage information.

This is basically the same image you can get from this NASA link, of Titan. However we don't have the image compression artefacts. We're also dimmer thanks to our hacky "shift down by 4", where we could be much smarter about the dynamic range - particularly since the higher values appear to be camera pixel glitches.
Putting the image into something like the gimp and playing with the threshold tools show we could be much better at using the entire dynamic range of the image.

And so onto the next thing: Since the 16 bit grayscale isn't handled well by Qt, and we want to do some more image processing on the raw data, then next time we'll look at replacing our Qt image handler with some OpenCV processing...

Cassini - Getting out the header information

Using Regular Expressions to Extract the tags

Rather than dump all the tags into a Hash, as we did with the Voyager data, we will use a regular expression to parse out the tag data from the image file.

In the case of a Vicar record then we can isolate the data value by looking for the occurrence of specific tag with a pattern like:

"one or more spaces" "the tag" "optional spaces" "=" "optional spaces"
In regexp speak this would be something like

QString label;
...
QString pattern;
pattern = "\\s+" + label + "\\s*\\=\\s*";

Where '\s' is regular expression speak for a whitespace character, "+" means "1 or more", and "*" means "0 or more". For the Qt code we have to put two backslashes due to the C string interpretation rules.

One place the above regular expression algorithm will fall down is actually on LBLSIZE, since this is at the start of the file and won't have a leading space. So we actually want the pattern formed to take account of this special case - we can do that with a simple flag marker:

bool start
...

    if (!start)
      pattern = "\\s+" + label;
    else
      pattern = label;

    pattern +="\\s*\\=\\s*";



We can put this into a Qt Regular Expression, and find a match in a byte array  with

QByteArray ba
...
QRegularExpression re;
...
  re.setPattern(pattern);
  QRegularExpressionMatch match = re.match(ba);

  if (match.hasMatch())
  {

...

And then we can locate the end of the match (i.e. the data value) with:
  int off = match.capturedEnd();

One minor point is that the QRegularExpression class is a Qt5 only - for Qt4 there's a QRegExp, which provides similar functionality but with a different API - I'm sticking to Qt5 for this.

At this point we extract the data value. This could be done with a more complex regular expression, but to keep it simple we'll just implement a character by character walk, keeping an eye out for quoting

  do
 {
  char ch = ba[off++];
  result.append(ch);

    if (ch == '\'')
      quoted = !quoted;
    if ((!quoted) && (ch == ' '))

      done = true;
  }while ((off < ba.size()) && (!done));

i.e. if we hit a whitespace outside a quote we're done, otherwise we just append data to the output.

The tags we're looking for

Since we're scanning for a specific set of tags, here's what we want to find:

 "LBLSIZE"
The size of the label header, in bytes. We can use this to locate the image data.

  "TYPE"
Type of data - this should always say 'IMAGE' for the files we're handling.

  "ORG"
The file Organisation; this tells us how the pixel data is arranged, and what the various size tags actually mean. In this case we're primarily looking for a value of 'BSQ' in this field, which means "Band Sequential".  This means the image data is arranged by per-line samples, arranged as a number of lines, and the lines are grouped into bands. The alternative arrangements can interleave the lines (BIL) or pixels (BIP). For now we can check for BSQ and only process files with this data layout.

  "FORMAT"
The format of pixel samples. We're going to deal with 'BYTE' (8bpp) and 'HALF' (signed 16bpp) only.

  "NL"
  "NS"
  "NB"

The count of Lines, samples and bands in the image. Essentially this gives us information we need to recover and correctly size the output image. Note that since we use the Record Size and prefix (RECSIZE & NBB) we don't actually need to use NS directly; more on this when we extract image data.

  "N1"
  "N2"
  "N3"

For BSQ Organised files these are equivalent to 'NS', 'NL' and 'NB' respectively. Since we only handle BSQ we use these interchangeably.

  "NBB"
Binary prefix bytes - this is used when we extract image data from the pixel area, by telling us how much of the line is binary prefix.

  "NLB"
Label area size, we remove this - more on this when we process the image.

  "RECSIZE"
Record size - this tells us about the underlying size of data chunks in the label and pixel data regions - more on this when we process the image.

  "MISSION_NAME"
  "INSTRUMENT_NAME"
  "IMAGE_TIME"
  "TARGET_NAME"
  "FILTER_NAME"
Image meta data we can use to classify the images.

And how we get them

This is fairly simple - we just walk a list of tags, extract the value using the regular expression parser (which for this example is in "getLabelValue()"), and drop it in a hash
i.e.
 static const char* fl_labels[] = {
  "LBLSIZE",
  "NBB",

...
 And then we have, from the previous code, something like:

QHash<QString, QString> _labels;
int  sz = sizeof(fl_labels)/sizeof(char*);
...
 for (unsigned int i=0; i < sz; i++)
  {
  QString value;
    value = getLabelValue(_data, fl_labels[i], i==0);
    _labels[key] = value;
  }

(Note the "i==0" which is special casing for the start flag in the regular expression function). So, once again, we do a file read to pull the data into a QByteArray with

QString nm
...
QFile fin;
...
  fin.setFileName(nm);
  if (!fin.open(QIODevice::ReadOnly))
    return false;

  _data = fin.readAll();
  fin.close();




and pass it through the tag parser above. Done!

Cassini - Getting the image data and something about the file format

the next few posts will cover images from the Cassini probe.

Getting the data

The main repository I'm grabbing data from is http://pds-imaging.jpl.nasa.gov/volumes/iss.html, and specifically the "Saturn EDR Data Sets" (EDR is "Experimental Data Record").

Note that the early versions of these volumes are primarily calibration data, and include Earth, Venus and Jupiter data sets. Spinning forward to a recent version then let's take
wget http://pds-imaging.jpl.nasa.gov/data/cassini/cassini_orbiter/coiss_2087.tar.gz
This gets us the 367M file coiss_2087.tar.gz, and this contains a bunch of Vicar format files. Vicar file formats date back to the 70's, and the Wikipedia Page contains some links to the data format, and converters.

Going into a bit more detail....

Vicar files and the detached Label file

The files under data/ come in pairs - one is the image file, and one is a detached header label, for example:
N1767153499_1.IMG
N1767153499_1.LBL
(The "'N" means this is from the Narrow Angle camera, and the leading 10 digit number is the spacecraft clock time at which the image was taken. The "_1" part is a revision label, which can be modified when the data is updated).

The Detached Label (.LBL) file should look familiar after parsing the Voyager object records - Similar to our Voyager header records this is basically a set of tag=value entries, with support for some specific per-Object namespaces. It also contains references to associated EDR file and the description of the file for the Vicar format as well (relative directory paths).

The files are described in the document/edrsis.txt file, which references the NASA specifications, but basically for the Vicar format image file then the format is as described in http://www-mipl.jpl.nasa.gov/PAG/public/vug/vug9.html.

So, when parsing into the Vicar file itself, we have the basic areas:
  • The VICAR Label
  • The Binary Label Header
  • The Image region
  • The End of Data Set label
Only the VICAR Label and Image region need to be in the file, and the other fields are optional.

The Image region contains lines of Pixel data, which may be prefixed with a Binary Prefix. This prefix is optional, but can be used to store per-line information.

The VICAR label must be parsed to get the image data format and layout. As with the Voyager files this is basically Key/Value pairs, however unlike Voyager they are not explicitly length delimited. In this case the key/value pairs are seperated by spaces. This has the complication that some quoted sequences may contain spaces, such as
  GAIN_MODE_ID='29 ELECTRONS PER DN'

For the files we're examining the field separators are actually multiple spaces, but to be slightly more robust we can look at the quoting.

Also the System Items in a VICAR header are in a fixed order, which should be enough to give us all the information we need to parse the image section. The first item should be a "LBLSIZE" entry, which will tell us the size of the label field, and we can go on and parse from there. However because we're going to be dealing with Venus express PDS files, which embed Vicar headers, then we won't be strict about this being the first item in the file (or the ordering).

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....

Saturday, 22 November 2014

Breaking Down The File

The Record Format

At the top level the IMQ files have a sequence of records, arranged in a basic "length/data" format.

So the file content is a series of data fields formatted as:

[L0] [L1][D0][D1]...[Dn]

Where L0 & L1 combine to give you a 16 bit length, and then D0..Dn are the target data.

Using some off the shelf Qt functions we can simply break out the entire file content  into a standard Qt byte holding structure, and then break out each record individually into a hash map.

Loading The File

Initially we can load the entire file into memory with a fragment like:
QFile fin;
QByteArray ba;


  fin.setFileName(nm);
  if (!fin.open(QIODevice::ReadOnly))
    return false;
...
  ba = fin.readAll();
  fin.close();

...
The QByteArray structure gives us an array of raw bytes with simple accessors and modifier hooks. So following this the "ba" contains all the file data. Since the files are relatively small (a couple of hundred K) this isn't a big deal.

Breaking Out The Fields

We can then break this down into into a set of independent fields, using QList to a manage a list of QByteArray objects, each of which holds a basic record from the input file, i.e.

  QList<QByteArray> _data;
  QByteArray ba;
  QByteArray bn;
  int rlength, length;

    bn = ba.left(2);
    ba.remove(0,2);

    length = (
              ((unsigned char)bn.at(1) << 8) |
              ((unsigned char)bn.at(0))
             );
    bn = ba.left(length);
    _data.append(bn);

    rlength = length;
    if ((rlength %2) != 0) {
      rlength++;
    }
    ba.remove(0,rlength);

In this piece of logic we:
  • Grab two bytes to get a length
  • Remove the length bytes from the start of the QByteArray
  • Get the given length  of data as a new QByteArray
  • Append the new QByteArray to the end of the QList
  • Remove the bytes from the main QByteArray, If the field length was odd, remove an extra byte
The extra byte is a side effect of the file format, which specifies that every record must contain an even number of bytes.
We should probably be more careful about the endianity of the 16 bit length value construction based on the host, but since this doesn't impact on our construction code this is left as "an exercise for the reader".

At this point we have a set of records, which break down into five distinct regions and for the Voyager images all follow the same layout:
  1. The image label - all the parameters that are associated with this image, such as the instrument used to capture, time of capture, etc as well as data are pointers (more on this later). This is a variable number of records, of which the last simply contains the string "END".
  2. The image histogram - Always two records, which combine to make a table of 256x32 bit integers, indicating the histogram of image elements.
  3. The Encoding histogram - Always three records, which combine to make a table of 511x32 bit integers. This is a set of offset/frequency values which are used to generate the binary tree for Huffman decompression (more on this later)
  4. Engineering table - Always one record, which contains "other" engineering data. For now I'm ignoring this field - we don't need it to decompress/view the image.
  5. The Image Object - This is 800 variable length records, each of which represents a single line of image data. The line is compressed and will extract to an 8 bit 800 pixel wide line, resulting in an 800x800x8bpp grey scale image.

How Many Values In the Encoding Histogram?

The decompression table is 511 entries, but the image only has 8bpp (256 values). At first glance this is a little odd.

This is because only the first pixel is an absolute value; the rest of the line is a set of "offset" values from the previous pixel - this leads to some efficient compression since the pixel-to-pixel differences tend to cluster and compress well, but the side effect is that the largest value swings are from -255 to +255 steps ( i.e. "0+255=255", or "255 -255=0") requiring more bytes to cover the range.

I'll go over this in a bit more detail when we actually come to uncompress the image.

Breaking down the header

The header is basically a set of entries each of which is of the form "Tag = Value".
There are only three exceptions to this general rule: 
  • END_OBJECT
  • END
  • Pure Comment records

The "END" is used to flag the end of the label header, and the start of the first image histogram field follows it.

Comment-only records start with "/*"and have no other data in them.

The  END_OBJECT is related to the "OBJECT" statement. OBJECT is used to qualify a particular set of associated data elements in the file, so for example the file we're using has an Object per histogram, and each Object describes the  specifics. Importantly objects may have overlapping tag values; i.e.:

[30]    "OBJECT = IMAGE_HISTOGRAM"    
[31]    " ITEMS = 256"
[32]    " ITEM_TYPE = VAX_INTEGER"
[33]    " ITEM_BITS                       = 32"
[34]    "END_OBJECT"
[35]    "OBJECT = ENCODING_HISTOGRAM"
[36]    " ITEMS = 511"
[37]    " ITEM_TYPE = VAX_INTEGER"
[38]    " ITEM_BITS = 32"
[39]    "END_OBJECT"

Actually since we know the item size and format is fixed you could just ignore these fields for now, however in practice I track the "current" object when parsing through and prefix the name into the data structure.

Data Locations

There's a couple of pointer fields in the original image data - these are fields which start with the "^" character. e.g.
"^IMAGE_HISTOGRAM = 56"
The pointers are offset by one from the array location (since the pointer values are "1" based), and can be used to locate the data structures by looking directly at the QList used to store the records.


Storing the header

Under Qt we can use the QHash template class as a simple lookup dictionary, so we can make a declaration like:
QHash<QString, QString> _labels;
This allows us to do a simple split at the "=" to separate the key and value, then we can just insert it in the hashtable with
_labels[key] = value;
And retrieve it with:
value = _labels[key];
So, for example, if we dump all the incoming header items into the hashtable we can then retrieve the name of the probe with
name = _labels["SPACECRAFT_NAME"];

Obviously when we insert these in the hash then the ordering is lost, and the pointer records are no longer useful.

Next up will be building the compression table and decompressing the image...

Things I'm glossing over

Removing comments, which start with "/*" from the records, using the .trimmed() and .simplified() methods to clean up the whitespace in the QByteArray entries and error handling throughout, none of which are particularly interesting...

Getting Raw Voyager Data

The main location I'm using for data is http://pds-rings.seti.org/archives/. This provides tar.gz files for a few missions.

To start with let's look at the "raw" Voyager 1 & 2 image archives. These are under VG_0xxx, and there's a "file by file" view  at http://pds-rings.seti.org/vol/VG_0xxx/.

Each compressed volume contains the same layout - at the top level the overview document AAREADME.TXT describes the file layout and content, and the  DOCUMENT/VOLINFO.TXT describes the data format.

Breaking down what appears where it's:
  • Uranus on 1-3
  • Subset of Saturn on 4-5
  • Subset of Jupiter on 6-8
  • Neptune 9-12
  • All of Jupiter 13-25
  • All of Saturn on 26-38
The volumes have a separate subdirectory for each astronomy target,  so one for the planet, one each for the major moons, etc.

For a starter I'll pull down C4394416.IMQ from VG_0005 as a reference file for the next bit of work, which will be to extract the basic file structure and header information.


So, what's this blog about?

So, what's this blog about?

Well, recently I read this story on Soylent News, about processing old Voyager probe data, and when reading up the background discovered the data archive of Voyager probe results, alongside Cassini, Venus Express and a bunch of other active probes. All the image data is open access, and available for download.

So, I started hacking on it, and this blog is the results of that hacking - it's a quick programmers guide on how to get at this data, and unpack it for processing and generating images.

As a disclaimer, there are much better ways to get image data if you're only interested in the post-processed versions: NASA and the ESA both provide image search engines which can be used to get images, and you're better off starting over at somewhere like the NASA Planetary Data System Archive. This set of pages is really for those of us who want to look behind the curtain at some of the starting data points and to get an idea of the processing done.

Updates will be irregular - really concentrating on when I have some new code or information on the probe data. Don't expect a daily update here.

And now, on with the show....