The computation of shortest paths is a key problem arising in a number of diverse application areas including geographic information systems, robotics, computer graphics, computer-aided design, medical computing and others. This has motivated the study and subsequent design of efficient algorithms for solving shortest path problems in different settings based on the geometric nature of the problem domain (e.g., two-dimensional (2-d), three-dimensional (3-d), surfaces, presence/absence of obstacles) and the cost function/metric (e.g., Euclidean, , link distance, weighted/unweighted, multi-criteria). In addition to its driver - the applications - the field has provided, and continues to do so, exciting challenges from a theoretical perspective. As a result, shortest path problems have become fundamental problems in areas of Computer Science such as Computational Geometry and Algorithmic Graph Theory.
The standard 3-d Euclidean shortest path problem of computing a shortest path between pair of points avoiding a set of polyhedral obstacles, denoted as the ESP3D problem, is known to be -hard even when the obstacles are parallel triangles in the space. It is not difficult to see that the number of combinatorially distinct shortest paths from a source point to a destination point may be exponential in the input size. Canny and Reif  used this to establish the -hardness of the ESP3D problem, for any metric, . In addition to this combinatorial hardness result, Bajaj  has provided an algebraic hardness argument that an exponential number of bits may be required. More recently, Mitchell and Sharir  gave -completeness proofs for the problem of computing Euclidean shortest paths among sets of stacked axis-aligned rectangles, and computing -shortest paths among disjoint balls. Given the -hardness of the ESP3D problem, work has focused on exploiting the geometric structure of the obstacles and/or on providing approximation algorithms. We will mention some of these approaches in Section 1.4.
In many applications, the Euclidean metric does not capture adequately the nature of the problem, for instance when the problem domain is not homogeneous. This motivates the weighted versions of the shortest path problem. For example in the 2-d case, consider triangulated regions where each triangle represents a particular terrain type such as water, rock, or forest. Here different weights capture the cost of traveling a Euclidean unit-length through each face. Incorporating weights makes the solution more difficult to obtain even in 2-d, but it does provide more realistic answers. It is known that light and other types of waves (e.g., seismic and sonic) travel along the shortest paths in heterogeneous media. Hence, algorithms solving the weighted shortest path problem (WSP3D) can be used for modeling wavefront propagation in such media. In the 3-d, a number of applications are non-homogeneous in nature and can be expressed using the weighted model. Next, we list some of such potential applications.
, seismic refraction and reflection methods are used based on measurements of the travel time of seismic waves refracted at the interfaces between subsurface layers of different densities. As such waves propagate along shortest paths and weighted shortest path algorithms may be used to produce more accurate and more efficient estimation of subsurface layer characteristics, e.g., the amount of oil contained in the subsurface. Another related application is the assessment of garbage dumps’ health. When a new garbage dump is built, sensors are placed at the bottom, and when the garbage dump starts to fill, waves from the top passing through the garbage to these sensors are used in order to determine the decomposition rate of the garbage .
Computation of 3-d shortest path have also been used to compute fastest routes for aircrafts between designated origin and destination points while avoiding hazardous, time-varying weather systems. Krozel et al.  investigate synthesizing weather avoidance routes in the transition airspace. Our weighted 3-d region model can be used to generalize that approach: instead of totally avoiding undesirable regions, one can assign penalty weights to them and then search for routes that minimize travel through such regions, while also avoiding unnecessarily long detours.
In medical applications simulation of sonic wavefront propagation is used when performing imaging methods as photoacoustic tomography or ultrasound imaging through heterogeneous tissue [12, 27]. In radiation therapy, domain features include densities of tissue, bone, organs, cavities, or risk to radiation exposure, and optimal radiation treatment planning takes this non-homogeneity into consideration.
The problem of time-optimum movement planning in 2-d and 3-d for a point robot that has bounded control velocity through a set of polygonal regions of given translational flow velocities has been studied by Reif and Sun . They state that this intriguing geometric problem has immediate applications to macro-scale motion planning for ships, submarines, and airplanes in the presence of significant flows of water or air. Also, it is a central motion planning problem for many of the meso-scale and micro-scale robots that have environments with significant flows that can affect their movement. They establish the computational hardness for the 3-d version of this problem by showing the problem to be PSPACE hard. They give a decision algorithm for the 2-d flow path problem, which has very high computational complexity, and they also design an efficient approximation algorithm with bounded error. The determination of the exact computational complexity of the 3-d flow path problem is posed as an open problem. Although, our weighted 3-d model does not apply directly to this setting, it can be used to construct initial approximations by assigning appropriate weights depending on the velocity and direction of the flows in different regions. In addition, the discretization scheme and the algorithmic techniques developed here can be used for solving the 3-d flow path problem.
1.2 Problem formulation
In this paper, we consider the following problem. Let be a connected 3-d domain consisting of tetrahedra with a positive real weight associated to each of them. The 3-d weighted shortest path problem (WSP3D) is to compute minimum-cost paths in from a fixed source vertex to all vertices of . The cost of a path in is defined as the weighted sum of the Euclidean lengths of the sub-paths within each crossed tetrahedron. We will describe and analyze an approximation algorithm for this problem that, for any real number , computes paths whose costs are at most times greater than the costs of the minimum cost paths. In Section 2, we describe our model in detail.
Note that the WSP3D problem can be viewed as a generalization of the ESP3D problem. Namely, given an instance of the ESP3D problem, one can find a large enough cube containing all the obstacles, tetrahedralize the free-space (i.e., exterior of the obstacles, but in the interior of the cube) and set equal weights to the resulting tetrahedra obtaining an instance of the WSP3D problem.
A key difference between Euclidean shortest path computation in 2-d and 3-d weighted domain is the -hardness already mentioned. Underlying this is the fact that, unlike in 2-d, Euclidean 3-d shortest paths are not discrete. Specifically, in 2-d, the edges of a shortest path (e.g., Euclidean shortest paths among obstacles in the plane) are edges of a graph, namely, the visibility graph of the obstacles including the source and the destination points. In contrast, in polyhedral 3-d domains, the bending points of shortest paths on obstacles may lie in the interior of the obstacles’ edges. Moreover, in weighted 3-d settings, bending points may even belong to the interior of the faces.
Furthermore, even in the case of weighted triangulated planar domains, the (weighted) shortest path may cross each of the cells times and may be composed of segments in total. Not only is the path complexity higher, but the computation of weighted shortest paths in 2-d turns out to be substantially more involved than in the Euclidean setting. In fact, there is not even an exact algorithm known, and the first approximation algorithm due to  had an time bound, where is the number of triangles in the subdivision. This problem has been actively researched since then, and currently the best known algorithm for the weighted region problem on planar subdivisions (as well as on polyhedral surfaces) runs in time . (Also, see  for a detailed literature review for the planar case.)
One of the classical tools of Computational Geometry is the Voronoi Diagram. This structure finds numerous applications (see e.g., ). It is also a key ingredient in several efficient shortest path algorithms. Researchers have studied these diagrams under several metrics (including Euclidean, Manhattan, weighted, additive, convex, abstract) and for different types of objects (including points, lines, curves, polygons), but somehow the computation of these diagrams in media with different densities (i.e., the refractive media) remained elusive. One of the main ingredients in solving the problem studied here is to compute (partial) additive Voronoi diagrams of points in refractive media. The generic techniques of Klein [15, 16] and Lê  do not apply in this case, as the bisecting surfaces do not satisfy the required conditions. In this paper, we make an important step towards the understanding and computation of these diagrams.
1.4 Previous related work
By now, shortest path problems in 2-d are fairly well understood. Efficient algorithms have been developed for many problem instances and surveys are readily available describing the state of the art in the field.
In 3-d, virtually all the work has been devoted to the ESP3D problem. Papadimitriou  suggested the first polynomial time approximation scheme for that problem. It runs in time, where is the number of bits of precision in the model of computation. Clarkson  provided an algorithm running in time, where is the ratio of the longest obstacle edge to the distance between the source and the target vertex, , and is the inverse Ackermann’s function.
Papadimitriou’s algorithm was revised and its analysis was refined by Choi et al.  under the bit complexity framework. Their algorithm runs roughly in time, where represents the time (or bit) complexity of multiplying -bit integers and . In , the same authors further developed their ideas and proposed a precision-sensitive algorithm for the ESP3D problem. In , Asano et al. proposed and studied a technique for computing approximate solutions to optimization problems and obtained another precision-sensitive approximation algorithm for the ESP3D problem with improved running time in terms of .
Har-Peled  proposed an algorithm that invokes Clarkson’s algorithm as a subroutine times to build a data structure for answering approximate shortest path queries from a fixed source in time. The data structure is constructed in roughly time. Agarwal et al.  considered the ESP3D problem for the case of convex obstacles and proposed an approximation algorithm running in time, where is the number of obstacles. In contrast to all other algorithms discussed here, the complexity of this algorithm does not depend on the geometric features of the obstacles. In the same paper, the authors describe a data structure for answering approximate shortest path queries from a fixed source in logarithmic time.
In the weighted (non-Euclidean) 3-d case no previous algorithms have been reported by other authors. In , we have announced and sketched a polynomial time approximation scheme for WSP3D problem that runs in time. The run-time improves to when all weights are equal. This algorithm can be used to efficiently solve the ESP3D problem. In this paper, we apply that approach, but develop the required details, apply new techniques, improve the complexity bounds, and provide a rigorous mathematical analysis.
1.5 Contributions of this paper
In this paper, we make several contributions to the fields of shortest path computations and the analysis of weighted 3-d regions model, as listed below.
We provide an approximation algorithm for solving the WSP3D problem in a polyhedral domain consisting of weighted tetrahedra. The algorithm computes approximate weighted shortest paths from a source vertex to all other vertices of in time, where is the user-specified approximation parameter and is a geometric parameter related to the aspect ratios of tetrahedra111See Lemma 3.1 for details on the value of .. The cost of the computed paths are within a factor of of the cost of the corresponding shortest paths.
As we stated above, the ESP3D problem, i.e., the unweighted version of this problem, is already -hard even when the obstacles are parallel triangles in the space . The time complexity of our algorithm, which is designed for the more general weighted setting, compares favorably even when applied in the Euclidean setting to the existing approximation algorithms.
Our detailed analysis, especially the results on additive Voronoi diagrams derived in Section 2, provides valuable insights into the understanding of Voronoi diagrams in heterogeneous media. This may open new avenues, for example, for designing an algorithm to compute discretized Voronoi diagrams in such settings.
Our approximation algorithms in 2-d have proven to be easily implementable and of practical value . Our algorithm for WSP3D presented here, in spite of being hard to analyze, essentially uses similar primitives, and thus has the potential to be implementable, practical, and applicable in different areas.
Our work provides further evidence that discretization is a powerful tool when solving shortest path-related problems in both Euclidean and weighted settings. We conjecture that the discretization methodology used here generalizes to any fixed dimension.
Furthermore, our discretization scheme is independent of the source vertex and can be used with no changes to approximate paths from other source vertices. This feature makes it applicable for solving all pairs shortest paths problem and for designing data structures for answering shortest path queries in weighted 3-d domains.
The complexity of our algorithm does not depend on the weights assigned to the tetrahedra composing , but it depends on their geometry. We analyze and evaluate that dependence in detail. Geometric dependencies arise also in Euclidean settings and in most of the previous papers. For example, in Clarkson , the running time of the algorithm depends on the ratio of the longest edge to the distance between the source and the target vertex. Applying known techniques (see e.g., ), such dependency can often be removed. Here, this would be possible provided that an upper bound on the number of segments on weighted shortest paths in 3-d is known. However, the increase in the dependency on in the time complexity that these techniques suffer from, which is of order , appears not to justify such an approach here. In our approach, the dependency on the geometry is proportional to the average of the reciprocal squared sinuses of the dihedral angles of the tetrahedra composing . Thus, when is large, many tetrahedra would have to be fairly degenerate so that this average to play a major role. We therefore conjecture that in typical applications, the approach presented here would work well.
1.6 Organization of the paper
In Section 2, we describe the model used throughout this paper, formulate the problem, present some properties of shortest paths in 3-d, and derive a key result on additive Voronoi diagrams in refractive media. In Section 3, we describe our discretization scheme which is a generalization of the 2-d scheme introduced in . In Section 4, we construct a weighted graph, estimate the number of its nodes and edges and prove that shortest paths in can be approximated by paths in the graph so constructed. In Section 5, we present our algorithm for the WSP3D problem. In Section 6, we conclude this paper.
2 Problem formulation and preliminaries
Let be a connected polyhedral domain in 3-d Euclidean space. We assume that is partitioned into tetrahedra , such that and the intersection of any pair of different tetrahedra is either empty or a common element (a face, an edge, or a vertex) on their boundaries. We call these tetrahedra cells. A positive weight is associated with each cell representing the cost of traveling in it. The cost of traveling along a boundary element of a cell is the minimum of the weights of the cells incident to that boundary element. We consider paths (connected rectifiable curves) that stay in . The cost of a path in is defined by , where denotes the Euclidean length of the intersection . Boundary edges and faces are assumed to be part of the cell from which they inherit their weight.
Given two distinct points and in , the shortest path problem in is to find a minimum cost path between and that stays in . We refer to the minimum cost paths as shortest paths. For a given approximation parameter , a path is an -approximation of the shortest path , if . Without loss of generality, we may assume that the points and are vertices of , since, if they are not, we can make them such by partitioning the cells where they belong. In this paper, we present an algorithm that, for a given source vertex and an approximation parameter , computes -approximate shortest paths from to all vertices of .
In this setting, it is well known 222The 2-d case was treated there, but the arguments readily apply to the 3-d model considered in this paper. that shortest paths are simple (non self-intersecting) and consist of a sequence of segments, whose endpoints are on the cell boundaries. The intersection of a shortest path with the interior of a cell, a face, or an edge is a set of disjoint segments. More precisely, each segment on a shortest path is of one of the following three types:
(i) cell-crossing – a segment that crosses a cell joining two points on its boundary;
(ii) face-using – a segment lying along a face of a cell;
(iii) edge-using – a segment along an edge of a cell.
We define linear paths to be paths consisting of cell-crossing, face-using, and edge-using segments exclusively. A linear path can be represented as the sequence of its segments or, equivalently, as the sequence of points , lying on the cell boundaries that are endpoints of these segments, i.e., , , and . The points that are not vertices of cells are called bending points of the path.
The local behavior of a shortest path around a bending point , lying in the interior of a face , is fully described by the directions of the two segments of the shortest path, and , that are incident to . The direction of each of these two segments is described by a pair of angles, which we denote by and , respectively. The in-angle is defined to be the acute angle between the direction normal to and the segment . Similarly, the out-angle is the acute angle between the normal and the segment . The angles and are the acute angles between the orthogonal projections of and with a reference direction in the plane containing the face , respectively (see Figure 1).
It is well known that when is a shortest path it is a linear path such that the angles
and are related by Snell’s law as follows:
Snell’s Law of Refraction: Let be a bending point on a geodesic path lying in the interior of a face of D. Let the segments preceding and succeeding in be and , respectively. Let the weights of the cells containing and be and , respectively. Then belongs to the plane containing and perpendicular to and the angles and are related by .
We refer to linear paths that are locally optimal (i.e., satisfy the Snell’s law) as geodesic paths. Hence, the shortest path between pair of vertices and is the geodesic path of smallest cost joining them. In the following we discuss some of the implications of Snell’s law on the local behavior of geodesic paths.
Hereafter, we denote by the ratio . Without loss of generality, we assume that , i.e., . Let be the acute or right angle for which . We refer to this angle as the critical angle for the face . From Snell’s law, it follows that . The case where deserves a special attention. In this case, must be a right angle and therefore the segment is a face-using segment. Furthermore, if the second endpoint of is in the interior of , then the segment following is inside the tetrahedron containing , and the out-angle at is equal to (see Figure 3 (b)). In summary, if is a face-using segment, then it is preceded and followed by segments lying in the cell with bigger weight and their corresponding in-angle and out-angles are equal to the critical angle .
In the next subsection, we study the properties of simple geodesic paths joining points in neighboring cells or in the same cell through a face-using segment. We define a function related to the cost of these geodesic paths and prove a number of properties that it possesses. These properties are essential to the design and the analysis of our algorithm.
2.2 Weighted distance function
Let be a plane in the three-dimensional Euclidean space. We denote the two half-spaces defined by by and and assume that positive weights and have been associated with them, respectively. We extend our model by assigning a weight to , so that if , and if . The latter case models the situation where the geodesic path joins two points in the same cell through a face-using segment on the boundary of that cell.
We refer to the half-spaces and as the lower and the upper half-space, respectively. Let be a point in the lower half-space at distance from and be a line parallel to in the upper half-space at distance from . Let be a Cartesian coordinate system such that the plane coincides with , has coordinates , and the line is described by .
We consider a point on and denote by the geodesic path between and . In this setting, the geodesic path is unique and thus coincides with the shortest path. We denote the cost of this path by , where is the -coordinate of . So, for fixed , can be viewed as a function defined for any real . We call the weighted distance function from to (Figure 2).
The local structure of the geodesic path is governed by Snell’s law. In the case where , the shortest path between and consists of two segments and , where the bending point is uniquely determined by Snell’s law (Figure 3 (a)). In the case where , the structure of the path is as follows. It is a single segment , provided that the angle between and the direction normal to the plane is smaller than or equal to the critical angle defined by . Or, if , it is in the plane perpendicular to containing and and consists of three segments , and , where the acute angles between the segments and , and the direction normal to the plane are equal to the critical angle , and the segment is in (Figure 3 (b)).
From these observations, it follows that, in all cases, weighted distance function can equivalently be defined by
In the case where , the minimum is achieved when , where is the unique solution in of the equation
where and . The latter leads to an algebraic equation of degree four and it is infeasible333Although the roots of a quartic can be expressed as a rational function of radicals over its coefficients, they are too complex to be analytically manipulated and used here. to evaluate explicitly.
The case where is easier, as in that case the geodesic path is either a straight line, or a three segment path as described above and illustrated in Figure 3 (b) and the function has an explicit representation, which is
where and . We refer to this case as the explicit case. In the next lemma we state and prove some useful properties of the function .
For a fixed , the weighted distance function has the following properties:
(a) It is continuous and differentiable.
(b) It is symmetric, i.e. .
(c) It is strictly increasing for .
(d) It is convex.
(e) It has asymptotes for as follows:
(e1) if then the asymptotes are ,
(e2) if then the asymptotes are ,
(e3) in the explicit case the asymptotes are ,
where is the critical angle.
Proof: In the explicit case (), all these properties follow straightforwardly from the explicit representation (3). So, we consider the case
From (1) and it follows that , where is the root of the equation (2). The root can be viewed as a function of , which by the implicit function theorem is continuous and differentiable. Hence property holds.
The property (b) follows from the observation that the value of the function is determined by the distance between the projections of the points and on , which is where is fixed.
To prove (c) we consider a point such that and denote by the corresponding root of (2). We have . Using the fact that the function is defined as the cost of the shortest path joining and we have
In order to prove (d), we show that for any three equidistant points on , i.e., such that , the second finite difference of the function is positive. We denote by and the bending points of the shortest paths from to and , respectively. Let be the middle point of the segment . Then, using the definition of and the convexity of the Euclidean distance function we obtain , which implies and (d).
Finally, we prove (e). Let us assume that . In this case, using Snell’s law we observe that when the bending point of the shortest path converges to the point (see Figure 2). Hence, we have . On the other hand
Combining these two limits we obtain and thus (e1) is valid for . The case where is symmetric.
In the case where we use Snell’s law and observe that the bending point of the shortest path converges to , that is . Then (e2) is established analogously to (e1).
2.3 Refractive Additive Voronoi diagram
Next we study Voronoi diagrams under the weighted distance metric defined above. Given a set of points in , called sites, and nonnegative real numbers , called additive weights, the additive Voronoi diagram for is a subdivision of space into regions , where ties are resolved in favor of the site with smaller additive weight. The regions are called Voronoi cells. In the classic case, has been defined as the Euclidean distance. Here, for , we use the weighted distance function .
Let and be two points in . We wish to study the intersection of the additive Voronoi diagram of and with with respect to the weighted distance function. Without loss of generality, we assume that and , where and are the additive weights assigned to and , respectively. We denote the intersection of the Voronoi cell of with by , or simply by when no ambiguity arises. We have
In Theorem 2.1, we will show that if and are at the same distance to , then the Voronoi cell restricted to the line has a very nice structure (i.e., it is an interval). Furthermore, in Remark 2.1, we will show that if and are not at the same distance to then Theorem 2.1 does not hold. In our algorithm, presented in Section 5, we use the information about the shape of in order to propagate approximate shortest paths in and it turns out that we need to only consider the sites that are restricted to be within a half-space of and at the same distance to . But, as we will see, this case in itself is mathematically challenging and provides valuable insights into the combinatorial structure of these diagrams.
The Voronoi cell is an interval on – possibly empty, finite or infinite.
Proof: First, consider the case when . We denote the set of points in such that by and observe that it is a half-plane perpendicular to . Therefore, the set of points on for which is either a single point, the whole line , or empty. Correspondingly, the Voronoi cell is either a half-line, empty, or the whole line and the theorem holds for .
So, we assume that . We consider the equation and claim that it cannot have more than two solutions. Before we prove that claim (Claim 2.1 below), we argue that it implies the theorem.
Assume that the equation has at most two solutions. If it does not have any or has just one solution, then the theorem follows straightforwardly. In the case where it has exactly two solutions, the cell has to be either a finite interval on , or a complement of a finite interval on . From the definition of the Voronoi cell and it follows that . We know that is either empty, a half-line, or the whole line. If it is either empty or a half-line then must be either empty or a finite interval and the theorem holds.
It remains to consider the case where is the whole line . We argue that can not be a complement to a finite interval. We have and therefore the line must be parallel to the half-plane . Furthermore, the plane containing is a perpendicular bisector of the segment and thus the point must have coordinates (Figure 2). In this setting, using Lemma 2.1(e), we observe that and have same asymptotes at infinity and thus . Therefore, the cell must be finite. The theorem follows.
Next we establish the validity of the claim used in the proof of Theorem 2.1.
The equation has at most two solutions.
We will prove the claim by showing that the function is unimodal, i.e., it has at most one local extremum. We establish this property in two steps. First, we prove a characterization property of a local extremum of (Proposition 2.1). Then, we show that there may be no more than one point possessing that property.
We focus our discussion on the case , since the other case is either simpler or can be treated analogously. We denote by and the bending points defining the shortest paths from and to , respectively. We assume that is oriented and denote by
the positive direction unit vector on. Furthermore, let and be the angles between vectors , and , respectively (see Figure 4).
These angles are completely defined by the angles and defining the corresponding shortest paths at the bending points and . Precisely, we have . Next, we prove that the angles and must be equal at any local extremum .
If is a local extremum of the function , then .
Proof: The proof is by contradiction. Let us assume that . We denote by the angle between vectors and . Similarly, denotes the angle between vectors and (Figure 4). The relation and Snell’s law readily imply that and , where . Thus, we have that . Without loss of generality, we assume that . Under these assumptions, we show the existence of two points on and on , such that:
By the assumption that is a local extremum of , it follows that, for any positive real number inside the interval , there are reals and , such that
On the other hand, if converges to zero, then converges to , and converges to . Therefore, the inequality implies that for a small enough , the inequalities
hold. From these inequalities, it follows that if we make the quadrilateral planar by rotation of the point around , then the obtained planar quadrilateral will be convex (Figure 5). Therefore we have
which proves (4) for and .
Next, we estimate the sum we show that there may be no more than one point possessing that property. We use (4) and obtain
On the other hand, by the definition of the weighted distance function (1), we have the inequalities
that contradict the previous strict inequality. Therefore, the angles and and consequently and must be equal.
Our next step is to show that there cannot be two points on satisfying Proposition 2.1. To do that, we study in more detail the relationship between the position of points and , and the angle . We observe that, for any fixed (the -coordinate of ), there is a one-to-one correspondence between the real numbers and the angles . That is, for a fixed , there is a one-to-one correspondence between the points on and the angle between the shortest path and the positive direction on . Hence, we may consider and study the well defined function . We prove the following:
The second mixed derivative of the function exists and is negative, i.e., .
Proof of Claim 2.1: We assume that Proposition 2.2 holds and will show that the function has at most one local extremum. Recall that the point has coordinates and let us denote the coordinates of the point by .
We first consider the case where . In this case, we observe that . In addition, the function is strictly monotone and, therefore, for any , the angels and are different. From Proposition 2.1 it follows that in this case the function has no local extremum.
Next, we consider the case . We assume, for the sake of contradiction, that has two local extrema, say and . By Proposition 2.1, and . We denote and