A Partial Differential Equation Obstacle Problem for the Level Set Approach to Visibility

TL;DR

This paper introduces a simplified PDE-based method for computing visibility sets from a point with obstacles represented as level sets, offering an alternative to ray tracing with demonstrated convergence and numerical examples.

Contribution

It proposes a new nonlinear obstacle PDE formulation for visibility problems, simplifying previous models and providing a convergent numerical scheme.

Findings

The PDE approach effectively computes visibility sets in 2D and 3D.

The numerical scheme converges under the Barles and Souganidis framework.

The method offers advantages over traditional ray tracing.

Abstract

In this article we consider the problem of finding the visibility set from a given point when the obstacles are represented as the level set of a given function. Although the visibility set can be computed efficiently by ray tracing, there are advantages to using a level set representation for the obstacles, and to characterizing the solution using a Partial Differential Equation (PDE). A nonlocal PDE formulation was proposed in Tsai et. al. (Journal of Computational Physics 199(1):260-290, 2004): in this article we propose a simpler PDE formulation, involving a nonlinear obstacle problem. We present a simple numerical scheme and show its convergence using the framework of Barles and Souganidis. Numerical examples in both two and three dimensions are presented.

Tables

Table 1. Table 1. Errors and order of convergence for Example 6.1 .

N	$h$	Error	Order
32	$1.29\text{×}{10}^{-01}$	$9.12\text{×}{10}^{-02}$	-
64	$6.35\text{×}{10}^{-02}$	$4.49\text{×}{10}^{-02}$	1.02
128	$3.15\text{×}{10}^{-02}$	$2.23\text{×}{10}^{-02}$	1.01
256	$1.57\text{×}{10}^{-02}$	$1.11\text{×}{10}^{-02}$	1.01
512	$7.83\text{×}{10}^{-03}$	$5.54\text{×}{10}^{-03}$	1.00
1024	$3.91\text{×}{10}^{-03}$	$2.76\text{×}{10}^{-03}$	1.00
2048	$1.95\text{×}{10}^{-03}$	$1.38\text{×}{10}^{-03}$	1.00
4096	$9.77\text{×}{10}^{-04}$	$6.91\text{×}{10}^{-04}$	1.00

Equations

$$\nabla \psi \cdot \frac{x-x^{\ast }}{|x-x^{\ast }|}=\min \left\{H(\psi -g)\nabla g\cdot \frac{x-x^{\ast }}{|x-x^{\ast }|},0\right\}$$

$$\min \{u(x)-g(x),(x-x^{\ast })\cdot \nabla u(x)\}=0$$

$$x\in S⟹t{x}^{\ast }+(1-t)x\in S,\text{ for all }t\in [0,1].$$

$$S_{\alpha }(u)$$

$$u(t{x}^{\ast }+(1-t)y)\le u(y),\text{ for all }y\in \Omega ,0\le t\le 1.$$

$$u(y)\le \alpha ⟹u(t{x}^{\ast }+(1-t)y)\le \alpha ,\text{ for all }y\in \Omega ,0\le t\le 1$$

$$S{S}^{-}(g)(x)=\sup \{v(x)\mid \text{ v is star-shaped with respect to x∗ and v≤g}\},$$

$$S{S}^{+}(g)(x)=\inf \{v(x)\mid \text{ v is star-shaped with respect to x∗ and v≥g}\}.$$

$$w(x)=\min \left\{g(y)\mid y=x+t(x-x^{\ast })\in \Omega ,t\ge 0\right\}$$

$$y\in \mathrm{argmin}\left\{g(y)\mid y=x+t(x-x^{\ast })\in \Omega ,t>0\right\}.$$

$$u(x)=\max \left\{g(y)\mid y=x^{\ast }+t(x-x^{\ast })\in \Omega ,t\in [0,1]\right\}.$$

$$y\in \mathrm{argmax}\left\{g(y)\mid y=x^{\ast }+t(x-x^{\ast })\in \Omega ,t\in [0,1)\right\}.$$

$$F(x,r,p,M)\le F(x,s,p,N)\text{whenever }r\le s\text{ and }N\le M.$$

$$u^{\ast }(x)$$

$$u_{\ast }(x)$$

$$F(x,u(x),\nabla \varphi (x))\le 0(\ge 0).$$

$$w(x)=\{\begin{array}{ll}-x-1& \text{if }x<-1,\\ 0& \text{if }-1\le x\le 0,\\ x+1& \text{if }x>0,\end{array}\text{and}u(x)=\{\begin{array}{ll}-x-1& \text{if }x<-2,\\ 1& \text{if }-2\le x\le 0,\\ x+1& \text{if }x>0.\end{array}$$

$$u(t{x}^{\ast }+(1-t)x)-u(x)=t\nabla u(x+s(x^{\ast }-x))\cdot (x^{\ast }-x).$$

$$\nabla u(x+s(x^{\ast }-x))\cdot (x^{\ast }-x)\le 0.$$

$$u(t{x}^{\ast }+(1-t)y)>u(y).$$

$$\nabla u(y+s(x^{\ast }-y))\cdot (x^{\ast }-y)>0$$

$$\nabla u(y+s(x^{\ast }-y))\cdot (y+s(x^{\ast }-y)-x^{\ast })\ge 0$$

$$⟺$$

$$u(x)\le u(x+h(x-x^{\ast }))$$

$$\varphi (x)\le \varphi (x+h(x-x^{\ast }))$$

$$\frac{\varphi (x+h(x-x^{\ast }))-\varphi (x)}{h}=\nabla \varphi (x+s(x^{\ast }-x))\cdot (x-x^{\ast })$$

$$\nabla \varphi (x)\cdot (x-x^{\ast })=(u(z)-u(y))(z-y)\cdot (x-x^{\ast })<0$$

$$\max \{w(x)-g(x),(x^{\ast }-x)\cdot \nabla w(x)\}=0$$

$$w(x)=\sup \{v(x)\mid \text{ v∈USC(Ω) is star-shaped with respect to x∗ and v≤g}\}.$$

$$\min \{u(x)-g(x),(x-x^{\ast })\cdot \nabla u(x)\}=0,$$

Full text

1. Introduction

In this article we consider the problem of finding the visibility set from a given viewpoint given a set of known obstacles using a Partial Differential Equation (PDE). In principle, the visibility set is simply given by ray tracing and there are numerous algorithms for solving the visibility problem using explicit representations of the obstacles [CT97, DDTP00, AS96a, AS96b].

Finding the visibility set plays a crucial role in numerous applications including rendering, visualization [HZ00], etching [SA97], surveillance, exploration [VTS14], navigation [LTC06], and inverse problems, to only name a few. Specifically, in [CT05] the level set framework [OS88, Set99] developed in [TCO*+*04] was extended to deal with the optimal placing of a single viewer or a group of viewers and $A$-to-$B$ optimal path planning, where optimality is measured in terms of the volume of the visible region. More recently, in [LT18] a convolutional neural network is proposed to determine the vantage points that maximize visibility in the context of surveillance and exploration, with the visibility sets of the training data being computed efficiently using the PDE formulation introduced in [TCO*+*04]. For applications which involve optimization of the viewpoint, the discontinuity of the visibility can make optimization more difficult. The advantage of using level set/PDE methods is the improved regularity of the solution.

It is clear then that the ultimate goal of the work is inverse problems involving visibility. As is the case with inverse problems, a better understanding of the forward problem is essential for better results of the more challenging inverse problem. In this work, we focus our attention in the forward problem and do not go further and study the inverse problem. We propose a simple formulation of the visibility problem - the visibility set is the subzero level set of the solution of a nonlinear obstacle problem.

In [TCO*+*04] the visibility problem is presented as a boundary value problem for a first order differential equation: the visibility set to a given viewpoint $x^{*}$ is given by $\{\psi(x)\geq 0\}$ where the function $\psi(x)$ is the solution of

[TABLE]

with $\psi(x^{*})=g(x^{*})$. Here $H(z)=\chi_{[0,\infty)}(z)$ is the characteristic function of $[0,\infty)$ and $g$ is a signed distance function to the obstacles, positive outside the obstacles and negative inside. Despite the complex nature of the operator in (1), in [KT08] the visibility function $\psi$ is shown to be the viscosity solution of an equivalent Hamilton-Jacobi type equation involving jump discontinuities in the Hamiltonian. A numerical scheme to solve this equation is presented and its convergence is established.

For our formulation, $g$ is still a signed distance function to the obstacles, but is instead negative outside and positive inside. Then the visibility set is given as $\{u(x)\leq 0\}$ where the function $u(x)$ solves the following nonlinear first order local PDE

[TABLE]

with $u(x^{*})=g(x^{*})$. This is not only considerably simpler than (1), but can also be generalized to allow multiple viewpoints as we will show. Moreover, each sublevel set of $u$ is in fact the visibility set of the corresponding superlevel set of $g$. Efficiencies then arise when the obstacles are given by the graph of a function (for example, heights of buildings). In this case, we can reduce the dimension of the problem, and compute the horizontal visibility set from a given height, using the level set representation. Similarly, if the $t$ superlevel set of $g$ represents the position of the obstacles at a certain time $t$ ($g$ can for instance be the solution of an Eikonal equation), then the PDE needs to be solved only once and the visibility set at any given time can be extracted from the corresponding sublevel set.

The paper is organized as follows. In section 2, we characterize visibility sets as star-shaped envelopes. In section 3 we derive the new visibility PDE and its generalization to multiple viewpoints. In section 4 we present the numerical scheme, while in section 5 we establish its convergence. Finally, in section 6 we present both two-dimensional and three-dimensional examples of visibility sets computed using the new proposed PDE.

2. Star shaped sets and functions

In this section we give an interpretation of the visibility set from a given point $x^{*}$ as the star-shaped envelope with respect to the point $x^{*}$. The definitions of star-shaped sets and envelopes are then extended to functions. Finally, we provide explicit formulas for the star-shaped envelopes of a function.

We star by recalling the definition of a star-shaped set.

Definition 2.1**.**

We say the set $S\subset\mathbb{R}^{n}$ is star-shaped with respect to $x^{*}$ if

[TABLE]

A simple example of a star-shaped set is a convex set. Indeed, convex sets are star-shaped with respect to every point inside. Moreover, intersections of convex sets are convex, but unions are not. As a consequence there is a natural (outer) convex envelope, but not an inner one. On the other hand, star-shaped sets are closed under both intersections and unions, which means one can define two star-shaped envelopes (with respect to $x^{*}$) for sets, the inner and outer envelopes.

Definition 2.2**.**

Given $S\subset\mathbb{R}^{n}$ and $x^{*}\in S$, the outer star-shaped envelope of $S$ with respect to $x^{*}$ is the intersection of all star-shaped sets with respect to $x^{*}$ which contain $S$. The inner star-shaped envelope of $S$ is the union of all star-shaped sets contained in it.

Looking at Figure 1, one immediately sees how the inner star-shaped envelope corresponds to the visibility subset of $S$ if there is an illumination source at the point $x^{*}$. The remainder of $S$ is the invisible part. One the other hand, the outer star-shaped envelope minus $S$ corresponds to the least amount of obstacles which would need to be removed so that all of $S$ is visible from $x^{*}$.

We now discuss star-shaped functions which are the main building block to characterize the visibility set as the solution of a nonlinear obstacle PDE. We write $S_{\alpha}(u)\equiv\{x\in\mathbb{R}^{n}\mid u(x)\leq\alpha\}$ for the $\alpha$-sublevel set of a function $u$ and let $\Omega$ be a star-shaped domain with respect to $x^{*}$.

Definition 2.3**.**

We say that the function $u:\Omega\to\mathbb{R}$ is star-shaped with respect to $x^{*}$ if

[TABLE]

for all $\alpha\in\mathbb{R}$.

Remark 2.1*.*

Notice the similarity to quasiconvex functions: $u$ is said to be quasiconvex if $S_{\alpha}(u)$ is convex for all $\alpha\in\mathbb{R}$.

We now characterize star-shaped functions with a zero-order condition.

Lemma 2.4**.**

A function $u:\Omega\to\mathbb{R}$ is star-shaped with respect to $x^{*}$ if and only if

[TABLE]

Proof.

By definition, $u$ is star-shaped with respect to $x^{*}$ if and only if for all $\alpha\in\mathbb{R}$, $S_{\alpha}(u)$ is star-shaped with respect to $x^{*}$. This is equivalent to the condition

[TABLE]

for all all $\alpha\in\mathbb{R}$, which in turn is equivalent to (2). ∎

We use this result to describe the monotonicity of a star-shaped function.

Proposition 2.5**.**

Let $u:\Omega\to\mathbb{R}$ be a function. Then $u$ is star-shaped with respect to $x^{*}$ if and only if $u$ is increasing along rays from $x^{*}$ to $x$. Moreover, if $u$ is star-shaped with respect to $x^{*}$, then $x^{*}$ is a global minimum of $u$.

Proof.

This follows immediately from Lemma 2.4. ∎

Next, in a similar way to star-shaped envelopes of a set, we define upper and lower star-shaped envelopes of a function $g$ with respect to a point $x^{*}$.

Definition 2.6**.**

Let $g\in C(\Omega)$ be bounded by below. The lower star-shaped envelope of $g$ with respect to $x^{*}$ is defined as

[TABLE]

while the upper star-shaped envelope is given by

[TABLE]

Remark 2.2*.*

We require that $g$ is bounded by below in order for the lower star-shaped envelope to be well defined since otherwise there would no star-shaped function with respect to $x^{*}$ bounded from above by $g$.

We finish this section by proving the following simple explicit solution formulas for the star-shaped envelope of a function.

Proposition 2.7**.**

Let $g\in C(\Omega)$ be bounded by below and define $w:\Omega\to\mathbb{R}$ to be given by

[TABLE]

with $w(x^{*})=\min_{x\in\Omega}g(x)$. Then $w$ is star-shaped with respect to $x^{*}$ and $w=SS^{-}(g)$.

Proof.

By the assumptions on $g$, $w$ is well-defined. By construction, $w$ is increasing along rays from $x^{*}$ to $x$, and so, by Proposition 2.5, $w$ is star-shaped with respect to $x^{*}$. Moreover, it is clear that $w\leq g$.

We want to show that $w=SS^{-}(g)$. Suppose by contradiction that it is not. This means that there exists a star-shaped with respect to $x^{*}$ function $v$ with $v\leq g$ such that $v(x)>w(x)$ for some $x\in\Omega$. Without loss of generality, assume that $x\neq x^{*}$. We have $w(x)<g(x)$ and that there exists $y\in\Omega$ such that

[TABLE]

Hence $v(x)>w(x)=w(y)=g(y)\geq v(y)$ and therefore $v$ is not increasing along the ray from $x^{*}$ to $x$. Finally, we invoke Proposition 2.5 to conclude that $v$ is not star-shaped with respect to $x^{*}$, which leads to the desired contradiction. ∎

Remark 2.3*.*

Intuitively, we can find $w$ by tracing the values from the boundary along rays to $x^{*}$ and taking the minimum of $g$ along the way.

Proposition 2.8**.**

Let $g\in C(\Omega)$ and define $u:\Omega\to\mathbb{R}$ to be given by

[TABLE]

Then $u$ is star-shaped with respect to $x^{*}$ and $u=SS^{+}(g)$.

Proof.

From the definition of $u$, it is clear that $u$ is increasing along rays from $x^{*}$ to $x$, and so, by Proposition 2.5, $u$ is star-shaped with respect to $x^{*}$. Moreover, $u\geq g$, again by definition of $u$.

We want to show that $u=SS^{+}(g)$. Suppose by contradiction that it is not. This means that there exists a star-shaped with respect to $x^{*}$ function $v$ with $v\geq g$ such that $v(x)<u(x)$ for some $x\in\Omega$. Without loss of generality assume $x\neq x^{*}$. We have $u(x)>g(x)$ and that there exists $y\in\Omega$ such that

[TABLE]

Hence $v(x)<u(x)=u(y)=g(y)\leq v(y)$, which means, just like in the proof of Proposition 2.7, that $v$ is not increasing along the ray from $x^{*}$ to $x$. Hence $v$ is not star-shaped with respect to $x^{*}$ according to Proposition 2.5 and we have obtained our contradiction. ∎

Remark 2.4*.*

This formula corresponds to the classic ray tracing algorithm to find the visibility set. We trace the values towards the boundary along rays from $x^{*}$ taking the maximum of $g$ along the way.

3. PDEs and Visibility

In this section, we present the new PDE formulation of visibility sets from a single viewpoint and its extension to multiple viewpoints. We start with a level set PDE interpretation for star-shaped envelopes. Given that the inner star-shaped envelope of a set corresponds to its visibility set, the PDE obtained computes the visibility set for each sublevel set. We then generalize it to multiple viewpoints.

3.1. Viscosity Solutions

Viscosity solutions [CIL92] provide the correct notion of weak solution to a class of degenerate elliptic PDEs which includes the PDEs considered here. We review it briefly here.

Let $S^{n}$ be the set of real symmetric $n\times n$ matrices, and take $N\leq M$ to denote the usual partial ordering on $S^{n}$, namely that $N-M$ is negative semi-definite.

Definition 3.1**.**

The operator $F(x,r,p,M):\Omega\times\mathbb{R}\times\mathbb{R}\times S^{n}\to\mathbb{R}$ is degenerate elliptic if

[TABLE]

Remark 3.1*.*

For brevity we use the notation $F[u](x)\equiv F(x,u(x),\operatorname{\nabla}u(x),D^{2}u(x))$.

Definition 3.2** (Upper and lower semi-continuous envelopes).**

The upper and lower semicontinuous envelopes of a function $u(x)$ are defined, respectively, by

[TABLE]

Definition 3.3** (Viscosity solutions).**

Let $F:\Omega\times\mathbb{R}\times\mathbb{R}^{n}$. We say the upper semi-continuous (lower semi-continuous) function $u:\Omega\to\mathbb{R}$ is a viscosity subsolution (supersolution) of $F[u]=0$ in $\Omega$ if for every $\phi\in C^{1}(\Omega)$, whenever $u-\phi$ has a local maximum (minimum) at $x\in\Omega$,

[TABLE]

Moreover, we say u is a viscosity solution of $F[u]=0$ if $u$ is both a viscosity sub- and supersolution.

Remark 3.2*.*

For brevity we use the notation $F[u](x)\equiv F(x,u(x),\operatorname{\nabla}u(x))$. In addition, when checking the definition of a viscosity solution we can limit ourselves to considering unique, strict, global maxima (minima) of $u-\phi$ with a value of zero at the extremum. See, for example, [Koi04, Prop 2.2].

3.2. Regularity

We briefly discuss the regularity of the star-shaped envelopes. We start by observing that the star-shaped functions need not be continuous.

Example 3.1**.**

In one dimension, the function $u(x)=1$ for $x\neq 0$ and $u(0)=0$ is star-shaped with respect to [math]. In two dimensions, take $A=\{(x,y)\in\mathbb{R}^{2}:xy=0\}$ and define $u$ as the characteristic function of the complement of $A$, i.e., $u(x)=0$ if $x\in A$ and $u(x)=1$ otherwise. Once again $u$ is star-shaped with respect to the origin, but it is not continuous. In fact, $u$ is lower semicontinuous.

As for the star-shaped envelopes of functions, the upper star-shaped envelope is continuous while the lower star-shaped envelope is only lower semicontinuous as it may be discontinuous at $x^{*}$.

Proposition 3.4**.**

Let $g\in C(\Omega)$ be bounded by below. Then $w=SS^{-}(g)$ is lower semicontinuous in $\Omega$ and continuous in $\Omega$ except at $x^{*}$, while $u=SS^{+}(g)$ is continuous in $\Omega$. If $x^{*}$ is a global minimum of $g$ then $w=SS^{-}$ is also continuous in $\Omega$.

Proof.

The proof follows from the solutions formulas (5) and (6) since we take the minimum and maximum of a continuous function $g$ along rays to and from $x^{*}$, respectively, as pointed out in Remarks 2.3 and 2.4. ∎

Example 3.2**.**

Let $g(x)=|x+1|$ and let $x^{*}=0$. Then the lower and upper star-shaped envelopes are given by

[TABLE]

A two-dimensional example is given in Figure 2: $g$ is given by the distance to two points and $x^{*}$ is chosen as a point on the $.3$ level set of $g$. The upper and lower envelopes are pictured. The lower one is discontinuous at $x^{*}$. Replacing $g$ with $\max(g,g(x^{*})$ leads to a function whose global minimum is attained at $x^{*}$ and therefore both star-shaped envelopes are continuous in this case. This is depicted in Figure 3.

3.3. A local PDE for star-shaped envelopes

We start by establishing a first order condition for star-shaped functions.

Proposition 3.5**.**

Suppose $u:\mathbb{R}^{n}\to\mathbb{R}$ is differentiable. Then $u$ is star-shaped with respect to $x^{*}$ if and only if $\operatorname{\nabla}u(x)\cdot(x-x^{*})\geq 0$ for all $x\in\mathbb{R}$.

Proof.

By the mean value theorem, given any $t\in[0,1]$ there exists $s\in(0,t)$ such that

[TABLE]

Thus if $u$ is star-shaped with respect to to $x^{*}$ we obtain using Lemma 2.4

[TABLE]

Taking the limit as $t\to 0$ leads to the desired inequality $\operatorname{\nabla}u(x)\cdot(x-x^{*})\geq 0$.

Now, suppose that $\operatorname{\nabla}u(x)\cdot(x-x^{*})\geq 0$ for all $x\in\mathbb{R}^{n}$. We argue by contradiction. Assume that $u$ is not star-shaped with respect to $x^{*}$. Hence, by Lemma 2.4, there are $y\in\Omega$ and $t\in[0,1]$ such that

[TABLE]

Then, again by the mean value theorem,

[TABLE]

for some $s\in(0,t)$. However, taking $x=y+s(x^{*}-y)$ in $\operatorname{\nabla}u(x)\cdot(x-x^{*})\geq 0$ leads to

[TABLE]

We have the desired contradiction and so the proof is complete. ∎

We are interested in star-shaped functions that may not be differentiable since typically the visibility set has corners. Thus we need to characterize star-shaped functions in the weak sense. We will do so using viscosity solutions. Throughout the rest of this section let $\Omega$ be a bounded star-shaped domain with respect to $x^{*}\in\mathbb{R}^{n}$ and denote its boundary by $\partial\Omega$.

Proposition 3.6**.**

Suppose $u\in USC(\Omega)$. Then $u$ is star-shaped with respect to $x^{*}$ if and only if $u$ is a viscosity subsolution of $\operatorname{\nabla}u(x)\cdot(x^{*}-x)=0$.

Proof.

Suppose $u$ is star-shaped with respect to $x^{*}$. Let $\phi\in C^{1}(\Omega)$ be such that $u-\phi$ has a local maximum at $x\in\Omega$. Without loss of generality assume that $\phi(x)=u(x)$. Then we have $u(y)\leq\phi(y)$ in a neighborhood of $x$. Since $u$ is star-shaped with respect to $x^{*}$, $u$ is increasing along arrays from $x^{*}$ to $x$ by Proposition 2.5 and therefore

[TABLE]

for $h>0$. Hence, given the choice of $\phi$,

[TABLE]

for $h>0$ sufficiently small. Since

[TABLE]

for $s\in(0,h)$, we obtain $\operatorname{\nabla}\phi(x)\cdot(x^{*}-x)\leq 0$ as $h\to 0$. This shows that $u$ is a viscosity subsolution of $\operatorname{\nabla}u(x)\cdot(x^{*}-x)=0$.

Suppose now that $u$ is a viscosity subsolution of $\operatorname{\nabla}u(x)\cdot(x^{*}-x)=0$ and that $u$ is not star-shaped with respect to $x^{*}$. Then there exists $y,z$ such that $u(y)\geq u(x^{*})$ and $u(y)>u(z)$ with $y$ lying on the line segment from $x^{*}$ to $z$. Without loss of generality assume that $y=\operatorname{{argmax}}_{x}u(x)$. We can then construct a linear $\phi\in C^{1}(\Omega)$ such that $\operatorname{\nabla}\phi=(u(z)-u(y))(z-y)$ and $u-\phi$ has a local maximum at $x$ with $(x-x^{*})\cdot(z-y)>0$. But then

[TABLE]

which contradicts our assumption. ∎

We can now finally write the PDEs for the upper and lower star-shaped envelopes of $g$ with respect to $x^{*}$.

Proposition 3.7**.**

Let $g\in C(\Omega)$ be bounded by below and let $w=SS^{-}(g)$ be the lower star-shaped envelope $w$ of $g$ with respect to $x^{*}$. Assume $w$ is continuous. Then $w$ is the viscosity solution of the obstacle problem

[TABLE]

along with boundary conditions $w=g$ on $\partial\Omega$.

Proof.

By Proposition 2.7, $w$ is star-shaped with respect to $x^{*}$ and therefore we can write

[TABLE]

Now, according to Proposition 3.6, $w$ is precisely the supremum of all subsolutions of the PDE and so the proof follows directly from Perron’s method. ∎

Proposition 3.8**.**

Let $g\in C(\Omega)$ and let $u=SS^{+}(g)$ be the upper star-shaped envelope (visibility) of $g$ with respect to $x^{*}$. Then $u$ is the viscosity solution of the obstacle problem

[TABLE]

along with $u(x^{*})=g(x^{*})$.

Proof.

By definition the upper star-shaped envelope is given by

[TABLE]

We observe that this is equivalent to

[TABLE]

and so $-u$ is the solution of

[TABLE]

by a similar reasoning to the one in Proposition 3.7. Now, since the equation can be rewritten as

[TABLE]

we conclude that $u$ is the solution of

[TABLE]

as desired. ∎

Remark 3.3*.*

Contrary to Proposition 3.7 we do not need to assume the continuity of the start-shaped envelope as this follows directly from the assumptions on $g$ (see Proposition 3.4).

3.4. Comparison Principle

An important property of elliptic equations, from which uniqueness follows, is the comparison principle that states that subsolutions lie below supersolutions. In addition, it also plays a crucial role when establishing the convergence of approximation schemes using the the theory of Barles and Souganidis [BS91], which we intend to do later on. However, in such setting, the comparison principle required is a strong comparison principle: The boundary conditions are satisfied in the viscosity sense. In general, such a comparison is only available when the solutions are continuous up to the boundary (see [CIL92] for more details).

We focus our attention in PDE (7) as it is the one we are most interested in: Its solution allows us to determine the visibility set. A strong comparison principle is however not satisfied: Notice that any function that is nonincreasing along any direction away from $x^{*}$ is a subsolution and thus given any supersolution we can always construct a subsolution that lies above it by adding a large enough constant. We can however circumvent this by requiring that the subsolution $u$ satisfies $u(x^{*})\leq g(x^{*})$. This will also prove to be enough to establish convergence of our numerical scheme. Intuitively, imposing that the subsolution $u$ satisfies $u(x^{*})\leq g(x^{*})$ guarantees that $u$ lies below $g$. This, together with the fact that the solution of (7) is the infimum of all supersolutions according to Perron’s method, is enough to reach the desired conclusion.

Proposition 3.9**.**

Let $g\in C(\Omega)$. Let $u$ be a viscosity subsolution of (7) such that $u(x^{*})\leq g(x^{*})$ and let $v$ be a supersolution of (7). Then $u\leq v$ in $\Omega$.

Proof.

Suppose there exists $x\in\Omega$ such that $u(x)>v(x)$ by contradiction. Since $v$ is a supersolution (7),

[TABLE]

where we used Perron’s characterization of $SS^{+}(g)$ and Proposition 2.8. In particular, we have $v(x)\geq g(x)$ and $v(x)\geq g(x^{*})$. Hence $u(x)>g(x)$ and $u(x)>~{}u(x^{*})$ by assumption on $u$. Therefore there exists a linear function $\phi$ such that $u-\phi$ has a local maximum at $y$ with $u(y)>g(y)$ and $(y-x^{*})\cdot\operatorname{\nabla}\phi(y)=u(x)-u(x^{*})>0$. We have derived a contradiction with the assumption that $u$ is a subsolution of (7) and the proof is complete. ∎

3.5. Visibility from multiple viewpoints

We are now interested in the visibility set from multiple viewpoints where a point is consider visible if is seen by at least one viewpoint. We follow the same ideas as before, starting by generalizing the definition of star-shaped set.

Definition 3.10**.**

We say that a set $S\subset\mathbb{R}^{n}$ is star-shaped with respect to $\{x^{*}_{1},\ldots,x^{*}_{r}\}$ if

[TABLE]

As in section 2 we can define star-shaped functions with respect to $\{x^{*}_{1},\ldots,x^{*}_{n}\}$ according to Definition 3.10. More importantly, Proposition 3.5 can be generalized.

Proposition 3.11**.**

Suppose $u:\mathbb{R}^{n}\to\mathbb{R}$ is differentiable and bounded by below. Then $u$ is star-shaped with respect to $\{x^{*}_{1},\ldots,x^{*}_{r}\}$ according to Definition 3.10 if and only if

[TABLE]

Remark 3.4*.*

The presence of the minimum in $t$ may not appear obvious at first, but it is crucial here. In order for a point $x$ to be visible from $x^{*}_{i}$ then $u$ must be increasing along the ray from $x^{*}_{i}$ to $x$ which guarantees that all the points along the ray will be in the visible set. Without the minimum in $t$ a point $x$ could be consider visible by first moving along a ray towards $x^{*}_{i}$ and then follow a different viewpoint.

In this case, $u=SS^{+}_{\{x^{*}_{1},\ldots,x^{*}_{r}\}}(g)$ is the solution of the following PDE

[TABLE]

Despite being a non-local PDE, unlike the previous ones with a single viewpoint $x^{*}$, we still have a fast solver available: The solution $u$ is given by

[TABLE]

where $u_{i}$ is the solution of

[TABLE]

This follows from noticing that

[TABLE]

In addition, the following solution formula is also available

[TABLE]

At this point, a natural question to ask is what other visibility definitions will lead to PDEs following the approach taken here. For instance, given $3$ viewpoints $\{x^{*}_{1},x^{*}_{2},x^{*}_{3}\}$ a point can be considered visible if (i) it is seen by at least two viewpoints; (ii) it is seen by $x^{*}_{1}$ or both $x^{*}_{2}$ and $x^{*}_{3}$ (iii) it is seen by all viewpoints. All these can be computed efficiently computed by first determining the visibility of each viewpoint and then taking the appropriate combination of maximums and minimums. However, determining a corresponding PDE following the approach considered here will in general fail. Consider for instance case (iii) where a point is visible if it is seen by all viewpoints. In this setting, one can no longer define star-shaped like sets as the visibility set may be disconnected (see Figure 4 for a simple example with two viewpoints). The fundamental difference is that it is no longer true that if a point $x$ is visible then all points along the ray from $x$ to the viewpoint are visible.

4. Convergent finite difference schemes

In this section, we discuss the numerical schemes used to solve the PDEs introduced in the previous section, focusing our attention on the visibility PDE (7). As we will see, the scheme proposed here is degenerate elliptic finite difference schemes for which there exists a well established convergence framework.

Before we begin we introduce some notation. For simplicity, we will assume we are working on the hypercube $[-1,1]^{n}\subseteq\mathbb{R}^{n}$. We write $x=(x_{1},\ldots,x_{n})\in[-1,1]^{n}$. The domain is discretized with a uniform grid, resulting in the following spatial resolution:

[TABLE]

where $N$ is the number of grid points used to discretize $[-1,1]$. We denote by $\Omega^{h}$ the computational domain which in our case reduces to $[-1,1]^{n}\cap h\mathbb{Z}^{n}$.

Our schemes are written as the operators $F^{h}[u]:C(\Omega^{h})\to C(\Omega^{h})$, where $C(\Omega^{h})$ is the set of grid functions $u:\Omega^{h}\to\mathbb{R}$. We assume they have the following form:

[TABLE]

where $u(\cdot)$ corresponds to the value of u at points in $\Omega^{h}$.

4.1. Finite difference for the visibility PDE

For simplicity, we present the scheme in the two-dimensional setting. The generalization to higher dimensions is straightforward.

Let $x\in\Omega^{h}$. Define the vector $v^{*}=x-x^{*}$ for $x\in\Omega^{h}$. If $x$ is one of the four grid points enclosing $x^{*}$, let $\tilde{x}=x^{*}$. Otherwise, let $\tilde{x}$ denote the intersection between the line that passes through $x$ and $x^{*}$ and a line segment formed by 2 of the 8 neighbors of $x$. Then $u_{v^{*}}(x)=(x-x^{*})\cdot\operatorname{\nabla}u(x)$ and its upwind approximation is given by

[TABLE]

where $\mathcal{I}_{h}$ is the piecewise linear Lagrange interpolant. The numerical scheme for the visibility PDE (7) is then given by

[TABLE]

4.2. Fast sweeping solver

We implement a fast sweeping solver to compute the solutions of (9). Solving the equation $F^{h}[u](x)=0$ for the reference variable, $u(x)$, leads to the update formula

[TABLE]

where $\tilde{x}$ was defined in the previous section. Since all characteristic are straight lines and flow away from $x^{*}$, the domain $\Omega^{h}$ only needs to be swept once (but in a very specific ordering). For simplicity, we will assume that $x^{*}\in\Omega^{h}$. Let $u_{i,j}$ denote the solution at $x_{i,j}$ where

[TABLE]

Set $(i^{*},j^{*})$ such that $x_{i^{*},j^{*}}\leq x^{*}<x_{i^{*}+1,j^{*}+1}$ (here the inequalities are interpreted component-wise). Set as well $u_{i^{*},j^{*}}=g(x^{*})$. The domain is then divided into four quadrants and each is swept in the following way:

• •

$i=i^{*},\ldots,N$, $j=j^{*},\ldots,1$ (sweeping bottom right square)

• •

$i=i^{*},\ldots,N$, $j=j^{*},\ldots,N$ (sweeping top right square)

• •

$i=i^{*}+1,\ldots,1$, $j=j^{*},\ldots,1$ (sweeping bottom left square)

• •

$i=i^{*}+1,\ldots,1$, $j=j^{*},\ldots,N$ (sweeping top left square)

Remark 4.1*.*

Alternatively, the solver can be seen as a direct consequence of the solution formula (6).

5. Convergence of Numerical Solutions

In this section we recall the notion of degenerate elliptic schemes and show that the solutions of the proposed numerical scheme converges to the solutions of (7) as the discretization parameter tends to zero. The standard framework used to establish convergence is that of Barles and Souganidis [BS91], which we state below. In particular, it guarantees that the solutions of any monotone, consistent, and stable scheme converge to the unique viscosity solution of the PDE.

5.1. Degenerate elliptic schemes

Consider the Dirichlet problem for the degenerate elliptic PDE, $F[u]=0$, and recall its corresponding finite difference formulation:

[TABLE]

where $h$ is the discretization parameter.

Definition 5.1**.**

$F_{h}[u]$* is a degenerate elliptic scheme if it is non-decreasing in each of its arguments.*

Remark 5.1*.*

Although the convergence theory in [BS91] is originally stated in terms of monotone approximation schemes (schemes with non-negative coefficients), ellipticity is an equivalent formulation for finite difference operators [Obe06].

Definition 5.2**.**

The finite difference operator $F_{h}[u]$ is consistent with $F[u]=0$ if for any smooth function $\phi$ and $x\in\Omega$,

[TABLE]

Definition 5.3**.**

The finite difference operator $F_{h}[u]$ is stable if there exists $M>0$ independent of $h$ such that if $F_{h}[u]=0$ then $\|u\|_{\infty}\leq M$.

Remark 5.2* (Interpolating to the entire domain).*

The convergence theory assumes that the approximation scheme and the grid function are defined on all of $\Omega$. Although the finite difference operator acts only on functions defined on $\Omega^{h}$, we can extend such functions to $\Omega^{h}$ via piecewise linear interpolation. In particular, performing piecewise linear interpolation maintains the ellipticity of the scheme, as well as all other relevant properties. Therefore, we can safely interchange $\Omega$ and $\Omega^{h}$ in the discussion of convergence without any loss of generality

5.2. Convergence of numerical approximations

Next we will state the theorem for convergence of approximation schemes, tailored to elliptic finite difference schemes, and demonstrate that the proposed scheme fits in the desired framework. In particular, we will show that the schemes are elliptic, consistent, and have stable solutions.

Proposition 5.4** (Convergence of approximation schemes [BS91]).**

Let $u$ denote the unique viscosity of the degenerate elliptic PDE $F[u]=0$ with Dirichlet boundary conditions for which there exists a strong comparison principle. For each $h$, let $u^{h}$ denote the solutions of $F^{h}[u]=0$, where the finite difference scheme $F_{h}[u]$ is a consistent, stable and elliptic scheme. Then $u^{h}\to u$ locally uniformly on $\Omega$ as $h\to 0$.

Now, we check that our scheme is consistent, stable and elliptic.

Lemma 5.5** (Consistency).**

The scheme is consistent.

Proof.

It is sufficient to show that

[TABLE]

which follows immediately from a Taylor expansion argument and the definition of $u^{h}_{v}$. ∎

Lemma 5.6** (Stability).**

Suppose $g$ is bounded. Then the scheme is stable.

Proof.

Since the solution $u^{h}$ satisfies

[TABLE]

it follows that $g(x^{*})\leq u^{h}(x)\leq\|g\|$. ∎

Lemma 5.7** (Ellipticity).**

The scheme is elliptic.

Proof.

The term $u(x)-g(x)$ is trivially elliptic. Since $\tilde{x}$ belongs to the line segment of two grid points, $\mathcal{I}_{h}u(\tilde{x})$ is a convex combination of neighboring grid values and therefore $u^{h}_{v}$ is elliptic. Hence $F^{h}$, being the minimum of two elliptic schemes, is also elliptic. ∎

Theorem 5.8**.**

Suppose $g\in C(\Omega)$. Then the solutions of (9) converge locally uniformly on $\Omega$ as $h\to 0$ to the unique viscosity solution of (7) along with $u(x^{*})=g(x^{*})$.

Proof.

By Lemmas 5.5, 5.6 and 5.7, we see that $F^{h}$ (9) is consistent, stable and elliptic.

In order to apply Proposition 5.4, our PDE must satisfy a strong comparison principle which is not the case here. However, inspecting the proof of Proposition 5.4 in [BS91], one notices that it is enough to show that $\overline{u}\leq\underline{u}$ in $\Omega$, where

[TABLE]

This will follow from Proposition 3.9 by proving that $\overline{u}(x^{*})\leq g(x^{*})$. Indeed, due to the continuity of $g$,

[TABLE]

∎

6. Numerical results

In this section we present the numerical results. We focus on solving equation (7) as it is the one we are mainly interested since the sublevel sets of its solution provide us with the visibility set. We present results both in two and three dimensions.

Example 6.1**.**

We start with a simple example to show a numerical convergence test. We take $x^{*}=(-1,-1)$ as the viewpoint and consider as the obstacle function the cone $g(x)=-\sqrt{x_{1}^{2}+x_{2}^{2}}$ and therefore each level-set of $g$ corresponds to a circle with center at the origin and a different radius. The exact solution was obtained using (6). The difference between the numerical solution and the exact solution in the $l_{\infty}$ norm is presented in Table 1, where we also confirm the expected first order convergence. Moreover, in Figure 5, we plot the level sets of both $g$ and the numerical solution, as well as their respective surface plots. As we can see, the visibility set is computed for each level set of $g$.

Example 6.2**.**

In this example we are interested in computing the visibility set where we have four different obstacles: two squares with centers $(-1.5,-0.2)$, $(0,0.3)$ and side lengths $0.5$, $1$ respectively and two circles with origins $(-0.3,1.5)$, $(-0.3,-1.4)$ and radius $0.5$. We achieve this by considering the obstacle function

[TABLE]

where

[TABLE]

and looking into the $0.5$ level set. We first solve (7) with $x^{*}=(-1.5,-1.4)$ and $x^{*}=(1.5,-0.3)$. We compute as well the visibility set with respect to $\{x^{*}_{1},x^{*}_{2}\}=\{(-1.5,-1.4),(1.5,-0.3)\}$ when a point is consider visible if it is seen by any of the viewpoints and by both viewpoints simultaneous. The former corresponds to the solution of (8). All results are displayed in Figure 6.

Example 6.3**.**

We consider a simple three-dimensional example where the obstacle function is given by

[TABLE]

where

[TABLE]

It can be interpreted as the visibility set of a $360^{\circ}$ camera in the middle of two buildings. The results are displayed in Figure 7.

Example 6.4**.**

The Stanford 3D Scanning Repository provides a dataset of 35947 distinct points that form what is known as the “Stanford Bunny”. In this example we considered it as the obstacle by taking the function $g$ to be the signed distance function to the dataset points. The results are displayed in Figure 8.

7. Conclusions

In this article, we described a new simpler PDE to compute the visibility set from a given viewpoint given a set of known obstacles. We proposed a finite difference numerical scheme to compute its solution and showed its convergence. We discuss the generalization of the result to multiple viewpoints and present a PDE for the visibility set where a point is visible if it is seen by at least one of many viewpoints. Numerical examples of different visibility sets computed as the solution of the new proposed PDE in both two and three dimensions are presented.

Acknowledgments

This material is based upon work supported by the Air Force Office of Scientific Research under award number FA9550-18-1-0167. The second author thanks the hospitality of the Mathematics and Statistics department of McGill University during its visit where the work for this paper was carried out.