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.

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.