I’ve developed a novel, automated image registration and georeferencing algorithm and written open source code that implements the algorithm. The algorithm is now mature enough for public use and I will soon (as of March 2018) be submitting an academic journal article that describes the algorithm and characterizes its error. When the article is accepted for publication, I will update with a link to it. In the meantime, this page provides a guide to using the algorithm. Ultimately, I hope to transition the project from Matlab to Python — Matlab is great for rapid prototyping, but neither Matlab nor many of its add-on toolkits (e.g., the Image Processing Toolkit) are free, though they are available at little to no cost at some academic institutions.
Step 1: Get the code
The code is available as a git repository hosted on bitbucket. To get the code from the command line (for Linux and OsX use the terminal; for Windows I suggest using git bash, https://git-scm.com/download/win) use the following command:
git clone https://bitbucket.org/mpatmudd/opengeoreg
This will create a directory called opengeoreg, and the code is in /opengeoreg/code/matlab.
Both the source code and main example script, georef_rudall15.m, are copiously commented. This example script is the best place to start. It is located at:
/opengeoreg/code/matlab/examples/rudall15/georef_rudall15.m
Step 2: Get the Input Files
To run the example, 15 digitized aerial images from Western Australia and a reference image are needed, as well as a “seed image” — one of the aerial images that has already been georeferenced. Download these images from here, then unzip them:
https://www.michaelholtonprice.com/rudall15.zip
The reference file, for example, is located at:
/rudall15/inputs/ref/AGRI_clip_rudall15.tif
Step 3: Run the Example in Matlab
Start Matlab, then issue the following command to add the OpenGeoReg code to the Matlab path (this is not permanent unless you save the path before exiting):
addpath(genpath(path to code))
For example, on a Windows machine I use the precise command is:
addpath(genpath(‘C:\Users\
Open the example script for editing by entering the following command at the Matlab prompt:
edit georef_rudall15.m
Change the base directory, baseDir, to point to the location of the downloaded files on your computer. For example, on a Windows machine I downloaded them to the following location, which is the current value in the example script:
baseDir = ‘C:\Users\mhp12\rudall15’;
Finally, run the example script by entering this command at the Matlab prompt:
georef_rudall15
This will perform the automated georeferencing and place the results in the geodatabase directory located at:
/rudall15/gdb
Details on the algorithm and the files that are created can be found in the comments in the source code.
If you have questions or encounter problems, please don’t hesitate to email me at michaelholtonprice@gmail.com.
Creating a Mosaic in ArcMap
Once the individual images are geo-referenced, a common next step is to create a single mosaic from them. This is surprisingly tricky to do for a couple reasons. First, there are legitimate technical challenges, such as achieving a good color balancing and generating sensible seamlines. Second, the software is finicky, poorly designed, and, worse, often buggy. I tried a variety of approaches and software packages and settled on a procedure in ArcMap that I outline here. It’s nonsense that mosaicking is this hard, but that’s how things are, at least for now.
(1) Create a new project in ArcMap. I called mine rudall15_auto.
(2) Connect to the folder with the images and where the project is saved (e.g., by using the Catalog)
(3) Add the images to the project
(4) Create a new file geodatabase in the project folder by right clicking it in the Catalog, going to New, then selecting File Geodatabase. I named my file geodatabase rudall15_auto.gdb.
(5) Create a new mosaic dataset within the file geodatabase by right clicking on the file geodatabase in the Catalog, going to New, then selecting Mosaic Dataset.
(6) In the dialog that comes up, give the mosaic dataset a name (I chose rudall15_auto in keeping with the previous pattern). Also specify the Coordinate System (for this example, WGS_1984_UTM_Zone_51S). Click OK.
(7) In the Catalog, right click on the mosaic dataset and select Properties. Under Defaults, make sure that Maximum Number of Rasters per Mosaic is greater than the number of images you are mosaicking. The default is 20. Click OK. [Editorial: The first time I mosaicked more than 20 images it took me almost a full day to figure out that this needs to be done.]
(8) In the Catalog, right click on the mosaic dataset and select Add Rasters. Change Input Data to Dataset. Navigate to the images and select them. Click OK.
(9) [If applicable] In the Catalog, right click on the mosaic dataset, go to Modify, then select Define NoData. Set the Bands for NoData Value to All Bands and, in the table, change the NoData value to 0 for All Bands. Click OK.
If you have a more complicated NoData situation, or if you have a mask that defines the active image area, things are much trickier. You will probably have to use ArcMap’s native GRID image format since ArcMap does not correctly implement no data with the geotiffs. The best help I can offer is python code I wrote to convert an image/mask pair in geotiff format to the GRID format:
https://www.michaelholtonprice.com/mask_and_create_grid_images.py
(10) [If you defined a NoData value in Step 9] In the Catalog, right click on the mosac dataset, go to Modify, then select Build Footprints. The defaults should be fine. Click OK.
(11) In the Catalog, right click on the mosaic dataset, go to Modify, then select Build Item Pyramids and Statistics. The defaults should be fine. Click OK. [Confusingly, Enhance -> Calculate Statistics does not create the statistics that are needed for the next steps.]
(12) In the Catalog, right click on the mosaic dataset, go to Enhance, then select Color Balance. You may want to choose different color balancing options, but I will describe what I did for the rudall15 example. For the Balance Method, select Histogram. For the Target Raster, select the reference image (AGRI_clip_rudall15.tif); you may need to connect the folder where it is located and/or navigate to its location. For the Stretch Type (under Pre-processing Options), select Adaptive. Click OK.
(13) In the Catalog, right click on the mosaic dataset, go to Enhance, then select Generate Seamlines. Again, you may want to choose different options, but I will describe what I did for the rudall15 example. For the Computation Method, select Radiometry. Scroll down to the Sliver Removal Options and use .5 and 100 for, respectively, the Minimum Thinness Ration and Maximum Sliver Size. Click OK.
This can be one of the most difficult and confusing points. When mosaicking large number of images especially, the seamline generation often will mysteriously fail. Alternatively, the seamline generation may finish, but the seamlines are not visible and cannot be removed from the Catalog (which implies that they were never really created, when in fact they were). I find it sometimes works to save, exit, and restart. You will know if this is successful by right clicking on the mosaic dataset and going to Remove. If Remove Seamlines is grayed out, you’ll have to do further de-bugging. Sometimes just repeating the whole process up to this stage works. It’s all very mysterious (and frustrating).
(14) In the Catalog, right click on the mosaic dataset and select Properties. Under Defaults -> Allowed Mosaic Methods deselect everything except Seamlines and None. Click OK.
(15) In the Catalog, right click on the mosaic dataset, go to Export, then select Raster to Different Format. In the ensuing dialog, change the format to TIff. Select the location and name of the output, mosaicked raster by editing the Output RasterDataset field. For example, I named the final mosaic rudall15_auto_radiometry.tif.