Fast Image Mosaicking using Multicore CPUs

In January 2007 I started working on the CZMIL project which is multiyear project involving the University of Southern Mississippi Department of Marine Science, Optech International, the US Army Corps of Engineers, and the Naval Oceanographic Office. The primary purpose of the CZMIL project is to develop a new Lidar collection system in support of the Joint Airborne Lidar Bathymetry Technical Center (JALBTCX) located at Stennis Airport in Kiln, MS. A secondary purpose of this project is to streamline the ongoing data collection and processing at JALBTCX. My work has been in streamlining the production of mosaics of digital camera images collected along with the Lidar data.

JALBTCX is a collaborative center for Lidar processing, research and development housed at the Stennis Airport. Government partners at JALBTCX include the US Army Corps of Engineers (USACE), the US Naval Meteorology and Oceanography Command, the National Oceanic and Atmospheric Administration (NOAA) and the Naval Oceanographic Office. JALBTCX conducts surveys using the Compact Hydrographic Airborne Rapid Total Survey (CHARTS) system which includes topographic and bathymetric lidar, a hyperspectral camera and a digital RGB camera.

Previous Mosaic Processing

The digital camera system in the CHARTS system is a DuncanTech DT4000 color camera. The system collects approximately 1 1600x1200 pixel color image per second as lidar data is collected. The original purpose for the camera was to assist in determining what features cause the responses in the lidar data. Seeing an image of the water or ground cover can quickly eliminate some confusion in lidar points.

Over time JALBTCX users started producing RGB orthomosaics to provide an overview of survey areas. These mosaics were produced using commercial software which required about 3 hours to produce a mosaic with about 600 input images. In addition there were size limits imposed which forced using coarse resolution elevation data as part of the process. Overall the quality of the mosaics was not great and the time was a processing bottleneck.

The problems with the previous mosaic solution were a result of performing orthorectification of 1 RGB image at a time with a subsequent mosaicking process. Each image processed resulted in the separate execution of external processes. The net effect was a fairly low CPU usage (usually under 25%) and slow processing.

A More Efficient Solution

One way to improve efficiency is to perform orthorectification and mosaicking as a single process, which reduces I/O requirements. Another issue is to avoid mapping output pixels multiple times. There is substantial overlap (~50%) between successive images and overlap between flight lines. When separate orthorectification is performed, each input pixel is mapped into the output space. This means that over half of the pixels produced in these separate orthorectified files will be overwritten when the mosaic is prepared.

The mosaic program uses a Voronoi image as an integral tool for improving efficiency. First the centers of each input image are mapped to their output UTM coordinates. From these centers a distance is calculated for each output pixel to the nearest image center. This is shown in the gray-scale image below. The image centers are located at the center of the dark circles. At the same time that the distance image is computed a Voronoi image is also computed. The Voronoi image (in color below) indicates which original image center is closest to a particular output pixel. In the Voronoi image below the mapping would start with pixels being mapped from the image identified by the yellow polygon. After about 50 pixels from the yellow image, pixels would be mapped from the blue image. Thus the Voronoi image tells the mosaic program which input image to use for each output pixel.

There are several improvements which result from using the Voronoi image as a guide in the mosaic process. As previously stated this allows mapping each output pixel exactly once, eliminating redundant calculcations. A second benefit is that the pixels mapped tend to be from the central parts of input images which improves the accuracy of the mapping. Finally it is fairly easy to divide the output image into tiles which can be processed in parallel taking advantage of the availability of multi-core CPUs.

Parallel Processes

The mosaic program uses 5 different types of processes to perform its work: a single master process, a DEM reader process, multiple RGB image reader processes, multiple mapper processes and a single writer process. In the initial version of the parallel mosaic program these processes are initiated by performing fork calls and executing specific functions. Due to Windows not directly supporting the fork system call, the software is being redesigned to use processes and threads, but the overall strategy is still much the same.

Master Process

The master process reads a parameter file which contains specific parameters such as camera boresight and lens values, names of other files, and a variety of processing parameters. The master reads a file containing the names and aircraft attitude of the RGB input images to be mosaicked. Then it computes the Voronoi image and determines the sequence for reading the input images and how many must be in memory to avoid deadlock. Then it starts the other processes and enters a loop where it services requests arriving over pipes.

DEM Reader

The DEM reader process reads a single digital elevation model image (32 bit floating point values) into a shared memory array. After it completes reading the file, it sends a message to the master announcing completion so that the master can allow the mapper processes to verify that the elevation data is available for processing.

RGB Readers

The RGB readers request from the master the image number of the next image to read. The master delays sending the reply until the shared memory region for the file is available for use. There are typically under 100 shared memory regions for RGB images though there can be many more images to read. The master knows when images are no longer needed and, thus, when their memory can be reused. When a reader completes reading an image, it sends a completion message to the master so that the mapper processes can access that image after it is ready. Eventully the master sends -1 image numbers to the readers to tell them to exit.

Mappers

There are multiple mapper processes - typically as many as the number of CPU cores in the system. Each mapper works on 1 tile at a time. A reasonable tile size is perhaps 512x512. The mappers request tile numbers from the master allowing the master to control the use of the memory for the tiles. After a mapper completes a tile, it sends a completion message to the master which ultimately allows the writer to write that tile to the output file. Eventually the master sends -1 tile number to the mappers to signal the end of processing.

Writer

There is a single writer process which writes the output file one tile at a time. The writer requests tile numbers from the master, writes the tiles and signals completion. Signaling completion of tile writing allows mappers to reuse tiles for different parts of the output image.

Modeling

The basic model used is a 3 dimensional affine transformation from pixel coordinates to UTM coordinates. The resampling uses the inverse of this transformation so that each output pixel can be mapped to an input pixel position. The output pixels have northing and easting values based on their location in the output file and elevation values from the DEM image. Typically the resolution of the DEM image is coarser than the output image, so bilinear interpolation is applied to elevation values. The northing, easting and elevation values are transformed into pixel coordinates and interpolation is applied to the input pixel values.

The mosaic program implements a radial lends distortion correction as part of the mapping. This correction is applied by precomputing a 2 dimensional array of scale factors which are multiplied times the pixel coordinates before adding the principal point values. The distortion is supplied as a set of coefficients for a 4th degree polynomial. As part of the precomputation the polynomial inverse is approximated by a 7th degree polynomial over the practical range of the polynomial (limited from 0 to 1000 pixels). The scale factor array is computed using this 7th degree inverse polynomial.

The user can select geographic (latitude/longitude) output coordinates if desired. In that case the actual mapping is still performed using UTM coordinates. The output geographic coordinates are converted to UTM prior to mapping. Since this conversion is non-trivial, every 10th row and column are converted using the PROJ4 library and interpolation is used to get accurate approximations during the mapping.

Results

The test area for this report is over Pickering, Ontario. The input consists of 608 RGB images. Most of the flight lines are approximately north-south or south-north, though some are east-west or west-east. A sample of the output image is shown below.

Performance

Performance testing was performed on an 8 core AMD Opteron computer. This system has 4 dual-core CPUs operating at 2.4 GHz. The system has 32 GB of RAM and runs a 64 bit version of Fedora Linux. The disk system uses SATA disks without RAID.  Linux allowed enabling and disabling CPUs as desired during testing. In addition Linux allowed flushing the disk cache so that each test properly reflected the time spent reading input images. Below you can see elapsed time and speedup graphs based on using 1 to 8 cores. In this test the speedup is nearly ideal through 4 cores. After 4 cores the system is probably limited by the single writer process which produced a compressed GeoTiff file. Eliminating compression would improve efficiency beyond 4 cores. Overall the time required is less than 1 minute for a task which had previously required 2 to 3 hours to complete.