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.