1. Download the Ideal data set.
2. List Ideal as a dataset inside your $ISISROOT/IsisPreferences

After you’ve done this, you are now able to run getsn and qnet successfully on the noproj’d imagery!

> getsn from=<LROC image
>  IDEAL/IDEAL/322147421.58037

This greatly simplifies the workflow I described for HiRISE in this post.

The first problem I’m going to write about has bugged me for a while. That problem is that each LROC-NA observation is actually 2 files and makes using ASP awkward. In the past I’ve been telling people to run perform 4 runs with stereo. Do the combinations of LE-LE, LE-RE, RE-LE, and RE-RE. However UofA had the same problem with HiRISE, which comes down in 20 different files. They had worked out an ISIS process chain that would noproj those files together among other things to make a single observation. I don’t know why such step doesn’t exist for LROC-NA but today I’ll show you what I’ve come up with.

If you try right now to noproj the RE cub file to the LE cube file you’ll find that the program errors out because an ideal camera model for LROC has not been defined. Definitions for ideal noproj cameras for all ISIS cameras are defined in a file at $ISIS3DATA/base/applications/ noprojInstruments003.pvl. Looking at that file we see that there are 5 elements that can be defined for the ideal camera model, TransY, ItransS, TransX, ItransL, and DetectorSamples. DetectorSamples is easy; it’s what the output image width should be in the final image.  The Trans* variables are measured in millimeters and define a focal plane offset from the original camera model we are using. I plan to noproj with the match file being the LE file. Thus Trans* references a shift from the original position of the LE sensors. The Itrans* variables are pixel conversion of the millimeter measurements. It’s used by ISIS to run the math the other direction. If Trans* and Itrans* variable don’t match, funny errors will occur where the CCDs noproj correctly but the triangulation in ASP is completely bonkers. Relating the two is easy, just use the pixel pitch. For LROC-NA that is 7 micrometers per pixel.

So how do we decide what to set the TransY and TransX values to be? If those values are left to zero, the LE image will be centered on the new image and the RE will be butted up beside but falling off the edge of the new image. A good initial guess would be to set TransY to be a shift half the CCD width. A better idea I thought to use was to look at the FK kernel and identify the angle differences between the two cameras and then use the instantaneous field of view measurement to convert to pixel offset between the two CCD origins. Use pixel pitch to convert to millimeters and then divide by two will be the shift that we want. Below are the final noproj settings I used for LROC.

 Group = "LUNAR RECONNAISSANCE ORBITER/NACL"
    TransY = 16.8833
    ItransS = -2411.9
    TransX = 0.6475
    ItransL = -92.5
    DetectorSamples = 10000
  End_Group

At this point we can noproj the imagery and then handmos them together. A naïve version would look something like the following.

> noproj from=originalRE.cub to=RE.cub match=originalLE.cub
> noproj from=originalLE.cub to=LE.cub match=origianlLE.cub
> cp LE.cub LERE_mosaic.cub
> handmos from=RE.cub mosaic=LERE_mosaic.cub

Alas, the LE and RE imagery does not perfectly align. In the HiRISE processing world we would use hijitreg to determine a mean pixel offset. There is no generic form of hijitreg that we can use for LROC-NA. There is the coreg application, but in all my tests this program failed to actually find any correct match points between the images. I’ve personally tried two successful methods. I’ve manually measured the offset using Qtie. Warning: Sending this images to Qtie requires twiddling with how serial numbers are made for ideal cameras. This is something I should document at some point as it would allow bundle adjusting ideal cameras like fully formed HiRISE and LROC images. My other solution was the example correlate program in Vision Workbench.  I did the following to make processing quick (5 minutes run time).

> crop from=LE.cub to=LE.centerline.cub sample=4900 nsamples=200
> crop from=RE.cub to=RE.centerline.cub sample=4900 nsamples=200
> parallel isis2std from={} to={.}.tif format=tiff ::: *centerline.cub
> correlate --h-corr-min -60 --h-corr-max 0 --v-corr-min 0 --v-corr-max 60 LE.centerline.tif RE.centerline.tif

That creates a disparity TIF file. The average of the valid pixels (third channel is 1) can then be used to solve for the mean translation. That translation can then be used during handmos by using the outsample and outline options. Ideally this would all be one program like hijitreg but that is for another time. Below is the final result.

Hijitreg actually does more than just solve for the translation between CCDs on HiRISE. It actually correlates a line of pixels between the CCDs in hope of determining the jitter of the spacecraft. I can do the same!

From the above plots, it doesn’t look like there is much jitter or the state of the data is not in form that I could determine. The only interesting bit here is that the offset in vertical direction changes with the line number. I think this might actually disparity due to terrain. The imagery I used for my testing was M123514622 and M123521405, which happened to be focused on the wall of Slipher crater. The NA cameras are angled 0.106 degrees off from each other in the vertical direction. Ideal stereo geometry would 30 degrees or 15 degrees, but 0.106 degrees should allow some disparity given the massive elevation change into the crater. I wouldn’t recommend using this for 3D reconstruction but it would explain the vertical offset signal. The horizontal signal has less amplitude but does seem like it might be seeing spacecraft jitter. However not enough to determine what the frequency is.

Anyways, making a noproj version of the LROC-NA observation is a huge boon for processing with Ames Stereo Pipeline. Using the options of affine epipolar alignment, no map projection, simple parabola subpixel, and it is possible to make a beautiful DEM/DTM in 2.5 hours. 70% of that time was just during triangulation because ISIS is single threaded. That would be faster with the application parallel_stereo (renamed from stereo_mpi).

“It does, but for pinhole camera models.”
“Why not satellites?”
“Well, they don’t really have epipolar lines but instead curves! Making rectification very difficult and not reversible.”

This has been my default response the last few years. But it has been a long time sense I first tried to understand this topic. Back then I thought the solution would require some difficult to implement polynomial mapping or rubbersheet transforms [1]. I’ve been with NASA a while now, maybe this time I can better grasp what is required.

We’ll I’ve gone back through the papers. Surprisingly epipolar rectification is still a popular subject to get published in the photogrammetric and remote sensing world. The new ones share a lot of my complaints in that proper resampling using the rigorous camera model is hard. However they show that for the height ranges that stereo modeling is interested in, the epipolar lines can be considered to be linear. Even more interesting, the papers note that all the epipolar lines are parallel [2][3]. They then attempt to solve for an affine transform for both the left and the right images that reduces Y disparity and minimizes image distortion. Solving for the best affine that doesn’t warp the image too much is the primary trick to this task and a few master’s theses have been written on the topic.

Some of the students like to introduce a new camera model, like the Parallel Projection model. I’d like to connect the dots and say that this is not a novel idea at all. In the computer science community this is known as the affine camera model. My favorite book has a whole chapter on the topic [4]. Hartley even goes into a few pages about solving for the affine fundamental matrix and then gives a teaser image about it’s application to stereo rectification for high altitude aerial imagery.

Math

Hartley never goes into how to do affine epipolar rectification in his book. I’d like to document my attempt at this topic as I believe it is a lot simpler than the remote sensing papers make it out to be. It is also simpler than some of the computer science papers make it to be [5].

The general plan of affine epipolar image rectification is best described in 3 steps.

1.)  Identify the slope of the epipolar lines and rotate them so they are horizontal. I’ll be using image correspondences to solve for the affine fundamental matrix, which will get me the epipolar lines and then rotation matrices.
2.)  Align and scale the Y axis of the right image to match the left image.
3.)  Align, scale, and solve for an appropriate skew to minimize the amount of searching a stereo algorithm would have to do in the X axis.

Solving for the fundamental matrix is an algorithm I took from [4]. An affine fundamental matrix differs from a regular fundamental matrix in the sense that its epipoles are at infinity. This makes the top left 2×2 area of the matrix zero. Only the remaining 5 elements contain any information. Solving for those elements we force the constraint that:

When you expand that equation out, you’ll see that the bottom right corner element is not multiplied by any of the image correspondence measurements. So we’ll solve for the other elements first and then we could back substitute for ‘e’ later. But I’m never going to use that element … so let’s forget ‘e’ exists. Solving for the other parameters is achieved by finding the eigen vector with the smallest singular value.

The last row and column of the affine fundamental matrix represent the epipoles of the left and right image.

Let’s use them to solve for the rotation matrices to flatten those epipolar lines. Here was my derivation where I attempted to avoid any trigonometric functions. In hindsight, this makes the algorithm harder to read. At least we avoided a CPU table lookup!

Now we are finally at step 2 where we need to find scaling and offset for the Y axis. This can be performed in a single blow using LAPACK least_squares routine if we just arrange the problem correctly.

Step 3 is solving for the scaling, skew, and offset for X. This too can be performed in a single blow using least_squares.

Again we’ll apply the X scaling and offset only to the right affine transform. However I thought that skew should be applied equally to both images so to minimize distortion. Splitting the skew evenly between both image transforms and how to apply the X scaling is achieved via:

Application to Mars Imagery

I wrote a C++ implementation of the above algorithm and applied it to the example stereo pairs ASP supplies for MOC, CTX, and HiRISE missions. I collect image correspondences for each stereo pair to work out what the search range needs to be in order to create a DEM. The area of this search range is proportional to how much work the stereo correlator would need to do. (Please ignore that in ASP, we use a hierarchal method that provides some optimization making processing time not linearly related to area of the search range.) I then compared the expected search range with alignment options of NONE, HOMOGRAPHY, and the new algorithm. Below are the results.

In the MOC case, homography alignment produced the same results as epipolar alignement. But once the images started to get bigger, it becomes clear that epipolar alignment is a winning strategy. The next test to be performed is to see if the stereo correlator still works on this imagery. We might be introducing skew or scaling errors that the correlator can’t recover from. I have yet to finish that part though.

Another fun note, epipolar rectification is how you produce anaglyphs! Here are few examples of what I was able to do with GIMP. I followed the guide linked here.

So … it looks like I need to put this in ASP 2.2 now. That’s what I’m going to be doing after this blog post.

References:

[1] Oh, Jaehong. Novel Approach to Epipolar Resampling of HRSI and Satellite Stereo Imagery-based Georeferencing of Aerial Images. Diss. The Ohio State University, 2011.
[2] Morgan, Michel, et al. “Epipolar resampling of space-borne linear array scanner scenes using parallel projection.” Photogrammetric engineering and remote sensing 72.11 (2006): 1255-1263.
[3] Wang, Mi, Fen Hu, and Jonathan Li. “Epipolar resampling of linear pushbroom satellite imagery by a new epipolarity model.” ISPRS Journal of Photogrammetry and Remote Sensing 66.3 (2011): 347-355.
[4] Hartley, Richard, and Andrew Zisserman. Multiple view geometry in computer vision. Vol. 2. Cambridge, 2000.
[5] Liansheng, Sui, Zhang Jiulong, and Cui Duwu. “Image rectification using affine epipolar geometric constraint.” Computer Science and Computational Technology, 2008. ISCSCT’08. International Symposium on. Vol. 2. IEEE, 2008.

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.

Previously I stopped the “Autotools applied build system for ISIS” from making the executables because we didn’t have a need for them in ASP. We just wanted to compile against ISIS’s camera models. Tonight however I wanted to use ISIS on some old RHEL5 machines we still have on the network at NASA Ames. Unfortunately, USGS stopped building ISIS for old systems like RHEL5! This blight caused me to now add in application building. Now BinaryBuilder produces ISIS applications and they operate anywhere the compilation does. You can also build “Autotools applied ISIS3″ by hand too.

There is one catch; AutotoolsForISIS doesn’t support any of ISIS’s GUIs. I’m lazy on Sunday and didn’t want to fiddle with the QISIS module or library thing they have. Qview, Qmos, Qtie and the like are thus not built. It also doesn’t build cnethist, hist, phohillier, or spkwriter due to linking issues I haven’t worked out yet. But the important stuff like spiceinit and all the *2isis and *cal applications are there and working.

Images