SPATIAL UNCERTAINTY IN LINE-SURFACE INTERSECTIONS WITH APPLICATIONS TO PHOTOGRAMMETRY

The fields of photogrammetry and computer vision routinely u se line-surface intersections to determine the point where a line intersects with a surface. The object coordinates of the intersection p oi t can be found using standard geometric and numeric algor ithms, however expressing the spatial uncertainty at the intersection poi t may be challenging, especially when the surface morpholo gy is complex. This paper describes an empirical method to characterize th e unknown spatial uncertainty at the intersection point by p ropagating random errors in the stochastic model using repeated random sampling methods. These methods accommodate complex surfa ce morphology and nonlinearities in the functional model, how ever the penalty is the resulting probability density funct ion associated with the intersection point may be non-Gaussian in nature. A formal hypothesis test is presented to show that straightfo rward statistical inference tools are available whether the data is Gaussian o r not. The hypothesis test determines whether the computed i n rsection point is consistent with an externally-derived known truth point. A numerical example demonstrates the approach in a ph otogrammetric setting with a single frame image and a gridded terrain eleva tion model. The results show that uncertainties produced by the proposed empirical method are intuitive and can be assessed with conv entional methods found in textbook hypothesis testing.


INTRODUCTION
Line-surface intersections are used in many fields to determine the 3D object coordinates of a point defined by the geometric intersection of a line with a surface.The technique is often used to determine the object coordinates of a point which appears in a single image.In the one hypothetical case where the line and surface are purposefully and intentionally defined to have no spatial error, then by logical extension, the intersection point also has no spatial error and no spatial uncertainty.However, in the other cases where the line and surface contain spatial uncertainty, then the intersection point must also contain spatial uncertainty.Characterizing this uncertainty is relatively straightforward when surface morphology is benign (e.g., planar), however uncertainty characterization becomes increasingly difficult in real-world scenarios where surface morphology turns complex and contains nonlinearities and discontinuities seen in high resolution urban models, for instance.
Line-surface intersections are also known as ray tracing in computer graphics and as single-ray backprojection by the photogrammetry community (Mikhail and Ackerman, 1976).Some also refer to line-surface intersections as collision detection, especially when a line intersects with a 3D point cloud (Klein and Zachmann, 2004).
The term uncertainty is used in this paper to allow for a broad description of the stochastic behavior of random variables, in particular random variables that are not necessarily normally distributed.In photogrammetry, covariance matrices are often used to convey the precision of random variables such as observations and parameters which are typically assumed to follow uni-or multi-normal distributions.Spatial uncertainty is conveniently conveyed in terms of the normal distribution, however recent publications show departures from this norm (Beekhuizen et al., 2011, Cuartero et al., 2010, Pollard et al., 2010).
The main contributions of this paper are to describe an empirical process to estimate the spatial uncertainty associated with an intersection point and to describe a hypothesis test to determine if the intersection point is consistent with an externally-defined known point.The first section describes the classical method of computing uncertainty using standard error propagation techniques and it touches on some of limitations of the classical technique.The second section describes an alternate method of uncertainty estimation using a Repeated Random Sampling (RRS) method.The Monte Carlo method (Press et al., 1988) is a commonly used RRS method which is used in this paper.The second section also describes a formal statistical hypothesis test and provides an extensive numerical example of the RRS method.

CLASSICAL METHOD OF UNCERTAINTY ESTIMATION
Classical error propagation is often used to describe the uncertainty at the intersection point by propagating random errors associated with the observables through a functional model to the line-surface intersection point.More generally, this method propagates the random stochastic properties of independent variables to dependent variables through a linear (or linearized) functional model (Mikhail and Ackerman, 1976) It is commonly assumed in the classical method that the random errors are normally distributed (i.e., Gaussian) and that systematic biases do not exist.
In photogrammetry the collinearity equations are often used as the functional model to relate a point on an image to its corresponding point in the object.In this case the collinearity equations define the spatial location of a line or ray in space, and Zs denotes the altitude of the object surface plane.Let the object point where the line intersects the plane be denoted by xI = [XI YI ZI ] T .Three equations express the linearized functional model that relates the dependent variables with the independent variables.The first two equations are the inverse collinearity equations and the third equation is the altitude of the plane in object space.The three equations are expressed as where the inverse collinearity equations are denoted by fX and fY , x and y denote the image coordinates, IO denotes a vector of sensor interior orientation parameters, and EO denotes a vector of sensor exterior orientation parameters (Mikhail et al., 2001) Let Σ denote the covariance matrix associated with the observables and let J denote the Jacobian matrix containing partial derivatives that relates dependent variables with independent variables in equation ( 1), or J = ∂(X I ,Y I ,Z I ) ∂(x,y,IO,EO,Zs) .Then the 3 by 3 covariance matrix associated with the intersection point, ΣII contains elements to generate a standard error ellipsoid.

Limitations of Classical Method
While the classical method of computing uncertainties is widely used, it has limitations, especially in complex terrain.In the context of line-surface intersections, the limitations are a consequence of the linearization of the functional model and the dispersion of random errors that must lie within this linearized space.Indeed, two relevant assumptions are (Mikhail and Ackerman, 1976): 1.If the functional model is nonlinear, then error propagation is valid for a region around the point that encompasses the dispersion of the random variables.
2. If y = g(x) denotes a continuous and differentiable function, then the inverse function x = h(y) must also exist and be continuous and differentiable as well.
In practice these assumptions are often satisfied by assuming the object space is planar.

Hypothesis Test
The classical hypothesis test to determine whether the observed intersection point xI is consistent with a known truth point, xµ consists of null and alternative hypotheses stated as where 0 is a 3 by 1 vector of zeros.The test statistic, T , is defined as where χ 2 r denotes the chi-square distribution with r = 3 degrees of freedom and The critical value is χ 2 r,1−α with a significance level α.Assuming a one-sided test, we reject the null hypothesis if T > χ 2 r,1−α and conclude the intersection point is inconsistent with the known point.

UNCERTAINTY ESTIMATION BY REPEATED RANDOM SAMPLING
This section proposes repeated random sampling as a flexible approach to accurately characterize the spatial uncertainty at a linesurface intersection point.Used for years in scientific and mathematical problems, repeated random sampling methods are applied here to yield the nature of intersection uncertainty in cases where the object surface may not be planar and where the a priori stochastic model associated with the observables is known.RRS methods expose the nature of intersection uncertainties through a computationally intensive process that produces an empirical cloud of points in the vicinity of the nominal intersection point by repeatedly perturbing input observables by pseudo-random numbers and executing the line-surface intersection.In practice the RRS method described here requires an expression or computational method that returns the X, Y, Z coordinates of the intersection point given the imaging variables, the surface variables, and a functional model that relates the two.A very simple example of this required expression appears in equation ( 1), however sophisticated computational ray tracing methods are often used in practice.Conceptually let the generic computational method be thought of as a callable software function called get intersection coordinates.This function would be called as [XI , YI , ZI ] = get intersection coordinates(x, y, IO, EO, surface) (6) Notice that the input elements to this function on the right hand side are observables (or perhaps perturbed observables), and the intersection coordinates are returned by the function on the left hand side.
RRS methods are flexible and can accommodate a wide range of input types.For instance, RRS methods support various userdefined a priori stochastic models and are not limited to multinormal distributions.Further, RRS methods accommodate object space surfaces defined in a variety of forms (e.g.planar, gridded 2-1 / 2 D elevation model, 3D point cloud, etc.) and they accommodate many interpolation methods (e.g.bilinear, bicubic, nearest neighbor, etc.).
The actual construction of the 3D point cloud is accomplished by computing one 3D point at a time until a total of N points are defined.N is a large positive integer that denotes the number of RRS trials.At each trial the original unperturbed observables are perturbed according to the a priori stochastic model and a linesurface intersection algorithm (Eq.6) returns the 3D object space coordinates of the intersection point.Figure 1 summarizes the process of generating the 3D point cloud.
The fidelity and spatial accuracy of RRS output point cloud is influenced by several factors to include the fidelity of the a priori input stochastic error model (e.g., Gaussian, exponential, etc.), the quality of the pseudorandom number generator, and the number of RRS trials executed.In practice, the quality of pseudorandom number generators is generally high, so the predominant factor influencing the accuracy of RRS output is the number of RRS trials executed.
The RRS method does not provide direct access to probabilities associated with the intersection point via a closed-form mathematical expression.The classical method, by contrast, offers complete access to probabilities associated with the intersection point via the multinormal distribution.Because the RRS method only delivers a 3D point cloud and does not deliver direct access to these probabilities, one is left to create an empirical density function using the 3D point cloud itself.One way to build this empirical density function is to generate a voxel space surrounding the point cloud and count the number of points that fall in each voxel.This realization of the empirical density function provides a framework where hypothesis testing can occur and is discussed later.

Size and Shape of Point Cloud
Several factors influence the size and shape of the the 3D point cloud generated by the RRS method.The factors fall into two separate categories: geometric factors and stochastic factors.The geometric factors include image collection geometry and terrain morphology, whereas the stochastic factors include a priori stochastic models and the nature of computational algorithms (deterministic or non-deterministic).Finally, terrain interpolation falls in both categories.The contribution from each factor is described next.
Geometry plays a large role in the resulting 3D point cloud.Conceptually, one can imagine a cone of uncertainty emanating from the sensor and projecting onto a surface with complex morphology.As image collection geometry changes from nadir perspectives to highly oblique perspectives, the effect of object space occlusions can lead to discontinuities in the projected uncertainties (Figure 2).Depending on the circumstance -ranging from benign, smooth and continuous surfaces to urban surfaces with sharp corners -these discontinuities may lead to bi-modal or even multi-modal point clouds.
Stochastic factors also affect the resulting point cloud.For instance, if the a priori error covariance of the observations mapped to image space follows a multinormal distribution, then the "cone" of uncertainties will be conically shaped for uncorrelated image coordinates or elliptically shaped for correlated image coordinates.By contrast, if the a priori distribution is uniformly distributed (i.e., a boxcar distribution), then the cone of uncertainty becomes more rectangular in shape.Separately, the nature of the computational algorithms that return the coordinates of the line-surface intersection may also affect the shape of the point cloud, especially if the algorithms are non-deterministic.Nondeterministic algorithms deliver a range of ouputs given a single input.The output from non-deterministic algorithms may contain random errors, systematic errors, or both.
Finally, surface interpolation methods fall into both the geometric and stochastic categories.In a geometric sense, the interpolation method can introduce extreme discontinuities in the resulting 3D point cloud when nearest neighbor interpolation methods are used, for instance.By contrast, a higher order interpolation method generally produces a smooth surface that may minimize the discontinuities in the resulting point cloud.From a stochastic perspective, surface interpolation may introduce both random and non-random errors into the resulting point cloud.The magnitude of interpolation errors is a function of terrain point spacing, terrain morphology, and terrain interpolation method (Hu et al., 2009)

Special Case
The classical and RRS methods produce nearly identical estimates of uncertainty in at least one special case.This special case occurs when the a priori stochastic model is multinormal and when the object space is planar.The uncertainty estimate is realized in the form of a 3 by 3 covariance matrix that follows the multinormal distribution.The reason the methods are only nearly identical and not exactly identical is because the classical method operates on a closedform expression of the a priori stochastic model that completely describes the theoretical probability density function (PDF) associated with the intersection point, whereas the RRS method attempts to closely approximate the theoretical PDF using a finite number of empirical samples.The difference between the two methods tends to decrease as the number of RRS samples increases.

Hypothesis Test
The null hypothesis under evaluation is whether the coordinates of an intersection point, xI = [XI YI ZI ] T , are consistent with the coordinates of a truth point, xµ = [Xµ Yµ Zµ] T .As in equation ( 3), the null and alternate hypotheses are stated as The objective of this test is no different than that of the classical method -to identify a density function and test statistic that are consistent with the null hypothesis.The primary challenge here is that the density function is unknown and must be created from empirical data points.Adding to the complexity, the density function is frequently non-Gaussian, potentially multi-modal in shape, and may contain random and non-random elements.
To meet these challenges and to accommodate the possibility of known a priori uncertainties in the truth point, xµ, an empirical density function is constructed on the simple difference metric d = xI − xµ.While other metrics could have been selected and may be viable, d is selected here because it is rooted in the null hypothesis statement and because it preserves the potentially complex shape and structure of the density function that originates from the line-surface intersection.Like the classical case, when the a priori stochastic model is Gaussian and the object space is benign (i.e., planar), the expected value of d is zero.However, in scenarios where the object space is complex, the density function of d can become irregular, causing expected values to differ from zero.
The actual construction of the density function is an exercise in density estimation from empirically observed data points.While other approaches to density estimation exist (Parzen, 1962, Mohamed et al., 2005), the appoach taken here is a relatively simple voxel method.The density function is constructed by placing the N observed data points into the voxel space.Let the voxel space contain M elements in each of the X, Y and Z directions.Let the integer indices at an arbitrary voxel be (i, j, k).
The outcome of this exercise is an expression -much like the analytical function in the classical case -whose total mass is 1 or Given access to the individual voxel densities, the proposed approach is to use standard elements of classical hypothesis testing as a guide with adaptations to accommodate the empirical nature of the problem and to accommodate the potentially complex shape and structure of the density function.
The p-value is a classical element that is reused here; it serves as a computed quantity to determine if the null hypothesis is rejected or not.To obtain the p-value let the integers (ip, jp, kp) denote the indices of the voxel containing the point Then the p-value , p, is computed as where A is the set of all indices (i, j, k) where the voxel density is less than or equal to the voxel density containing point d or The p-value is small when d is inconsistent with the null hypothesis and large when d is consistent with the null hypothesis.Let α denote the user-defined significance level of the test.The decision rule is to reject the null hypothesis when p < α.
The actual construction of the density function associated with d is obtained by augmenting the steps used to obtain the RRS point cloud described in Figure 1 with additional steps to accommodate truth point uncertainty and to place the observed data points into a voxel space.When it exists, the truth point uncertainty inflates the geometric size of ρ d .Let the origin of the density function lie at a point x I nominal defined by the original, unperturbed observations.Also, in the context of the N RRS trials, let perturbed Parameter value and standard deviation focal length = 100mm ± 0.01mm Table 1: Frame image parameters used in example values be denoted with an elevated tilde (e.g., xI denotes the perturbed intersection point).Then the empirical data points that make up the density function are defined by three steps which are repeated at each trial: 1) perform the ray intersection using perturbed observations to obtain xI , 2) add a perturbation to xµ according to the truth covariance matrix Σµµ to obtain xµ, and 3) obtain a single datapoint d = xI − xµ.The N data points obtained in this manner are then placed in the voxel space and the density function about d is computed.Figure 3 summarizes these steps.

Example
The theory above is applied to a real-world single-image example in this section.The data consist of a synthetic gridded elevation model and a synthetic image (Figures 4  and 5) to test whether the coordinates of an intersection point are consistent with the coordinates of an externally-defined ground truth point.Linear units are in meters and angular units are in decimal degrees.For simplicity it is assumed that all observables are normally distributed and uncorrelated.To simplify the example, the gridded elevation data is assumed to have no absolute or relative horizontal error (i.e., σ = 0).By contrast the absolute elevation uncertainty at each node of the grid is σ = ±1 meter with no correlation between nodes.Bilinear interpolation is used to define terrain elevations between nodes.The elements of the standard frame image are described in  Figure 5 illustrates the notional geometry of the problem as well as the location of 100 intersection points relative to the gridded elevation model.These intersection points are computed using the procedure outlined in Figure 1.A bimodal distribution of the RRS point cloud is evident in Figure 5 where one mode lies near the highest elevation in the elevation model and the second mode lies at the lower elevation beside it.
Given the object-space points in the RRS point cloud, the next step is to compute the density of each voxel using the procedure outlined in Figure 3.This voxel space is illustrated in Figure 6 where an isosurface is rendered on the density function constructed from the entire set of N = 100, 000 points.As expected, Figures 5 and 6 exhibit a similar bi-modal distribution since they are produced from the same RRS simulation.
From a qualitative point of view, we reject the null hypothesis, Ho, because the point d lies outside the 1 − α = 0.95 confidence surface; we would accept Ho if d lied inside the surface (Figure 6).From a quantitative perspective, the null hypothesis is rejected because the p-value , p = 0.002, is less than the signifi-cance level α = 0.05.In summary, Therefore, according to this proposed methodology, since p < α, the ground truth point is considered inconsistent with intersection point and we reject Ho in this case.

DISCUSSION
While the proposed method appears to correctly handle challenges associated with non-Gaussian data, it suffers from long computation times required to create an RRS point cloud containing thousands or hundreds of thousands of points.In the example it took neary 24 hours of computation time on a modern desktop computer using a crude 4-threaded parallel processing scheme.
Given today's computing power it may be difficult to generate the RRS point cloud and compute the intersection uncertainty in near real-time.The computational task would be difficult even under ideal circumstances with highly-optimized multi-core and multi-threaded computational resources.The computational burden grows substantially when one imagines computing intersection uncertainties for every pixel on an image.

CONCLUSIONS AND FUTURE
A repeated random sampling method has been proposed to characterize the spatial uncertainty at a point where a line intersects a surface.The classical method of characterizing spatial uncertainty is restricted to cases where the object surface is planar.The RRS method, on the other hand, makes no assumptions concerning surface morphology.The RRS method does not provide direct access to a closed-form expression of the probability density function; consequently an empirical density function is created from the 3D data points produced by the RRS method.The empirical density function is realized through a simple technique using voxels.A formal hypothesis test is designed to determine whether an externally-defined point and its associated uncertainty is consistent with the intersection point and its uncertainty.An empirical density function is constructed from the RRS process and a test statistic is used to determine if the null hypothesis is rejected or accepted.
Future efforts will focus on techniques to improve and enhance the fidelity of the RRS empirical density function.While the voxel method described in this paper is intuitive, other sophisticated methods of density estimation may result in greater fidelity of the density function while at the same time requiring fewer RRS trials.Likewise, methods such as stratified sampling seek to exhaustively interrogate the sample space using fewer trials (Wikipedia, 2010).Also, computation time may be improved by using an intelligent convergence test to determine when the density function has reached a steady state and does not benefit from additional RRS trials.

Figure 2 :
Figure 2: A cone of uncertainty originates at the sensor and projects onto the object.

Figure 1 :
Figure 1: Procedure to obtain empirical 3D point cloud using RRS method

Figure 3 :Figure 4 :
Figure 3: Procedure to obtain the empirical density function associated with the null hypothesis 1

Figure 5 :Figure 6 :
Figure5: Synthetic data set illustrating notional geometry between image and object.The 3D point cloud in object space is produced using a synthetic image and terrain data described in the text.The symbol ⊗ denotes the nominal location of the intersection point in the image and raster elevation model, while points in the RRS point cloud are denoted by "•".For visual clarity, the first 100 of the N = 100, 000 points are plotted.The shape and bi-modal structure of the point cloud closely resembles the density function in Figure6.
Table1.A total of N = 100, 000 trials are used to describe the density function surrounding x I nominal .The coordinates of the nominal intersection point are computed using the nominal values appearing in Table1and Figure4to produce x I nominal ≈ [32.78 30.00 5.78] T .The known, externally-derived truth point is defined with coordinates xµ = [30.0029.00 4.00] T with a covariance matrix equal to I3.Consequently, the difference vector for this sample data is d = xI − xµ ≈ [2.78 1.00 1.78] T .Voxels are 0.5 meters on a side.