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.