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.

Rendering the Moon from AMC

Ames Stereo Pipeline is currently a candidate in the running for NASA’s Software of the Year award. We needed a pretty graphic and decided that making a cool and possibly realistic rendering of Moon would fit the bill. This is a little more difficult than simply hill shading because the Moon has a specular component to it. Hill shading can be interpreted as being only the diffuse component of the phong model. An interesting example of the Moon’s specular compoent is this picture taken with a Hasselblad during Apollo 17.

Below, are videos of my results where the Sun’s projected coordinates sweep from 90 W longitude to 90 E. Both these views are showing map projected imagery, thus this isn’t a true perspective shot. The difference between these videos is the input observer’s altitude above the surface. Lower altitude and more of the specular component can be seen.

I’m using nothing but Apollo Metric imagery for this example. The DEM source was our product for LMMP. The Albedo source was the Apollo Metric Albedo map that Dr. Ara Nefian produced and will eventually be in NASA’s PDS. The photometric model was the Lunar-Lambertian model as described by McEwen’s paper. Shadows were not rendered because that seemed harder than I could accomplish in 24 hours.

Processing Antarctica

I’ve been sick all last week. That hasn’t stopped me from trying to process World View imagery in bulk on NASA’s Pleiades supercomputer. Right now I’m just trying to characterize how big of a challenge it is to process this large satellite data on a limited memory system for an upcoming proposal. I’m not pulling out all the tricks we have to insure that all parts of the image correlate. Still that hasn’t stopped ASP from producing this interesting elevation model of a section of Antarctica’s coastline, just off of Ross Island. Supposedly Marble Point Heliport is in this picture (QGIS told me it was the blue dot at the bottom of the coastline).

I’m using homography alignment, auto search range, parabola subpixel, and no hole filling. The output DEMs were rasterized at 5 meters per pixel. The crosses or fiducials in the image are posted 5 km apart. This represents a composite of 10 pairs of WV01 stereo imagery from 2009 to 2011 and no bundle adjustment or registration has been applied. The image itself is just a render in QGIS where the colorized DEM has had a hillshade render of the same DEM overlayed at 75% transparency.

I haven’t investigated why more of the mountains didn’t come out. When it looks like a whole elevation contour has been dropped, that’s likely because auto search range didn’t guess correctly. When it looks like a side of the mountain didn’t resolve, that’s likely because there was shadow or highlight saturation in the image. Possibly it could also be that ASP couldn’t correlate correctly on such a steep slope.

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.

Winging a DEM for a mission using World View 1

The group I work for at NASA has a big robot that likes to drive in a quarry at speed. Doing this is risky as we could easily put the robot in a position to hurt itself or hurt others. One of things we do to mitigate the risk is by having a prior DEM of the test area. The path planning software can then use the DEM to determine where it is and what terrain is too difficult to cross.

Since ASP recently gained the ability to process Digital Globe and GeoEye imagery (more about that in a later post), I was given a request to make a DEM from some World View 1 imagery they purchased. The location was Basalt Hills, a quarry at the south end of the San Luis Reservoir. To process this imagery with any speed, it is required to map project the imagery on some prior DEM. My choices were SRTM or NED. In my runs, both DEMs have problems. SRTM has holes in the middle of it that needed to be filled so ASP would work correctly. NED had linear jumps in it that ASP couldn’t entirely reverse in its math.

I ended up using SRTM as a seed to create my final DEM of the quarry. If you haven’t seen this, the process looks like the following commands below in ASP 2.0+. What’s happening is that ASP uses an RPC map projection to overlay the imagery over SRTM. When it comes time for ASP to triangulate, it reverses math it used to map project, and then in the case of Digital Globe it will triangulate using the full camera model. Another thing worth noting is that ASP needs control over how the interpolation is performed when doing RPC map projection. This forces us not to use the GDAL utilities during this step and instead use our own custom utility.

parallel rpc_mapproject --tr 0.5 \
      --t_srs'"+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs"' \
      filled_srtm_dem.tif {} {.}.XML {.}.srtm.crop.tif ::: left.TIF right.TIF
stereo left.srtm.crop.tif right.srtm.crop.tif left.XML right.XML \
      r1c1_srtm_crop/r1c1_srtm_crop filled_srtm_dem.tif

Afterwards we got a pretty spiffy result that definitely shows more detail than the prior DEM sources. Unfortunately the result was shifted from the NED DEM source that my crew had previously been using. This ideally would be fixed by bundle adjusting the World View camera locations. It was clearly needed as most of our projected rays only came within 3 meters of each other. Unfortunately ASP doesn’t have that implemented.

EDIT: If I had paid closer attention to my data I would have noticed that a large part of the differences I was seeing between my DEM and USGS’s NED was because the NED data uses a vertical datum. My ASP DEM are referenced against the WGS84 ellipsoid. NED data is referenced against WGS84 plus the NAVD88. This would account for a large part of the 30 meter error I was seeing. (11.19.12)

My “I’m-single-with-nothing-going-on-tonight” solution was the Point Cloud Library. It has promising iterative closest point (ICP) implementations inside it and will eventually have the normal distribution transform algorithm in it. It also has the benefit of having its libraries designed with some forethought compared to the hideous symbol mess that is OpenCV.

PCL's pcd_viewer looking at the quarry.

I achieved ICP with PCL by converted my PC (point cloud) file from ASP into a PCL PCD file [1]. I also converted the NED DEM into a PCD file [2]. I then subsampled my ASP point cloud file to something more manageable by PCL’s all-in-memory tactics [3]. Then I performed ICP to solve for the translation offset I had between the two clouds [4]. My offset ended up being about a 40 meter shift in the north and vertical direction. I then applied this translation back to the ASP PC file [5] so that the DEM and DRG could be re-rendered together using point2dem like normal.

I wrote this code in the middle of the night using a lot of C++ because I’m that guy. Here’s the code I used just for reference in the event that it might help someone. Likely some of the stuff I performed could have been done in Python using GDAL.

1. convert_pc_to_pcd.cc
2. convert_dem_to_pcd.cc
3. pcl_random_subsample.cc
4. pcl_icp_align.cc
5. apply_pc_offset.cc

After rendering a new DEM of the shifted point cloud, I used MDenoise to clean up the DEM a bit. This tool is well documented at its own site (http://personalpages.manchester.ac.uk/staff/neil.mitchell/mdenoise/).

I’ve also been learning some QGIS. Here are some screen shots where you can see the improved difference map between NED and my result after ICP. Generally this whole process was very easy. It leaves me to believe that with some polish this could make a nice automated way to build DEMs and register them against a trusted source. Ideally bundle adjustment would be performed, but I have a hunch that the satellite positioning for Earth targets is so good that very little shape distortion has happen in our DEM triangulations. I hope this has been of interest to some of you out there!

Difference map between the USGS NED map and ASP's WV01 result.

Creating Control Networks and Bundle Adjusting with ISIS3

Bundle Adjustment is the process of back solving for a camera’s trajectory and pose. This process needs to be performed for most satellites images at some point or another because there is always an error in the camera location. Satellite position and velocity is usually found though radio communications via the Deep Space Network via methods such as measuring the antennae direction, time of flight, and doppler effects on the signal. The spacecraft’s pose is usually made from outward facing cameras called star-trackers. All of this is surprisingly accurate considering that it’s a measurement made from at least one planet away but it doesn’t meet the demand of photogrammetrists who wish to register images to meter level precision.

In this article, I’ll run an example of producing a control network and running jigsaw with USGS’s ISIS3 software. I’ll be playing today with two CTX images (P02_001918_1735_XI_06S076W and P02_001984_1735_XI_06S076W.IMG) that looked at the southwest end of Candor Chasma on Mars. Everything that is practiced here will equally apply to other missions with some additional number fiddling.

You probably don’t need this reminder, but before you can do anything, these files need to be ingested into ISIS. I’m going to also take this time to radiometric calibrate and attach spice data. The parallel you see in my examples is GNU Parallel; its not usually installed on systems by default. I strongly recommend that everyone gets it and learns to use it as it is a time saving utility.

parallel mroctx2isis from={} to={.}.cub ::: *.IMG
parallel spiceinit from={} ::: *.cub
parallel ctxcal from={} to={.}.cal.cub ::: *.cub

Now that you have beautiful images of Mars, lets break into the new stuff. We are going to start by building a control network. Control Networks are databases of image measurements that are used during a bundle adjustment. It defines a location in 3D space called Control Points and the pixel locations for which that point projects into, called Control Measures. Before we go too far, we need to build some metadata so that the control network code knows how the images overlap.

parallel footprintinit from={} ::: *cal.cub
echo *cal.cub | xargs –n1 echo > cube.lis
findimageoverlaps from=cube.lis overlaplist=overlap.lis

In the commands above, we have created 2 files, cube.lis and overlap.lis, that we’ll be repeatedly using. An interesting side note, footprintinit has created a vector layer inside each of the cube files that shows the lat-lon outline of the image. If one is so inclined, that vector layer can be extracted with an “isis2gml label=Footprint” call. That gml can then be rendered with the gdal_rasterize tool or can be converted to KML with the ogr2ogr tool.

Since most of the pretty NASA cameras are linescan, we are trying to bundle adjust a trajectory and thus need many control points. About 20-30 points are required. Ideally these control points would be distributed evenly across the surface of each image. ISIS has provided the autoseed command to help with that.

autoseed fromlist=cube.lis overlaplist=overlap.lis onet=control.net deffile=autoseed.def networkid=ctx pointid=???? description=mars

The settings of how autoseed works is defined in the definitions file, autoseed.def. I haven’t given you this; so let’s take a look into what should be inside that file.

Group = PolygonSeederAlgorithm
      Name = Grid
      MinimumThickness = 0.1
      MinimumArea = 1
      XSpacing = 8000
      YSpacing = 8000
End_Group

The minimum thickness defines the minimum ratio between the sides of the region that can have points applied to it. A choice of 1 would define a square and anything less defines thinner and thinner rectangles. The minimum area argument defines the minimum square meters that must be in an overlap region. The last two are the spacing in meters between control points. I played with those two values so that I could get about 50 control points out of autoseed. Having more control points just makes for more work later on in this process.

After the autoseed command, we finally have a control network that we can view with ISIS’s qnet utility. Run qnet in the terminal and a window should pop up. You’ll then have to click ‘open’ and select the cube.lis file and then the control.net file that we created earlier. In the control network navigator window, select the drop down menu so that you can select ‘cubes’. Highlight the names that show up on the left side and then press the ‘view cubes’ button.

You can now see the location of the control points in the two images that we have been playing with. However the alignment between the control measures is imperfect at this point. We can see this visually by requesting that qnet show us a single control point. In the control network navigator window, select the drop down menu and select points. Then in the left side, double click point ‘0001’. A new window should have popped up called ‘Qnet Tool’. You should click on the ‘geom’ option in the bottom right of the window. The two pictures of this window show the control point projected into the reference image (left) and then the second image on right.

You can click the right image or use the arrow keys to reposition the control measure. You can also click the play button on the bottom left of the window so that reference image keeps flipping between images. I prefer to operate with that window playing as it clearly shows the misalignment between measures. An example is show left if you click on the picture.

We could at this point fix these 50 or so control points by hand using qnet. There is instead a better option. ISIS’s pointreg is an automatic control measure registration tool that tends to get about 75% of the points correct. The command to use it looks like the following:

pointreg fromlist=cube.lis cnet=control.net onet=control_pointreg.net deffile=autoRegTemplate.def

Again all the settings are in the definition file. Here are the contents of 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

If you are a user of Ames Stereo Pipeline, these settings should look pretty similar. The search chip defines the search range for which pointreg will look for matching imagery. The pattern chip is simply the kernel size of the matching template. You will likely have to redefine the search range when you are working with new imagery. Use qnet to get an idea for what your pixel shifts are search ranges should be. A search range of 75 pixels in all directions is what happened to work for me with these specific CTX images.

With those settings, I was able to register 38/47 existing control points! The reset I’ll have to register manually in qnet. Using qnet to quickly register points is a bit of a fine art that I don’t think I can describe here. Maybe when I have free time I could make a video.

After you cleaned up the last 9 control points in qnet, we should have a complete control network in the control_pointreg.net file. We are ready to run jigsaw and update the camera pointing. Here’s the command:

jigsaw fromlist=cube.lis update=yes twist=no radius=yes cnet=control_pointreg.net onet=control_ba.net

From my jigsaw’s terminal output, I can see that the starting projection error for this control network was 58.9 pixels on average. This is the sigma0 value under iteration 1. After jigsaw converged by rotating the cameras’ pose and by moving the control points, the sigma0 dropped to 1.2 pixels. This is quite an improvement that should help in DTM quality. If you happen to mess up and write a camera solution to your input files that is incorrect, you can simply run spiceinit to revert your changes.

In the next version of Ames Stereo Pipeline (ver 1.0.6 or 2.0), we’ll probably be providing the ability to render the triangulation error of a DTM. Triangulation error is simply how close the projected rays of the image came together. It is one of many measurements that can be used to judge the quality of a DTM. I’ve gone ahead and rendered DTMs that use both the jigsaw’d and non versions of these CTX images. On the upright right is their hillshaded colormap output. Visually, there’s not a noticeable difference between the two. However the height range has changed drastically. The orignal data produced height ranges between 0.6 km and 4.9 km however the bundle adjusted data produces a range between -8.9 km and -4.4 km. The colormap figure I’ve produced uses two different color scales for those DTMs just to simply show that the DTM hasn’t pivoted any. The only change was a height drop.

I’ve also produce colormap output of the triangulation error. Here we can really see that jigsaw has made a difference for good. The color green represents a triangulation error of 1 meter, blue 0.1 meter, and yellow 10 meters. From the figure on the left, it’s clear to show that for the most part bundle adjustment improved every pixel in the image.

I’m sorry that this has been a run on article. Writing this was also an experiment for me. I hope I’ve shown you how to use ISIS’s control network tools and I’ve managed to show myself that fixed ground control points in jigsaw seem to be required. I have very little trust in the absolute height values in the DTM. I think their relative measurements are correct but I was definitely not expecting the 10 km drop in height between non and bundle adjusted solutions.