Automatic production of HiRISE DEMs

After scraping Dr. Shane Byrne’s website, I realized that there are only 3110 completed HiRISE stereo pairs. These stereo pairs seem to require about 150 GB of intermediate storage to make their way through ISIS and ASP. Afterwards the 1 m /px projected results come to a size of about 200 MB each, or 607 GB total. I have 3 TB free on my home server, so I decided to see how many I could process unattended. It’s interesting to me to see the ASP failure cases and it would give me something to write about every 2 weeks or so. I also just really like pretty pictures.

I’m using the latest ASP development branch to do this processing. It has some critical improvements over ASP 2.0. We managed to fix the deadlock problem that was happening in our parallel tile cache system. We also fixed the segfault that sometimes happens in point2dem. Finally we fixed a bug in our RANSAC implementation that has gone unnoticed since 2006. All of these fixes were critical for making sure the code wouldn’t kill itself. We also killed our random seed so that users could always exactly reproduce their results. All ASP sessions should reproduce prior ASP sessions. These improvements will be available in the next release of ASP, along with RPC support.

Here’s my script. I’ll be clearing the backup archives made by pixz (parallel xz) manually every now and then once I’m sure I no longer have interest in the stereo pair. Also, all of this data is only processed with parabola subpixel. Bayes EM will always improve the result and never harm anything. Since my goal was just to process a lot of HiRISE and see where things broke, I decided to forego this extra refinement step. Anyways, at the resolution I’m sharing with all of you, Parabola still looks pretty sexy. Finally, at no point am I specifying the search range for these stereo pairs. They’re all being processed completely automatically.

#!/bin/bash

# expects ./hirise2dem_v2.sh <left prefix> <right prefix> <user name> <comment>
set -e
set -x

left_prefix=$1
right_prefix=$2
username=$3
comment=$4

hiedr2mosaic.py -h > /dev/null || { echo >&2 "I require ASP to be installed Aborting."; exit 1; }

function write_w_conv_to_height {
    # write_w_conv_to_height
    percent=$1
    min=$2
    max=$3
    string=$4
    log=$5
    height=$(echo "$min + $percent*($max - $min)" | bc | awk '{ printf "%.1f\n", $1}')
    echo $height","$string","$height >> $log
}

function download_hirise_single_img {
    phase_prefix=${1:0:3}
    orbit_prefix=${1:4:4}
    obsv_prefix=${1:0:15}
    if [ -e "$1" ]
    then
        echo "  Found " $1 >> ${left_prefix}.log
    else
        echo "  Downloading " http://hirise-pds.lpl.arizona.edu/PDS/EDR/${phase_prefix}/ORB_${orbit_prefix}00_${orbit_prefix}99/$obsv_prefix/$1 >> ${left_prefix}.log
        wget -O $1 http://hirise-pds.lpl.arizona.edu/PDS/EDR/${phase_prefix}/ORB_${orbit_prefix}00_${orbit_prefix}99/$obsv_prefix/$1
    fi
}

function download_hirise_cube {
    if [ -e "$1_IMGBackup.tpxz" ]
    then
        echo "  Found " $1_IMGBackup.tpxz >> ${left_prefix}.log
        xz -d -k -c $1_IMGBackup.tpxz | tar xv
        hiedr2mosaic.py $1*IMG
        rm $1*IMG
    else
        for i in {0..9}; do
            download_hirise_single_img $1_RED${i}_0.IMG
            download_hirise_single_img $1_RED${i}_1.IMG
        done
        hiedr2mosaic.py $1*IMG

        # Compress and backup the IMGs. This will reduce them to 1/3rd
        # their size.
        tar cvf $1_IMGBackup.tar $1*IMG
        pixz $1_IMGBackup.tar
        rm $1*IMG
    fi
}

# Get the files
mkdir -p $left_prefix
echo ${username} >> ${left_prefix}.log
echo ${comment} >> ${left_prefix}.log
cd $left_prefix
echo "Start: " `date` > ${left_prefix}.log
if [ ! -e "${left_prefix}_RED.mos_hijitreged.norm.cub" ]; then
    download_hirise_cube $left_prefix
    echo "Left Download Fin: " `date` >> ${left_prefix}.log
else
    echo "  Found " ${left_prefix}_RED.mos_hijitreged.norm.cub >> ${left_prefix}.log
fi
if [ ! -e "${right_prefix}_RED.mos_hijitreged.norm.cub" ]; then
    download_hirise_cube $right_prefix
    echo "Right Download Fin: " `date` >> ${left_prefix}.log
else
    echo "  Found " ${right_prefix}_RED.mos_hijitreged.norm.cub >> ${left_prefix}.log
fi

# Work out a nice latitude to project with
caminfo from= ${left_prefix}_RED.mos_hijitreged.norm.cub to=caminfo.txt
echo "Preprocessing Fin: " `date` >> ${left_prefix}.log
latitude=$(grep CenterLatitude caminfo.txt | awk -F "=" '{print $2}' | awk '{if ($1 > 0 ) {print int($1+0.5)} else {print int($1-0.5)}}')
srs_string="+proj=eqc +lat_ts=${latitude} +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=3396000 +b=3396000 +units=m +no_defs"

# # Process the stereo
stereo ${left_prefix}_RED.mos_hijitreged.norm.cub ${right_prefix}_RED.mos_hijitreged.norm.cub ${left_prefix}/${left_prefix} --stereo ../stereo.default
echo "Stereo Fin: " `date` >> ${left_prefix}.log
point2dem ${left_prefix}/${left_prefix}-PC.tif --orthoimage ${left_prefix}/${left_prefix}-L.tif --tr 2 --t_srs "${srs_string}" --nodata -32767
echo "Point2dem Fin: " `date` >> ${left_prefix}.log
mv ${left_prefix}/${left_prefix}-D{RG,EM}.tif .

# Compress the cube files
pixz ${left_prefix}_RED.mos_hijitreged.norm.cub
echo "Compressed Left Fin: " `date` >> ${left_prefix}.log
pixz ${right_prefix}_RED.mos_hijitreged.norm.cub
echo "Compressed Right Fin: " `date` >> ${right_prefix}.log

# Compress the old Stereo file
tar cvf ${left_prefix}.tar ${left_prefix}
pixz ${left_prefix}.tar
echo "Final Compression Fin: " `date` >> ${left_prefix}.log

# Generate colormap and LUT for QGIS
gdalinfo -stats ${left_prefix}-DEM.tif > tmp.txt
maximum=$(grep MAXIMUM tmp.txt | awk -F "=" '{n=$2/100; printf("%i\n",int(n+0.5))}' | awk '{n=$1*100.0; print n}')
minimum=$(grep MINIMUM tmp.txt | awk -F "=" '{n=$2/100; printf("%i\n",int(n-0.5))}' | awk '{n=$1*100.0; print n}')
rm tmp.txt

echo "# QGIS Generated Color Map Export File" > ${left_prefix}_LUT.txt
echo "INTERPOLATION:INTERPOLATED" >> ${left_prefix}_LUT.txt
write_w_conv_to_height 0.0 $minimum $maximum "255,120,255,255" ${left_prefix}_LUT.txt
write_w_conv_to_height 0.2 $minimum $maximum "120,120,255,255" ${left_prefix}_LUT.txt
write_w_conv_to_height 0.4 $minimum $maximum "120,255,255,255" ${left_prefix}_LUT.txt
write_w_conv_to_height 0.6 $minimum $maximum "120,255,120,255" ${left_prefix}_LUT.txt
write_w_conv_to_height 0.7 $minimum $maximum "255,255,120,255" ${left_prefix}_LUT.txt
write_w_conv_to_height 0.9 $minimum $maximum "255,120,120,255" ${left_prefix}_LUT.txt
write_w_conv_to_height 1.0 $minimum $maximum "255,255,255,255" ${left_prefix}_LUT.txt

# Render hillshade
gdaldem hillshade ${left_prefix}-DEM.tif ${left_prefix}_HILLSHADE.tif
colormap --lut /home/zmoratto/raid/asp_test_data/Mars/HiRISE/LMMP_color_medium.lut -s ${left_prefix}_HILLSHADE.tif --min $minimum --max $maximum ${left_prefix}-DEM.tif

# Render MOLA regional crop and create a difference map
/home/zmoratto/raid/asp_test_data/Mars/HiRISE/tools/generate_mola_secondary ${left_prefix}-DEM.tif

On top of producing a colormap and orthoimage, I decided to render a difference map against MOLA. I’m using the gridded MOLA product that is shipped in ISIS. I then bicubicly interpolate it to 1/4th the resolution of my HiRISE DEM. I also render an over-crop to make a hillshaded render of MOLA for the background of all of my pictures. I think it shows how little information is available in MOLA at the resolution I’m rendering at. The code I used to do this is available below and is written against Vision Workbench.

using namespace vw;
using namespace vw::cartography;
namespace fs = boost::filesystem;

double taste_nodata_value( std::string const& filename ) {
  boost::scoped_ptr rsrc( DiskImageResource::open( filename ) );
  return rsrc->nodata_read();
}

int main( int argc, char **argv ) {

  // expect ./generate_mola_secondary <dem>

  std::string input_dem( argv[1] );
  std::string isis_prefix( getenv("ISIS3DATA") );

  // Load georeference information
  GeoReference input_georef, mola_georef;
  read_georeference( input_georef, input_dem );
  read_georeference( mola_georef, isis_prefix+"/base/dems/molaMarsPlanetaryRadius0005.cub" );
  DiskImageView<float> input_image( input_dem );
  DiskImageView<float> mola_image( isis_prefix+"/base/dems/molaMarsPlanetaryRadius0005.cub" );

  // Decide on an input box to rasterize (preferrably a 30% pad around the input dem)
  BBox2i output_bbox = bounding_box( input_image );
  output_bbox.expand( 0.15 * ( input_image.cols() + input_image.rows() ) );

  // Down sample by 2 since there is very little content here
  ImageViewRef<float> output_mola = resample( crop( geo_transform( mola_image, mola_georef, input_georef,
                                                                   PeriodicEdgeExtension(), BicubicInterpolation() ), output_bbox), 0.25 );
  GeoReference output_georef = resample( crop( input_georef, output_bbox ), 0.25 );

  // Render output MOLA DEM with same projection as input
  {
    boost::shared_ptr rsrc( new DiskImageResourceGDAL( fs::basename( input_dem ) + "_MOLA.tif", output_mola.format() ) );
    write_georeference( *rsrc, output_georef );
    block_write_image( *rsrc, output_mola,
                       TerminalProgressCallback("asp", "MOLA Context:") );
  }

  // Render difference map where nodata is -1.
  ImageViewRef< PixelMask<float> > input_mask =
    create_mask( input_image, taste_nodata_value( input_dem ) );
  ImageViewRef< float > difference_map =
    apply_mask(abs(output_mola - crop(subsample(crop(edge_extend(input_mask, ValueEdgeExtension >( PixelMask() )), output_bbox),4),bounding_box(output_mola))),-1);
  {
    boost::shared_ptr rsrc( new DiskImageResourceGDAL( fs::basename( input_dem ) + "_MOLADiff.tif", difference_map.format() ) );
    rsrc->set_nodata_write( -1 );
    write_georeference( *rsrc, output_georef );
    block_write_image( *rsrc, difference_map,
                       TerminalProgressCallback("asp", "MOLA Difference:") );
  }

  return 0;
}

Unfortunately for this post, I didn’t correctly record the processing times since I was still debugging my epic shell script above. In the future, I hope to give better estimates of time. All of this processing is happening on my personal server since it has finally gotten cold enough that I want my heaters on in San Jose. That server has a single AMD FX 8120 processor with 8 cores and a generous 16 GB of DD3 1600 RAM. There’s also three 2 TB hard drives in a software RAID5 that seem to give me pretty epic read bandwidth. This machine cost me $1200 when I built it 6 months ago and is cheaper than my government supplied laptop. I think the most important bit is that the FX 8120 supports the wide AVX instruction set. I compiled all the code with GCC 4.7 so that I could take advantage of this.

Here’s ESP_011277_1825 and ESP_011910_1825, requested by user cweitz. This first pair also shows my first error. The search range solved inside D_sub wasn’t large enough to see the crater floor. This caused a trickle down error where the integer correlation took many days because its search range was defined by white noise and still didn’t contain values that matched the actual disparity. Since it couldn’t zero in, it did wide searches at every pyramid level. Ick.

ESP_014167_1300 and ESP_021947_1300, requested by mcewen. Flat enough that everything just worked perfectly. This pair finished very quickly (about 1.5 days).

ESP_018181_1730 and ESP_018893_1730, requested by rbeyer. Another pair that was pretty flat and just worked.

ESP_018959_1730 and ESP_019025_1730, requested by rbeyer. Again, another perfect stereo pair.

Seems like major improvements could be had if we focus on improving results that happen in the generation of D_sub. Since D_sub only takes a few seconds to calculate, having a large search range on that step wouldn’t be too bad and would fix ESP_011277_1825’s problem. Then I just need to focus on filtering D_sub so noise doesn’t propagate to the integer correlator who works on full resolution.

Bundle Adjusting HiRISE

Last week, I showed off a method for processing LRO-NAC. Now I’m going to show an even more difficult process with HiRISE. Each observation of LRO-NAC was 2 CCDs and it made for a lot of click work. In the case of HiRISE, it has a whopping 10 CCDs. This will make bundle adjustment very tricky if we treated all 20 files as individual cameras. To get around this problem we’ll deploy a new feature in Jigsaw (the Observation option) and we’ll have to modify the HiEDR2Mosaic script that is released with Ames Stereo Pipeline.

HiRISE Preparation

I recently found out that there is a nifty UofA site from Shane Byrne that details all the stereo pairs captured by HiRISE. It is available at this link. I then stalked the user ‘mcewen’ and process a stereo pair he selected. If I was a nice guy I would process stuff selected by user ‘rbeyer’, but he images boring places. Loser! Anyways, for this demo I’ll be processing ESP_013660_1475 and ESP_013950_1475. UofA says it’s a “Gullied 35 Kilometer Diameter Impact Crater in Promethei Terra”. Whatever, I just want 3D.

At this point I would normally tell you to run HiEDR2Mosaic blindly. Unfortunately that won’t work because that script will attempt to project the outer CCDs into the RED5 or RED4’s frame of reference. Any projection is bad because it requires using the spacecraft’s ephemeris, which we don’t trust and we haven’t corrected yet. It also won’t work because noproj will drop some observation serial or whatnots. To fix this, we want to stop HiEDR2Mosaic right before it does noproj. Then we’ll perform our jigsaw and then afterwards we want to resume the process of HiEDR2Mosaic. Doing that required the modification that I checked into Github here. Just download the ‘py.in’ file and rename it to ‘py’. The ‘in’ suffix just means that our build system is going to burn the ASP version into the script at compile time.

Now the run looks like the following:

download all ESP_013660_1475 IMGs
download all ESP_013950_1475 IMGs
hiedr2mosaic.py --stop-at-no-proj ESP_013660*IMG
hiedr2mosaic.py --stop-at-no-proj ESP_013950*IMG

Bundle Adjusting

No magic here, it’s same process as usual for creating a control network. I just picked a special XSpacing of 200 meters and YSpacing of 3 km for Autoseed. This was to make sure that there would be control points on all 10 CCDs. I was guessing that a single HiRISE CCD had a swath of ~500 meters. I also made the MinimumThickness a very small number (something like .0000001) since the HiRISE CCDs are very thin strips. After Autoseed, you then proceed to manually clean up the control network and it will take a very long while. Then you should perform a couple of jigsaw runs to identify control points that have mistakes and correct them in qnet. However, for HiRISE all your jigsaw runs should use the option, observations=yes. This says that images with the same image serial should be treated as the same camera or observation. So CCDs 0 through 9 will be treated as a single camera. Without this flag, jigsaw will never converge.

You might be asking why I didn’t use this for LRO-NAC, that camera also had multiple CCDs. That’s because in the case of LRO-NAC, the CCDs are in two separate optical housing whose position and angle changes noticeably with the thermal cycling of the spacecraft. On HiRISE, all the CCDs are on the same optical plane and they all use the same optics. It is not noticeable at all. If you use the observation option for LRO-NAC, you’ll find that it is impossible to get a sigma0 value under 2 pixels.

In my last post, I tried an idea where I didn’t use any ground control points and tried to make jigsaw auto register to the default ISIS DEM, which is usually the best altimeter data available. I wanted to try that again, unfortunately HiRISE doesn’t have a big enough footprint to cover much detail in MOLA so I added 2 CTX images to the jigsaw problem. Those images were B10_013660_1473_XN_32S256W and B11_013950_1473_XN_32S256W. I made their control network separately and then used cnetmerge to add them to the HiRISE network. I then proceeded to match a bunch of the control points between the two imagers. This would have been easier if I had just processed the images from the beginning together.

Another problem I noticed from the last post was that with every jigsaw run, all control points would start from the ‘apriori’ position. I was updating their radius with cnetnewradii, but their latitude and longitude kept being locked back to their original position. I wanted this to be like ICP, so I modified my script so that after every jigsaw run, ‘apriori’ latitude and longitude would be replaced by their ‘adjusted’ solution. Below is that code.

Adjusted2Apriori.py:

#!/usr/bin/env python
# Expect ./adjusted2apriori.py <input pvl> <output pvl>
import sys

inf = open(sys.argv[1], 'r')
outf = open(sys.argv[2], 'w')

delayed_write = False
delayed_lines = []

for line in inf:
    # Search for Apriori X
    if 'AprioriX' in line:
        delayed_write = True

    if 'AdjustedX' in line:
        delayed_lines.insert(5,line.replace('Adjusted','Apriori'))
    if 'AdjustedY' in line:
        delayed_lines.insert(6,line.replace('Adjusted','Apriori'))
    if 'AdjustedZ' in line:
        delayed_lines.insert(7,line.replace('Adjusted','Apriori'))
        delayed_write = False
        for dline in delayed_lines:
            outf.write( dline )
        delayed_lines = []

    if not delayed_write:
        outf.write( line )
    else:
        if 'AprioriX' in line or 'AprioriY' in line or 'AprioriZ' in line:
            pass
        else:
            delayed_lines.append(line)

bundleadjust.sh:

#!/bin/bash

input_control=control_comb_pointreg_const.net
radius_source=/home/zmoratto/raid/isis/isis3.4.1_ubuntu1204/data/base/dems/molaMarsPlanetaryRadius0005.cub

cp $input_control control_loop.net

for i in `seq 1 50`; do
    echo Iteration $i

    # Convert point's apriori position to be adjusted position
    cnetbin2pvl from= control_loop.net to= control_loop.pvl
    ./adjusted2apriori.py control_loop.pvl control_loop2.pvl
    cnetpvl2bin from= control_loop2.pvl to= control_loop.net

    cnetnewradii cnet= control_loop.net onet= output.net model= $radius_source getlatlon= apriori
    mv output.net control_loop.net

    jigsaw fromlist=cube.lis radius=yes twist=yes cnet= control_loop.net  onet= output.net update=yes spsolve= position camsolve= velocities observations= yes maxits= 100
    mv output.net control_loop.net

    #Gathering statistics for user monitoring
    grep Sigma0: bundleout.txt
    list_length=`wc -l < bundleout.txt`
    interesting_part=`grep -n "POINTS DETAIL" bundleout.txt | awk -F ":" '{print $1}'`
    tail -n $(expr ${list_length} - ${interesting_part}) bundleout.txt | grep RADIUS --color=no | awk -F " " '{print $4}' | awk '{sum+=$1; sumsq+=$1*$1;} END {print "stdev = " sqrt(sumsq/NR - (sum/NR)**2) " meters";}'
    tail -n $(expr ${list_length} - ${interesting_part}) bundleout.txt | grep RADIUS --color=no | awk -F " " '{print $4}' | awk '{mean+=$1} END {print "mean = " mean/NR " meters";}'
done

Running jigsaw in this case was just running my script.

./bundleadjust

Unfortunately I didn’t see the standard deviation of the radius correction reducing ever. So possibly my whole idea of a jigsaw without control points is flawed, or I need a better height source. I really want to go back and redo LRO-NAC now.

After we finish our bundle adjustment, we now finish our HiEDR2Mosaic to create a single image of a HiRISE observation. This looks like the following with the modified script. Notice specifically that I’m giving the script the histitch cube files this time instead of the IMG files.

./hiedr2mosiac.py --resume-at-no-proj ESP_013660*histitch.cub
./hiedr2mosiac.py --resume-at-no-proj ESP_013950*histitch.cub

Processing in ASP

I ran the stereo command like this and I used the default stereo.options file. That just means everything is in full auto and I’m only using parabola subpixel. Thus, my output DEM looks little ugly up close.

stereo ESP_013660_1475_RED.mos_hijitreged.norm.cub ESP_013950_1475_RED.mos_hijitreged.norm.cub HiRISE/HiRISE

Then I used the following point2dem option. Notice I’m using a custom latitude of scale so that the crater of interest will be circular in the projection.

point2dem --t_srs "+proj=eqc +lat_ts=-32 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=3396000 +b=3396000 +units=m +no_defs" --orthoimage HiRISE-L.tif HiRISE-PC.tif --tr 2

Results

Here are the difference maps between MOLA and a non Bundle Adjusted ASP HiRISE DEM (Raw), a Bundle Adjusted DEM with no GCPs, a Bundle Adjusted DEM with no GCPs but with the addition of the CTX imagery. We can see that the bundle adjustment helps definitely; it removes a 200-meter error.

Adding the CTX imagery definitely helped reduce the error against MOLA. However if we look at a DEM that can be created from the additional 2 CTX images we processed, we’ll see that there are still large pockets of error around the rim of the crater. When I flip back and forth between the MOLA hillshade and the CTX hillshade, I think I can definitely see a shift where the CTX DEM is too low and slightly shrunk.

So possibly my whole ‘no-gcps’ idea might be bunk. For the case of Mars, it is really easy to go get ground control points against the THEMIS mosaic or a processed HRSC observation. You’ll just want to chain THEMIS Mosaic to CTX and then to HiRISE since there is such a large resolution change. I’ll leave it as a homework assignment for someone to work out the exact commands you need to run. At least now you should understand how to bundle adjust HiRISE.

New CTX Layer in Google Earth

The Google community silently released a few new features for their Mars mode in Google Earth. MER-B, Opportunity, now has an updated traverse path thanks to a fellow at the Unmanned Spaceflight forum along with a new base map that Ross created. However I’m really excited about a new CTX Global Map that is available. Below is a screen shot:

From this high view you can see that that CTX hasn’t imaged all of Mars. This is expected. CTX isn’t trying to image all of Mars, it is simply the context imager for HiRISE. Meaning that CTX tends to roll tape only when HiRISE is. Never the less, the imagery is still beautiful and, in my opinion, shows more detail that the default base layer from an HRSC composite.

This mosaic was created by simply downloading all CTX data from NASA’s PDS servers and then calibrating and map projected the imagery with USGS’s ISIS software. The composite was then made with in house software from our team at NASA Ames Research Center. This is a Vision Workbench combination plus our Plate Filesystem. It is the same software we used to do a global MOC-NA and HiRISE mosaic for Microsoft’s World Wide Telescope. Unfortunately, I believe those servers have bit the dust. Regardless, if you have free time, I encourage you to check out the CTX Mosaic and Opportunity traverse in Google Earth. Click the planet, select “Mars”, and then in the bottom left select “CTX Mosaic” under “Global Maps”.

ASP v2 Status Update

If you watch your email or ISIS’s website like a hawk, you’ll notice ISIS 3.4.0 just dropped this week. So where’s Ames Stereo Pipeline in all of this?

The current plan is to get Ames Stereo Pipeline 2.0 out before the Planetary Data Workshop on June 25th in Flagstaff. At this meeting I hope to be giving a tutorial and showing off all the cool stuff the community and the team have been developing for ASP.

http://astrogeology.usgs.gov/groups/Planetary-Data-Workshop

I currently have two limiting factors that I hope will get resolved in the next 2 weeks:

  • ASP and Binary Builder need to agree with the changes in ISIS 3.4.0.
  • The NASA Ames Lawyers need to re-license Vision Workbench under the Apache 2 license. ASP has already been re-licensed to Apache 2, however one of its prime dependencies is Vision Workbench. This is a legal block against releasing binaries.

Thank you all for your patience and for doing neat projects with Ames Stereo Pipeline. I enjoy seeing the pretty pictures at science conferences.

Things to look forward in this next release are:

  • Digital Globe Image support.
  • Faster and Memory efficient integer correlator.
  • GDAL projection interface for point2dem.
  • Parameter settings from command line and configuration file.

Making well registered DEMs with ISIS and Ames Stereo Pipeline

Measurements of the Satellite’s position and pose are often not perfect. Those minor mistakes can cause gross offsets when producing DEMs from cameras with such long focal lengths. Correcting this is achieved with bundle adjustment. I’ve discussed this before, but I’d like to expand my previous article by adding ground control points.

Ground Control points are special measurements in an image that have well known geodetic position. On earth, this might be achieved with a large target (an expedition’s tent) that can be observed with a satellite and has a GPS beacon. We don’t currently have such an analog on the Moon or Mars. The best bet is to instead compare an image against a good world map and then pick out individual features. For Mars, I think the best map to register against is the THEMIS global Day IR 100 m map. For the Moon, I would pick the 100 m LRO-WAC DTM and Image mosaic.

In the following sections, I’m going to run through an example for how to register a pair of stereo CTX observations to the THEMIS global map. For those who are interested, my example images are of Hrad Vallis and use P17_007738_2162_XN_36N218W and P18_008028_2149_XN_34N218W.

Getting a THEMIS Global Map Tile into ISIS

The 11th version of the THEMIS Day IR map is available from the following ASU link:

http://www.mars.asu.edu/data/thm_dir_100m/

Unfortunately the data you download is just a simple PGM file that mentions nothing about how it was georeferenced. That information can only be obtained through careful reading from the website. You however just need to know that the data is in Simple Cylindrical projection, it’s Geocentric, it’s sampled at 592.75 PPD, and that the files are indexed by their lower left corner.

I’m going to leaving it up to you to read the documentation to see what the following commands in ISIS do. That documentation is available here from USGS’s website. The tile I’m processing is “lat30_lon120.pgm” and I choose it specifically because it contains my example CTX images. I had to use an intermediate GDAL command because std2isis couldn’t read a pgm file.

gdal_translate -of GTiff lat30_lon120.pgm lat30_lon120.tif
std2isis from=lat30_lon120.tif to=lat30_lon120.cub
maptemplate map=map.template projection=SIMPLECYLINDRICAL resopt=ppd resolution=592.75 clon=0 targetname=Mars targopt=user lattype=planetocentric eqradius=3396190.0
maplab sample=0 line=0 coordinates=latlon lat=60 lon=120 from=lat30_lon120.cub map=map.template

Creating Image to Image Measurements

I’ve discussed this process before, and I’m not going to re-describe it. (I’m quite lazy). However, the following is a quick recap of the commands I used to do it.

Autoseed.def:

Group = PolygonSeederAlgorithm
      Name = Grid
      MinimumThickness = 0.01
      MinimumArea = 1
      XSpacing = 4000
      YSpacing = 4000
End_Group

AutoRegTemplate.def:

Object = AutoRegistration
   Group = Algorithm
     Name         = MaximumCorrelation
     Tolerance    = 0.7
   EndGroup

   Group = PatternChip
     Samples = 19
     Lines   = 19
     MinimumZScore = 1.5
     ValidPercent = 80
   EndGroup

   Group = SearchChip
     Samples = 75
     Lines   = 75
   EndGroup
 EndObject

The ISIS Commands:

parallel spiceinit from={} ::: *cub
parallel footprintinit from={} ::: *cub
echo *cub | xargs -n1 echo > cube.lis
findimageoverlaps from=cube.lis overlaplist=overlap.lis
autoseed fromlist=cube.lis overlaplist=overlap.lis onet=control.net deffile=autoseed.def networkid=ctx pointid=???? description=mars
pointreg fromlist=cube.lis cnet=control.net onet=control_pointreg.net deffile=autoRegTemplate.def
qnet

Then you’ll need to do clean up in ISIS’s Qnet after you’ve completed the above commands. In my experience, pointreg will only correlate maybe 75% of your control measures correctly. To clean these up quickly, just go into qnet and select that you want to filter your control points to show only those that have been “ignored”. Then go and save. I had an excess of control points, and decided to just delete a lot of my ignore points. You can see a distribution of my control points in the following Qnet screenshot.

Creating Ground Control Measurements

Finally, we are breaking new ground here. We are going to measure some ground control points against the THEMIS map tile we downloaded earlier. You could do this in Qnet, however in practice I find this segfaults a lot (I’m using ISIS 3.3.1). So instead we are going to do this process with Qtie.

Open Qtie by simply just typing “qtie” into your terminal. Then click open button or command “O”. It will ask you for a basemap and you should select your THEMIS tile that is in cube format. It will then ask you for a “cube to tie to base”, select your first CTX image of the stereo pair. When it asks you for a control network, click cancel. If you actually give it your pointreg’d control network, it will go into an infinite loop because that control network has no GCPs.

You’ll then want to use the ‘tie’ tool to create a bunch of GCP. If you don’t remember, you have to right click on your CTX image to create a new measurement. You’ll want to capture at least 3 GCPs if you are working with a frame camera. Do about 10 evenly distributed GCPs if you have a linescan camera. CTX is a linescan camera, so have fun clicking.

When matching a low-resolution image to a high-resolution image, I find it very helpful to turn on the affine alignment option in the “Tie Point Tool”. That is the “Geom” radio button on the right side of the window. I also find it very helpful to have the left window constantly flipping between the two input images at a high rate. You can achieve this by clicking the play button in the lower left and then lowering the number next to the play button. You’ll also want to zoom in on your basemap.

On the left is an example of a match I found. This process is difficult and takes some patience. You’ll want to look at your match from different scales just to make sure that you haven’t aligned to some pareidolia.

Once you finish, you oddly have to click the diskette in the “Tie Point Tool” to save your control network. I’ve saved mine as “gcp.net”. If you were selecting ground control points against the LRO-WAC Global Image mosaic (for registering images of the moon), you might want to take a detour here and use the “cnetnewradii” command. That way you can have your GCPs use radius measurements from the LRO-WAC Global DTM. Otherwise your radius measurements will come from ISIS’s default sources of an interpolated LOLA or MOLA for Mars.

We now want to merge our GCP control network to our control network we created with pointreg. This is done with the following:

cnetmerge inputtype=cnets cnet=control_pointreg.net cnet2=gcp.net onet=control_gcp.net networkid=CTX_with_gcp description="against themis 100m"
qnet

Back in qnet, you’ll want to tie the other CTX image to your new GCPs. You can do this by open each GCP point up by double clicking it, and then selecting “Add Measure(s) to Point” in the Qnet Tool window. A screenshot is shown left. You’ll then need to manually align the control measures. Be sure not to move your reference control measure from the first CTX image. You’ll also want to change your GCP’s PointType from “Fixed” to “Constrained”. This will allow your bundle adjustment session to move your GCP to reduce its reprojection error. However it will be rubber banded to the measurement you made in Qtie. This is helpful because your basemap source will usually have mistakes.

After having registered all your GCPs to the other image, and also having set their type to “constrained”, you’ll want to do one last thing. Select all your ground control points in the Control Network Navigator and then press the “Set Apriori/Sigmas” button. I have a screenshot below. The sigmas are your uncertainty measurement for your GCP and basically tell bundle adjustment just how much it can ignore or trust your GCP measurements. Since our measurements were all made from the same source, we’re using the same sigma for all our GCPs. In my opinion a good sigma value is two times the pixel resolution of your basemap. So in this example, I’m putting the GCPs as being accurate to 200 meters in longitude, latitude, and radius directions.

Bundle Adjusting

Alright, it’s time to bundle adjust. This almost exactly the same command we ran before, however we can turn on a new option we couldn’t before:

jigsaw fromlist=cube.lis radius=yes twist=no cnet=control_gcp.net onet=control_jigsaw.net spsolve=position

This time we can solve for the spacecraft’s position since we have GCPs! Do a test run without updating the spice. If it converges to a sigma0 less than 2 px, then you can go ahead an add the “update=yes” option. If you are working with LRO-NAC images, you can likely converge to less than 1 px sigma0. If you unfortunately didn’t converge or your sigma0 has a high value, you likely have a messed up measurement in your control network. This probably came from pointreg, where it managed to match to the incorrect repeated texture. You can identify the mistaken control point by looking at the output residuals.csv file that was created by jigsaw. Just go an manually align the control points in Qnet that show a high residual error. The error should be obvious in Qnet. If it is not, then likely jigsaw fitted to the outliers. If that happened, the control points with low residual errors will have the obvious misregistration in Qnet.

My final DEM results

Creating the above measurements and performing jigsaw took me about 2 hours. Creating the DEMs in Ames Stereo Pipeline took another 3 hours. Here are my results with and without my jigsaw solution rendered in Google Earth. I’ve overlayed a hillshaded render of my DEMs at 50% opacity on top of the THEMIS Map I registered to. My DEMs also look a little crummy since I just used parabola subpixel for speed in processing.

Bundle Adjusted DEM

Non Bundle Adjusted DEM

Interesting GDAL Usage Examples

Despite first appearances, GDAL is way more than an image I/O library. I routinely use it resize images, transcode images, and reproject images. But did you know that it can also draw contour lines, write kml, and orthorectify RPC images? Even if you did, I think it’s good to document here because I keep forgetting.

Render Contour Lines

Publishing DEM’s created by ASP is a tricky because it is easiest visualize height change through a color gradient. That loses its usefulness when publishing to a conference where they print in gray scale. A snazzier option would be a hill-shaded render of a DEM with contour lines plotted on top. You can see an example on right and the commands below:

gdaldem hillshade <DEM> output.tif -z 0.000012
gdal_contour -i 50 <DEM> contour
gdal_rasterize -burn 0 -l contour contour output.tif

The “-z” value on gdaldem hillshade requires some playing around to figure out what is best for your image. The “-i” value on gdal_contour is the interval between lines in units of height defined by the DEM (for ASP products, it’s meters).

KML Outline of a Cube File’s footprint

Rendering of KML outputEver wanted to provide context for your ISIS cube file image on Google Earth? These commands allow you visualize the image on right. The first two lines are ISIS commands.

footprintinit from=<ISIS Cube File>
isis2gml from=<ISIS Cube File> to=output.gml label=footprint
ogr2ogr -f KML output.kml output.gml

Orthorectify an RPC image on a DEM

This one doesn’t apply to planetary imagery, but most Earth images tend to have an RPC camera model attached to them. Download a DEM from USGS’s NED, and then with GDAL you can drape your source imagery over the DEM.

 gdalwarp -of GTiff -co TILED=yes -co COMPRESS=lzw -r cubicspline \
    -rpc -to RPC_DEM=<your DEM> <RPC image> output.tif -dstnodata 0