Thoughts on improving ASP Stereo

Ames Stereo Pipeline’s (ASP) current integer correlator leaves a bit to be desired. Currently it does poorly in scenes with aggressively changing slopes. It is also a coin flip if it finishes in an hour or several days. So I’ve been working on researching a new correlator by reading, implementing, and applying to satellite imagery a select few of the top performers from the Middlebury stereo competition. I started this a long time ago with PatchMatch and I never gave a good conclusion. Now I will summarize by experiences and give a short introduction into the current solution I’m pursuing.

Algorithm Shoot Out!

Semi Global Matching: [1] This is a world recognized well performing stereo algorithm. I don’t need to say its graces. The cons in my opinion are that it uses a lot of memory and that it is only applicable to 1-D searching. For ASP we like to have 2-D searching solution, or optical flow, to handle flaws in the user’s input data and because some users have actual used us for the creation of velocity maps. We might have been to get around the inaccuracies in our users data and the horrors of linescan cameras by calculating a local epipolar vector for each pixel after a bundle adjustment. But I believe we wouldn’t catch the vertical CCD shifts and jitter seen in HiRISE and World View satellites. As for the memory problem, there have been derivative SGM algorithms to fix this problem, but I didn’t evaluate them.

PatchMatch: [2] I really love the idea of starting with a uniform noise guess for the disparity and then propagating lowest cost scores to the neighbors. There were a couple downsides to this algorithm for satellite processing. 1. The cost metric of absolute differencing intensities and gradients performed much worse than an NCC cost metric in the arctic. 2. The run time was horrible because each pixel evaluation didn’t reuse previous comparison used by neighboring pixels. 3. Their slanted window needed to be adapted to support slants in the vertical direction as well as the horizontal for our optical flow demands. I couldn’t find a formulation that would stop the algorithm from cheating by defining the window as 90 degrees from the image geometry. In other words, the PatchMatch algorithm kept finding out that the correlation score was minimal if you define the kernel as having no area.

Despite all of this, a plain jane implementation of PatchMatch using NCC and non-slanted windows performs the same as a brute force dense evaluation of a template window across all disparity values. This also means that places were brute force search fails, so would PatchMatch. But, maybe for extremely large search ranges, PatchMatch might be worth its processing time. I will keep this in the back of mind forever.

PatchMatch with Huber Regularization: [3] This is a neat idea that is built on top of Steinbruecker and Thomas Pock’s “Large Displacement Optical Flow Computation without Warping” [4]. (Seriously though, Thomas Pock hit a gold mine with lets apply a regularity term to everything in computer vision and show an improvement.) I eventually learned how to implement primal dual convex optimization using Handa’s guide [5]. I realize now that everything I need to know is in Heise’s paper [3], but it took me a long time to understand that. But I never implement exactly what the paper described. They wanted a smoothness constraint applied to both the disparity and the normal vector used to define the correlation kernel. Since I couldn’t define a slanted correlation kernel that worked both in horizontal and vertical directions as seen in PatchMatch, I just dropped this feature. Meaning I only implemented a smoothness constraint with the disparity. Implementing this becomes a parameter tuning hell. I could sometimes get this algorithm to produce a reasonable looking disparity. But if I ran it for a few more iterations, it would then proceed to turn slopes into constant disparity values until it hit a color gradient in the input image. So it became a very difficult question for me of, at what point in the iterations do I get a good result? How do I know if this pretty result is actually a valid measurement and not something the smoothness constraint glued together because it managed to out weight the correlation metric?

In the image I provided above, you can see a slight clustering or stair-casing of the disparity as the smoothness constraint wants disparity values to match their neighbors. Also, random noise spikes would appear and neither the total variance (TV) term or the data term would remove them. They are stable minimas. I wonder if a TVL1 smoothnss term would be better than a TVHuber.

As Rigid As Possible Stereo under Second Order Smoothness Priors: [6] This paper again repeats the idea seen in PatchMatch Huber regularization of having a data term, a regularization term, and theta that with increasing iterations forces the two terms to converge. What I thought was interesting here was their data term. Instead of matching templates between the images for each pixel, instead break the image into quadratic surfaces and then refine the quadratic surfaces. This is incredibly fast evaluating even when using a derivative free Nelder Mead simplex algorithm. Like several orders of magnitude faster. Unfortunately this algorithm has several cons again. 1. They wanted to use the cost metric seen in PatchMatch that again doesn’t work for the satellite imagery of the arctic that I have been evaluating. 2. The data term is incredibly sensitive to its initial seed. If you can’t find a surface that is close to the correct result, the Nelder Mead algorithm will walk away. 3. This algorithm with a smoothness prior is again a parameter tuning hell. I’m not sure that what I tune up for my images will work equally well for the planetary scientists as well as the polar scientists.

Fast Cost-Volume Filtering for Visual Correspondence and Beyond: [7] This is an improvement algorithm to the KAIST paper about Adaptive Support Weights. [8] (Hooray KAIST! Send us more of your grad students.)  They say hey, this is actually a bilateral filter that Yoon is talking about.  They also recently read a paper about performing a fast approximate of the bilateral filter by using a ‘guided’ filter. In the end this is similar to a brute force search except now there is fancy per pixel weighting for each kernel based on image color. This algorithm is easy to implement but fails to scale to larger search regions just like brute force search. Yes this can be applied in a pyramidal fashion but I think in the next section that I’ve hit on a better algorithm. I wouldn’t count this algorithm out all together though. I think it has benefit as a refinement algorithm to the disparity, specifically in cases of urban environments with hard disparity transitions.

What am I pursuing now?

Our users have long known that they could get better results in ASP by first map projecting their input imagery on a prior DEM source like SRTM or MOLA. This reduces the search range. But it also warps the input imagery so that from the perspective of the correlator, the imagery doesn’t have slopes anymore. The downside is that this requires a lot of work on the behalf of the user. They must run a bunch more commands and must also find a prior elevation source. This prior elevation source may or may not correctly register with their new satellite imagery.

My coworker Oleg hit upon an idea of instead using a lower resolution disparity, smoothing it, and then using that disparity to warp the right image to the left before running the final correlation. It’s like map projecting, except with out the maps, camera models, and prior existing elevation source. I’ve been playing with it and made a pyramidal version of this idea. Each layer of the pyramid takes the previous disparity, smooths it, and the warps the right image to the left. Here is an example of a disparity produced with this method up against current ASP correlator’s result. I have single thread rough prototype variant and an in-progress parallel variant I’m working on.

Looks pretty good right? There are some blemishes still that I hope to correct. Surprisingly the parallel implementation of this iterated warping correlator is 2x faster than our current pyramid correlator. Another surprising feature is that the runtime for this mapping algorithm is mostly constant despite the image content. For consecutive pyramid levels, we’ll always be searching a fixed square region, whereas the original ASP pyramid correlator will need to continually adapt to terrain it sees. Once I finish tuning this new algorithm I’ll write another post on exactly why this is the case. There is also a bit of a black art for smoothing the disparity that is used for remapping the right image.

Conclusion

I’m pretty excited again about finding a better correlator for ASP. I still have concerns about how this iterative mapping algorithm will handle occlusions. I also found out that our idea is not completely new. My friend Randy Sargent has been peddling this idea for a while [9]. He even implemented it for the Microscopic Imager (MI) on board the Mars Exploration Rovers. I didn’t even know that software existed! But they used homography matrices, while our ‘new’ idea is using a continuous function. In the end, I hope some of you find my diving into stereo research papers useful. I learned about a lot of cool ideas. Unfortunately very few of them scale to satellite imagery.

Reference

[1] Hirschmuller, Heiko. “Accurate and efficient stereo processing by semi-global matching and mutual information.” Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on. Vol. 2. IEEE, 2005.
[2] Bleyer, Michael, Christoph Rhemann, and Carsten Rother. “PatchMatch Stereo-Stereo Matching with Slanted Support Windows.” BMVC. Vol. 11. 2011.
[3] Heise, Philipp, et al. “PM-Huber: PatchMatch with Huber Regularization for Stereo Matching.” Computer Vision (ICCV), 2013 IEEE International Conference on. IEEE, 2013.
[4] Steinbrucker, Frank, Thomas Pock, and Daniel Cremers. “Large displacement optical flow computation withoutwarping.” Computer Vision, 2009 IEEE 12th International Conference on. IEEE, 2009.
[5] Handa, Ankur, et al. Applications of Legendre-Fenchel transformation to computer vision problems. Vol. 45. Tech. Rep. DTR11-7, Department of Computing at Imperial College London, 2011.
[6] Zhang, Chi, et al. “As-Rigid-As-Possible Stereo under Second Order Smoothness Priors.” Computer Vision–ECCV 2014. Springer International Publishing, 2014. 112-126.
[7] Rhemann, Christoph, et al. “Fast cost-volume filtering for visual correspondence and beyond.” Computer Vision and Pattern Recognition (CVPR), 2011 IEEE Conference on. IEEE, 2011.
[8] Yoon, Kuk-Jin, and In So Kweon. “Adaptive support-weight approach for correspondence search.” IEEE Transactions on Pattern Analysis and Machine Intelligence 28.4 (2006): 650-656.
[9] Sargent, Randy, et al. “The Ames MER microscopic imager toolkit.” Aerospace Conference, 2005 IEEE. IEEE, 2005.

Jigsawing ISIS Ideal Cameras

Ideal cameras that I’m talking about are images that have been noproj’d together, such as LRO-NAC and HiRISE observations. Previously I thought this was impossible because Qnet required unique serial numbers for all observations and getsn on a noproj’d image just returned an ‘unknown’. Turns out ISIS can totally do this. You just need to make a couple modifications.

  1. Download the Ideal data set.
  2. List Ideal as a dataset inside your $ISISROOT/IsisPreferences

After you’ve done this, you are now able to run getsn and qnet successfully on the noproj’d imagery!

> getsn from=<LROC image
>  IDEAL/IDEAL/322147421.58037

This greatly simplifies the workflow I described for HiRISE in this post.

Automatic Production of HiRISE DEMs 4

Server is still running. ASP has changed during this run. We introduced a new IP filtering technique, MPI Parabola was sped up, and added a faster point2dem. I don’t remember when these changes happened to my server while processing these images. I also added a timeout for stereo pairs that take longer than 3 days to process. This is to stop a bad pair from just dominating the machine for 2 weeks. During this run, I also changed the settings so that two stereo pairs are always being processed at the same time. This is to keep my CPU pegged.

HiRISE Processing Log
Left Image Time What Happened
ESP_011543_1665 N/A Timed out
ESP_011563_2200 Lost Too narrow search
ESP_011582_1730 Lost Too narrow search
ESP_011592_1595 1d 10h I dunno
ESP_011608_1525 3h Too narrow search
ESP_011661_1410 N/A Timed out
ESP_011662_1750 1d 5h Pretty good
ESP_011672_1395 1d 12h Awesome
ESP_011675_1470 20h Too narrow search
ESP_011676_1700 1d 10h Awesome
ESP_011677_1655 2h Too narrow search
ESP_011683_1540 2d I dunno
ESP_011688_1760 N/A Timed out
ESP_011701_1790 11h Okay
ESP_011715_1800 15h Pretty good
ESP_011717_1910 2d 12h Pretty good
ESP_011720_1835 21h I dunno
ESP_011722_1460 N/A Timed out
ESP_011728_0985 N/A Timed out
ESP_011740_1830 1d 7h Too narrow search
ESP_011743_1725 5h Pretty good
ESP_011761_1485 13h Too narrow search
ESP_011765_1780 N/A Failed to download
ESP_011767_1640 4h Pretty good
ESP_011772_1790 2h Too narrow search
ESP_011773_1460 21h Bad search range
ESP_011780_1515 7h Awesome
ESP_011785_1875 4h Too narrow search
ESP_011818_1505 2h Awesome

Images

Automatic Production of HiRISE DEMs 3

It’s still kinda cold in San Jose so I’m still processing HiRISE DEMs. Please look at this post to see how I’m creating these elevation models.

HiRISE Processing Log
Left Image Time What Happened
ESP_011377_1835 1h 30m Good
ESP_011384_2185 7h 20m Good
ESP_011385_1695 13h 10m Goodish
ESP_011386_2065 6h 0m Goodish
ESP_011392_1370 2h 0m Good
ESP_011397_1535 2.5 weeks Correlated Noise
ESP_011399_2045 N/A Failed to make cube files
ESP_011406_0945 4h 0m Too small search range
ESP_011417_1755 12h 40m Good
ESP_011425_1775 3h 20m Good
ESP_011440_1335 3h 20m Good
ESP_011443_1380 18h 10m Too small of search range
ESP_011447_0950 1h 40m Good
ESP_011494_1535 2h 30m Good
ESP_011504_1700 1 week Not sure – Maybe poor texture?
ESP_011523_1695 7h 20m Not sure

Images

Prettiest DEM Award of the bunch goes to user ‘jwray’ for picking out ESP_011494_1535 and ESP_011639_1535. Supposedly this is “Sirenum crater floor LTLD with Al-OH”. I don’t know what that means.

Automatic Production of HiRISE DEMs 2

I’m still heating my apartment with HiRISE processing. For details how I’m creating this, please look at this previous post. Let’s look at the recent results!

ESP_011265_1560 and ESP_011331_1560, requested by timparker. This stereo pair was a nightmare. D_sub said part of the south part of the image could be correlated even though only one of the images saw that part. Then in full resolution integer correlation, the algorithms went into an aggressive search for a match that didn’t exist. The end result meant that stereo took 15 days to complete! Yikes!

ESP_011287_2165 and ESP_011564_2165, requested by mcewen. Stereo finished in 5 hours.

ESP_011290_1800 and ESP_011778_1800, requested by mcewen. Stereo finished in 12 hours. It looks like stereo mis-guessed the correlation parameters for the crater chain, thus the long processing time.

ESP_011293_1710 and ESP_019218_1710, requested by cweitz. Stereo finished in 5 days! This might be because there is a lot of disparity happening this image. This is one of those problems where local homography transforms or epipolar alignment would make a major contribution in speed.

ESP_011310_1395 and ESP_011811_1395, requested by jwray. Stereo finished in 2 days. I’m disappointed that the mesa’s edges didn’t come out.

ESP_011314_1585 and ESP_011595_1585, requested by mcewen. Interest point detection failed in preprocessing. ASP just gave up and never completed anything. It is not clear to me why this didn’t work.

ESP_011319_2100 and ESP_011464_2100, requested by ltornabene. Stereo finished in 2 hours.

ESP_011339_1665 and ESP_011972_1665, requested by ltornabene. Stereo finished in 9 hours.

ESP_011350_0945 and ESP_011351_0945, requested by cjhansen. Stereo finished in 2 hours.

ESP_011365_1365 and ESP_011642_1365, requested by alxla. Stereo finished in 11 hours.

ESP_011371_1730 and ESP_011516_1730, requested by rbeyer. Stereo finished in 15 hours.

ESP_011372_1730 and ESP_012295_1730, requested by mcewen. Stereo finished in 9 hours.

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.