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.

Creating a registered LRO-NAC DEM without GCPs

I’ve previously written a tutorial on performing a bundle adjustment, jigsaw in ISIS lingo, using CTX. Let’s try something a little bit more complex with LRO-NAC. If you are not previously aware, LRO-NAC is actually two separate optical assemblies that are joined by the spacecraft bus. Ideally this spacecraft hardware would be solid as a rock and wouldn’t bend. Unfortunately that’s not the case, so what would have been a bundle adjustment between 2 cameras is now 4 in the case of LRO-NAC.

LRO-NAC Data Preparation

Picking a stereo pair has always been the most difficult part of the job for me. I honestly don’t know how to find stereo pairs. So instead I’m just going to recreate something ASU did in SOCET SET. If you are not aware, ASU has processed about 100 DEMs with LRO-NAC. I’m going to recreate their Moore F Crater DEM using M125713813 and M125720601. You can download the raw data via PDS.

I picked this stereo pair blindly and it ended up being extremely difficult to register correctly. In order to help out, I added a fifth image to the bundle adjustment to give another orbit’s worth of spacecraft ephemeris into the measurements. That image was M110383422LE and I picked it because it overlapped 2 of my previous images and had the same lighting conditions as my previous images.

Once you have downloaded all the images, prep the files in ISIS with the following:

parallel "lronac2isis from={} to={.}.cub; lronaccal from={.}.cub to={.}.cal.cub; spiceinit from={.}.cal.cub; rm {.}.cub" ::: *IMG

* The ‘parallel’ command is GNU Parallel, a non-standard system utility that I love dearly.

Attaching a custom Elevation Source

From a previous article I wrote about, we saw that map projecting against the WAC Global DTM gave much better results than working against the sparse LOLA that is the default for map projecting by in ISIS. So I’m going to do that now. Once we have a successful jigsaw result, then we’ll map-project. If you forgot how to attach the WAC Global DTM, here are the commands I used. My image pair was at 37.3N and 185E. So the tile WAC_GLD100_E300N2250_100M.IMG from this link was what I needed.

pds2isis from=WAC_GLD100_E300N2250_100M.IMG to=WAC_GLD100_E300N2250_100M.cub
algebra from=WAC_GLD100_E300N2250_100M.cub to=WAC_GLD100_E300N2250_100M.lvl1.cub operator=unary A=1 C=1737400
demprep from=WAC_GLD100_E300N2250_100M.lvl1.cub to=WAC_GLD100_E300N2250_100M.lvl2.cub
rm WAC_GLD100_E300N2250_100M.cub WAC_GLD100_E300N2250_100M.lvl1.cub

parallel spiceinit from={} shape=USER model=$PATH_TO/WAC_GLD100_E300N2250_100M.lvl2.cub ::: *cal.cub

Making the control network

Making the control network is pretty straightforward and is just time consuming. I performed the commands I previously outlined here on this blog.

Autoseed.def:

Group = PolygonSeederAlgorithm
      Name = Grid
      MinimumThickness = 0.01
      MinimumArea = 1
      XSpacing = 600
      YSpacing = 1000
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 = 120
     Lines   = 100
   EndGroup
 EndObject

Commands

parallel footprintinit from={} ::: *cal.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=lronac pointid=???? description=moon
pointreg fromlist=cube.lis cnet=control.net onet=control_pointreg.net deffile=autoRegTemplate.def
qnet (editting and fixing control_pointreg.net)

I had to spend a lot of time in Qnet for this stereo pair. About half of my automatic matches didn’t work. (Possibly my search range was too small?) In order to correct this, I had to then manually edit the output control network in Qnet. I then filter down all the control points to the ones marked ‘Ignored’. Then I just manually align the measurements. You just need to get the measures close to the right spot by clicking on the same feature. Afterwards, you can then click the ‘Register’ button to get subpixel accurate alignment.

Autoseed and AutoReg also have a bad habit of only matching LE images to LE, and RE images to RE. It is a good idea to run along the edge of an LE image and try to find common points that are seen in both RE images. These are going to be the few control points that have more than 2 control measures. I also had to do this between my 5th image and the other 2 images that it had overlap with. This problem was likely because the overlap between these 3 images didn’t meet the requirement of MinimumThickness in Autoseed.def.

Normally at this point I would then recommend to you to start picking out ground control points. I tried that. Then I went a little insane. I’m better now. Instead we’re going to do crazy stuff in the next section.

The reason we’re not using GCPs, is that the only source to align against is the WAC Image mosaic. I tried this and I could easily get a bundle adjustment result that aligned well with the image mosaic. However the output DEM would always show a massive error between my height results and the results of LOLA and the WAC Global DTM. After adding the 5th image and clamping the trust of my GCP to 50 meters in radius, I got a result that matched LOLA and WAC DTM in an okay manner but didn’t really align at all to the WAC image mosaic. Things made sense once I compared the image mosaic to a hillshade of WAC DTM.

WAC’s Image mosaic doesn’t agree with WAC’s Global DTM! Ghaaaa! Madness! I was so incredibly disappointed. I didn’t know what to do at this point, so I just dropped the ground control points.

Bundle Adjusting

Alrighty then, time to start bundle adjusting with jigsaw. Likely on the first run things are going to blow up because bad measurements still exist in the control network. My solution to this was to first have jigsaw solve for only a few camera parameters and then check out the output control network for features that had the highest error. In Qnet you can filter control measures that have errors higher than a certain amount of pixels. After I fixed those measurements manually, I resave the control network as my input network and then redo the jigsaw.

Each time I redo jigsaw, I let jigsaw solve for more parameters. My pattern was:

jigsaw fromlist=cube.lis radius=no twist=no cnet= control_pointreg_gcp.net onet=control_jigsaw.net update=no
qnet
jigsaw fromlist=cube.lis radius=no twist=yes cnet= control_pointreg.net onet=control_jigsaw.net observations=no update=no
qnet
jigsaw fromlist=cube.lis radius=yes twist=yes cnet= control_pointreg.net onet=control_jigsaw.net observations=no update=no
qnet
jigsaw fromlist=cube.lis radius=yes twist=yes cnet= control_pointreg.net onet=control_jigsaw.net observations=no update=no camsolve=velocities
qnet

After iterating many times, I eventually produced a control network and jigsaw result that had a sigma0 of .34 pixels. However, if we’re to write this camera solution to file and then make a DEM, we would see that I have worse alignment with LOLA than if we hadn’t bundle adjusted at all.

So how can we get a good bundle adjustment result that aligns well with LOLA? One idea would be to pick GCPs directly against LOLA using their gridded data product. I know folks in the USGS who have done exactly this. However I don’t even see how they can determine correct alignment, so I’m going to ‘nix that idea. Another idea would be to have Jigsaw not solve for radius (radius=no), but instead use only the radius values provided by LOLA. This will keep problem in the ball park, but the resulting jigsaw solution will have a sigma0 around 15 pixels! Not solving for radius essentially ties our hands behind our back in the case of LRO-NAC when our radius source is so incredibly low res.

What we need instead is to simply constrain / rubber-band our Control Points to have radius values that are similar to LOLA. Unfortunately, how do we determine what correct Lat Long to sample the LOLA radius for each control measure? My solution was to not define that. I did two things: (1) I changed all my control measures from ‘free’ to ‘constrained’ with a sigma of 200,200,50. This can be performed with cneteditor. (2) I then wrote and ran the following bash script.

#/bin/env bash

input_control=control_pointreg_all_constrained.net
radius_source=/Users/zmoratto/Data/Moon/LROWAC/DTM/WAC_GLD100_E300N2250_100M.lvl2.cub

cp $input_control control_loop.net

for i in {0..100}; do
    echo Iteration $i
    cnetnewradii cnet= control_loop.net onet= output.net model= $radius_source getlatlon= adjusted
    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
    mv output.net control_loop.net

done

Your mind should be blown at this point. I’ve essentially written an iterative closest point algorithm using the ISIS utilities jigsaw and cnetnewradii. This process is only able to work because we have a perfect control network that we vetted with the previous jigsaw runs. This also works because our LOLA data here is not flat, and has plenty of texture. Finally, LRO-NAC was close to the correct solution. Imagine images with horrible attitude data, like Apollo, would just blow up with this technique.

On each iteration, the sigma0 shouldn’t change. Or if it changes, it should be extremely small. That’s because a translation or scale change of the control points and camera location has no effect on the reprojection error. What should be reducing on each iteration is the standard deviation of the radius correction for each control point as listed in bundleout.txt. The sparse 3D model of the problem in Jigsaw is slowly fitting itself to the LOLA DEM but constrained only by the camera geometry. The standard deviation of the radius correction will never go to zero and will eventually plateau. This is because there will always be some error against LOLA since large portions of the measurements are interpolations between orbits.

Disclaimer: I’ve only tried this technique for this LRO-NAC stereo pair. I don’t know if it will work for others. However this was the only technique I could come up with that would produce results that best aligned with LOLA. I couldn’t get anything good to come out of using WAC Image Mosaic GCPs.

Map projecting for Stereo

Now that we have an updated ephemeris, we’re ready to perform a map projection to speed up stereo. Be sure to not accidentally spiceinit your files, which would wipe your updated ephemeris. Notice that I’m map projecting everything at a low resolution. This is because I’m in a hurry for this demo. In a production environment for LRO-NAC, we would likely map project at 1 MPP or 0.5 MPP. This step and ASP’s triangulation are usually the longest part of this entire process.

parallel cam2map pixres=mpp resolution=5 from={} to={/.}.map.cub ::: ../*cal.cub

Running ASP

Nothing special here, I’m just processing the LE-LE, RE-RE, and LE-RE combination. RE-LE stereo pair wouldn’t resolve anything at 5 m/px. Afterwards, I gdalbuildvrt the outputs together. The massive lines separating the individual DEMs that make up this stereo pair may concern you. However that effect will diminish if I had map projected my inputs at a higher resolution. The reason being that, correlation is much like convolution, in that we erode the input image by a half kernel size along the image border. That border gets perceivably smaller if the image is of higher resolution.

parallel -n 2 stereo M125713813{1}.cal.map.cub M125720601{2}.cal.map.cub {1}{2}/{1}{2} ::: LE LE RE RE LE RE
parallel -n2 point2dem -r moon --t_srs \"+proj=eqc +lat_ts=37 +lat_0=0 +lon_0=180 +x_0=0 +y_0=0 +a=1737400 +b=1737400 +units=m +no_defs\" --tr 5 --orthoimage {1}{2}/{1}{2}-L.tif {1}{2}/{1}{2}-PC.tif --nodata -32767 ::: LE LE RE RE LE RE
gdalbuildvrt aspba_drg.vrt */*-DRG.tif
gdalbuildvrt aspba_dem.vrt */*-DEM.tif
parallel gdal_translate -co TILED=yes -co COMPRESS=LZW -of GTiff {} {.}.tif ::: *.vrt

Conclusion and Comparison

Now the moment you’ve been waiting for. The point where I have to prove I’m not crazy. I have three versions of this stereo pair. There is the result produced by ASU using SOCET SET and then there’s ASP with and without Bundle Adjustment. The look of ASP’s result will be a little shoddy because I map projected my input at 5 m/px which prevent me from get the most detail out. It did however save me a lot processing time while I experimented this method out.

Here are the hillshades for those 3 datasets overlaid a hillshade of the LRO-WAC Global DTM which is 100 m/px. There’s not much to tell, but you can at least see that everyone seems to have decent horizontal alignment.

Here are the difference maps between those 3 datasets and the LRO-WAC Global DTM.

Here are the difference maps between those 3 datasets and the LOLA Gridded Data Product, which I resampled to 100 m/px. I think this shows straight forward the improvement bundle adjustment made to the output DEM.

Just for giggles, I also did a difference between the ASU DEM and the bundle adjusted ASP DEM. They’re pretty close to each other. There seems to just be a simple tilt between the two or that the ASP DEM is just shifted to a slightly lower latitude. Who’s more correct is anyone’s guess for now.

 

Ideally at this point I would then compare the ASU DEM and the BA’d ASP DEM against the raw LOLA shot points. This is where curiosity got me and where I move on to the next personal project. Maybe someone else would be interested in a more rigorous review.

Also, I’m not sure if the 5th image is required.

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