LatLon Notebooks

I finally got around to checking my mail at work and found some LatLon notebooks! They’re all island themed with topo maps on their covers. The one pictured is my favorite since it features Mars and the elevation data was from CTX data processed by ASP. Yay!

Thank you Aitor Garcia. I’m now the engineer with the best stationery at NASA ARC.

PC Align

Last month we had a new release of Ames Stereo Pipeline, version 2.3! We’ve performed a lot of bug fixing and implementing new features. But my two new prized features in ASP are pc_align and lronac2mosaic. Today I’d like to only introduce pc_align, a utility for registering DEMs, LIDAR points, and ASP point clouds to each other. All you have to do is specify an input, a reference, and an approximate estimate for how bad you think the misplacement is.

Does that sound like magic? Under the hood, pc_align is performing an implementation of the iterative closest point algorithm (ICP). Specifically we are using internally libpointmatcher library from ETH. What ICP does is iteratively attempt to match every point of the input with the nearest neighbor in the reference point cloud. ICP then solves for a transform that would globally reduce the distances between the current set of matches. Then it repeats itself and performs a new round of matching input to reference and then again solves for another global step. Repeat, repeat, repeat until we no longer see any improvements in the sum of distances between matches.

This means that failure cases for PC_align are when the reference set is too coarse to describe the features seen in the input. An example would be having a HiRISE DEM and then only having a single orbit of MOLA that intersects. A single line of MOLA does nothing to constrain the DEM about the axis of the shot line. The 300 meter post spacing might also not be detailed enough to constrain a DEM that is looking at small features such as dunes on a mostly flat plane. What is required is a large feature that both the DEM and the LIDAR source can resolve.

CTX to HRSC example

Enough about how it works. Examples! I’m going to be uncreative and just process CTX and Gale Crater because they’re fast to process and easy to find. I’ve processed the CTX images P21_009149_1752_XI_04S222W, P21_009294_1752_XI_04S222W, P22_009650_1772_XI_02S222W, and P22_009716_1773_XI_02S223W using stereo options “–alignment affineepipolar –subpixel-mode 1”. This means correlation happens in 15 minutes and then triangulation takes an hour because ISIS subroutines are not thread safe. I’ve plotted these two DEMs on top of a DLR HRSC product, H1927_0000_DT4.IMG. It is important to note that this version of the DLR DEM is referenced against the Mars ellipsoid and not the Aeroid. If your data is referenced against the Aeroid, you’ll need to use dem_geoid to temporarily remove it for processing with pc_align who only understand ellipsoidal datums.

In the above picture, the two ASP created DEMs stick out like sore thumbs and are misplaced by some 200 meters. This is due to pointing information for MRO and subsequently CTX being imperfect (how much can you ask for anyway?). You could run jigsaw and that is the gold standard solution, but that takes a lot of manual effort. So instead let’s use pc_align with the commands shown next.

> pc_align --max-displacement 200 H1927_0000_DT4.tif P22-PC.tif \
        --save-transformed-source-points  -o P22_align/P22_align \
> point2dem --t_srs "+proj=sinu +lon_0=138 +x_0=0 +y_0=0 \
        +a=3396000 +b=3396000 +units=m +no_defs" \
        --nodata -32767 P22_align/P22_align-trans_source.tif

There are 3 important observations to make from the command line above. (1) We are using the max displacement option to set the upper bound of how bad we think we are, 200 meters. (2) I’m feeding the ASP PC file as the input source instead of ASP’s DEM. This is because in (3) with the save transformed source points we’ll be writing out another PC file. PC align can only export PC files so we always have to perform another round of point2dem. It is possible to run stereo -> point2dem -> PC Align -> point2dem, but it means you are unneccesarily resampling your data once. Using the PC file directly from stereo saves us from potential aliasing and removes a point2dem call.

Here’s the final result where both CTX DEMs are plotted on top of the HRSC DEM. Everything looks really good except for that left edge. This might be because the DEM was rendered with incorrect geometry and the output DEM is subtly warped from a perfect solution.

Another cool feature that pc_align author Oleg Alexandrov added was recording the beginning and ending matching errors in CSV files. They’re found with the names <prefix>-{beg,end}_errors.csv. You can load those up in QGIS and plot theirs errors to visualize that pc_align uniformly reduced matching error across the map. (Thanks Ross for showing me how to do this!)

CTX to MOLA

Quite a few MOLA shots can be found inside the CTX footprints. Above is a plot of all MOLA PEDR data for Gale crater. I was able to download this information in CSV format from the MOLA PEDR Query tool from Washington University St. Louis. Conveniently PC_align can read CSV files, just not in the format provided by this tool. PC_align is expecting the data to be in format long, lat, elevation against ellipsoid. What is provide is lat, long, elevation against aeroid, and then radius. So I had to manually edit the CSV in Excel to be in the correct order and create my elevation values by subtracting 3396190 (this number was wrong in first draft) from the radius column. The other added bit of information needed with CSV files is that you’ll need to define the datum to use. If you don’t, pc_align will assume you’re using WGS84.

> pc_align --max-displacement 200 P22-PC.tif mola.csv\
        -o P22_mola/P22_mola --datum D_MARS --save-inv-trans \
> point2dem --t_srs "+proj=sinu +lon_0=138 +x_0=0 +y_0=0 \
        +a=3396000 +b=3396000 +units=m +no_defs" \
        --nodata -32767 P22_mola/P22_mola-trans_reference.tif

Two things to notice in these commands, the inputs are backwards from before and I’m saving the inverse transform.  You can keep things in the same order as when I was aligning to HRSC,  it is just that things will run very slowly. For performance reasons, the denser source should be considered the reference and then you must request the reference to be transformed to the source. You’ll likely routinely be using this inverse form with LIDAR sources.

In the end I was able to reduce alignment error for my CTX DEMs from being over 50 meters to being less than 15 meters against MOLA and from over 100 meter to 40 meters error against HRSC. A result I’m quite happy with for a single night processing at home. You can see my final composited MOLA registered CTX DEMs on the left. The ASP team will have more information about pc_align in LPSC abstract form next year. We also hope that you try out pc_align and find it worth regular use in your research.

Update:

I goofed in the MOLA example! Using D_MARS implies a datum that is a sphere with 3396190 meter radius. I subtracted the wrong number from MOLA’s radius measurement before (the value 3396000). That probably had some effect on the registration result shown in the pictures, but this mistake is smaller than the shot spacing of MOLA. Meaning the horizontal registration is fine, but my output DTMs are 190 meters higher than they should have been. FYI, D_MOON implies a datum that is a sphere with radius 1737400 meters.

Mass CTX Processing

A few weeks back, Ross Beyer presented my blog posts on autonomous HiRISE DEM processing at the HiRISE team meeting in Tucson, AZ. This brought about a question of could this be performed for CTX. Of course! ASP can be applied to do bulk processing for any of the missions it supports. Earth or any place ISIS and CSpice have defined coordinate system for. Just put in some safeties about run time into the processing scripts because ASP still occasionally goes mad and eats a whole bunch of processing time for no output. (Were working on it!)

Processing CTX stereo pairs however is in fact a little more difficult compared to the HiRISE processing I was doing before. HiRISE lists all their stereo pairs from their website or Dr. Shane Byrne’s website. There’s no equivalent for CTX. Luckily for me, some folks at UofA wrote PairendipityCTX (Chris Schaller?). They provided Ross and I a detailed report of overlapping files and other statistics. I cut out everything but the filenames for my own use and you can get a copy of the list with its 1,542 stereo pairs here.

Another difference was how these two missions stored their data on PDS. I can look at a HiRISE filename and work out its download path in PDS with no trouble. CTX on the other hand seems to have arbitrary volume breaks, thus the download URL is not predictable. My solution is a bad solution, but a quick solution. I wrote a python script that scraped PDS’s servers and identified all the CTX images it has and what their URLs are. I then just ‘grep’ for the URL in the processing scripts. Here’s the resulting text file that lists all of the 50,708 CTX images that existed at the time of my scraping. This is a mean trick because my script can make HTTP requests much faster than a human can. In a sense, I’m DOS’ing the PDS servers. So please copy my results rather than my methods.

Processing scripts

Previously with my autonomous HiRISE processing efforts, I just wrote a Bash script and then spawned it multiple times on my server using GNU parallel. That’s great, but looking around my office I saw a lot of unused computers that I’d like to do my bidding. What I wanted was a job management system like PBS, which is what they use on the super computer. But PDS is a little too heavy and sometimes cost money. So I instead found two alternative projects that seemed to fit the bill, Gearman and Celery. Gearman is the minimalist. Celery required a database backing and multiple ports open with each slave worker. I decided to use Gearman since it seemed simpler to learn.

My code for this little project is available on Github in my Mars3DGearman project. I hope to eventually add HiRISE support. Here’s how it plays out. All machines make a working directory that contains the CTX stereo list, the CTX url lookup list, and the folders DEM and DRG. One machine is designated the server, in my case, my workstation at home. It starts the only instance of ctx_processor.py with the help of a backing ‘gearmand’ executable (gearman daemon). All of the slaves then SSH back to my server and forward 4730, the port used by Gearman for communication. Then all the slave machines can start one or more instances of ctx_worker.py.

The part I haven’t worked out is how to send home the completed results to my main server. The ctx_worker script produces a DEM and orthophoto and then just dumps it locally into the DEM and DRG folder. Gearman does allow sending binary strings back the main server, but I’m betting a 100 MB long string would break something. I chose a different route. Since I already have all the slaves SSH’ing back to my main server, I decided to simply rsync all the slaves’ DEM and DRG folder back home. I don’t have to re-enter my password as I’ve enabled SSH ControlMaster which re-uses previous connections. For now, I just put that rsync in a watch command that fires every 2 hours. Alternatively it could be inside the ctx_worker script. A better bet would be to use SSH keys.

Another worthwhile improvement compared to my HiRISE processing scripts is the inclusion of a timeout method for each step. When it comes to CTX, if the correlation doesn’t finish in 2 hours, I’m not interested. This timeout is achieved through the run_cmd and process_timeout functions in the ctx_worker script. The Internet helped me a lot in figuring out how to make that a reality. Thanks Internet!

Results

These last few days I’ve roped 4 machines into doing my bidding. So far they’ve produced ~260 DEMs. 5 new DEMs completed just while I was writing this article. There are still some hiccups in the process. But when the stars align, I seem to produce over 50 new DEMs every day. I’m not going to show you all of them as that felt like a lot of work to post on to the blog. Instead I’m just going to show off a couple screenshots of some interesting places in Google Mars. The color ramp is a little funky because someday I need to learn to reference everything against the Mars Aeroid and not just the sphere.

Not everything looks as great as those screenshots. Here are some places that failed to correlate. I’m not sure what went wrong and I unfortunately don’t have the time to investigate.

In conclusion, this is just another proof of concept of what is possible. I hope that someday someone else will attempt this and do a better job than me. ASP is not perfect, but it can achieve a lot of processing on its own that could be beneficial to the scientific community.

Related

Shean, D. E., et al. “MRO CTX Stereo Image Processing and Preliminary DEM Quality Assessment.“ Lunar and Planetary Institute Science Conference Abstracts. Vol. 42. 2011.

David Shean, MSSS, and Larry Edwards actually already attempted this once before! In the above abstract you can see that they went above and beyond my little demo and processed 1180 stereo pairs. They also developed some initial steps for registering the output DEMs to MOLA and plotted the relationship convergence angle has on the outcome of a stereo pair.

Moar Processes

Oleg slipped a new feature in ASP 2.1 that I didn’t call too much attention to. That is the “–left-image-crop-win” option which causes the code to process only a subsection of the stereo pair as defined in the coordinates of the left image. Pair this ability with GDAL’s buildvrt to composite tiles and you have yourself an interesting multiprocess variant of ASP. In the next release we’ll be providing a script called “stereo_mpi” that does all of this for you and across multiple computers if you have an MPI environment setup (such as on a cluster or supercomputer).

Our code was already multithreaded and could bring a single computer to its knees with ease. Being a multiprocess application allows us to take over multiple computers. It also allows us to speed up sections of the code that are not thread-safe through parallelization. That is because processes don’t share memory across each other like threads do. Each process gets their own copy of the non-thread-safe ISIS and CSpice libraries and can thus run them simultaneously. However this also means that our image cache system is not shared among the processes. I haven’t noticed this to be too much of a hit in performance.

I still have an account on NASA’s Pleiades, so I decided to create a 3D model of Aeolis Planum using CTX imagery and 3 methods now available in the ASP code repo. Those options are the traditional stereo command using one node, stereo_mpi using one node, and finally stereo_mpi using two nodes. Here are the results:

Aeolis Planum processing runs on Westmere Nodes
Command Walltime CPU Time Mem Used
stereo 00:44:40 08:46:57 1.71 GB
stereo_mpi –mpi 1 00:32:11 11:28:44 5.77 GB
stereo_mpi –mpi 2 00:21:46 10:10:22 5.52 GB

The stereo_mpi command is faster in walltime compared to traditional stereo command entirely because it can parallel process the triangulation step. Unfortunately not every step of ASP can be performed with multiple processes due to interdependencies of the tiles. Here’s a quick handy chart for which steps can be multiprocessed or multithreaded. (Well … we could make the processes actually communicate with each other via MPI but … that would be hard).

ASP Steps: Multiprocess or Multithreaded
PPRC CORR RFNE FLTR TRI
Multithread x x x x DG/RPC sessions only
Multiprocess x x x

Just for reference, here’s my VWRC file I used for all 3 runs and the PBS job script for the 2 node example. All runs were performed with Bayes EM subpixel and homography pre-alignment.

[general]
default_num_threads = 24
write_pool_size = 15
system_cache_size = 200000000
#PBS -S /bin/bash
#PBS -W group_list=#####
#PBS -q normal
#PBS -l select=2:model=wes
#PBS -l walltime=1:30:00
#PBS -m e

# Load MPI so we have the MPIEXEC command for the Stereo_MPI script
module load mpi-mvapich2/1.4.1/intel

cd /u/zmoratto/nobackup/Mars/CTX/Aeolis_Planum_SE
stereo_mpi P02_002002_1738_XI_06S208W.cal.cub P03_002279_1737_XI_06S208W.cal.cub mpi2/mpi2 --mpiexec 2 --processe
s 16 --threads-multi 4 --threads-single 24

Come to think of it, I was probably cheating the traditional stereo by leaving the write pool size set to 15.

Update 2/4/13

I also tried this same experiment with the HiRISE stereo pair of Hrad Vallis that we ship in our binaries. Unfortunately the single node runs didn’t finish in 8 hours and were shut off. Again, this is homography alignment plus Bayes EM subpixel. Everything would have finished much sooner if I used parabola.

HiRISE Hrad Vallis processing runs on Westmere Nodes
Command Walltime CPU Time Mem Used Completed
stereo 08:00:24 106:31:38 2.59 GB Nope
stereo_mpi –mpi 1 08:00:28 181:55:00 10.0 GB Nope
stereo_mpi –mpi 6 02:18:19 221:41:52 44.9 GB Yep

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

Planetary Data Workshop

2 Weeks ago I had the pleasure of going to Flagstaff AZ for the Planetary Data Workshop. It was a conference for scientists to express their needs and for engineer types to discuss their solutions. As can be guessed from the name, our topics were the dispersal, the processing, and the tools used for scientific imagery of the planets in our solar system (primarily from NASA’s robotic missions).

I was at the conference to discuss my software product, Ames Stereo Pipeline. I gave two talks and eventually they’ll be available on YouTube. I also gave an hour long tutorial that is now online. It’s the video above. I’m not sure how interesting it is to watch but it was a lot of fun performing the tutorial. I now realize how much of a nerd I sound like. I ended up throwing away my prepared HiRISE and LRO-NAC imagery. The crowd that attended seemed more interested in the mass processing of CTX imagery. I, unfortunately, did not have any CTX data on my laptop. Instead, I asked Fred Calef from JPL for a CTX stereo pair that he wanted processed. To my benefit, ASP v2 processed it autonomously without hassle during the demo! My laptop managed to chunk through the data in 15 minutes. I spent most of the tutorial just talking about the ancillary files and what users can look into to see if their output will turn out all right.

Shooting from the hip for a tutorial could have bitten me pretty badly. But I think ASP really has improved a lot and is ready for mass production environments. I’m trying to push for one for earth polar imagery but there are many more datasets that could have this same treatment. I think that idea became clear to the 30 guests who attended my tutorial. We’ve had an uptick in downloads and I hope that means I’ll be seeing some cool 3D models in the future.

Sidenote:  I found out that Jay Laura from Penn State has a blog going called Spatially Unadjusted. He’s a GIS guy who also uses ASP (Aww yeah). He presented a poster on his experience of using ASP v1.0.5 for LRO-NAC imagery.

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