Open Access
Add to BibSonomyAbstract
A distortion mapping and computational image unwarping method based on a network interpolation that uses radial basis functions is presented. The method is applied to correct distortion in an off-axis head-worn display (HWD) presenting up to 23% highly asymmetric distortion over a 27°x21° field of view. A 10−5 mm absolute error of the mapping function over the field of view was achieved. The unwarping efficacy was assessed using the image-rendering feature of optical design software. Correlation coefficients between unwarped images seen through the HWD and the original images, as well as edge superimposition results, are presented. In an experiment, images are prewarped using radial basis functions for a recently built, off-axis HWD with a 20° diagonal field of view in a 4:3 ratio. Real-time video is generated by a custom application with 2 ms added latency and is demonstrated.
©2012 Optical Society of America
1. Introduction
In areas where images are used for data extraction and interpretation, it is often necessary to correct the distortion present in the images or, at least, obtain accurate knowledge of the relationship between any point in the real object and its counterpart in the image. In the context of photogrammetry and remote sensing, photographic data acquisition devices need to be calibrated in order to obtain quantitative information about 3D-objects [1]. In computer and machine vision, inferring precise 3D information about the real world from computer images requires accurate calibration of a camera’s extrinsic and intrinsic parameters [2]. In virtual reality, distortion correction of the projected images is needed for reliable use of equipment [3]. In the context of diagnostic imaging, medical images (i.e. X-ray, endoscopic) also undergo distortion and require distortion-correction [4]. Finally, unwarping algorithms are commonly applied to interferogram processing in optical testing.
Optical distortion introduced by imaging systems includes the Seidel aberration and higher order terms [5], but also decentering and tilt effects such as keystone [6]and sigmoidal distortion [4]. Design constraints, such as projection on planar surfaces, or non-constrained distortion in compact head-worn displays (HWD) also create the need for distortion compensation [3,7]. Thus, an adequate model of the distortion created by the imaging system (or other physical processes) is critical to ensure usability. A pair of points consisting of an object point and its corresponding distorted point in the image is called a correspondence, of which a set is illustrated in Fig. 1 . The search for a distortion model of a given system generally aims to model a set of correspondences over a desired field of view (FOV).
In this paper, we first review known approaches to distortion mapping. In Section 3, we detail a method for radial basis function (RBF) based distortion mapping and apply it to the distortion of an off-axis HWD. Results of image unwarping based on an image-rendering feature of optical design software are presented along with a quantitative analysis of the method. In Section 4, the RBF distortion mapping method is applied to a compact eyeglass format off-axis HWD, yielding an accuracy of less than 1 RMS pixel displacement error, and motivates the use of RBFs for distortion mapping of as-built optics. Finally, in Section 5, we demonstrate real-time distortion mapping for laboratory hardware.
2. Overview of distortion mapping methods
Three main distortion-mapping methods can be found in the literature: physical-model based mapping methods, global methods, and local methods [4]. Physical-model based methods model the physical phenomena occurring during the imaging process, namely optical distortion produced by the imaging system and sigmoidal distortion in radiographic images. Optical distortion is typically decomposed into two components—pure radial distortion with respect to the optical center for rotationally symmetric systems and non-radial distortion, also called tangential distortion, for rotationally nonsymmetric systems. Radial component mappings include polynomial expansions used in conventional optical aberration theory [5], computer vision, and photogrammetry [2,3,7,8]. More specific mappings exist to model high distortion amounts where polynomial mapping is not sufficiently accurate, such as the “fish-eye transform” (FET) model [9] and the FOV model [10] for fish-eye lenses and computer vision wide-angle, omnidirectional cameras [11]. Non-radial distortion models comprising both radial and tangential components are not yet frequent in the literature [3] and are either used to enhance mapping accuracy by accounting for imperfect rotational symmetry of imaging systems, or more rarely, when rotationally nonsymmetric optical systems are considered such as in [12] and [13].One of the few tangential distortion models found in the literature is a power series expansion of two contributions referred to as “decentering distortion” and “thin prism distortion”. The expansion to a2nd-order power series in the perfect Cartesian coordinates is frequently used in the photogrammetry or computer vision and graphics areas [8,14–16] and is based on past studies of tangential distortion causes in optical systems [17,18].
Global-mapping methods are proposed as an alternative to physical-model based methods to increase mapping accuracy when the distortion type is not known a priori or when it is difficult to match to a single physical distortion model. Global methods consist of fitting the system distortion with an analytic function comprised of non-physical parameters. Polynomials can be used in global analytic methods [4,11]. In the remote sensing or camera calibration fields, rational functions are frequently studied and efficiently used for distortion mapping [19, 20, 11]. Several special cases were found to be suitable for a specific application: rational polynomial functions, up to 3rd-order, for compensation of computer vision cameras [20], cubic rational polynomials for modeling complex cameras used in remote sensing applications [19], and quadratic rational function based-distortion mapping for wide-angle and catadioptric optical systems [11].
One accuracy limitation of physical-model based and global analytic methods is their inability to model local distortion deformations; for instance, noise in the correspondences measurements or local distortions intrinsic to the system, such as misalignment effects. Alternative models are local mapping methods, which use sets of local patches and interpolation algorithms widely used in the image warping field to model image deformation to arbitrary shapes. Specifically, scattered data interpolation techniques are used to define interpolation functions over a finite set of irregularly located data points, efficiently constraining correspondences while interpolating in-between points [21].
The local distortion mapping technique that uses basis function networks is investigated in this paper. The mapping function is constructed as a linear combination of basis functions (commonly radial-symmetric). Examples of RBFs are thin plate splines [22],Gaussians [4,23], or multiquadrics [21]. Basis function network mappings are not limited to a specific type of distortion and are reported as presenting better distortion correction accuracy when compared to global analytic methods either in the context of TV cameras calibration [23] or remote sensing [25] and perform as well as global polynomials when local distortions are present, as for X-ray distortion correction [4,22].
Depending on the application, the distortion mapping process connects distorted image coordinates to perfect input coordinates or vice-versa. Once a model is chosen, its parameters (e.g. physical parameters, global analytic models parameters, basis functions weights) are determined by solving the linear system composed of the equations linking correspondences pairs (objects points and distorted image points) through the model. Advantages and shortcomings of a model depend on the considered application, such as required mapping accuracy or degree of distortion unwarping.
3. Principle of the RBF-based distortion mapping method and its application to a computational unwarping task
The distortion mapping method based on RBFs will be illustrated using the optical design software (CODE V®) model of an off-axis eyeglass display [12]. A 2D-layout and distortion grid over the FOV is shown in Fig. 2 . Computational correction of the HWD distortion consists of prewarping images to display them on a microdisplay in a way that counteracts the distortion of the optical system. In other words, an image on a microdisplay is prewarped, or pre-distorted, such that when viewed through the optical system, will appear unwarped. The distortion correction can be decomposed into two steps: 1) search for a warping function that approximates the distortion function of the HWD system when traced from visual space (see black rectangular grid in Fig. 2(b) that represents a grid in visual space) to the microdisplay (see the reddish distorted grid in Fig. 2(b); each image displayed on the microdisplay will need to be prewarped to this distortion map)and 2) using this approximation function, prewarp the images to display. These two steps are presented next, followed by an assessment of the unwarping process as further detailed in section 3.3.
Fig. 2 (a) 2D layout of the starting design for the off-axis HWD under study [12].(b) CODE V distortion grid computation.
3.1Optical system distortion mapping
The distortion pattern of the HWD system shown in Fig. 2(b) shows a large (i.e. 23%) tangential distortion component. Per the light reverse travel law, imprinting any prewarped image following this distortion pattern on the microdisplay would result in the observation of the undistorted pattern by the observer's eye in visual space. Using the distortion grid feature output as data, the desired distortion mapping is the one linking object points of the system and their distorted counterparts. The FOV to be considered when running the distortion grid feature should be chosen so that the distorted pattern fits exactly in the microdisplay under study. However, such an exact definition of the working FOV will actually exhibit low mapping accuracy at the edges due to inherent edge effects occurring in the interpolation processes. Thus, a 25%-wider FOV is used to increase the overall mapping accuracy in this step of the process.
For nequal212 points in visual space (evenly positioned on a rectangular grid satisfying the 25% oversized FOV requirement), distorted point coordinates were computed by the Distortion Grid option in CODE V to give the correspondences needed for the determination of the mapping parameters. Normalization to the horizontal FOV maximum value was performed for computation convenience. Using the local RBF-based mapping method with N multiquadric basis functions [21], each correspondence between a point on the regular grid M’j(xj’,yj’) and its distorted counterpart Mj(xj,yj) is written as follows:
where, Ri represents the ith basis function centered at (xbasis_center_i, ybasis_center_i), with basis function weights αxy,i; pm(xj, yj) is a polynomial of degree m that assures polynomial precision of degree m; j is an integer that runs from 1 to N; and χ is a scaling factor to be chosen heuristically based on the specific application. Figure 3 shows a graphical representation of the algorithm.Fig. 3 A visual representation of the local RBF-based distortion mapping method. An RBF is located at each correspondence, j (bottom). Using the correspondences and Eq. (1), a weighted amplitude is given to each RBF (middle) to create the top interpolated surface. The interpolated surface (top) represents all the possible x-coordinate values in the distorted image for any pixel (x,y) in the original image. A separate, yet similar figure must be made for the distorted y-coordinates using Eq. (2).
Here, the number of basis functions was taken to be equal to the number of correspondences (n = N).The basis centers were placed at the original grid points location. The characteristic radii (ri)i = 1..n of all basis functions were taken to be constant and equal to the minimum spacing between 2 original points on the grid. The multiquadric exponent,µ, was set equal to −2, as in [21]. The 2N-equations linear system formed by Eqs. (1) and (2) repeated for the N correspondences was solved by least squares minimization to find the basis functions weights (αxi)i = 1..n and (αyi)i = 1..n. Once learned, those weights completely define the RBF-based mapping of the HWD distortion. The RBF distortion mapping is compared to the distortion grid output in Fig. 4 yielding accuracies higher than 3 μm across a high density of test points within the effective FOV, which we found to be more than 100 times better than that obtained with the next best model discussed by Weng [8].
Fig. 4 Comparison of data and RBF modeled distorted points showing normalized x coordinates along the horizontal axis and the normalized y coordinate along the vertical axis: (a) Set of N correspondences only; (b) Set of Ntest = 412 test point including correspondences. The resulting RMS and maximum absolute error between data points and RBF model-computed points is 4 × 10−17 mm and 4 × 10−15 mm, respectively, considering only correspondences; 4 × 10−15 mm and 0.058 mm, respectively, considering Ntest = 1012 test points; 10−5 mm and 0.003 mm considering all Ntest test points, but only over the effective FOV.
3.2 Image Prewarping Using RBF-based distortion mapping
An image prewarping algorithm based on the previous RBF-based distortion mapping function was developed in MATLAB, considering images with resolution (Nyx Nx) (rows by columns) to be displayed with the HWD from [12]. Two algorithms were tested, one faster (method 1), but suffering from input image pixel loss, and the other more accurate (method 2), using oversampling to prevent such pixel loss. For both methods 1 and 2, the prewarping of an input image consists of building a storage matrix of size (Ny_storage x Nx_storage), from the (Ny x Nx) input image pixels and forming the final prewarped image from it.
For method 1, Nx_storage equals Nx and Ny_storage equals Ny; for method 2, Nx_storageequals3Nx-2 and Ny_storageequals3Ny-2. The 3-oversampling ratio in method 2 was heuristically determined to fulfill the total pixel retrieval condition for the specific distortion under study. The (Ny x Nx) input image pixels are evenly positioned in this storage matrix; identical to the input image for method 1, and re-indexed in [1, 3Ny-2] x [1, 3Nx-2] intervals for method 2. Before warping them using the RBF-based mapping, the (Ny x Nx) pixels in the storage matrix are re-indexed into real coordinates through affine operations [24]in order for the warping process to be performed in the same domain as the RBF mapping, with the central image pixel located at (0,0) and the extreme pixels abscissa and ordinate equal to ± 0.8 and ± 0.8simage respectively, with simage being the input image aspect ratio (simage = 3/4), and the 0.8 factor related to the use of the 25% oversized FOV during the mapping weights determination. All of the (Ny x Nx) pixels with new indices are then warped following the RBF-based mapping of Eqs. (1) and (2) with the previously learned weights, (αxi)i = 1..n and (αyi)i = 1..n. They are then re-indexed back to the storage matrix indices domain [1, Ny_storage ] x [1, Nx_storage], also through affine transformations defined so that the extreme warped pixels exactly correspond to the storage matrix edge pixels (there are extreme rows if the aspect ratio of the distorted pattern formed by the warped pixels is smaller than the storage matrix aspect ratio, extreme columns if it is larger).
After those operations, the storage matrix contains a number of pixels smaller than (Ny x Nx) for method 1, and equal to (Ny x Nx) for method 2 (by definition of this method). The precise number of pixels in the storage matrix for method 1 depends on the specific distortion mapping. In this case, the (Ny x Nx) input pixels are not moved to (Ny x Nx) distinct locations due to the existence of pixel overlaps. This occurs when initially adjacent pixels are moved closer than 1 pixel apart by the distortion function, and thus, have identical distorted coordinates. Similar to these pixel overlaps, holes occur between pixels that are moved further than 1 pixel apart. These holes are corrected with a filling algorithm. After the hole-filling correction step, the storage matrix contains the desired hole free Ny_storagex Nx_storage image, which has been prewarped following the desired distortion mapping. Figure 5 shows, as an example, an input grayscale image (a grid) with SVGA resolution (i.e. Ny x Nx equal 600 x 800) and its corresponding prewarped image for both methods. Uniform gray pixels were added around the effective distorted image during the prewarping process for further assessment convenience of the unwarped images. Smoother prewarped images in the case of method 2 may be observed in Fig. 5.
Fig. 5 Images showing RBF-based distortion mapping. (a) Input Image (SVGA resolution), (b) prewarped image (method 1), (c) prewarped image (method 2).
3.3 Computational Distortion Unwrapping Simulation
To assess the distortion unwarping resulting from the previous RBF-based distortion mapping and corresponding image prewarping, the visualization process through the HWD was simulated in CODE V using the 2D IMS option for the direct HWD model (i.e. microdisplay in object space) that computes the appearance of a 2D input object imaged through an optical system. The input object consisted of an (Ny x Nx) equal 600x800 prewarped image filling the display FOV (18.5 mm-diagonal) and translated to account for the reverse system (the visual space for the reverse system is the object space) chief ray offset in the prewarped image. The resulting rendered images from methods 1 and 2 are shown in Fig. 6 .
Fig. 6 Resulting unwarped images (a) Original images. (b) Unwarped images (method 1). (c) Unwarped images (method 2).
The black areas around the rendered images were added by the 2D IMS option and do not correspond to any real features. The unwarping efficacy of the initial image pattern (encompassed inside the gray border) can be visually observed. Figure 7 compares initial images (downsampled) and final unwarped images (cropped around the effective area) for improved visualization; the main visual difference between methods 1 and 2 is the greater smoothness of edges in method 2 compared to method 1.
Fig. 7 Side-by-side comparison with original images. (a) Input images after downsampling and cropping to unwarped image size. (b) Unwarped cropped images (method 1). (c) Unwarped cropped images (method 2).
Direct pixel-by-pixel comparison between the original image and the RBF-unwarped version is a challenging quantitative metric due to differences in their respective generation. The unwarped image, computed from the original image, experienced an averaging process during the prewarping step in the case of method 2, and the 2D IMS option affected it by both modifying its pixel gray level caused by the system vignetting and resampling it. Thus, the original image to be used for a cross-comparison needs to be downsampled to this new resolution and modified to have the same gray level. In addition, such a pixel-to-pixel comparison should be considered as a “case study”, as in the case of the computational distortion of a visual instrument such as the presented HWD.
The unwarped images (A) were quantitatively compared to the original images (B) using the correlation coefficient (CC) given by
Where N is the number of pixels, Ai and Bi are the pixel values of the ith pixel, and µA,B are the mean values of A and B. A value of CC equal 1 denotes perfect correlation between two images. The unwarped images resulting from the 2D IMS option were cropped along the uniform gray border enclosing the original image pattern and corrected from the system vignetting effect using the rendered image of a uniform white rectangle. The original image was downsampled in MATLAB to the cropped unwarped image resolution using the imresize function.The obtained values for the correlation coefficient are listed in Table 1 . The results for the four tested images indicate a higher similarity to input images for images obtained with method 2. The absolute values of the CC were found to be lower for images with a larger proportion of finer edges (“Grid 1” and “Grid 2”), which can be explained by the higher impact of all processing steps on the high frequency features location and gray level computation. Unwarped “natural” images (“Landscape 1” and “Landscape 2”) present CC-values close to the perfect 1 value.
Table 1. Correlation Coefficients between CODE V 2D IMS UnwarpedImages and Original Images
As an additional quantitative assessment of the distortion unwarping process efficacy, an edge detection algorithm (Canny method [26]) was applied in the case of the “Landscape 1” and “Landscape 2” images to both input and unwarped images to compare the location and shape of their edges. The results were superimposed as two components (red for the original image edges, green for the unwarped image edges) of an RGB image, so as to observe coincidence (in yellow) and differences of the edges. The edge detection results are similar between methods 1 and 2, and the RBG images for method 2 are presented in Fig. 8 . Overall superposition of the input and unwarped images edges can be observed. Some edges present a 1-pixel shift that might be attributed to the different methods of generation of the input and unwarped images, affecting the edge detection algorithm precision. A maximum edge location difference of 3pixels was found for the “Landscape 1” and of 2 pixels for the “Landscape 2” images.
Fig. 8 Edge comparison between original and unwarped images using method 2 for(a) Landscape 1 and (b) Landscape 2.
4. Prewarping an image for an off-axis HWD
In order to detail the accuracy of the RBF-based distortion mapping, an image was prewarped for the off-axis, eyeglasses-format HWD shown in Fig. 9 , which was designed and built by Optical Research Associates (ORA) (now Synopsys Inc.) [27]. The HWD is an off-axis design consisting of tilted and decentered elements. The distortion map for this HWD was provided in the form of a table of correspondences that was created by painstakingly mapping each individual pixel’s origin and destination in ray trace software in order to establish a gold standard for comparison. We will show that, using an RBF-based distortion mapping, we are able to emulate (with < 1pixel RMS displacement error) the distortion of the HWD by sampling only a small fraction of the actual correspondences and generating the remaining correspondences using the RBF algorithm.
Fig. 9 An optical see-through, eyeglasses-format, head worn display [Courtesy of Optical Research Associates].
The correspondences coordinates are the source of the RBF centers, and thus, serve the same purpose as the data from the Distortion Grid function in the optical design software used in the previous section. Note that the Distortion Grid function only generates a small sample of correspondences from discrete points in the FOV. From the entire set of correspondences established to get the gold standard, we select a small subset of them and aim to predict the entire set with the RBF mapping function. As previously mentioned, there are inherent edge effects that occur as a result of RBF interpolation, so one must cleverly select the subset of correspondences.
Two separate selection algorithms were implemented and combined to create the full subset. The first algorithm chooses only points around the edge of the image and in the corners. This is crucial because, as explained in section 3, there are inherent edge effects due to the RBF interpolation process. The second algorithm uniformly selects points throughout the body of the image. The combination of these two selection algorithms provides a sufficient subset of correspondences for the RBF distortion process. Using this subset of correspondences, the RBF weighting coefficients are generated using the process previously outlined and a full table of correspondences is created using those weights together with Eq. (1), (2), and (3). This one-time process takes under 5 seconds on a standard PC. The resulting image can be seen in Fig. 10(c) . For comparison, the image created by the original (gold standard) correspondence table is shown in Fig. 10(b). As mentioned earlier, χ in Eq. (3) is chosen heuristically per application. Simulations showed that χ equal 24 yields the most accurate results, as quantified in Fig. 11(b) .Results also show that the process is not sensitive to the exact value of χ around 24.
Fig. 10 (a) Original input image, (b) Prewarped image using the exact (gold standard) correspondence table, (c) Prewarped image using an RBF-generated correspondence table,(d) Prewarped image showing the location of the sampled correspondences.
Fig. 11 (a) Plot showing RMS pixel displacement error for the RBF distortion method for subsets containing different amounts of correspondences. There is no accuracy gain for sampling more than 120 points. (b) Plot showing RMS pixel displacement error for the RBF distortion method for various values of the parameter, χ, in Eq. (3). The optimal value is 24 and the sensitivity to a small change is low.
To compare the two prewarped images, the correspondence table generated by the RBF process was cross-checked, pixel-by-pixel, with the gold standard correspondence table. For each pixel, the straight-line distance was calculated from the pixel’s theoretical destination to its actual destination (ideally, that distance would be zero). That distance was then squared and a running sum of these squared distances was computed for all pixels. To complete the RMS calculation, the squared-distances-sum was divided by the number of pixels, and then a square root was taken. For a subset of 121 correspondences, a 0.805 RMS pixel displacement error was obtained. Figure 11(a) shows results for simulations using subsets of different numbers of correspondences. It can be seen that the RMS pixel displacement error decreases as the number of correspondences increases, until ~120 correspondences where there is no gain in accuracy from adding correspondences. As more correspondences are added, the points become more densely packed and reach a point at which the basis functions experience a high degree of overlap, thus decreasing the effectiveness of each, or, more technically, the RBF interpolant becomes computationally unstable for higher numbers of correspondences [28].
This result plays a role of significance for designers who only need to provide an 11 x 11 table of distorted points together with their final design specification. This result, however, plays a role of even larger significance for as-built optical systems. The optical design is often not available for a given system, so if one were required to know the distortion characteristics of an optical system that knowledge would need to be gained experimentally. We showed that it is possible to map out a complex distortion over the full FOV of an optical system by sampling the distortion at a small subset of points over the FOV. Thus, using the RBF distortion mapping process, many fewer measurements need to be taken to completely model the distortion of an as-built optical system.
5. Real-time distortion mapping for use with HWDs
For a HWD displaying the large degree of distortion, such as the one shown in section 4, real-time prewarping is necessary to avoid requiring users to view only content that has been specifically prepared for use with such a HWD; without real-time prewarping, watching prerecorded and processed video would be, essentially, the only application. The key component to any real-time computational process is computation speed, which is why a GPU accelerated application, which utilizes parallel processing techniques, is required. A GPU accelerated prewarping application was built for use with the HWD showing the distortion shown in Fig. 10.
Before delving into the implementation details, it is important to have a basic understanding of how computer video rendering works (i.e. how pixels are moved from storage or capture data to the stream leaving the video card). Video data is commonly stored and read as YUV encoded pixels. Since the hardware that drives computer monitors expects RGB pixel values, the YUV encoded pixels need to be re-encoded as RGB pixels, or transcoded. The final stage is to write the pixels out to the video hardware's frame buffer and let the hardware render the video stream to the monitor.
The past twenty years have seen great advances towards making this process fast enough to make latencies insignificant for standard video sizes. Various encoding formats were developed, shrinking data storage requirements while adding another transcoding stage. Multi-threaded operating environments provided multiple processing lines for distributing the workload, first with multiple processors and later using multiple core processors. Data transfer architectures were developed to transfer data from devices to memory and back again without taking CPU time. Some of the largest gains in the past decade have come from video hardware development and the rise of both the parallel processing capabilities and the ability to program the rendering pipeline in the GPU.
Our algorithm takes advantage of many of these advances while using features that have been available on video cards for at least six years. We used the video library libVLC, which has been optimized for rendering video. It allows us to interrupt the process before the pixels are sent to the computer monitor and enables us to write those pixels to the texture memory of the video hardware.
To ensure generality of the application, the only required inputs are a video source and a correspondence table describing the distortion. The basic approach to implementing the correspondence table is to create a distortion mesh with texture coordinates, to then render a scene or video (sequence of scenes) directly to texture memory on the video card and, finally, to bind that texture to the warped mesh. Rendering the scene into texture memory has the added benefit of handling the color-sampling problem that can occur when minifying textures and allows the user to choose between a bilinear color filter and a (Manhattan) distance filter. The correspondence table is only read once during initialization, and the memory it occupies can be freed immediately after it has been processed.
There are a number of approaches that one could use to leverage OpenGL accelerated hardware on the video board to implement the correspondence table, but the fastest approach uses a framebuffer object (FBO) to render the scene or video directly to texture memory. FBOs work better in this case than pixel buffer objects (PBO), as one has to explicitly copy to texture with the PBOs (glCopyTexSubImage2D). It should be obvious that using the GPU to handle the color sampling problem beats generating a pixel by pixel drawing of the distortion while preserving as much detail as possible from the original scene. Another benefit of using a FBO is that one can create arbitrarily sized regions of texture memory as opposed to being restricted to rendering into a window with dimensions that are powers of two, which can cause a loss of information if the scene needs to be minified to satisfy the power-of-two constraint.
The correspondence table is used to generate a mesh with the appropriate texture coordinates. We then sort the table according to which eye the correspondence pair relates and according to the original image space coordinates. We use this sorted table to iterate, at regular intervals, through the original image space, lookup the corresponding pixel location in the final image (distorted) space, and add that location to the mesh. This is a one-time only process that occurs during initialization and lasts ~5 seconds for a ~800 x 600 image size.
During the rendering loop, we render a scene or video to the FBO and then make a call to the display list for the distortion mesh. The loop runs as fast as the CPU can render the scene to texture memory, as the number of triangles involved is low in the final rendering stage. If one were streaming a video to the warp display, we would only need to stream into a FBO that is the size of the original video instead of generating a pixel-by-pixel distortion in the final window. The prewarped video plays at the same frame rate as the input video source. The only system requirement of this application is a video card that supports FBOs (which the majority of video cards manufactured since 2005 support). A screenshot and associated video of the finished application can be seen in Fig. 12 .
Fig. 12 Custom real-time distortion correction application. While the application itself is designed for use with any generic distortion map, the images shown correspond to the ORA HWD detailed above. The input video is split into the separate outputs video feeds for the Left and Right eyes with the correct distortion applied. The right eye output is flipped vertically because it was necessary for the micro display in the HWD to be inverted to comply with constraints of the electronics. A video clip of the application is presented (Media 1).
Computation speed is a key factor for real-time distortion correction. With our algorithm, for 1280x700 video, only ~2ms are added to the normal workload. For NTSC video, a new frame is presented approximately every 30 ms and for 60 Hz video, a new frame is presented every 15 ms. For each of these video frame rates, the added latency due to the distortion correction (~2 ms) is negligible. Also, as stated previously, the correspondence table is loaded into memory only at initialization and is cleared from memory thereafter, resulting in a computationally efficient process.
6. Conclusion and discussion
A distortion mapping method based on the use of multiquadric RBFs was presented and applied to a computational distortion correction task of an off-axis HWD computer model. The novelty of this method resides in its ability to map complex off-axis optical distortion with knowledge of distortion at only an 11 x 11 grid points to an accuracy of less than one pixel RMS pixel displacement error. It also enables the determination of the distortion of an as-built optical system with only a few measurements (i.e. 121 points). This ability is significant for a system with off-axis distortion, as many measurements typically need to be taken throughout the FOV to accurately describe the distortion. With the RBF method, the number of measurements required is minimized, also decreasing measurement-induced errors.
The corresponding mapping coupled to a filling algorithm was used to prewarp digital images to display with the HWD by a final method that is fast (i.e. seconds at initialization and ms in real-time runs) and accurate. Also, the efficacy of the computational prewarping process was evaluated using a rendering feature integrated in optical design software. Quantitative assessment of the resulting distortion corrected images provided correlation coefficients between the initial and final images higher than 0.77 for grid images with exclusively high-frequency content and greater than 0.99 for “natural” images. A prewarping case study was performed on actual hardware, the ORA HWD. Using 121 “real” correspondences, the RBF-based distortion method produced an image with an RMS pixel displacement error of 0.8 pixels. Finally, using a GPU accelerated application, real-time prewarping of input signal was demonstrated in ~2ms, making optical hardware with large and non-symmetric distortions, such as off-axis eye-glass format displays, available for use in a variety of real-time applications including Augmented Reality.
Acknowledgments
This work was supported by the National Science Foundation GOALI grant ECCS-1002179, The NYSTAR Foundation C050070, a Carl Zeiss Fellowship to Aaron Bauer, a Provost Multidisciplinary Fund grant, and AUGMATE Corporation. We thank Paul Glenn of Bauer Associates for creating the gold standard lookup table for the ORA HWD against which we compared our method. We thank Synopsys Corporation for the loan of the off-axis HWD employed in the experiments, as well as for the student license of CODE VTM.
References and links
1. W. Faig, “Calibration of close-range photogrammetric systems: Mathematical formulation,” Photogramm. Eng. Remote Sensing 41, 1479–1486 (1975).
2. R. Y. Tsai, “A versatile camera calibration technique for high-accuracy 3D machine vision metrology using off-the-shelf TV cameras and lenses,” IEEE J. Robot. Autom. 3(4), 323–344 (1987). [CrossRef]
3. W. Robinett and J.P. Rolland, “A computational model for the stereoscopic optics of a head mounted-display,” Presence (Camb. Mass.) 1, 45–62 (1992).
4. P. Cerveri, C. Forlani, A. Pedotti, and G. Ferrigno, “Hierarchical radial basis function networks and local polynomial un-warping for X-ray image intensifier distortion correction: a comparison with global techniques,” Med. Biol. Eng. Comput. 41(2), 151–163 (2003).
5. W. T. Welford, Aberrations of Optical Systems (CRC Press, 1986).
6. J. M. Sasian, “Image plane tilt in optical systems,” Opt. Eng. 31(3), 527 (1992). [CrossRef]
7. B. A. Watson and L. F. Hodges, “Using texture maps to correct for optical distortion in head-mounted-displays,” in Proceedings of the Virtual Reality Annual Int. Symposium, 172–178 (1995).
8. J. Weng, P. Cohen, and M. Herniou, “Camera calibration with distortion models and accuracy evaluations,” IEEE Trans. Pattern Anal. Mach. Intell. 14(10), 965–980 (1992). [CrossRef]
9. A. Basu, S. Licardie, A. Basu, and S. Licardie, “Alternative models for fish-eye lenses,” Pattern Recognit. Lett. 16(4), 433–441 (1995).
10. F. Devernay, O. Faugeras, F. Devernay, and O. Faugeras, “Straight lines have to be straight,” Mach. Vis. Appl. 13(1), 14–24 (2001).
11. D. Clause and A. W. Fitzgibbon, “A rational function lens distortion model for general cameras,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 213–219 (2005).
12. J. P. Rolland, “Wide-angle, off-axis, see-through head-mounted display,” Opt. Eng. 39(7), 1760–1767 (2000). [CrossRef]
13. K. Fuerschbach, J. P. Rolland, and K. P. Thompson, “A new family of optical systems employing φ-polynomial surfaces,” Opt. Express 19(22), 21919–21928 (2011). [CrossRef] [PubMed]
14. C. Slama, Manual of Photogrammetry, (Amer. Soc. of Photogrammetry, 1980).
15. J. Heikkilä and O. Silven, “A four-step camera calibration procedure with implicit image correction,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 1106–1112 (1997).
16. J. Wang, F. Shi, J. Zhang, and Y. Liu, “A new calibration model of camera lens distortion,” Pattern Recognit. 41(2), 607–615 (2008). [CrossRef]
17. A. E. Conrady, “Decentered lens systems,” Monthly Notices of The Royal Astr Society 39, 384–390 (1919).
18. D. C. Brown, “Decentering distortion of lenses,” Photogramm. Eng. Remote Sensing 32, 444–462 (1966).
19. R. I. Hartley and T. Saxena, “The cubic rational polynomial camera model,”in Proceedings of Defense Advanced Research Projects Agency Image Understanding Workshop, 649–653 (1997).
20. G. Q. Wei and S. D. Ma, “A Complete two-plane camera calibration method and experimental comparisons,” in Proceedings of IEEE 4th International Conference on Computer Vision, 439–446 (1993).
21. D. Ruprecht, H. Muller, D. Ruprecht, and H. Muller, “Image warping with scattered data interpolation,” IEEE Comput. Graph. Appl. 15(2), 37–43 (1995).
22. C. C. Leung, P. C. K. Kwok, K. Y. Zee, and F. H. Y. Chan, “B-spline interpolation for bend intra-oral radiographs,” Comput. Biol. Med. 37(11), 1565–1571 (2007). [CrossRef] [PubMed]
23. P. Cerveri, S. Ferrari, and N. A. Borghese, “Calibration of TV cameras through RBF networks,” Proc. SPIE 3165, 312–318 (1997). [CrossRef]
24. G. E. Martin, Transformation Geometry: An Introduction to Symmetry (Springer, 1982).
25. D. N. Fogel, “Image Rectification with Radial Basis Functions: Applications to RS/GIS Data Integration,” in Proceedings of the Third International Conference on Integrating GIS and Environmental Modeling,(Sante Fe, 1996).http://www.ncgia.ucsb.edu/conf/SANTA_FE_CD-ROM/sf_papers/fogel_david/santafe.html
26. X. Zhu, R. M. Rangayyan, and A. L. Ellis, Digital Image Processing for Ophthalmology: Detection of the Optic Nerve Head (Morgan & Claypool, 2011), Chap 3.
27. J. P. McGuire Jr., “Next-generation head-mounted display,” Proc. SPIE 7618, 761804, 761804-8 (2010). [CrossRef]
28. B. Fornberg, E. Larsson, and N. Flyer, “Stable computations with Gaussian radial basis functions,” SIAM J. Sci. Comput. 33(2), 869–892 (2011). [CrossRef]
