\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
Tenth MSU Conference on Differential Equations and Computational
Simulations. \newline
\emph{Electronic Journal of Differential Equations},
Conference 23 (2016),  pp. 47--57.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2016 Texas State University.}
\vspace{9mm}}

\begin{document} \setcounter{page}{47}
\title[\hfilneg EJDE-2016/Conf/23 \hfil Digital elevation modeling]
{Digital elevation modeling via curvature interpolation for LiDAR data}

\author[H. Kim, J. L. Willers, S. Kim \hfil EJDE-2016/conf/23 \hfilneg]
{Hwamog Kim, Jeffrey L. Willers, Seongjai Kim}

\address{Hwamog Kim \newline
 Department of Mathematics and Statistics,
 Mississippi State University, \newline
 Mississippi State, MS 39762, USA}
\email{hk404@msstate.edu}

\address{Jeffrey L. Willers \newline
 USDA-ARS,  Genetics and Precision Agriculture Research, \newline
 Mississippi State, MS 39762, USA}
\email{jeffrey.willers@ars.usda.gov}

\address{Seongjai Kim \newline
 Department of Mathematics and Statistics,
 Mississippi State University, \newline
 Mississippi State, MS 39762, USA}
\email{skim@math.msstate.edu}

\thanks{Published March 21, 2016}
\subjclass[2010]{65M06, 62H35, 65D05}
\keywords{Digital elevation model; curvature interpolation method (CIM);
\hfil\break\indent surface reconstruction; point cloud data}

\begin{abstract}
 Digital elevation model (DEM) is a three-dimensional (3D) representation
 of a terrain's surface -- for a planet (including Earth), moon, or asteroid
 -- created from point cloud data which measure terrain elevation.
 Its modeling requires surface reconstruction for the scattered data,
 which is an ill-posed problem and most computational algorithms become
 overly expensive as the number of sample points increases.
 This article studies an effective partial differential equation
 (PDE)-based algorithm, called the \emph{curvature interpolation method} (CIM).
 The new method iteratively utilizes curvature information,
 estimated from an intermediate surface, to construct a \emph{reliable}
 image surface that contains all of the data points.
 The CIM is applied for DEM for point cloud data acquired by
 light detection and ranging (LiDAR) technology.
 It converges to a piecewise smooth image,
 requiring $\mathcal{O}(N)$ operations independently of the number
 of sample points, where $N$ is the number of grid points.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\allowdisplaybreaks

%\def\curvatureq#1{\nabla\cdot\Big(\frac{\nabla#1}{|\nabla#1|^q}\Big)}
%\def\curvature#1{\nabla\cdot\Big(\frac{\nabla#1}{|\nabla#1|}\Big)}
%\def\curvaturew#1{\nabla\cdot\Big(\frac{\nabla#1}{|\nabla#1|^{\omega}}\Big)}
%\def\curvaturenu#1{\nabla\cdot\Big(\frac{\nabla#1}{|\nabla#1|^{\nu}}\Big)}



\section{Introduction} \label{Introduction}

Point clouds are gained by scanning \emph{three-dimensional} (3D) objects
using various measuring techniques.
The point cloud represents the set of points, each of which is defined
by $(x,y,z)$ coordinates.
Point clouds are used for many purposes, including 3D
\emph{computer-aided design} (CAD) modeling for manufactured parts
(reverse engineering), meteorology/quality inspection, visualization,
animation, mass customization applications, and geosciences
\cite{Glassner:95,Watson:92}.
In applications, point clouds are usually converted to polygon mesh
or triangle mesh models, NURBS surface models, or CAD models through a
process commonly referred to as \emph{surface reconstruction}.

There are many techniques for surface reconstruction.
Some approaches build a network of triangles over the existing vertices
of the point cloud (Delaunay triangulation,
marching triangles \cite{HiltonETAL:96}, and
ball-pivoting \cite{BernardiniETAL:99}),
while other approaches convert the point cloud into a volumetric distance
field and reconstruct the implicit surface through a marching
cubes algorithm \cite{HoppeETAL:92}.
However, in practice, the most cited are polynomial
interpolation such as nearest-neighbor, linear, and cubic methods
\cite{LehmannGonnerSpitzer:99,Turkowski:90} due to their simplicity;
these methods are easy to implement, but offer only low-quality results.
Inverse-distance methods \cite{Shepard:68} are also used, although they
are computationally expensive and become impractical as the number of
samples increases.
Another common interpolation model for point cloud data is the method
of thin-plate splines, which is based on radial basis functions
\cite{TurkOBrien:99};
the method is hard to be practical due to a high computational
complexity.
See \cite{ChenSuter:98,FaulPowel:99} for efforts for the
reduction of computational complexity of the method.
In \emph{light detection and ranging} (LiDAR) data processing,
a frequently-used interpolation algorithm is the
\emph{inverse-distance weighting} (IDW) method \cite{Shepard:68}.

Surface (image) reconstruction from point cloud (arbitrarily sampled)
data is a challenging problem particularly when no constraint is imposed
on their locations.
The problem is ill-posed, first of all, and numerical methods solving
its optimization formulation become overly expensive as the number of
sample points increases \cite{ArigovindanETAL:05,BeatsonETAL:00}.
Furthermore, it is often the case that the constructed image is not an
interpolator but an approximator, i.e., the reconstructed surface may
not include all the data values.

In this article, the authors are interested in a novel PDE-based method
called the \emph{curvature interpolation method} (CIM), for surface
reconstruction of terrain elevation data acquired by
LiDAR technology.
The new method utilizes a curvature information which is estimated
from an intermediate surface of the point cloud data and plays a role of
driving force for the construction of a reliable image surface.
It is often the case that the constructed image surface does not contain all
of the data values due to the estimated curvature.
However, the misfit can be corrected by a recursive application of the
CIM.
The CIM is first studied for image zooming by one of the authors;
see \cite{ChaLeeKim:14_RCIM} and \cite{KimChaKim:11_CIM}.
It has been verified that the CIM shows a \emph{minimum} oscillatory
behavior, and yet it results in piecewise smooth images containing all
the data values.

The article is organized as follows.
The next section briefly reviews the CIM and its recursive application,
as preliminaries.
Section~\ref{DEM} presents our digital elevation modeling strategies
including an effective method for the reduction of Moir\'e effect
inheritent in LiDAR data.
Section~\ref{Numerical} gives a set of numerical experiments, showing
effectiveness of the suggested algorithm.
At the end of this section, we summarize and
conclude our experiments.

\section{Preliminaries} \label{Preliminaries}

This section presents a brief review for
the curvature interpolation method (CIM)
studied for image zooming \cite{ChaLeeKim:14_RCIM,KimChaKim:11_CIM}.

\subsection{Curvature interpolation method (CIM)} \label{ss_CIM}

Image zooming is a processing task to enlarge images by applying
interpolation methods.
The CIM was a PDE-based model; it begins with a selection of a
curvature-related term which is to be estimated from the low resolution
(LR) image, interpolated to the high resolution (HR) image grid, and
incorporated as a driving force for the construction of HR image.
PDE-based models that employ the (mean) curvature itself as the
smoothing operator (e.g., the total variation (TV) model
\cite{RudinOsherFatemi:92})
are known to have a tendency to converge to a piecewise constant image
\cite{ChanShen:05,MarquinaOsher:00}.
Such a phenomenon is called  \emph{staircasing}.
Thus the curvature would better be replaced by a more effective and
convenient operator $\mathcal{K}$.
In \cite{KimChaKim:11_CIM}, the authors adopted the following
\emph{gradient-weighted} (GW) curvature 
\begin{equation}
 \mathcal{K}(u)= -|\nabla{u}|
\nabla\cdot\Big(\frac{\nabla u}{|\nabla u|}\Big).
\label{eqn_K}
\end{equation}

Let $\Omega$ and $\widetilde{\Omega}$ be the original LR image domain and its
$\alpha$-times magnified HR image domain, $\alpha>1$, respectively.
Let $\widehat{\Omega}^0$ denote the set of pixel points in $\widetilde{\Omega}$ which can be
expressed as $\alpha\mathbf{p}$, where $\mathbf{p}$ is a pixel point in $\Omega$.
Then the CIM \cite{KimChaKim:11_CIM} can be outlined as follows.
\begin{itemize}
\item[(I)]
Compute the GW curvature of the given LR image $v^0$:
\begin{equation}
 \mathcal{K}=\mathcal{K}(v^0) \quad \text{on }\Omega.
\label{eqn_Ku0}
\end{equation}

\item[(II)]
Interpolate $\mathcal{K}$ to obtain $\widehat{\mathcal{K}}$ on $\widetilde{\Omega}$.

\item[(III)]
Solve, for $u$ on $\widetilde{\Omega}$, the following constrained problem
\begin{equation}
 \mathcal{K}(u)=\frac{1}{\alpha^2}\widehat{\mathcal{K}} \quad u|_{\widehat{\Omega}^0}=v^0.
\label{eqn_Curvature_Solve}
\end{equation}
\end{itemize}

In the above algorithm, the GW curvature measured from the LR image is
interpolated and incorporated as an explicit driving force for the same
GW curvature model on $\widetilde{\Omega}$.
The driving force would help the model construct the HR image more
effectively, enforcing the resulting image to satisfy the given curvature
profile.


However, because the involved interpolation operations, the constructed
surface (the solution of \eqref{eqn_Curvature_Solve}) may have image
values different from those in the corresponding LR grid points.
A natural remedy for this drawback is to update image values iteratively
by utilizing the difference between the LR image and the last updated
image in the LR grid (misfit).

\subsection{CIM with iterative refinement} \label{ss_IR-CIM}

When the CIM is applied for image zooming as in \S\ref{ss_CIM},
the curvature of the LR image is first computed and then interpolated
for an approximation of the curvature of the HR image.
Such an approximated curvature shows a reasonable accuracy so that the
CIM in each iteration can produce a reliable correction term to update
the image surface.
In image zooming, the given image may be viewed as an LR approximation
of the target HR image.
However, the case is quite different for the surface reconstruction
for point cloud data, because the data loci are nonuniform and
it is hard to estimate the curvature.
Thus we first have to introduce an efficient scheme to estimate the
surface and its curvature.
As a strategy, we will construct an \emph{intermediate surface}, from
which useful curvature information would be obtained.

Let $\Omega$ be the image domain and $\Omega^0$ the set of data pixels
where image values are initialized as $v^0$.
Our new surface reconstruction algorithm to be presented below involves
three major steps: the construction of an intermediate surface,
the curvature evaluation and smoothing, and surface reconstruction.
When the reconstructed surface does not contain all of the prescribed
image values ($v^0$), the difference can be corrected by applying
the procedure iteratively.
The following outlines our reconstruction algorithm.
\begin{itemize}
\item[] Initialize $\mathbf{u}_0=0$, on the image domain $\Omega$
\item[] Select the tolerance $\tau>0$
\item[] For $k=1,2,\cdots$
\begin{itemize}
\item[(i)] Compute the misfit on $\Omega^0$: 
\[
  \mathbf{r}_{k-1}=\mathbf{v}^0-\mathbf{u}_{k-1}
\]
 \item[(ii)]
 If $\|\mathbf{r}_{k-1}\|<\tau$, stop

 \item[(iii)]  Construct an intermediate surface $\phi_k$:
\begin{equation} 
\begin{gathered}
 -\Delta \phi_k=0, \quad \mathbf{x}\in\Omega\setminus\Omega^0 \\
 \phi_k=\mathbf{r}_{k-1}, \quad \mathbf{x}\in\Omega^0 \\
 \nabla\phi_k\cdot\mathbf{n}=0, \quad \mathbf{x}\in\partial\Omega
 \end{gathered} \label{eqn_R-CIM_PC}
\end{equation}

 \item[(iv)]  Evaluate $A_k$ and $K_k$ from $\phi_k$:
\[  
K_k=A_k\boldsymbol{\phi}_k \approx \mathcal{K}(\phi_k)
\]
\item[(v)]  Smoothen $K_k$ to get $\widetilde{K}_k$

\item[(vi)]  Construct the correction surface $\mathbf{w}_k$:
\[ 
A_k\mathbf{w}_k=\widetilde{K}_k 
\]
\item[(vii)] Update:
\[
 \mathbf{u}_k=\mathbf{u}_{k-1}+\mathbf{w}_k
\]
\end{itemize}
\end{itemize}


Here $\mathbf{v}^0$ is the vector representation of the sampled data $v^0$,
$\mathbf{r}_{k-1}$ denotes the misfit defined on the sample points $\Omega^0$,
$\phi_k$ is an intermediate surface, a solution of an interior point value
problem of the Laplace equation,
$\mathbf{n}$ denotes the unit outward normal defined on the image boundary
$\partial\Omega$, and $\|\cdot\|$ is either $L^2$ norm or the maximum norm.
The equation $\nabla\phi_k\cdot\mathbf{n}=0$ is called the no-flux boundary
condition.
The interior point value problem in
\eqref{eqn_R-CIM_PC}(iii) is incorporated for the construction of an
intermediate surface, due to its simplicity.
See \cite{CordellHakranKimETAL:Arbitrarily} for effective
computational methods including finite difference schemes,
a smoothing strategy, and algebraic solvers.

In this article, we apply \eqref{eqn_R-CIM_PC} for digital elevation
modeling for LiDAR point cloud data; we call the algorithm the
\emph{CIM with iterative refinement} (IR-CIM).


\section{Digital Elevation Modeling} \label{DEM}

\subsection{LiDAR data acquisition} \label{ss_LiDAR_Data}


For the previous decade or so, the \emph{light detection and ranging} (LiDAR)
technology has grown in popularity, meeting the need and replacing
conventional surveying techniques which are time-consuming and
labor-intensive
\cite{ChoiETAL:08,HosoiETAL:10,McKinionETAL:10,MitasovaETAL:05}.
LiDAR data are acquired in a form of point cloud; helicopter or
fixed-wing LiDAR systems scan the surface below the aircraft,
collecting reflected light signals in a scan rate of 50,000 to 100,000
pulses per second, achieving high accuracies up to 5cm, for most cases.
Due to their high resolution and rich information content,
LiDAR data are utilized for a wide range of applications with different
requirements in terms of resolution, accuracy, and surface representation.
Applications include digital elevation model (DEM) topography, flood
risk mapping, watershed analysis, coastal erosion analysis, landslides,
tree canopy analysis, transmission line mapping, and urban applications
\cite{Cary:09}.
LiDAR is an active remote sensing technique, analogous to radar, but
using laser light.

\begin{figure}[htb]
\begin{center}
 \includegraphics[width=0.27\textwidth]{fig1a} % lidar_acquisition.png}
 \includegraphics[width=0.32\textwidth]{fig1b} % Aircraft_Trajectory_May_12.png
 \includegraphics[width=0.38\textwidth]{fig1c} \\ %Aircraft_Scan_Coverage.png
 (a) \hfil (b)\hfil (c)\
\end{center}
\caption{LiDAR data acquisition:
(a) a schematic illustration of data collection,
(b) the aircraft trajectory for a field survey over Mississippi farms near
 Mississippi State University, May 12, 2011, and
(c) one of the LiDAR scan coverages.}
\label{fig_LiDAR_Acquisition}
\end{figure}

See Figure~\ref{fig_LiDAR_Acquisition}, which depicts a schematic
illustration of LiDAR data acquisition, along with the aircraft trajectory
and the LiDAR scan coverage for a field survey over Mississippi farms
near Mississippi State University (MSU), May 12, 2011.
LiDAR instruments built on an aircraft measure the round-trip time for a
pulse of laser energy to travel between the sensor and a target.
This incident pulse of energy
(usually with a near-infrared wavelength for vegetation studies)
reflects off of canopy (branches, leaves) and ground
surfaces and back to the instrument where it is collected by a telescope.
The travel time of the pulse, from initiation until it returns to the
sensor, provides a distance or range from the instrument to the object.
The acquired information is then transformed, with the aid of a minimum
of four GPS satellites, to obtain a 3D position fix.
The individual data points are collected to form a set of point cloud data.
As one can see from Figure~\ref{fig_LiDAR_Acquisition}(c),
the data are collected with the scan strips overlapped,
in order to densely cover the scan area by cloud points.


\begin{figure}[htb]
\begin{center}
 \includegraphics[width=0.49\textwidth]{fig2a} % Carr_Area_Data_Gap_m.png}
 \includegraphics[width=0.49\textwidth]{fig2b} \\ %Carr_Area_Elevation_Surface_m.png}
 (a) \hfil  (b)
\end{center}
\caption{LiDAR soil survey over Mississippi farms near MSU:
(a) a missing gap in the point cloud data and
(b) the elevation surface processed by the IDW built in ArcMap.}
\label{fig_DATA_Gap_DEM}
\end{figure}


\subsection{Moir\'e patterns} \label{ss_Moire_Patterns}

DEMs derived from LiDAR techniques are growing in popularity as a tool
for use in soil survey, in particular.
This form of remotely sensed elevation data can serve multiple purposes,
not the least of which is the visual interpretation of landforms and
soil parent materials.

However, as one can see from Figure~\ref{fig_DATA_Gap_DEM}(a), the
data may involve missing gaps; for which it is necessary to either
perform extra rounds of data acquisition or incorporate a very effective
interpolation algorithm for the reconstruction of reliable elevation
surfaces.
Figure~\ref{fig_DATA_Gap_DEM}(b) depicts a elevation surface processed
by the inverse-distance weighting (IDW) method built in ArcMap.
(ArcMap is one of state-of-the-art geospatial processing programs.)
In the figure, white regions are related to trees and buildings and the
ground elevation decreases in the SE direction.
The parallel features running in the SSW-NNE direction are processing
artifacts involved during the surface reconstruction of the data which
have been acquired through scan lines in the same direction with
multiple overlapped scan strips.
Thus the parallel features are kinds of Moir\'e interference patterns.
Contours and topographic parameters (gradients, curvatures) derived
from such elevation surfaces \emph{must} include a noisy pattern; for most
applications, they have been further processed (filtered, smoothed),
often manually, to make them suitable \cite{MitasovaETAL:05}.


\subsection{Correction strategy for Moir\'e effect} \label{ss_Correction}

As aforementioned, it is a common practice in LiDAR data acquisition
that the data are collected with the scan strips overlapped.
However, due to various technical reasons, data sets obtained from
different scans covering the same overlapped area may have misfits
in elevation values.
When these data sets are used without an appropriate correction, the
resulting surface may involve serious artifacts including oscillatory
patterns.

\begin{figure}[htb]
\begin{center}
 \includegraphics[width=0.47\textwidth]{fig3} %overlap_scans3.png}
\end{center}
\caption{A schematic illustration of LiDAR scan coverages.
The dash-dot lines (in blue) indicate the center of the scan strips
and the red bullets in the overlapped areas are check points.}
\label{fig_LiDAR_Schematic}
\end{figure}

The elevation misfits can be eliminated using local elevation averages
that are obtained from each of the data sets and measured over the
overlapped scan areas.
For a simple presentation, we first assume that the scan coverages
overlap maximum twice as shown in Figure~\ref{fig_LiDAR_Schematic}.
In the figure, $S_i$, $i=1,2,3$, denote the scan strips, the dash-dot
lines (in blue) indicate the center of the scan strips, and $C_i$
are points on the center lines.
The check points $P_{12}$ and $P_{23}$ represent centroids of overlapped
areas between the scans $S_1$ and $S_2$ and between $S_2$ and $S_3$,
respectively.
The data correction begins with partitioning the overlapped areas
and setting check points, the two rows of red points shown in
Figure~\ref{fig_LiDAR_Schematic}.
Then the average elevation values on each of overlapped scan strips are
computed in a vicinity of the check points, using the data values at
pixels that are scanned and assigned from both sides of the adjacent
scan strips.
If the average elevations obtained from the two different scan strips
are different at the check points, the data on the two adjacent strips
can be adjusted \emph{piecewise linearly} so that the resulting data sets
have the same average elevations at the check points.

\begin{figure}[htb]
\begin{center}
 \includegraphics[width=0.47\textwidth]{fig4} % overlap_scans3_Cor.png}
\end{center}
\vskip-1mm
\caption{A schematic illustration for data correction through local
elevation averages.}
\label{fig_LiDAR_Dat_Cor}
\end{figure}

In Figure~\ref{fig_LiDAR_Dat_Cor}, we assume that the average elevation
at $P_{12}$ computed from $S_1$ ($E_{12}$) is larger than that from $S_2$
($E_{21}$).
Then the data elevation values in $S_1$ and $S_2$ can be corrected
by adjusting the data values on half-strips linearly.
For example, assuming that $\overrightarrow{C_1C_2}$ is parallel
to the $y$-axis, the elevation values can be adjusted as follows.
\begin{equation}\label{eqn_S2_Correction}
\begin{gathered}
 (x,y,z)\in S_1 \mapsto (x,y,z'), \\[1mm]
 (x,y,z)\in S_2 \mapsto (x,y,z''),
\end{gathered}
\end{equation}
where
\begin{gather*}
 z' = z+\frac{(E_{21}-E_{12})/2}{p_{12,y}-c_{1,y}}\,(y-c_{1,y}),\\
 z'' = z+\frac{(E_{21}-E_{12})/2}{c_{2,y}-p_{12,y}}\,(c_{2,y}-y).
\end{gather*}
Here $c_{1,y},\,c_{2,y},\,p_{12,y}$ are the second coordinates
of $C_1,\,C_2,\,P_{12}$, respectively.
Note that when $y=p_{12,y}$,
$$
 z'=z+(E_{21}-E_{12})/2, \quad z''=z+(E_{21}-E_{12})/2;
$$
the corrections have the same magnitude from the both sides.

The algorithm for the elimination of Moir\'e patterns should also
incorporate the following concerns.


\noindent$\bullet$  \emph{Array of check points}:
The misfits, the differences of the average elevations on overlapped areas,
may differ when measured from different regions.
Thus we have introduced a line of check points in each of
overlapped areas, aligned parallel to the $x$ direction, as shown in
Figures~\ref{fig_LiDAR_Schematic}.
The Moire effects can be reduced effectively by applying piecewise bilinear
functions, as explained earlier in this subsection.

\noindent$\bullet$ \emph{Multiple overlaps}:
The scan coverages may overlap more than twice for some regions;
in practice, it is often the case that some regions are covered three times.
In this case, the array of check points should include more rows in the cross
direction of the scan flight (the $y$ direction),
in order to represent the misfits more appropriately.
For example, let the data overlap $m$ times in the scan strip $S_k$.
Then the misfit correction function (MCF) for the data in $S_k$ can be
defined to be a polynomial of degree at most $m-1$ in the $y$ direction.
The MCF is still piecewise linear in the $x$ direction.

\noindent$\bullet$ \emph{Accuracy of MCFs}:
It is occasionally the case that data loci are distributed with a largely
varying density; the data may involve missing gaps.
Thus, local elevation averages obtained from raw data may not accurately
represent the real elevation averages for some regions and therefore
the MCFs can be erroneous.
In order to overcome the difficulty, we may consider the following strategy.
\begin{itemize}
\item[(i)]
Construct the surfaces for each of the scan strips
using the IR-CIM \eqref{eqn_R-CIM_PC}.

\item[(ii)]
Compute the local elevation averages using the elevation values from
the reconstructed surfaces.

\item[(iii)]
Construct MCFs for each of overlapped strips.

\item[(iv)]
Apply the MCFs to the reconstructed surfaces.

\item[(v)]
Blend the corrected surfaces for a single surface covering the whole scan area.

\item[(vi)]
Improve the blended surface by applying the IR-CIM.
\end{itemize}
The final step is not expensive computationally, because
the blended surface would be a good approximation of the final result.

\section{Numerical experiments} \label{Numerical}

This section gives numerical examples to show effectiveness of the IR-CIM
and the Moir\'e-pattern reduction algorithm applied for DEM modeling.
The 12 strips of LiDAR point cloud data acquired over
Mississippi farms near Mississippi State University, as shown in
Figure~\ref{fig_LiDAR_Acquisition}(c), are utilized for the modeling.
The data set includes approximately 37 million points counted including
multiple arrivals.


\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.96\textwidth]{fig5} % Surface_LiDAR_DATA_All_s.png}
\end{center}
\caption{Final image surface in 96 million pixels,
covering a region of 3km$\times$2km square
with a quarter-meter resolution.}
\label{fig_LiDAR_DATA_All}
\end{figure}

Figure~\ref{fig_LiDAR_DATA_All} depicts the final elevation image in 96
million pixels, which covers a region of 3Km$\times$2Km square with a
quarter-meter resolution.
In this study, the Moir\'e effect is reduced using
the elevation averages measured from the data, applying the formulas
in \eqref{eqn_S2_Correction}.
About 17\% of the pixel values of the image can be assigned directly
from the elevation values in the data and other pixel values are computed
using the IR-CIM \eqref{eqn_R-CIM_PC}.
The algorithm (written in a combination of C and C++) converges in 3
iterations, taking 67.2 seconds on a desktop computer for the whole
processing.

\begin{figure}[htb]
\begin{center}
 \includegraphics[width=0.47\textwidth]{fig6a} % Moire_effect.png
 \includegraphics[width=0.47\textwidth]{fig6b} \\ %Moire_effect_fixed.png
 (a)\hfil (b) 
\end{center}
\caption{Moir\'e effect and its correction.}
\label{fig_Moire_Effect}
\end{figure}


Figure~\ref{fig_Moire_Effect} shows image surfaces, shown in a prism view
of 1cm resolution in elevation, covering the region marked with a square
(in red) in Figure~\ref{fig_LiDAR_DATA_All}.
When the data values are used without correction of Moir\'e effect,
the reconstructed surface involves oscillatory patterns as in
Figure~\ref{fig_Moire_Effect}(a).
On the other hand, our Moir\'e-effect reduction algorithm
has eliminated the oscillatory patterns effectively as depicted
in Figure~\ref{fig_Moire_Effect}(b).
In the last image, Moir\'e effect is hardly observable even shown
in 1cm resolution in elevation.


\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.47\textwidth]{fig7} % moire_correction_2d.png
\end{center}
\caption{Correction of elevation values, shown on a vertical line segment.
The solid curve (in black) represents the corrected elevation values
over an overlapped region of two scan strips.}
\label{fig_moire_correction_2d}
\end{figure}

To further investigate effectiveness of the suggested algorithm,
we select a vertical segment in the mid of the square-marked region in
Figure~\ref{fig_LiDAR_DATA_All} and compare elevation values.
In Figure~\ref{fig_moire_correction_2d}, the dotted and dashed curves
indicate elevation values obtained from separate scan strips of an
overlapped region, while the solid curve (in black) represents the
corrected elevation values over the overlapped region of the two scan strips.
As one can see from the figure, the elevation values are corrected
using a locally linear correction function effectively, to result in
a continuous and reliable image surface.

\subsection*{Conclusions} %\label{Conclusions}

Surface reconstruction from point cloud data is a challenging problem
particularly when no constraint is imposed on their locations, as for
LiDAR data.
This article has applied an effective PDE-based algorithm, called
the curvature interpolation method (CIM), for a set of LiDAR data
acquired from a field survey over Mississippi farms near
Mississippi State University.
For the reduction of oscillatory patterns appearing over overlapped
regions of scan strips, we have introduced an effective misfit correction
function (MCF) which is piecewise linear.
A recursive application of the CIM with the suggested MCF has converged
in 3-4 iterations to produce reliable and piecewise smooth image
surfaces that introduce no Moir\'e effect, taking less than a second per
million pixels on common desktop computers.

\subsection*{Acknowledgments}

 H. Kim and S. Kim were supported in part by NSF grant
DMS-1228337.


\begin{thebibliography}{00}

\bibitem{ArigovindanETAL:05}
M.~Arigovindan, M.~{S\"uhling}, P.~Hunziker,  M.~Unser;
 \emph{ Variational image reconstruction from arbitrarily spaced samples: {A} fast
 multiresolution spline solution}, IEEE Trans. Image Process., 14 (2005),
 pp.~450--460.

\bibitem{BeatsonETAL:00}
R.~Beatson, W.~A. Light,  S.~Billings; \emph{Fast solutions of radial
 basis function interpolation equations: {Domain} decomposistion methods},
 SIAM J. Sci. Comput., 22 (2000), pp.~1717--1740.

\bibitem{BernardiniETAL:99}
F.~Bernardini, J.~Mittleman, H.~Rushmeier, C.~Silva, G.~Taubin;
 \emph{ The ball-pivoting algorithm for surface reconstruction}, IEEE Transactions on
 Visualization and Computer Graphics, 5 (1999), pp.~349--359.

\bibitem{Cary:09}
T.~Cary; \emph{New research reveals current and future trends in {LiDAR}
 applications}, Industry Insights, (2009), pp.~8--9.

\bibitem{ChaLeeKim:14_RCIM}
Y.~Cha, G.~Y. Lee,  S.~Kim; \emph{Image zooming by curvature
 interpolation and iterative refinement}, SIAM J. Imaging Sciences, 7 (2014),
 pp.~1284--1308.

\bibitem{ChanShen:05}
T.~Chan, J.~Shen; \emph{Image Processing and Analysis}, SIAM,
 Philadelphia, 2005.

\bibitem{ChenSuter:98}
F.~Chen, D.~Suter; \emph{Using fast multipole method for accelerating
 the evaluation of splines}, IEEE Comput. Sci. Eng., 5 (1998), pp.~24--31.

\bibitem{ChoiETAL:08}
Y.-W. Choi, Y.-W. Jang, H.-J. Lee,  G.-S. Cho;
 \emph{Three-dimensional
 {LiDAR} data classifying to extract road point in urban area}, IEEE
 Geoscience and Remote Sensing Letters, 5 (2008), pp.~725--729.

\bibitem{CordellHakranKimETAL:Arbitrarily}
W.~Cordell, H.~Kim, J.~Willers,  S.~Kim;
 \emph{Image reconstruction for
 arbitrarily spaced data using curvature interpolation}, in Proceedings of The
 2013 International Conference on Image Processing, Computer Vision, and
 Pattern Recognition, 2013, pp.~506--512.

\bibitem{FaulPowel:99}
A.~Faul and M.~Powel; 
\emph{Proof of convergence of an iterative technique
 for thin plate spline interpolation in two dimensions}, Adv. Comput. Math.,
 11 (1999), pp.~183--192.

\bibitem{Glassner:95}
A.~Glassner; \emph{Principles of Digital Image Synthesis}, Morgan
 Kaufmann, San Mateo, CA, 1995.

\bibitem{HiltonETAL:96}
A.~Hilton, A.~J. Stoddart, J.~Illingworth, T.~Wandeatt;
 \emph{Marching triangles: {Range} image fusion for complex object modeling}, 
in Proceedings of International Conference on Image Processing, vol.~1, 1996, pp.~381--384.

\bibitem{HoppeETAL:92}
H.~Hoppe, T.~DeRose, T.~Duchamp, J.~McDonald, W.~Stuetzle;
 \emph{ Surface reconstruction from unorganized points}, in SIGGRAPH'92: Proceedings
 of the 19th Annual Conference on Computer Graphics and Interactive
 Techniques, New York, NY, USA, 1992, ACM, pp.~71--78.

\bibitem{HosoiETAL:10}
F.~Hosoi, Y.~Nakai, K.~Omasa;
 \emph{Estimation and error analysis of  woody canopy leaf area density profiles 
using {3-D} airborne and ground-based scanning {Lidar} remote-sensing techniques}, 
IEEE Trans. on Geoscience and  Remote Sensing, 48 (2010), pp.~2215--2223.

\bibitem{KimChaKim:11_CIM}
H.~Kim, Y.~Cha, S.~Kim;
\emph{Curvature interpolation method for image
 zooming}, IEEE Trans. Image Process., 20 (2011), pp.~1895--1903.

\bibitem{LehmannGonnerSpitzer:99}
T.~Lehmann, C.~{G\"onner}, K.~Spitzer;
 \emph{Survey: {Interpolation}
 methods in medical image processing}, IEEE Trans. Medical Imaging, 18 (1999),
 pp.~1049--1075.

\bibitem{MarquinaOsher:00}
A.~Marquina, S.~Osher;
 \emph{Explicit algorithms for a new time
 dependent model based on level set motion for nonlinear deblurring and noise
 removal}, SIAM J. Sci. Comput., 22 (2000), pp.~387--405.

\bibitem{McKinionETAL:10}
J.~McKinion, J.~Willers, J.~Jenkins;
 \emph{Comparing high density  {LIDAR} and medium resolution {GPS} generated 
elevation data for predicting  yield stability}, Computers and Electronics in Agriculture, 74 (2010),
 pp.~244--249.

\bibitem{MitasovaETAL:05}
H.~Mitasova, L.~Mitas, R.~S. Harmon;
 \emph{Simultaneous spline approximation and topographic analysis for {Lidar} 
elevation data in  open-source {GIS}}, IEEE Geoscience and Remote Sensing Letters, 
2 (2005),  pp.~375--379.

\bibitem{RudinOsherFatemi:92}
L.~Rudin, S.~Osher, E.~Fatemi;
 \emph{Nonlinear total variation based
 noise removal algorithms}, Physica D, 60 (1992), pp.~259--268.

\bibitem{Shepard:68}
D.~Shepard;
 \emph{A two-dimensional interpolation function for
 irregularly-spaced data}, in Proceedings of the 1968 ACM National Conference,
 1968, pp.~517--524.

\bibitem{TurkOBrien:99}
G.~Turk J. O\'Brien;
 \emph{Shape transformation using variational
 implicit functions}, in Proc. ACM SIGGRAPH, Los Angeles, CA, 1999,
 pp.~335--342.

\bibitem{Turkowski:90}
K.~Turkowski;
 \emph{Filters for common resampling tasks}, in Graphics Gems
 I, A.~S. Glassner, ed., Academic Press Professional, Inc., San Diego, CA,
 USA, 1990, pp.~147--165.

\bibitem{Watson:92}
D.~Watson;
 \emph{Contouring: A Guide to the Analysis and Display of
 Spatial Data}, Pergamon, New York, 1992.

\end{thebibliography}

\end{document}