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.


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

comment=$4 -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
    height=$(echo "$min + $percent*($max - $min)" | bc | awk '{ printf "%.1f\n", $1}')
    echo $height","$string","$height >> $log

function download_hirise_single_img {
    if [ -e "$1" ]
        echo "  Found " $1 >> ${left_prefix}.log
        echo "  Downloading "${phase_prefix}/ORB_${orbit_prefix}00_${orbit_prefix}99/$obsv_prefix/$1 >> ${left_prefix}.log
        wget -O $1${phase_prefix}/ORB_${orbit_prefix}00_${orbit_prefix}99/$obsv_prefix/$1

function download_hirise_cube {
    if [ -e "$1_IMGBackup.tpxz" ]
        echo "  Found " $1_IMGBackup.tpxz >> ${left_prefix}.log
        xz -d -k -c $1_IMGBackup.tpxz | tar xv $1*IMG
        rm $1*IMG
        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 $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

# 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
    echo "  Found " ${left_prefix}_RED.mos_hijitreged.norm.cub >> ${left_prefix}.log
if [ ! -e "${right_prefix}_RED.mos_hijitreged.norm.cub" ]; then
    download_hirise_cube $right_prefix
    echo "Right Download Fin: " `date` >> ${left_prefix}.log
    echo "  Found " ${right_prefix}_RED.mos_hijitreged.norm.cub >> ${left_prefix}.log

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

3 thoughts on “Automatic production of HiRISE DEMs

  1. This is awesome. I would like to work on a conversion of your stuff to Python scripts (just a personal preference, but also easier for me to manage CPU processes).
    One question: When you say “perfect” do you imply that it approaches the standard HiRISE DEM quality or even improves on it? I have heard the opinion about ASP that it’s nice for creating “quickview-DEMs”, but when using DEM’s for calculations (I’d be calculating solar irradiations) it might be better to wait for a ‘proper’ hand-made DEM? What’s your take on that?

  2. Can you comment a bit more about your MOLA diff? Is it really true that you get up to 400 m (height?) difference to MOLA in most of your products? Can this be true?

  3. Feel free to take/convert whatever you like. I’m currently working on a rewrite myself using Python and Gearman. Ross Beyer challenged me to process CTX in bulk. I’ll write up all of that once I finished.

    I said “perfect” in the sense that ASP managed to correlate the whole image. ASP can deduce better quality than shown here (and the sequential blog posts) with the better Bayes EM subpixel method turned on. It just takes more CPU time than I’m willing to spend for engineering tests.

    As for its quality compared to the handmade DEMs, I can’t really say. I haven’t done a comparison. I imagine we’re in the ballpark with the NGATE option in Socet SET in terms of features resolved. Where we differ is the property you pointed out, the high error in difference to MOLA. When producing DEMs with ASP, we don’t modify or try to refine the stored camera location and pointing. We just use the camera data raw as if it is truth. Unfortunately for missions to other planets, the spacecraft ephemeris available is quite bad. Good enough to target an imager by, but fairly poor for mosaicking and 3D reconstruction. This is where the human operator at the Socet SET station triumphs. Usually the operator will refine camera geometry by hand via matching to prior maps or other data sources. With that ability they are able to insure that the DEM is placed in the correct spot and that it mostly agrees with MOLA.

    As for calculating solar irradiation, which I assume is something along the lines of using the normals in a phong model, I bet you would be okay. Subtract off the average error against MOLA and then manually fix the horizontal against the THEMSIS mosaic, and you are probably in the ballpark. The local variations in our DEM should be pretty close to reality. I believe we just have a large translation error. Though, I do recommend testing this idea on an area where you can both produce an ASP made DEM and download a copy of a UofA made DEM. If you would like a GeoTIFF DEM file for any of the HiRISE results I’ve shown on the blog, please email me directly at my Z.M.MORATTO NASA GOV address. I keep a copy since I haven’t ran out of hard drive space yet.

Leave a Reply

Your email address will not be published. Required fields are marked *