MODELLING QUASI-ORTHOGRAPHIC CAPTURES FOR SURFACE IMAGING

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).

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,

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

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

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

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

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}$.

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$.

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.

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

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 ).

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.

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).

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.

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]$.

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.

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$.

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 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$.

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.

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.

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).

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.

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.

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

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.

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

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

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.