\documentclass[letterpaper,12pt]{article} \usepackage{geometry} \usepackage{float} %\geometry{verbose,letterpaper,tmargin=2cm,bmargin=2cm,lmargin=2cm,rmargin=2cm} \usepackage{color} \usepackage{xcolor} \usepackage[utf8]{inputenc} \usepackage{graphicx} \usepackage{subfigure} \usepackage{CJKutf8} \usepackage{amsfonts,amssymb} \usepackage{amsmath} \usepackage{amsthm} \usepackage[T1]{fontenc} \usepackage[utf8]{inputenc} \usepackage{lmodern} \usepackage{tikz} \usepackage{setspace} % math theorem \newtheorem{fact}{Fact} \newtheorem{lemma}{Lemma} \newtheorem{theorem}[lemma]{Theorem} \newtheorem{definition}[lemma]{Definition} \newtheorem{assumption}[lemma]{Assumption} \newtheorem{corollary}[lemma]{Corollary} \newtheorem{proposition}[lemma]{Proposition} \newtheorem{exercise}[lemma]{Exercise} \newtheorem{claim}[lemma]{Claim} \newtheorem{remark}[lemma]{Remark} \newtheorem{prob}{Problem} \newtheorem{conjecture}{Conjecture} % computer science % \usepackage[linesnumbered,boxed,ruled,vlined,longend]{algorithm2e} \usepackage[linesnumbered,boxed,ruled,longend]{algorithm2e} % cite and link \usepackage{cite} \usepackage{url} \usepackage[final]{hyperref} % adds hyper links inside the generated pdf file \hypersetup{ colorlinks=true, % false: boxed links; true: colored links linkcolor=blue, % color of internal links citecolor=blue, % color of links to bibliography filecolor=magenta, % color of filelinks urlcolor=blue } \title{Ruppert's Delaunay Refinement Algorithm} \author{Tong, Shen. \footnote{Email: 385989829@qq.com} Wei, Cao.\footnote{Email: fatboy\_cw@163.com } Yifei, Jin \footnote{Email: bluewould@yeah.net} } \begin{document} \maketitle \begin{abstract} Delaunay triangulation maximizes the smallest angle among all possible triangulations of a given input and hence is a powerful discretization tool. However, Delaunay triangulation can have arbitrarily small angles depending on the input configuration. Thus, Delaunay refinement algorithms which iteratively insert additional points were developed to remedy the problem. Delaunay refinement method is arguably the most popular due to its theoretical guarantee and performance in practice. The first step of a Delaunay refinement algorithm is the construction of a constrained or conforming Delaunay triangulation of the input domain. This initial Delaunay triangulation is likely to have bad elements. Delaunay refinement algorithm then iteratively adds new points to the domain to improve the quality of the mesh and to ensure that the mesh conforms to the boundary of the input domain. The points inserted by the Delaunay refinement are \emph{Steiner points}. A sequential Delaunay refinement algorithm typically adds one new vertex at each iteration. Each new vertex is chosen from a set of candidates -- the circumcenters of bad triangles and the mid-points of input segments. Chew \cite{Chew} showed that Delaunay refinement can be used to produce quality-guaranteed in two dimensions. Ruppert extended the technique for computing not only quality-guaranteed but also size-optimal triangulations. In our experiment report, we implement the Ruppert algorithm. And give a friendly interactive application. \end{abstract} \section{Algorithm} \label{sec:01} \subsection{Background} \label{sec:01:01} In two dimensions, the input domain $\Omega$ is usually represented as \emph{planar straight line graph} (PSLG) -- proper planar drawing in which each edge is mapped to a straight line segment between its two endpoint \cite{Ruppert}. The segments express the boundaries of $\Omega$ and the endpoints are the vertices of $\Omega$. \begin{definition}[diametral circle] The diametral circle of a segment is the circle whose diameter is the segment. Equivalently, the diametral circle of a segment is also the smallest circle that encloses the segment. \end{definition} \begin{definition}[encroach] A segment is said to be encroached if a vertex lies strictly inside its diametral circle, or if the segment does not appear as an edge of the triangulation. \end{definition} \begin{definition}[Constrained Delaunay triangulation] constrained Delaunay triangulation is a generalization of the Delaunay triangulation that forces certain required segments into the triangulation. \footnote{Note: Because a Delaunay triangulation is almost always unique, often a constrained Delaunay triangulation contains edges that do not satisfy the Delaunay condition. Thus a constrained Delaunay triangulation often is not a Delaunay triangulation itself.} \end{definition} \begin{definition}[Radius-edge Ratio] Radius-edge ration of a triangle is the ratio of its circumradius to the length of its shortest side. \end{definition} \begin{definition}{skinny triangle} Skinny trinangle, alternately named bad triangle, is a triangle whose radius-edge ratio is larger than a specified constant $\beta$.\footnote{in this experiment report, we define $\beta=\sqrt{2}$} \end{definition} Note, in the refinement triangulation, we demand all triangles' angle larger than some constant $\alpha$. It is easy to prove, in two dimensions, an upper bound $\beta$ on the radius-edge ration implies a lower bound $\alpha = \arcsin(1/2\beta)$ on the smallest angle and vice versa. \subsection{Ruppert's Algorithm} \label{sec:01:02} The pseudocode Ruppert's algorithm is listed in \ref{algo:01}. As an example , the figure \ref{fig:01} give a explanation. The algorithm input an PSLG graph. At first step, we do a Delaunay triangulation of the input vertices. But there are some segments of subsegments are not strongly Delaunay, just like the figure \ref{fig:subfig:c}, which means these segments do not satisfies \emph{empty circle property}\footnote{Two sites $p_i$ and $p_j$ are connected by an edge in the Delaunay triangulation, if and only if there is an empty circle passing through $p_i$ and $p_j$ .}. And some vertices lie strictly inside a diametral circle of some segments, just like figure \ref{fig:subfig:d}. In order to include the segments, second step, we add its mid-points. Because there are some segment locates outside the interval region, at third step, we just delete them. Then at fourth step, we process the skinny triangle just like the figure \ref{fig:subfig:e}. we add the circumcenter. If the circumcenter is not cause any encroach, we accept it and re-triangulate the graph. If not we delete it and add the encroach segment. \footnote{Here note that: the circumcenter is deleted, in fact the segment is just encroach a virtual point. The reason is guarantee the following circumcenter locate in the internal region which we prove later}. Through iterations, we could get final refinement Delaunay triangulation. \begin{figure}[!htb] \centering \subfigure[A sample input PSLG]{ \label{fig:subfig:a} \includegraphics[width=0.25\textwidth]{fig01.png}} \subfigure[Delaunay triangulation of the input vertices, note that some input segment are missing]{ \label{fig:subfig:b} \includegraphics[width=0.25\textwidth]{fig02.png}} \subfigure[Vertex insertion restores the missing segment]{ \label{fig:subfig:c} \includegraphics[width=0.25\textwidth]{fig03.png}} \subfigure[Encroached subsegment is split]{ \label{fig:subfig:d} \includegraphics[width=0.25\textwidth]{fig04.png}} \subfigure[Delete the outer segments]{ \label{fig:subfig:e} \includegraphics[width=0.25\textwidth]{fig05.png}} \subfigure[Shinny triangle's circumcenter is inserted]{ \label{fig:subfig:f} \includegraphics[width=0.25\textwidth]{fig06.png}} \subfigure[This circumcenter encroaches upon a segment, and is rejected]{ \label{fig:subfig:g} \includegraphics[width=0.25\textwidth]{fig07.png}} \subfigure[Although the circumcenter was rejected, the segment it encroached upon is still marked for bisection]{ \label{fig:subfig:h} \includegraphics[width=0.25\textwidth]{fig08.png}} \subfigure[Final mesh]{ \label{fig:subfig:j} \includegraphics[width=0.25\textwidth]{fig09.png}} \caption{The flow of Ruppert's Algorithm} \label{fig:01} \end{figure} \begin{algorithm}[!] \caption{Ruppert's refinement triangulation algorithm} \label{algo:01} \SetKwFunction{DelaunayTriangulation}{DelaunayTriangulation} \SetKwFunction{AddMidOfEncroach}{AddMidOfEncroach} \SetKwFunction{DelExternalSegment}{DelExternalSegment} \SetKwInOut{Input}{input} \SetKwInOut{Output}{output} \Input {points, segment, threshold.} \Output{refinement triangulation algorithm.} $T$: = \DelaunayTriangulation(points) \; $Q$: = the set of encroached segments and skinny triangles \; $T, Q$: = \AddMidOfEncroach(T,Q) \; $T$: = \DelExternalSegment(T) \; \While{$Q$ is not empty:}{ \eIf{$Q$ contains a segment s:} {insert the midpoint of s into $T$ \;} ( $Q$ contains poor quality triangle $t$:) { \eIf{the circumcenter of $t$ encroaches a segments $s$: } {add $s$ to $Q$ \;} {insert the circumcenter of $t$ into $T$ \;} } update $Q$ \; } return $T$ \; \end{algorithm} \subsection{Algorithm Analysis} \label{sec:01:03} \subsubsection{Circumcenter lies in boundaries} \label{sec:01:03:01} The algorithm is not difficult, but their some question remained to answer. First is whether $v$ always locate in boundaries or not. If not, the algorithm won't convergence. Fortunately, we can prove all circumcenter of triangles lies in boundaries. \begin{theorem} Let $T$ be a segment-bounded Delaunay triangulation (hence, any edge of $T$ that belongs to only one triangle is a subsegment). Suppose that $T $ has no encroached subsegments. Let $v$ be the circumcenter of some triangle $t$ of $T$. Then $v$ lies in $T$.\cite{Richard} \end{theorem} \begin{proof} Suppose for the sake of contradiction that $v$ lies outside $T$. Let $c$ be the centroid of $t$. $c$ clearly lies inside $T$. Because the triangulation is segment-bounded, the line segment $cv$ must cross some subsegment $s$, as Figure \ref{fig:thm1} illustrates. Because $cv$ is entirely contained in the interior of the circumcircle of $t$, the circumcircle must contain a portion of $s$ ; but the Delaunay property requires that the circumcircle be empty, so the circumcircle cannot contain the endpoints of $s$. Say that a point is inside $s$ if it is on the same side of $s$ as $c$, and outside $s$ if it is on the same side of $s$ as $v$ . Because the center $v$ of the circumcircle of $t$ is outsides, the portion of the circumcircle that lies strictly insides (the bold arc in the illustration) is entirely enclosed by the diametral circle of $s$. The vertices of $t$ lie upont’s circumcircle and are (not strictly) insides. Up to two of the vertices oftmay be the endpoints of $s$, but at least one vertex oftmust lie strictly inside the diametral circle of $s$. But $T$ has no encroached subsegments by assumption; the result follows by contradiction. \end{proof} \begin{figure}[!htb] \centering \includegraphics[width=0.6\textwidth]{figa.png} \caption{: If the circumcenter v of a triangle t lies outside the triangulation, then some subsegment s is encroached.} \label{fig:thm1} \end{figure} \subsubsection{Convergence} \label{sec:01:03:02} The claim that Ruppert's algorithm produces nicely graded meshes is based on the fact that the spacing of vertices at any loxation in the mesh is within a constant factor of the sparsest possible spacing. To formalize the idea of "sparsest possible spacing'', Ruppert introduces a function called the \emph{local feature size} \begin{definition}[local feature size] local feature size of each point $x \in \mathcal{R}^2 $ , denoted by $lfs_{\Omega}(x)$, is the radius of the smallest disk centered at $x$ that touch two different features which is the vertices and boundary segments of $\Omega$. \end{definition} Figure \ref{fig:b} illustrates the notion by giving examples of such disks for a variety of points. \begin{figure}[!htb] \centering \includegraphics[width=0.8\textwidth]{figb.png} \caption{the radius of each disk illustrated is the local feature size of the point at its center} \label{fig:b} \end{figure} It is not difficult to prove the lemma below: \begin{theorem} For any PSLG X, and any two points $u$ and $v$ in the plane, \begin{equation*} lfs(v) \leq lfs(u) +|uv| \end{equation*} \end{theorem} \begin{proof} Since $lfs(u)$ represent a disk centered at $u$ and intersects two non-incident features of $X$. So the disk whose center is at $v$ and radius is $lsf(u) + |uv|$ also intersect at least two features since it contains disk $lfs(u)$. Hence $ lfs(v) \leq lfs(u) +|uv|$. \end{proof} Ruppert prove that radius-edge ratio larger than $\sqrt{2}$ and any two incident segments in the input PSLG are separated by an angle of $60^{\circ}$ or greater, that the algorithm produces meshes that are nicely graded and size-optimal. \begin{theorem} \label{thm terminate} Let $lfs_{\min}$ be the shortest distance between two non-incident entities (vertices or segments) of the input PSLG suppose that any two incident segments are separated by an angle of at least $60^{\circ}$, and a triangle is considered to be skinny if its circumradius -to-shortest edge ratio is larger than $B$, where $B \geq \sqrt{2}$ Ruppert's algorithm will terminate, with no triangulation edge shorter than $lfs_{min}$, where $lfs_{min}$ is the $\min_{u} {lfs(u)}$, where $u$ is chosen from among the input vertex. \end{theorem} The proof of the theorem is not difficult but tedious. You can refer reference \cite{Richard}. Since no edge shorter than $lfs_{min}$, and the region is limited. Therefore the algorithm must terminate. \subsubsection{Size-Optimality} \label{sec:01:03:03} Theorem \ref{thm terminate} guarantees that no edge of the final mesh is smaller than $lfs_{\min}$, but it is not satisfying spatially graded mesh. What following is a proof that each edge of the output mesh has length proportional to the local feature sizes of its endpoints. Here omit the proof detail, just give the result. \begin{theorem} For any vertex $v$ of the output mesh, the distance to its nearest neighbor $w$ is at least $\frac{lfs(v)}{D_s+1}$, where $D_s$ is a constant larger than 1. \end{theorem} \subsection{Implement trick} \label{sec:01:04} The algorithm has several implement trick, we list some important below: \begin{itemize} \item At first we count the encroach segment. The segment is not strongly Delaunay , i.e. the segment is not used by Delaunay triangulation. There must be a vertex on or inside its diametral circle. This observation is important because it unifies the theoretical treatment of missing subsegments and encroached subsegments that are not missing. \item A subsegment may be tested for encroachment by inspecting only those vertices that appear directly opposite the subsegment in a triangle. Consider Figure \ref{fig:c}. Both of the vertices (v and w ) opposite the segment $s$ lie outside the diametral circle of $s$. Because the mesh is constrained Delaunay, each triangle's circumcircle is empty (on its side of $s$), and therefore the diametral circle of s is empty. The Delaunay refinement process may cause other subsegments to become encroached. According to the observation above, we don not check all segment each time but just test each of the edges that appear opposite the vertex in some triangle is enough. \item Each time we insert a new circumcenter, we need know the point locate in which triangle. So we need a point a location algorithm to remain the information. But in this algorithm, we already know the Delaunay triangulation, so there exist a simple point location algorithm. Just long the segment between circumcenter and orthocenter. Search all the triangle along the segment, we could find the circumcenter location. \begin{figure}[!htb] \centering \includegraphics[width= 0.9\textwidth]{figc.png} \caption{if the apices(here, v and w) of the triangles that contain a subsegment $s$ are outside the diametral circle of $s$, then no vertex lies in the diametral circle of $s$, because the triangles are Delaunay} \label{fig:c} \end{figure} \end{itemize} \section{Code Structure} \label{sec:02} The code of the this project is written by \emph{Python 2.7} based on \emph{PyOpenGL} and \emph{PyQt5}. Only the initial delaunay triangulation is implemented with package \emph{triangle} provided y CMU. The code structure of main algorithm part is associated in following way: \begin{itemize} \item cgalgo.py \begin{itemize} \item Sgn(x):Check the sign of a float number, the eps is setted by $1e12$. \item TurnLeft(v):return the $90$ degree left turn vector of input vector. \item Cross(a,b, c):return the multiplication cross of vector $\vec{ba}$ and $\vec{ca}$. \item GetDistance(u,v):return the Euclidean distance of point $u$ and $v$. \item GetIntersection(lu, lv):return the intersection of line $lu$ and $lv$. \item GetBisector(u, v):return bisector of line $uv$ \item GetCircleCenter(a, b, c):return circumcenter of triangle $abc$. \item GetSegIntersection(a, b, c, d):return the intersection of segment $ab$ and $cd$ \item InTriangle(p, a, b, c):check whether point $p$ is inside of triangle $abc$. \item InCircle(p, a, b, x):checht whether the point $x$ is inside of circumcircle of triangle $pab$. \end{itemize} \item Ruper.py \begin{itemize} \item InitializeDelaunay: return the initial delaunay triangulation implemented with third part package \emph{triangle}. \item TrianglePointLocation: find the circumcenter of a skinny triangle by flipping triangle sequentially. \item IsEncroached:check whether a segment is a encroached segment. \item IsSkinny:check whether a triangle is a skinny triangle. \item UpdateDelaunay: Do swap test and maintain some necessary information when the triangulation changed. \item RecoverTriangles: recover the triangulation if a circumcenter is failed to insert into triangulation. \item SplitSegment:split a segment by inserting a midpoint of corresponding segment. \item InsertCircleCenter:insert a circumcenter of a skinny triangle. \item EliminateSegment:Main loop of segments processing. \item EliminateAngle:Main loop of triangles processing. \item RemoveOutSide:Remove the triangles outside the planar graph. \end{itemize} \end{itemize} \section{Software instruction} \label{sec:03} This part you can see our homepage \url{http://caow13.github.io/Rupers_Algorithm/}. \section{Experiment} \label{sec:04} We list some comparison about origin Delaunay triangulation and Ruppert's refinement algorithm. Just like figures \ref{fig:cmp1} and \ref{fig:cmp2} \begin{figure}[!htb] \centering \subfigure[Holes- Delaunay]{ \label{fig:aa} \includegraphics[width=0.4\textwidth]{figaa.png}} \subfigure[Holes- Ruppert]{ \label{fig:ab} \includegraphics[width=0.4\textwidth]{figab.png}} \subfigure[cmu- Delaunay]{ \label{fig:ba} \includegraphics[width=0.4\textwidth]{figba.png}} \subfigure[cmu- Ruppert]{ \label{fig:bb} \includegraphics[width=0.4\textwidth]{figbb.png}} \subfigure[double\_hex- Delaunay]{ \label{fig:ca} \includegraphics[width=0.4\textwidth]{figca.png}} \subfigure[double\_hex- Ruppert]{ \label{fig:cb} \includegraphics[width=0.4\textwidth]{figcb.png}} \caption{experiment results(1)} \label{fig:cmp1} \end{figure} \begin{figure}[!htb] \centering \subfigure[face- Delaunay]{ \label{fig:da} \includegraphics[width=0.4\textwidth]{figda.png}} \subfigure[face- Ruppert]{ \label{fig:db} \includegraphics[width=0.4\textwidth]{figdb.png}} \subfigure[key- Delaunay]{ \label{fig:ea} \includegraphics[width=0.4\textwidth]{figea.png}} \subfigure[key- Ruppert]{ \label{fig:eb} \includegraphics[width=0.4\textwidth]{figeb.png}} \caption{experiment results(2)} \label{fig:cmp2} \end{figure} % todo \section{Conclusion and Future work} \label{sec:05} The Ruppert's algorithm is an fundamental work about refinement Delaunay triangulation. But the algorithm time complexity is difficult to analysis. We already know there are some algorithm can achieve $O(m+n\log(n))$ where $m$ is the extra points we add. The further work you can see reference \cite{Optimal} % bibliography \begin{thebibliography}{99} \bibitem {Optimal} Har-Peled, Sariel, and Alper \"{U}ng\"{o}r. "A time-optimal Delaunay refinement algorithm in two dimensions." Proceedings of the twenty-first annual symposium on Computational geometry. ACM, 2005. \bibitem {Chew} L.P.Chew. Constrained Delaunay triangular meshes. Technical Report TR-89-983, Dept. Comput. Sci. Cornell Univ., Ithaca, NY, Apr. 1898 \bibitem{Ruppert} J. Ruppert/ A new and simple algorithm for quality 2-dimensional mesh generation. In \emph{Proc. 4th ACM-SIAM Sympos. Discrete Algorithm, pages 83-92, 1993} \bibitem{Richard} Lecture Notes on Delaunay Mesh Generation, Jonathan Richard Shewchuk, September 20, 1999 \end{thebibliography} \end{document}