# MODELLING QUASI-ORTHOGRAPHIC CAPTURES FOR SURFACE IMAGING

###### Abstract

Surveillance and surveying are two important applications of empirical
research. A major part of terrain modelling is supported by photographic
surveys which are used for capturing expansive natural surfaces using a wide
range of sensors- visual, infrared, ultrasonic, radio, etc. A natural surface
is non-smooth, unpredictable and fast-varying, and it is difficult to capture
all features and reconstruct them accurately. An orthographic image of a
surface provides a detailed holistic view capturing its relevant features. In
a perfect orthographic reconstruction, images must be captured normal to each
point on the surface which is practically impossible. In this paper, a
detailed analysis of the constraints on imaging distance is also provided. A
novel method is formulated to determine an approximate orthographic region on
a surface surrounding the point of focus and additionally, some methods for
approximating the orthographic boundary for faster computation is also
proposed. The approximation methods have been compared in terms of
computational efficiency and accuracy.

Index Terms—  Digital Elevation Map, Field of View, Surface Imaging,
Orthography, Curvatures, Imaging Surface

##  1 Introduction

Visual representation of topographical data or Digital Terrain Modelling has a
wide-spread application, be it Google Maps, Navigation, Geological surveys,
agriculture, disaster management, Astronomical studies and mapping of extra
terrestrial surfaces, Archaeological studies, Sociological studies or even
deciding national policies. The capturing technology has progressed a lot over
the past decade and nowadays a mixture of visual, radio, infrared, laser and
radar sensors are used to capture terrain information [3]. However, natural
terrains are highly uneven and for vast surfaces, when very small features are
the focus of studies, it is easy to lose them while reconstruction or even
while capturing. An orthographic projection or image of a surface is a
holistic representation of all its features. To construct an idealistic
orthographic projection of the whole surface, captures at all points need to
be taken separately. This approach is highly impractical, implausible and
incomprehensible due to resource limitations. Thus an approximation of local
orthographic region needs to be formulated for practical purposes and methods
to determine capture points for the coverage of the surface is also needed.

The physical parameters of an imaging system is guided by the optics and the
sensors of the device[2], among which, Field of View (FOV) and Working
Distance (WD) are important parameters to be considered for orthographic
imaging. Usually, an object-space telecentric lens[1] with the entrance pupil
at infinty, is used for eliminating the perception of depth and creating
orthographic images. A more commonly used technology is orthophotography.
Orthophotographs are commonly used in geographic information systems (GIS) as
it is ’map-accurate’. A digital elevation model (DEM) is often required to
create accurate orthophotos. An aerial image or a satellite image of a terrain
or surface which is geometrically corrected or orthorectified such that the
photo or image is essentially an orthographic projection of the terrain. An
orthophotograph can be used to measure distances accurately because it is an
almost accurate depiction of the Earth’s surface being adjusted for
topographic relief, perspective distortion and camera tilt.[4]

In this paper, a mathematical approach is taken to generate orthographic
captures of surfaces. In section 2, the derivation and analysis of imaging
surfaces is given and the mathematical bounds on working distance has been
formulated. Section 3 deals with the formulation of approximate orthography
and in section 4 the algorithms for computation of orthographic bounds for
curves and regions for surfaces have been provided. In section 5, several
methods for approximating orthographic boundary for faster computation are
proposed and compared and in section 6, the contributions of this paper have
been summarized along with potential future extensions.

(a) Gray-scale DEM

(b) Generated surface plot

Fig. 1: a. A gray-scale DEM denoting a terrain. (Source: Creating Heightfields
and Details on Terrain RAW Files, wiki.secondlife.com). b. Surface plot
generated in Matlab.

##  2 IMAGING SURFACE

###  2.1 Generating Surface Topographies

A digital terrain elevation map (DTEM) is a digital image of a terrain or a
topography map where the pixel intensity at a point gives the relative
elevation of the point. There are various ways of allocating pixel values
(Gray-scale or RGB colormap) in a DTEM. The images that were used for the
purpose of this thesis are Gray-scale DTEMs, where the elevation of a point or
location in the digital map can range from $0$ to $255$, the white intensity
pixels denoting the highest elevation points and the black intensity pixels
denoting the lowest intensity points. Fig. 1(a) is an example of a gray-scale
digital elevation map.

The obtained image of the height-map is first smoothed out because rough
surfaces with abrupt changes creates difficulty in further processing. The
double precision matrix $I$ is used as a topographical surface for later
processing. The surface generated in Matlab for the DTEM in figure 1(a) is
shown in figure 1(b).

###  2.2 Formulation

Let the working distance be denoted as $d$, i.e., it is assumed that to
capture a point $P(x,y,z)$ on the surface $S$ (given by the bi-variate
function $f(x,y)$), the camera needs to be placed at
$P^{\prime}(x^{\prime},y^{\prime},z^{\prime})$ at a height $d$ along the
normal to the surface at $P$. So,

$z=f(x,y).$ (1)

Then the surface normal at point $P(x,y,z)$ is given as

$\vec{n}=\left[\frac{\partial f(x,y)}{\partial x},\frac{\partial
f(x,y)}{\partial y},-1\right].$ (2)

If variables $p$ and $q$ are defined as

$p=\frac{\partial f(x,y)}{\partial x}\ \ and\ \ q=\frac{\partial
f(x,y)}{\partial y},$ (3)

then the surface normal can be written as $[p,q,-1]$. The unit normal vector
at $P$ is

$\displaystyle\hat{n}$
$\displaystyle=\frac{\vec{n}}{|{\vec{n}}|}=\frac{\vec{n}}{\sqrt{p^{2}+q^{2}+1}}$
(4)
$\displaystyle=\left[\frac{p}{\sqrt{p^{2}+q^{2}+1}},\frac{q}{\sqrt{p^{2}+q^{2}+1}},\frac{-1}{\sqrt{p^{2}+q^{2}+1}}\right].$

Now, using Eq. 4 as derived above, the corresponding imaging point,
$P^{\prime}(x^{\prime},y^{\prime},z^{\prime})$, at a height $d$ from point $P$
and along the unit surface normal $\hat{n}$ is given by

$\displaystyle\vec{P^{\prime}}$ $\displaystyle=\vec{P}+d\cdot\hat{n}$ (5)
$\displaystyle=\vec{P}+\left[\frac{d\cdot
p}{\sqrt{p^{2}+q^{2}+1}},\frac{d\cdot
q}{\sqrt{p^{2}+q^{2}+1}},\frac{-d}{\sqrt{p^{2}+q^{2}+1}}\right].$

Therefore the co-ordinates of point
$P^{\prime}(x^{\prime},y^{\prime},z^{\prime})$ can be derived using Eq. 5. The
imaging surface $S^{\prime}$ at imaging distance $d$ can be parameterized in
terms of $x$ and $y$ as shown in Eq. 6. Here $p=\frac{\partial
f(x,y)}{\partial x}$ and $q=\frac{\partial f(x,y)}{\partial y}$.

$\vec{S^{\prime}}=\left[\begin{array}[]{c}{x+\frac{d\cdot
p}{\sqrt{p^{2}+q^{2}+1}}}\\\ {y+\frac{d\cdot q}{\sqrt{p^{2}+q^{2}+1}}}\\\
{f(x,y)-\frac{d}{\sqrt{p^{2}+q^{2}+1}}}\end{array}\right]$ (6)

The surface plots in Fig. 2(f) demonstrate the imaging surfaces for the
surface $S$ given by $f(x,y)=cos(x)+cos(y)$ in the range $(-5\leq x\leq 5)$
and $(-5\leq y\leq 5)$ calculated and plotted at different values of $d$.

(a) Surface normal plot

(b) Contour plot

(c) $d=0.5$

(d) $d=0.5$

(e) $d=1$

(f) $d=1$

Fig. 2: The surface $S$ given by $f(x,y)=cos(x)+cos(y)$ and the imaging
surfaces $S^{\prime}$ plotted in Matlab

.

In case the surface $S$ is represented a double precision matrix ($I$) as
shown in section 2.1, then instead of calculating mathematical gradients
($\frac{\partial f(x,y)}{\partial x}$ and $\frac{\partial f(x,y)}{\partial
y}$) for finding surface normals, the numerical gradients can be calculated as
an approximation.

###  2.3 Analysis of Imaging Curves

If at any point on $S^{\prime}$ is below the surface $S$, then those points
are inaccessible and hence cannot be used as imaging points, i.e. if
$z^{\prime}<f(x^{\prime},y^{\prime})$,
$P^{\prime}(x^{\prime},y^{\prime},z^{\prime})$ cannot be an imaging point for
$P(x,y,z)$ at height $d$. Visualizing the variation of the imaging surface
with imaging height $d$ is difficult for bi-variate functions. Hence the
following analysis is done for imaging curves and virtual bounds on $d$ is
derived.

Let $C$ be the target curve given by $y=f(x)$. The curve can be parametrized
by vector $\vec{P}$ as

$\vec{P}(x)=\left[\begin{array}[]{c}{x}\\\
{y}\end{array}\right]=\left[\begin{array}[]{c}{x}\\\
{f(x)}\end{array}\right].$ (7)

Proceeding similarly as in section 2.2, the coordinates($x^{\prime}$ and
$y^{\prime}$) of imaging curve $C^{\prime}$ located at distance $d$ can be
parametrized in terms of $x$ (ref Eq.8 ).

$\displaystyle x^{\prime}$ $\displaystyle=x-\frac{d\cdot
f^{\prime}(x)}{\sqrt{1+(f^{\prime}(x))^{2}}}$ (8) $\displaystyle\text{and}\ \
\ y^{\prime}$ $\displaystyle=f(x)+\frac{d}{\sqrt{1+(f^{\prime}(x))^{2}}}.$

If a point $P^{\prime}(x^{\prime},y^{\prime})$ on curve $C^{\prime}$ is such
that it satisfies $y^{\prime}(x)<f(x^{\prime})$, then it lies below the curve
$C$($f(x)$) and hence it cannot be accepted as a valid imaging point. In other
terms for a $d$ to be valid, curves $C$ and $C^{\prime}$ should not intersect
at any point. This gives a mathematical bound for imaging height $d$ -

  * •

$d>0$

  * •

$d<D$ such that $\forall\ d\geq D,\ \exists$ some $x$ in $dom(f)$, s.t.
$y^{\prime}(x)<f(x^{\prime})$, where $x^{\prime}$ and $y^{\prime}$ are as
given in Eq.8.

The mathematical upper bound $D$ depends on the curvature or nature of the
function $f(x)$ and also the imaging range, i.e., the range of values of $x$
that is to be imaged. D can be calculated numerically by solving Eq.9 and
applying the bisection algorithm.

$\displaystyle y^{\prime}(x)=f(x^{\prime})$ (9) $\displaystyle or,$
$\displaystyle f(x^{\prime})=f(x)+\frac{d}{\sqrt{1+(f^{\prime}(x))^{2}}}$

For $d<D$, the above equation will have no solution and for $d\geq D$, the
above equation will have one or more solution(s).

Algorithm 1 Finding upper bound $D$ by Bisection

1:Set lower bound of $d$, $L=0$.

2:Set a very large upper bound $U$ such that by replacing $d=U$ in Eq.9, it
has a solution.

3:Set $D=(L+U)/2$.

4:Replace $d=D$ in Eq.9 and look for a solution.

5:if solution exists then

6: Set $U=D$

7:else

8: Set $L=D$

9:Repeat steps 5-8 until convergence criterion is met.

For some smooth functions, there may not be any upper limit on $d$ (i.e.
$D=\infty$). For those functions, Eq. 9 does not have a solution for any
$d>0$, i.e. $y^{\prime}(x)>f(x^{\prime})$ for all positive values of $d$ and
all $x$ in $dom(f)$. However, the practical upper bound depends on the
limitations of resolution of the capturing device and also the concerned
application.

(a) $f(x)=sin(x)$

(b) $f(x)=x^{2}$

(c) $f(x)=\exp{\sqrt{|x|}}$

Fig. 3: a. As $d$ increases, $C^{\prime}$ moves further away from $C$, thus
the upper bound on $d$, $D=\infty$. b. As $d$ increases, $C^{\prime}$ moves
further away from $C$. For lower values of $d$, $C^{\prime}$s do not intersect
$C$ and therefore are valid imaging curves. But for bigger values of $d$, they
intersect $C$ and hence are invalid. Here the upper bound on $d$, $D\approx
2.6$. c. Here $C$ is non-smooth and it is non-differentiable at $x=0$. In this
case, $C^{\prime}$ generated for any $d>0$ is invalid, as it always intersects
$C$, i.e. there are always points around $0$ which generate invalid imaging
points. Also, as $d$ increases, the number of invalid imaging points
increases.

The imaging surfaces for some curves have been computed and plotted in Matlab
as shown in figure 3. The target curves($C$) are shown in blue along with the
imaging curves($C^{\prime}$) for different values of $d$. It is interesting to
note that all non-smooth curves do not generate invalid imaging curves. For
example, $f(x)=|mx|$ is non-differentiable at $x=0$, but generates valid
imaging curves for all $m\in(0,1]$.

##  3 MODELLING ORTHOGRAPHIC IMAGING

###  3.1 Definition and Assumptions

A practical approximation of orthography is to consider a very
small($\epsilon$) angular field of view(FOV) and the points on the surface
within this $\epsilon$-FOV to be roughly orthographic and the boundary of the
region consisting of such points is the orthographic boundary. Here $\epsilon$
is a very small angle ($\sim\ 10^{\circ}-20^{\circ}$). Also, for orthographic
imaging of a surface point $P$, the imaging point must be along the surface
normal at $P$ and at a height $d$. The region in a capture which lies within
the orthographic boundary is defined as the $\epsilon$-Orthographic Image.

###  3.2 Circular Case

For a circle $C$ (Fig. 4a), that radius to a point on the circumference is
always orthogonal to the tangent at that point. Consequently, the centre $O$
of the circle satisfies the properties of a valid imaging point for any point
$P$ on the circle. So the imaging point at $O$ can be used to capture a length
of $2\pi R$ or the entire circumference as shown in figure below. The total
number of captures required to cover the entire circle is
$\lceil\frac{2\pi}{\epsilon}\rceil$. Here imaging height is $d=R$.

For an imaging point located at an eccentric point $Q$ (Fig. 4b) at distance
$x$ from $O$, the normal from only two diagonally opposite points on the
circumference passes through it. In this case the total length of $C$ that can
be imaged can be proved to be $2\epsilon R$.

(a) Imaging point at centre $O$

(b) Eccentric imaging point at $Q$

Fig. 4: The case of circle as the target curve.

###  3.3 Derivation of $\epsilon$-Orthography

####  3.3.1 Bounds for Curves

Let us consider a curve $C$ (Fig.5) given by a univariate function $f(x)$. The
tangent($\vec{T}$) and normal($\vec{N}$) vectors at point $P(x,f(x))$ are
given as

$\vec{T}(x)=\left[\begin{array}[]{c}{1}\\\ {f^{\prime}(x)}\end{array}\right]\
\ \ \ \vec{N}(x)=\left[\begin{array}[]{c}{-f^{\prime}(x)}\\\
{1}\end{array}\right].$ (10)

Let point $P^{\prime}(x^{\prime},f(x^{\prime}))$ be situated at a small
distance $\Delta x$ to the left of $x$. Let $p=f^{\prime}(x)$ and
$p^{\prime}=f^{\prime}(x^{\prime}).$If $\Delta x$ is very small then
$f(x^{\prime})$ and $f^{\prime}(x^{\prime})$ can be approximated as

$\displaystyle f(x^{\prime})$ $\displaystyle=f(x-\Delta x)\approx f(x)-\Delta
x\cdot f^{\prime}(x)$ (11) $\displaystyle f^{\prime}(x^{\prime})$
$\displaystyle=f^{\prime}(x-\Delta x)\approx f^{\prime}(x)-\Delta x\cdot
f^{\prime\prime}(x),$

therefore,

$\displaystyle p^{\prime}$ $\displaystyle=p+\Delta p$ (12)
$\displaystyle\Delta p$ $\displaystyle\approx-\Delta x\cdot
f^{\prime\prime}(x)$

Fig. 5: Illustrating the variables for derivation.

Tangent $\vec{T^{\prime}}$ and normal $\vec{N^{\prime}}$ are constructed at
$P^{\prime}$. The normals $\vec{N}$ and $\vec{N^{\prime}}$ intersect at $Q$ at
an angle $\phi$. Therefore,

$\displaystyle cos(\phi)$
$\displaystyle=\frac{\vec{N}\cdot\vec{N^{\prime}}}{|\vec{N}|\cdot|\vec{N^{\prime}}|}$
(13) $\displaystyle=\frac{1}{\sqrt{p^{2}+1}\sqrt{p^{\prime
2}+1}}\cdot\left[\begin{array}[]{c}{-p}\\\
{1}\end{array}\right]\cdot\left[\begin{array}[]{c}{-p^{\prime}}\\\
{1}\end{array}\right]$
$\displaystyle=\frac{pp^{\prime}+1}{\sqrt{p^{2}+1}\sqrt{p^{\prime 2}+1}}.$

So,

$\phi=cos^{-1}\Big{(}\frac{pp^{\prime}+1}{\sqrt{p^{2}+1}\sqrt{p^{\prime
2}+1}}\Big{)}.$ (14)

Also, if the line joining $P^{\prime}$ and the imaging point $O$ intersect
$\overline{OP}$ at angle $\theta$, as $\Delta x$ is much smaller compared to
$d$,

$\displaystyle tan(\theta)=\frac{\Delta x}{d}$ (15) $\displaystyle or,$
$\displaystyle\theta=tan^{-1}\big{(}\frac{\Delta x}{d}\big{)}.$

$\epsilon$-orthographic bounds are dependent on both the FOV and the curvature
at the concerned point. For a point $P^{\prime}$ on $C$ to lie within the
$\epsilon$-orthographic region for capturing point $O$ at a height $d$ from
the point $P$, it must satisfy- 1. $\theta\leq\epsilon$, and 2.
$\phi\leq\epsilon$.

Condition 1 is necessary so that the point $P^{\prime}$ lies within the
$\epsilon$-FOV. Condition 2 is required because as curvature of $C$ increases
around $P$, although a point close to it may remain within the $\epsilon$-FOV
bound, the high curvature causes very small region around $P$ to be
approximately orthographic. With reference to Fig. 5, if curvature at $P$
increases, $OP$ and $OP^{\prime}$ may differ very much and then $P$ and
$P^{\prime}$ cannot be considered in the same orthographic region.

####  3.3.2 Boundary for Surfaces

The derivation is very similar to that for the curves. A surface $S$ is given
by a bi-variate function, $z=f(x,y)$. Given a central point $P(x,y,z)$ on the
surface, an imaging height of $d$ and useful FOV $\epsilon$, the goal is to
find the orthographic boundary surrounding $P$, or in other words, the area
around $P$ which can be considered as an approximate orthographic image.

From Eqs. 2 and 3, the surface normal at point $P(x,y,z)$ is $[p,q,-1]$, where
$p=\frac{\partial f(x,y)}{\partial x}$ and $q=\frac{\partial f(x,y)}{\partial
y}$. The Hessian matrix of $f(x,y)$ is

$H=\left[\begin{array}[]{cc}{\frac{\partial^{2}f(x,y)}{\partial
x^{2}}}&{\frac{\partial^{2}f(x,y)}{\partial x\partial y}}\\\
{\frac{\partial^{2}f(x,y)}{\partial y\partial
x}}&{\frac{\partial^{2}f(x,y)}{\partial y^{2}}}\end{array}\right].$ (16)

Let us take a point $P^{\prime}(x^{\prime},y^{\prime},z^{\prime})$ very close
to $P$ such that

$\displaystyle x^{\prime}$ $\displaystyle=x+\Delta x$ (17) $\displaystyle
y^{\prime}$ $\displaystyle=y+\Delta y.$

As $\Delta x$ and $\Delta y$ are very small quantities, the change in the
surface normal vector is also small. So the surface normal at $P^{\prime}$,
$\vec{N^{\prime}}=[p^{\prime},q^{\prime},-1]$ and it can be approximated as

$\displaystyle\left[\begin{array}[]{cc}{\Delta p}\\\ {\Delta
q}\end{array}\right]=H\cdot\left[\begin{array}[]{cc}{\Delta x}\\\ {\Delta
y}\end{array}\right]$ (18) $\displaystyle p^{\prime}=p+\Delta p$
$\displaystyle q^{\prime}=q+\Delta q$

Similarly as Eqs. 13 and 15, $\phi$ and $\theta$ are calculated as

$\displaystyle cos(\phi)$ $\displaystyle=\hat{N}.\hat{N^{\prime}}$ (19)
$\displaystyle=\frac{\vec{N}\cdot\vec{N^{\prime}}}{|\vec{N}||\vec{N^{\prime}}|}$
$\displaystyle=\frac{1}{\sqrt{p^{2}+q^{2}+1}\sqrt{p^{\prime 2}+q^{\prime
2}+1}}\left[\begin{array}[]{c}{p}\\\ {q}\\\
{-1}\end{array}\right]\cdot\left[\begin{array}[]{c}{p^{\prime}}\\\
{q^{\prime}}\\\ {-1}\end{array}\right]$
$\displaystyle=\frac{pp^{\prime}+qq^{\prime}+1}{\sqrt{p^{2}+q^{2}+1}\sqrt{p^{\prime
2}+q^{\prime 2}+1}}.$

So,

$\phi=cos^{-1}\Big{(}\frac{pp^{\prime}+qq^{\prime}+1}{\sqrt{p^{2}+q^{2}+1}\sqrt{p^{\prime
2}+q^{\prime 2}+1}}\Big{)}$ (20)

and

$\displaystyle tan(\theta)=\frac{\sqrt{\Delta x^{2}+\Delta y^{2}}}{d}$ (21)
$\displaystyle or,$ $\displaystyle\theta=tan^{-1}\Big{(}\frac{\sqrt{\Delta
x^{2}+\Delta y^{2}}}{d}\Big{)}.$

Now, if $P^{\prime}$ belongs to the orthographic region around point $P$ for
an imaging height $d$, then both $\theta\leq\epsilon$ and $\phi\leq\epsilon$.
In all the aforementioned derivations, the gradient components at $P^{\prime}$
have been calculated by approximations for fast computation. Where
computational capability is not an issue, the actual gradients can be
calculated by differentiating the function $f(x,y)$.

##  4 IMPLEMENTATION

###  4.1 Algorithms

The algorithm for computing the orthographic bounds for a smooth curve $f(x)$
at a central point $P_{0}(x_{0},f(x_{0}))$, for $\epsilon$ angular FOV,
resolution $dx$ and imaging height $d$ is stated.

Algorithm 2 Finding Orthographic Bounds of a Curve

1:procedure Left Orthographic Bound($x_{0},\epsilon,d,dx$)

2: Find $p_{0}=f^{\prime}(x_{0})$

3: Set $x=x_{0}-dx$

4: Find $p=f^{\prime}(x)$ and calculate $\phi$ using Eq.14

5: Set $\Delta x=|x-x_{0}|$

6: Calculate $\theta$ using Eq.15

7: while $\phi\leq\epsilon$ and $\theta\leq\epsilon$ do

8: Set $x=x-dx$

9: Find $p=f^{\prime}(x)$ and calculate $\phi^{\prime}$ using Eq.14

10: $\phi=\phi^{\prime}$’

11: Set $\Delta x=|x-x_{0}|$

12: Calculate $\theta^{\prime}$ using Eq.15

13: $\theta=\theta^{\prime}$

14: return $x=x+dx$

15:procedure RightOrthographic Bound($x_{0},\epsilon,d,dx$)

16: Find $p_{0}=f^{\prime}(x_{0})$

17: Set $x=x_{0}+dx$

18: Find $p=f^{\prime}(x)$ and calculate $\phi$ using Eq.14

19: Set $\Delta x=|x-x_{0}|$

20: Calculate $\theta$ using Eq.15

21: while $\phi\leq\epsilon$ and $\theta\leq\epsilon$ do

22: Set $x=x+dx$

23: Find $p=f^{\prime}(x)$ and calculate $\phi^{\prime}$ using Eq.14

24: $\phi=\phi^{\prime}$’

25: Set $\Delta x=|x-x_{0}|$

26: Calculate $\theta^{\prime}$ using Eq.15

27: $\theta=\theta^{\prime}$

28: return $x=x-dx$

(a) Convex part

(b) Concave part

Fig. 6: The orthographic bounds for convex and concave parts of a sine curve
for increasing $d$.

Algorithm 2 has been implemented in Matlab and the bounds have been calculated
and plotted for a convex curve and a concave curve as shown in Fig. 6. It is
to be noted that although $d$ keeps increasing, the bounds do not spread after
a point. If the bounds were a function of $\epsilon$ only, then with increase
in $d$, they would have spread apart indefinitely which is not a true
characteristic of orthography. This shows that $\epsilon$-orthography not only
depends on the FOV but also the curvature.

The $\epsilon$-orthographic boundary can be numerically calculated using
Algorithm 3 for a surface $S$ ($z=f(x,y)$) at a central point
$P_{0}(x_{0},y_{0},z_{0})$, for $\epsilon$ angular FOV, resolutions $dx$ and
$dy$, and imaging height $d$.

Algorithm 3 Finding Orthographic Boundary of a Surface

1:Find surface normal components $p_{0}$ and $q_{0}$ at $P_{0}$.

2:Create an empty point set $P$ for storing the eligible points inside the
orthographic region.

3:Append $P_{0}$ to $P$

4:Compute the number of $x$ or $y$ co-ordinates in the grid.
$n_{x}=(r_{H}-r{L})/dx=R_{x}/dx$ and $n_{y}=R_{y}/dy$.

5:Set $s=max(n_{x},n_{y})$

6:Set $buff=0$

7:for $n=1:s$ do

8: Vector $out=PairGen(n)$

9: Set $count=0$

10: for $j=1:length(out)$ do

11: $x_{1}=x_{0}+dx\cdot out(1,j)$

12: $y_{1}=y_{0}+dy\cdot out(2,j)$

13: Compute surface normal vector at $P_{1}(x_{1},y_{1},z_{1})$

14: Calculate $\phi$ using Eq.20

15: Calculate $\theta$ using Eq.21

16: if $\theta\leq\epsilon\ \ and\ \ \phi\leq\epsilon$ then

17: Append $P_{1}(x_{1},y_{1},z_{1})$ to

18: $count=count+1$

19: if $count=0$ then

20: $buff=buff+1$

21: if $buff>3$ then

22: Break from loop

(a) $d=2,\ (x_{0},y_{0})=(0,0)$

(b) $d=2,\ (x_{0},y_{0})=(0,0)$

(c) $d=2,\ (x_{0},y_{0})=(0,-1)$

(d) $d=2,\ (x_{0},y_{0})=(0,-1)$

(e) $d=2,\ (x_{0},y_{0})=(-1,-1)$

(f) $d=2,\ (x_{0},y_{0})=(-1,-1)$

Fig. 7: Orthographic regions drawn on curve $f(x,y)=cos^{2}(x)+cos^{2}(y)$
shown in white. The figures on right show the boundary shape. The central
point $(x_{0},y_{0})$ is plotted in red. (Here $\epsilon=10^{\circ}$)

The function PairGen generates a vector of all pairs of integers $n_{1}$ and
$n_{2}$ such that $|n_{1}|+|n_{2}|=n$, i.e. all co-ordinates located at
absolute distance $n$. Algorithm 3 has been implemented in Matlab to generate
$\epsilon$-orthographic regions on smooth curves (Fig. 7). This algorithm is
also valid for non-smooth curves, for which numerical gradients can be
calculated at non-differentiable points.

###  4.2 Special Surfaces

Conjecture: Points on surfaces of constant Gaussian curvature [6] form
$\epsilon$-orthographic regions of same area for constant imaging height $d$.
The upper bound on $d$ depends on the nature (parameters) of such surfaces.

Surfaces of constant curvatures can be classified into the following three
classes-

1) Zero Curvature Surfaces:  A surface with Gaussian curvature($\kappa$) equal
to zero at all points is a plane. For a plane, which is inherently
orthographic, the calculated region is thus of same shape as the FOV. As the
FOV is considered circular in all our calculations, the orthographic region is
thus circular for a planar surface, the radius of which depends on the imaging
height as given by Eq.22. Thus the problem of finding optimal capture points
is reduced to a Circle Packing problem.

2) Positive Curvature Surfaces:  A surface with equal positive Gaussian
curvature($\kappa$) at all points is sphere. Using the definition of
$\epsilon$-orthography, it can be shown that for a sphere, the orthographic
regions are also circular and of constant radii, dependent on the imaging
height $d$ (Fig.8a–b). This is due to the fact that the two principal
curvatures ($\kappa_{1}$ and $\kappa_{2}$) at any point on the sphere are
equal and constant. However, unlike a plane, sphere is not inherently
orthographic but behaves like one.

(a) 3D view of the surface

(b) 2D top view of the surface

(c) 3D view of the surface

(d) 2D top view of the surface

Fig. 8: a–b. $\epsilon$-Orthographic regions plotted on a sphere- a surface of
constant positive Gaussian curvature. c–d. $\epsilon$-Orthographic regions
plotted on a pseudosphere- a surface of constant negative Gaussian curvature
(large regions are plotted for demonstration).

3) Negative Curvature Surfaces:  A surface with equal negative Gaussian
curvature($\kappa$) at all points is a pseudosphere [5]. Unlike the other two
cases, for a pseudosphere, the orthographic boundaries are not circular, and
the limit on imaging height $d$ is dependent on the radius $a$ (Fig.8c–d).
Among the two principal curvatures ($\kappa_{1}$ and $\kappa_{2}$) calculated
at any point on the surface, one is positive and the other is negative. As we
move along the surface, from the flat region to the narrow region, the
magnitude of the positive principal curvature increases and the negative
principal decreases such that the product of the two (Gaussian curvature)
remains the same. As a consequence of this property, it has been empirically
observed that the size of the $\epsilon$-orthographic regions remain the same,
although the shape may vary (figure 8c–d).

###  4.3 Limitations

The calculation of orthographic regions is computationally expensive,
specially for higher resolutions. Consequently, finding the overlap between
two regions is also expensive. Also, the entire orthographic region needs to
be calculated to find the boundary. For non-smooth surfaces, gradients cannot
be calculated at non-differentiable points and for natural surfaces, with fast
varying curvatures, orthographic boundaries are difficult to calculate.
Moreover, the exact boundaries cannot be used to calculate optimal capture
points, where the boundaries need to be computed at multiple points
simultaneously, which is repetitive and slow.

##  5 APPROXIMATION OF ORTHOGRAPHIC BOUNDARY

###  5.1 Approaches

As pointed out in the limitations (4.3), calculating the exact orthographic
region and hence the boundary is not practical because it is computationally
expensive and thus time consuming. However, instead of considering exact
boundaries, they can be approximated to some regular shapes for faster
boundary computation as well as calculating overlap between regions as
explored in the following approaches.

(a) Polygonal

(b) Elliptical

(c) Circular-I

Fig. 9: a. Boundary points detected for $N=8$. Here $\theta=45^{\circ}$.
Actual boundary $C$ shown in bold. b. Actual boundary $C$ and approximated
elliptical boundary $B$. c. Actual boundary $C$ shown in blue and approximated
circular boundary $B$ shown in red.

1) Polygonal Approximation:  In this approach, $N$ points are calculated in
$N$ different directions from the central point. By setting $\Delta x$ and
$\Delta y$ in Eq.20 according to the direction of calculation, and by
calculating $\theta$ and $\phi$ using Eqs.20 and 21, and checking at each step
to see whether they maintain the $\epsilon$ constraint, the boundary point in
the concerned direction can be computed. The directions in which the boundary
points must be calculated, should be at equal angles to each other at the
central point $(x_{0},y_{0})$, i.e. the directions should be at an angle
$\theta=\frac{360^{\circ}}{N}$ from each other. The $N$ boundary points are
joined to form the $N$-polygonal boundary.

Fig.9(a) shows an example of polygonal approximation of the boundary. Here the
boundary points $P_{i}$’s are evaluated in 8 directions centered at
$O(x_{0},y_{0})$. Obviously, larger the number of directions taken, better
will be the approximate boundary. This approach may lead to both over-
estimation and under-estimation of boundary depending on the convexity of the
boundary curve.

2) Elliptical Approximation:  This is an extension or further approximation of
the polygonal approximation approach but here only even sided polygons are
considered. Boundary points are calculated in $N$ different equiangularly
spaced directions. Hence, we get $N$ boundary points $P_{i},\ \ i=1,2,...,N$.
Now, the distances between the diagonally opposite boundary points is
calculated, and thus we have $N/2$ diagonals ($d_{i}$). The maximum and
minimum diagonals are considered, $d_{max}=max(d_{i})$ and
$d_{min}=min(d_{i})$. The boundary is approximated as an ellipse with the
major axis as $d_{max}$ and the minor axis as $d_{min}$ and the major axis is
aligned along the longest diagonal.

In Fig.9(b), the maximum length diagonal is $\overline{P_{4}P_{8}}$ and the
minimum length diagonal is $\overline{P_{2}P_{6}}$. The major axis of the
constructed ellipse($B$) is $\overline{P_{4}P_{8}}$ and the minor axis is of
same length as $\overline{P_{2}P_{6}}$. It is to be noted that central point
$O(x_{0},y_{0})$ is not the centre of the ellipse.

3) Circular Approximation–I:  This is a further simplification of the
elliptical approximation. In this case, the orthographic boundary is
approximated as a circle. Boundary points are calculated in $N$ different
equiangularly spaced directions as discussed in polygonal case. Hence, we get
$N$ boundary points $P_{i},\ \ i=1,2,...,N$. Now, the distances of the
boundary points from the central point $(x_{0},y_{0})$ is calculated, and thus
we have $N$ distances ($d_{i}$). The average of all the distance lengths is
calculated, $d_{avg}=\Sigma d_{i}$.The boundary is approximated as a circle
centered at $(x_{0},y_{0})$ and of radius $R=d_{avg}$. In the illustrated
figure 9(c), $d_{i}=\overline{OP_{i}}$, and the average of all 8
$\overline{OP_{i}}$’s is calculated. The boundary circle $B$ is constructed
with centre at $O(x_{0},y_{0})$ and radius equal to the average of $OP_{i}$’s.

4) Circular Approximation–II:  This approach is the simplest and most
intuitive approach among the ones discussed. It is based on the observation
that on a smooth surface, the orthographic regions are small at places of high
Absolute Gaussian Curvature [6] and relatively larger at places where it is
low.

If the orthographic region or boundary is estimated by a circle, an inverse
relation between the radius of the boundary and the curvature of the central
point can be formulated. Also, for a planar surface or a zero curvature
surface, the orthographic region is circular with radius

$R=d\cdot tan(\epsilon),$ (22)

where $d$ is the imaging distance and $\epsilon$ is the useful FOV as
discussed in the derivation of $\epsilon$-orthography. For any point on a non-
planar surface having an absolute curvature $|K|\geq 0$, the boundary will
shrink from this circle. Therefore, if the region boundary is approximated by
a circle of radius $r$, $r\leq R$.

Now, let us consider a surface $S$ and its Gaussian curvature ($K$). Given the
bounds of the surface, the maximum absolute curvature is calculated.

$K_{max}=\max_{x,y}|K(x,y)|,\ \ \ \ (x,y)\ \in\ Bound(S).$ (23)

Fix a ratio ($m$) between the largest radius possible $(r_{max}=R)$ for the
points of least absolute curvature and the least radius possible for the point
having curvature $K_{max}$. So, $m=\frac{r_{max}}{r_{min}}$. Therefore, the
radius of the approximated circular boundary can be expressed as a function of
point $(x,y)$ as

$r(x,y)=R-\frac{|K(x,y)|}{K_{max}}\cdot R\cdot(1-\frac{1}{m}).$ (24)

The value of $m$ can be tuned by experimental observations.

###  5.2 Comparison and Analysis

The four different approaches for approximating the orthographic boundary can
be compared in terms of accuracy of approximation and computation time. For
Approach 1 (Polygonal Approximation), 2 (Elliptical Approximation) and 3
(Circular Approximation - I), the computational time depends on the number of
directions $N$ or the number of boundary points used for approximation, which
increases with increase in $N$. However in Approach 2, the diagonals have to
be compared and the equation of ellipse has to be calculated, which is the
most time consuming among the four. For Approach 1 and 2, calculation of
overlap among the regions is most complicated because of their polygonal and
elliptical shape respectively, whereas, in Approach 3 and 4, it is much easier
due to their circular shape. The accuracy of approximation decrease from 1 to
4, polygonal being the most accurate and the second circular approximation
being the crudest. However, Approach 4 is the most intuitive and easiest to
compute and can be used for further optimization but has the highest error,
specially for natural or fast-varying surfaces.

##  6 CONCLUSION AND FUTURE WORK

Orthographic imaging is a crucial tool for terrain survey or terrain mapping.
Although technological improvements have been made widely in the devices to
capture visual or other sensor data of a surface, a proper and efficient
algorithm for reconstructing the surface topography and creating an
orthographic projection of the terrain is lacking. This paper aims to study
and analyze this problem and propose novel methods for solving it. A technique
for generating topographical surface from elevation maps has been proposed. A
detailed study of imaging surfaces and bounds on imaging height has been
presented. The effects of imaging height and angular field of view for
capturing orthographic views have been formulated and analyzed in detail. A
novel method for calculating orthographic boundaries have been proposed and
demonstrated. Different methods of approximating the orthographic boundaries
have been proposed and compared.

As a future extension of this study, better approximations of orthographic
boundaries should be explored, and in case of which, how and by how much the
results are affected must be studied. Methods for computing orthographic views
by combining visual data with that from devices not capturing visual data can
be explored and incorporated. Faster algorithms for real-time computation of
orthographic boundaries should be explored for large-scale computation
problems. The approximation can be used for computing the optimal number of
capture points to cover the whole surface with minimum overlap among them.

## References

  * [1] Telecentric lenses: Basic Information and Working principles.  Opto Engineering, 2008.
  * [2] Gregory Hollows.  Imaging Optics Fundamentals.  Proceedings of the 17th ACM SIGSPATIAL international conference on advances in geographic information systems, pages 1–16, 2014.
  * [3] Steven Scheding, Jeff Leal, Mark Bishop, and Salah Sukkarieh.  Terrain Mapping in Real-Time: Sensors and Algorithms.  Geospatial Information and Agriculture, 2001.
  * [4] Gary S. Smith.  Digital Orthophotography and GIS.  ESRI Conference.
  * [5] H Steinhaus.  Mathematical Snapshots, $3^{rd}$ ed. New York: Dover, 1999.
  * [6] Eric W. Weisstein.  Gaussian Curvature.  Mathworld– A Wolfram Web Resource.