Merge branch 'master' of https://github.com/libigl/libigl

Former-commit-id: d33ef6b8fc88bcb54d9ebf70213bd710144f2c34

README.md

@@ -7,9 +7,7 @@
 > Get started with:
 >
 ```bash
-git clone https://github.com/libigl/libigl.git
-cd libigl
-git submodule update --init --recursive
+git clone --recursive https://github.com/libigl/libigl.git
 ```
 
 libigl is a simple C++ geometry processing library. We have a wide
@@ -137,6 +135,25 @@ We hope to fix this, or at least identify which functions are safe (many of
 them probably work just fine). This requires setting up unit testing, which is
 a major _todo_ for our development.
 
+## Git Submodules 
+Libigl uses git submodules for its _optional_ dependencies,
+in particular, those needed by the OpenGL viewer to run the examples in the
+[tutorial](tutorial/tutorial.html). Git submodules allow use to treat clones of
+other libraries as sub-directories within ours while separating our commits.
+Read the [documentation](http://git-scm.com/docs/git-submodule) for a detailed
+explanation, but essentially our libigl repo stores a hash for each of its
+subrepos containing which version to update to. When a change is introduced in
+a dependencies repo we can incorporate that change by pulling in our sub-repo
+and updating (i.e.  committing) that change to the hash.
+
+When pulling new changes to libigl it's also a good idea to update changes to
+subrepos:
+
+```bash
+git pull
+git submodule update -- recursive
+```
+
 ## How to contribute
 
 If you are interested in joining development, please fork the repository and
@@ -193,7 +210,8 @@ Libigl is a group endeavor led by [Alec
 Jacobson](http://www.cs.columbia.edu/~jacobson/) and [Daniele
 Panozzo](http://www.inf.ethz.ch/personal/dpanozzo/). Please [contact
 us](mailto:alecjacobson@gmail.com,daniele.panozzo@gmail.com) if you have
-questions or comments. We are happy to get feedback!
+questions or comments. For troubleshooting, please post an
+[issue](https://github.com/libigl/libigl/issues) on github.
 
 If you're using libigl in your projects, quickly [drop us a
 note](mailto:alecjacobson@gmail.com,daniele.panozzo@gmail.com). Tell us who you

documentation/active-set.tex

@@ -1,326 +0,0 @@
-\documentclass[12pt]{diary}
-\title{Active set solver for quadratic programming}
-\author{Alec Jacobson}
-\date{18 September 2013}
-
-\renewcommand{\A}{\mat{A}}
-\renewcommand{\Q}{\mat{Q}}
-\newcommand{\RR}{\mat{R}}
-\newcommand{\Aeq}{\mat{A}_\textrm{eq}}
-\newcommand{\Aieq}{\mat{A}_\textrm{ieq}}
-\newcommand{\beq}{\vc{b}_\textrm{eq}}
-\newcommand{\bieq}{\vc{b}_\textrm{ieq}}
-\newcommand{\lx}{\Bell}
-\newcommand{\ux}{\vc{u}}
-\newcommand{\lameq} {\lambda_\textrm{eq}}
-\newcommand{\lamieq}{\lambda_\textrm{ieq}}
-\newcommand{\lamlu}  {\lambda_\textrm{lu}}
-
-\begin{document}
-
-\begin{pullout}
-\footnotesize
-\emph{Disclaimer: This document rewrites and paraphrases the ideas in
-\url{http://www.math.uh.edu/~rohop/fall_06/Chapter3.pdf},
-\url{http://www.math.uh.edu/~rohop/fall_06/Chapter2.pdf}, and
-\url{http://www.cs.cornell.edu/courses/cs322/2007sp/notes/qr.pdf}. Mostly this
-is to put everything in the same place and use notation compatible with the
-\textsc{libigl} implementation. The ideas and descriptions, however, are not novel.}
-\end{pullout}
-
-Quadratic programming problems (QPs) can be written in general as:
-\begin{align}
-\argmin \limits_\z &
-  \z^\transpose \A \z + \z^\transpose \b + \text{ constant}\\
-\text{subject to } & \Aieq \z ≤ \bieq,
-\end{align}
-where $\z \in \R^n$ is a vector of unknowns,  $\A \in \R^{n \times n}$ is a (in
-our case sparse) matrix of quadratic coefficients, $\b \in \R^n$ is a vector of
-linear coefficients, $\Aieq \in \R^{m_\text{ieq} \times n}$ is a matrix (also
-sparse) linear inequality coefficients and $\bieq \in \R^{m_\text{ieq}}$ is a
-vector of corresponding right-hand sides. Each row in $\Aieq \z ≤ \bieq$
-corresponds to a single linear inequality constraint.
-
-Though representable by the linear inequality constraints above---linear
-\emph{equality} constraints, constant bounds, and constant fixed values appear
-so often that we can write a more practical form
-\begin{align}
-\argmin \limits_\z &
-  \z^\transpose \A \z + \z^\transpose \b + \text{ constant}\\
-\text{subject to } & \z_\text{known} = \y,\\
-                   & \Aeq \z = \beq,\\
-                   & \Aieq \z ≤ \bieq,\\
-                   & \z ≥ \lx,\\
-                   & \z ≤ \ux,
-\end{align}
-where $\z_\text{known} \in \R^{n_\text{known}}$ is a subvector of our unknowns $\z$ which
-are known or fixed to obtain corresponding values $\y \in \R^{n_\text{known}}$,
-$\Aeq \in \R^{m_\text{eq} \times n}$ and $\beq \in \R^{m_\text{eq} \times n}$
-are linear \emph{equality} coefficients and right-hand sides respectively, and
-$\lx, \ux \in \R^n$ are vectors of constant lower and upper bound constraints.
-
-\todo{This notation is unfortunate. Too many bold A's and too many capitals.}
-
-This description exactly matches the prototype used by the
-\texttt{igl::active\_set()} function.
-
-The active set method works by iteratively treating a subset (some rows) of the
-inequality constraints as equality constraints. These are called the ``active
-set'' of constraints. So at any given iterations $i$ we might have a new
-problem:
-\begin{align}
-\argmin \limits_\z &
-  \z^\transpose \A \z + \z^\transpose \b + \text{ constant}\\
-\text{subject to } & \z_\text{known}^i = \y^i,\\
-                   & \Aeq^i \z = \beq^i,
-\end{align}
-where the active rows from 
-$\lx ≤ \z ≤ \ux$ and $\Aieq \z ≤ \bieq$  have been appended into
-$\z_\text{known}^i = \y^i$ and $\Aeq^i \z = \beq^i$ respectively.
-
-This may be optimized by solving a sparse linear system, resulting in the
-current solution $\z^i$. For equality constraint we can also find a
-corresponding Lagrange multiplier value. The active set method works by adding
-to the active set all linear inequality constraints which are violated by the
-previous solution $\z^{i-1}$ before this solve and then after the solve
-removing from the active set any constraints with negative Lagrange multiplier
-values. Let's declare that $\lameq \in \R^{m_\text{eq}}$, $\lamieq \in
-\R^{m_\text{ieq}}$, and $\lamlu \in \R^{n}$ are the Lagrange multipliers
-corresponding to the linear equality, linear inequality and constant bound
-constraints respectively. Then the abridged active set method proceeds as
-follows:
-
-\begin{lstlisting}[keywordstyle=,mathescape]
-while not converged
-  add to active set all rows where $\Aieq \z > \bieq$, $\z < \lx$ or $\z > \ux$
-    $\Aeq^i,\beq^i \leftarrow \Aeq,\beq + \text{active rows of} \Aieq,\bieq$
-    $\z_\text{known}^i,\y^i \leftarrow \z_\text{known},\y + \text{active indices and values of} \lx,\ux$
-  solve problem treating active constraints as equality constraints $\rightarrow \z,\lamieq,\lamlu$
-  remove from active set all rows with $\lamieq < 0$ or $\lamlu < 0$
-end
-\end{lstlisting}
-
-The fixed values constraints of $\z_\text{known}^i = \y^i$ may be enforced by
-substituting $\z_\text{known}^i$ for $\y^i$ in the energy directly during the
-solve.  Corresponding Lagrange multiplier values $\lambda_\text{known}^i$ can
-be recovered after we've found the rest of $\z$.
-
-The linear equality constraints $\Aeq^i \z = \beq^i$ are a little trickier. If
-the rows of 
-$\Aeq^i \in \R^{(<m_\text{eq}+ m_\text{ieq}) \times n}$ are linearly
-independent then it is straightforward how to build a Lagrangian which enforces
-each constraint. This results in solving a system roughly of the form:
-\begin{equation}
-\left(
-\begin{array}{cc}
-\A      & {\Aeq^i}^\transpose \\
-\Aeq^i  & \mat{0}
-\end{array}
-\right)
-\left(
-\begin{array}{l}
-\z\\\lambda^i_\text{eq}
-\end{array}
-\right)
-=
-\left(
-\begin{array}{l}
--\onehalf\b\\
--\beq^i
-\end{array}
-\right)
-\end{equation}
-\todo{Double check. Could be missing a sign or factor of 2 here.}
-
-If the rows of $\Aeq^i$ are linearly dependent then the system matrix above
-will be singular. Because we may always assume that the constraints are not
-contradictory, this system can still be solved. But it will take some care.
-Some linear solvers, e.g.\ \textsc{MATLAB}'s, seem to deal with these OK.
-\textsc{Eigen}'s does not.
-
-Without loss of generality, let us assume that there are no inequality
-constraints and it's the rows of $\Aeq$ which might be linearly dependent.
-Let's also assume we have no fixed or known values. Then the linear system
-above corresponds to the following optimization problem:
-\begin{align}
-\argmin \limits_\z &
-  \z^\transpose \A \z + \z^\transpose \b + \text{ constant}\\
-\text{subject to } & \Aeq \z = \beq.
-\end{align}
-For the sake of cleaner notation, let $m = m_\text{eq}$ so that $\Aeq \in \R^{m
-\times n}$ and $\beq \in \R^m$.
-
-We can construct the null space of the constraints by computing a QR
-decomposition of $\Aeq^\transpose$:
-\begin{equation}
-\Aeq^\transpose \P = \Q \RR =
-\left(\begin{array}{cc}
-\Q_1 & \Q_2
-\end{array}\right)
-\left(\begin{array}{c}
-\RR\\
-\mat{0}
-\end{array}\right)=
-\Q_1 \RR
-\end  {equation}
-where $\P \in \R^{m \times m}$ is a sparse permutation matrix, $\Q \in \R^{n
-\times n}$ is orthonormal, $\RR \in \R^{n \times m}$ is upper triangular. Let
-$r$ be the row rank of $\Aeq$---the number of linearly independent rows---then
-we split $\Q$ and $\RR$ into $\Q_1 \in \R^{n \times r}$, $\Q_2 \in \R^{n \times
-n-r}$, and $\R_1 \in \RR^{r \times m}$.
-
-Notice that 
-\begin{align}
-\Aeq \z &= \beq \\
-\P \P^\transpose \Aeq \z &= \\
-\P \left(\Aeq^\transpose \P\right)^\transpose \z &= \\
-\P\left(\Q\RR\right)^\transpose \z &= \\
-\P \RR^\transpose \Q^\transpose \z &= \\
-\P \RR^\transpose \Q \z &= \\
-\P \left(\RR_1^\transpose \mat{0}\right)
-  \left(\begin{array}{c}
-  \Q_1^\transpose\\
-  \Q_2^\transpose
-  \end{array}\right) \z &=\\
-\P \RR_1^\transpose \Q_1^\transpose \z + \P \0 \Q_2^\transpose \z &= .
-\end{align}
-Here we see that $\Q_1^\transpose \z$ affects this equality but
-$\Q_2^\transpose \z$ does not. Thus
-we say that $\Q_2$ forms a basis that spans the null space of $\Aeq$. That is,
-$\Aeq \Q_2 \w_2 =
-\vc{0},
-\forall \w_2 \in \R^{n-r}$. Let's write 
-\begin{equation}
-\z = \Q_1 \w_1 + \Q_2 \w_2.
-\end{equation}
-We're only interested in solutions $\z$ which satisfy our equality constraints
-$\Aeq \z = \beq$. If we plug in the above then we get
-\begin{align}
-\P \Aeq \left(\Q_1 \w_1 + \Q_2 \w_2\right) &= \beq \\
-\P \Aeq \Q_1 \w_1 &= \beq \\
-\P \left(\RR^\transpose \quad
-\mat{0}\right)\left(\begin{array}{c}\Q_1^\transpose\\
-\Q_2^\transpose\end{array}\right) \Q_1 \w_1 &= \beq \\
-\P \RR^\transpose \Q_1^\transpose \Q_1 \w_1 &= \beq \\
-\P \RR^\transpose \w_1 &= \beq \\
-\w_1 &= {\RR^\transpose}^{-1} \P^\transpose \beq.
-\end{align}
-% P RT Q1T z = Beq
-% RT Q1T z = PT Beq
-% Q1T z = RT \ PT Beq
-% w1 = RT \ PT * Beq
-% z = Q2 w2 + Q1 * w1 
-% z = Q2 w2 + Q1 * RT \ P * Beq
-% z = Q2 w2 + lambda0
-So then
-\begin{align}
-z &= \Q_1 {\RR^\transpose}^{-1} \P^\transpose \beq  + \Q_2 \w_2 \\
-z &= \Q_2 \w_2 + \overline{\w}_1\\
-\end{align}
-where $\overline{\w}_1$ is just collecting the known terms that do not depend on the
-unknowns $\w_2$.
-
-Now, by construction, $\Q_2 \w_2 + \overline{\w}_1$ spans (all) possible solutions
-for $\z$  that satisfy $\Aeq \z = \beq$
-% Let:
-% z = Q2 w2 + Q1 * w1 
-% z = Q2 w2 + Q1 * RT \ P * Beq
-% z = Q2 w2 + lambda0
-% Q2 w2 + lambda0 spans all z such that Aeq z = Beq, we have essentially
-% "factored out" the constraints. We can then simply make this substitution
-% into our original optimization problem, leaving an \emph{unconstrained}
-% quadratic energy optimization
-\begin{align}
-\argmin \limits_{\w_2} & (\Q_2 \w_2 + \overline{\w}_1)^\transpose \A (\Q_2 \w_2 + \overline{\w}_1)  +
-(\Q_2 \w_2 + \overline{\w}_1)^\transpose \b + \textrm{ constant}\\
-&\text{\emph{or equivalently}} \notag\\
-\argmin \limits_{\w_2} & 
-\w_2^\transpose \Q_2^\transpose \A \Q_2 + \w_2^\transpose \left(2\Q_2^\transpose \A
-\Q_2 \overline{\w}_1 + \Q_2^\transpose \b\right) + \textrm{ constant}
-\end{align}
-which is solved with a linear system solve
-\begin{equation}
-\left(\Q_2^\transpose \A \Q_2\right) \w_2 = 
-  -2\Q_2^\transpose \A \Q_2 \overline{\w}_1 + \Q_2^\transpose \b
-\end{equation}
-\todo{I've almost surely lost a sign or factor of 2 here.}
-If $\A$ is semi-positive definite then $\Q_2^\transpose\A\Q_2$ will also be
-semi-positive definite. Thus Cholesky factorization can be used. If $\A$ is
-sparse and $\Aeq$ is sparse and a good column pivoting sparse QR library is
-used then $\Q_2$ will be sparse and so will $\Q_2^\transpose\A\Q_2$.
-
-
-After $\w_2$ is known we can compute the solution $\z$ as 
-\begin{equation}
-\z = \Q_2 \w_2 + \overline{\w}_1.
-\end{equation}
-
-Recovering the Lagrange multipliers $\lameq \in \R^{m}$ corresponding to $\Aeq
-\z = \b$ is a bit trickier. The equation is:
-\begin{align}
-\left(\A\quad \Aeq^\transpose\right)
-\left(\begin{array}{c}
-\z\\
-\lameq
-\end{array}\right) &= -\onehalf\b \\
-\Aeq^\transpose \lameq &= - \A \z - \onehalf\b \\
-\Q_1^\transpose \Aeq^\transpose \lameq &= - \Q_1^\transpose \A \z - \onehalf\Q_1^\transpose \b \\
-% AT P = Q_1 R
-% AT = Q_1 R PT
-%
-% Q1T AT
-% Q1T Q_1 R PT
-% R PT
-& \textrm{\emph{Recall that $\Aeq^\transpose\P = \Q_1\RR$ so then $
-\Q_1^\transpose \Aeq^\transpose = \RR \P^\transpose$}}\notag\\
-\RR \P^\transpose \lameq &= - \Q_1^\transpose \A \z - \onehalf\Q_1^\transpose \b \\
-\lameq &= \P
-\RR^{-1}
-\left(-\Q_1^\transpose \A \z - \onehalf\Q_1^\transpose \b\right).
-\end{align}
-
-\alec{If $r<m$ then I think it is enough to define the permutation matrix as
-rectangular $\P \in \R^{m \times r}$.}
-
-We assumes that there were now known or fixed values. If there are, special
-care only needs to be taken to adjust $\Aeq \z = \beq$ to account for
-substituting $\z_\text{known}$ with $\y$. 
-
-\hr
-\clearpage
-
-The following table describes the translation between the entities described in
-this document and those used in \texttt{igl/active\_set} and
-\texttt{igl/min\_quad\_with\_fixed}.
-\begin{lstlisting}[keywordstyle=,mathescape]
-$\A$: A
-$\b$: B
-$\z$: Z
-$\Aeq$: Aeq
-$\Aieq$: Aieq
-$\beq$: Beq
-$\bieq$: Bieq
-$\lx$: lx
-$\ux$: ux
-$\lamieq$: ~Lambda_Aieq_i
-$\lamlu$: ~Lambda_known_i
-$\P$: AeqTE
-$\Q$: AeqTQ
-$\RR$: AeqTR1
-$\RR^\transpose$: AeqTR1T
-$\Q_1$: AeqTQ1
-$\Q_1^\transpose$: AeqTQ1T
-$\Q_2$: AeqTQ2
-$\Q_2^\transpose$: AeqTQ2T
-$\w_2$: lambda
-$\overline{\w}_1$: lambda_0
-\end{lstlisting}
-
-
-% Resources
-%  http://www.cs.cornell.edu/courses/cs322/2007sp/notes/qr.pdf
-%  http://www.math.uh.edu/~rohop/fall_06/Chapter2.pdf
-%  http://www.math.uh.edu/~rohop/fall_06/Chapter3.pdf
-
-
-\end{document}

documentation/implemented-papers.tex

@@ -1,50 +0,0 @@
-\documentclass[12pt]{diary}
-\immediate\write18{bibtex \jobname}
-
-\title{Papers implemented in \textsc{libigl}}
-\author{Alec Jacobson}
-\date{last revised 22 April 2014}
-
-\begin{document}
-This document serves as a companion reference to better list the references to
-scientific articles implemented within \textsc{libigl}. It will no doubt be
-incomplete.
-
-\paragraph{\texttt{cotmatrix}, \texttt{massmatrix}}
-build discrete operators on triangle and tetrahedral meshes. 
-\cite{Pinkall:1993:CDM,meyer03ddo,Jacobson:THESIS:2013}. For tet meshes, we no
-longer use the ``by the book'' FEM construction \`a la \cite{Sharf:2007vv},
-rather a purely geometric approach \cite{Barth:1994,Xu:1999}.
-
-\paragraph{\texttt{harmonic}} solves a Laplace problem (equivalently
-minimizes the Dirichlet energy) with some simple boundary conditions
-\cite{HarmonicCoodinates07}. There's also an option to solve
-``higher order Laplace problems'' (bi-Laplace, tri-Laplace, etc.)
-\cite{Botsch:2004:AIF,sorkine04lsm,Jacobson:MixedFEM:2010}.
-
-\paragraph{\texttt{bbw/}} implements ``bounded biharmonic
-weights'' \cite{Jacobson:BBW:2011}.
-
-\paragraph{\texttt{svd3x3/arap}} is a generalized implementation
-for solving ``as-rigid-as-possible'' (ARAP) mesh deformation or parameterization
-problems \cite{ARAP_modeling:2007,Liu:2008:ALA,Chao:2010:ASG}.
-
-\paragraph{\texttt{svd3x3/arap\_dof}} implements ``FAST'',
-which is simultaneously a reduced form of ARAP and a method for automatically
-choosing skinning transformations \cite{Jacobson:FAST:2012}.
-
-\paragraph{\texttt{dqs}} implements ``Dual quaternion skinning''
-\cite{Kavan:2008:GSW}.
-
-\paragraph{\texttt{lbs}} implements ``linear blend skinning'', also known as
-``skeletal subspace deformation'', or ``enveloping''. This technique is often
-attributed to \cite{Magnenat-Thalmann:1988:JLD}.
-
-\paragraph{\texttt{winding\_number}} implements ``Generalized Winding Numbers''
-\cite{Jacobson:WN:2013}
-
-\bibliographystyle{acmsiggraph}
-\bibliography{references} 
-
-\end{document}
-__END__

documentation/references.bib.REMOVED.git-id

@@ -1 +0,0 @@
-346cd2c8206ee2c83e7aef50adb832e483320d09

include/igl/AABB.h

@@ -484,7 +484,7 @@ inline std::vector<int> igl::AABB<DerivedV,DIM>::find(
   if(is_leaf())
   {
     // Initialize to some value > -epsilon
-    Scalar a1=1,a2=1,a3=1,a4=1;
+    Scalar a1=0,a2=0,a3=0,a4=0;
     switch(DIM)
     {
       case 3:
@@ -511,7 +511,6 @@ inline std::vector<int> igl::AABB<DerivedV,DIM>::find(
           // Hack for now to keep templates simple. If becomes bottleneck
           // consider using std::enable_if_t 
           const Vector2S q2 = q.head(2);
-          Scalar a0 = doublearea_single(V1,V2,V3);
           a1 = doublearea_single(V1,V2,q2);
           a2 = doublearea_single(V2,V3,q2);
           a3 = doublearea_single(V3,V1,q2);
@@ -519,6 +518,12 @@ inline std::vector<int> igl::AABB<DerivedV,DIM>::find(
         }
       default:assert(false);
     }
+    // Normalization is important for correcting sign
+    Scalar sum = a1+a2+a3+a4;
+    a1 /= sum;
+    a2 /= sum;
+    a3 /= sum;
+    a4 /= sum;
     if(
         a1>=-epsilon && 
         a2>=-epsilon && 
@@ -624,7 +629,7 @@ igl::AABB<DerivedV,DIM>::squared_distance(
   using namespace std;
   using namespace igl;
   Scalar sqr_d = min_sqr_d;
-  assert(DIM == 3 && "Code has only been tested for DIM == 3");
+  //assert(DIM == 3 && "Code has only been tested for DIM == 3");
   assert((Ele.cols() == 3 || Ele.cols() == 2 || Ele.cols() == 1)
     && "Code has only been tested for simplex sizes 3,2,1");
 
@@ -1007,9 +1012,14 @@ inline void igl::AABB<DerivedV,DIM>::leaf_squared_distance(
       Ele(m_primitive,1) != Ele(m_primitive,2) && 
       Ele(m_primitive,2) != Ele(m_primitive,0))
   {
-    const RowVectorDIMS v10 = (V.row(Ele(m_primitive,1))- V.row(Ele(m_primitive,0)));
-    const RowVectorDIMS v20 = (V.row(Ele(m_primitive,2))- V.row(Ele(m_primitive,0)));
-    const RowVectorDIMS n = v10.cross(v20);
+    assert(DIM == 3 && "Only implemented for 3D triangles");
+    typedef Eigen::Matrix<Scalar,1,3> RowVector3S;
+    // can't be const because of annoying DIM template
+    RowVector3S v10(0,0,0);
+    v10.head(DIM) = (V.row(Ele(m_primitive,1))- V.row(Ele(m_primitive,0)));
+    RowVector3S v20(0,0,0);
+    v20.head(DIM) = (V.row(Ele(m_primitive,2))- V.row(Ele(m_primitive,0)));
+    const RowVectorDIMS n = (v10.cross(v20)).head(DIM);
     Scalar n_norm = n.norm();
     if(n_norm > 0)
     {

include/igl/barycenter.cpp

@@ -39,4 +39,5 @@ template void igl::barycenter<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::M
 template void igl::barycenter<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::barycenter<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template void igl::barycenter<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::barycenter<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&);
 #endif

include/igl/boolean/mesh_boolean.cpp

@@ -10,6 +10,7 @@
 #include <igl/cgal/peel_outer_hull_layers.h>
 #include <igl/cgal/remesh_self_intersections.h>
 #include <igl/remove_unreferenced.h>
+#include <igl/mod.h>
 #include <igl/unique_simplices.h>
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 #include <iostream>
@@ -199,9 +200,11 @@ IGL_INLINE void igl::boolean::mesh_boolean(
 #endif
   Matrix<bool,Dynamic,1> from_A(CF.rows());
   // peel layers keeping track of odd and even flips
-  Matrix<bool,Dynamic,1> odd;
+  VectorXi I;
   Matrix<bool,Dynamic,1> flip;
-  peel_outer_hull_layers(EV,CF,CN,odd,flip);
+  peel_outer_hull_layers(EV,CF,CN,I,flip);
+  // 0 is "first" iteration, so it's odd
+  Array<bool,Dynamic,1> odd = igl::mod(I,2).array()==0;
 
 #ifdef IGL_MESH_BOOLEAN_DEBUG
   cout<<"categorize..."<<endl;

include/igl/cgal/peel_outer_hull_layers.cpp

@@ -21,13 +21,13 @@ template <
   typename DerivedV,
   typename DerivedF,
   typename DerivedN,
-  typename Derivedodd,
+  typename DerivedI,
   typename Derivedflip>
 IGL_INLINE size_t igl::cgal::peel_outer_hull_layers(
   const Eigen::PlainObjectBase<DerivedV > & V,
   const Eigen::PlainObjectBase<DerivedF > & F,
   const Eigen::PlainObjectBase<DerivedN > & N,
-  Eigen::PlainObjectBase<Derivedodd > & odd,
+  Eigen::PlainObjectBase<DerivedI> & I,
   Eigen::PlainObjectBase<Derivedflip > & flip)
 {
   using namespace Eigen;
@@ -52,12 +52,11 @@ IGL_INLINE size_t igl::cgal::peel_outer_hull_layers(
   // keep track of iteration parity and whether flipped in hull
   MatrixXF Fr = F;
   MatrixXN Nr = N;
-  odd.resize(m,1);
+  I.resize(m,1);
   flip.resize(m,1);
   // Keep track of index map
   MatrixXI IM = MatrixXI::LinSpaced(m,0,m-1);
   // This is O(n * layers)
-  bool odd_iter = true;
   MatrixXI P(m,1);
   Index iter = 0;
   while(Fr.size() > 0)
@@ -82,7 +81,7 @@ IGL_INLINE size_t igl::cgal::peel_outer_hull_layers(
     vector<bool> in_outer(Fr.rows(),false);
     for(Index g = 0;g<Jo.rows();g++)
     {
-      odd(IM(Jo(g))) = odd_iter;
+      I(IM(Jo(g))) = iter;
       P(IM(Jo(g))) = iter;
       in_outer[Jo(g)] = true;
       flip(IM(Jo(g))) = flipr(Jo(g));
@@ -108,7 +107,6 @@ IGL_INLINE size_t igl::cgal::peel_outer_hull_layers(
         }
       }
     }
-    odd_iter = !odd_iter;
     iter++;
   }
   return iter;
@@ -117,23 +115,23 @@ IGL_INLINE size_t igl::cgal::peel_outer_hull_layers(
 template <
   typename DerivedV,
   typename DerivedF,
-  typename Derivedodd,
+  typename DerivedI,
   typename Derivedflip>
 IGL_INLINE size_t igl::cgal::peel_outer_hull_layers(
   const Eigen::PlainObjectBase<DerivedV > & V,
   const Eigen::PlainObjectBase<DerivedF > & F,
-  Eigen::PlainObjectBase<Derivedodd > & odd,
+  Eigen::PlainObjectBase<DerivedI > & I,
   Eigen::PlainObjectBase<Derivedflip > & flip)
 {
   using namespace std;
   Eigen::Matrix<typename DerivedV::Scalar,DerivedF::RowsAtCompileTime,3> N;
   per_face_normals(V,F,N);
-  return peel_outer_hull_layers(V,F,N,odd,flip);
+  return peel_outer_hull_layers(V,F,N,I,flip);
 }
 
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template size_t igl::cgal::peel_outer_hull_layers<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<bool, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
-template size_t igl::cgal::peel_outer_hull_layers<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<bool, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
+template size_t igl::cgal::peel_outer_hull_layers<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template size_t igl::cgal::peel_outer_hull_layers<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/cgal/peel_outer_hull_layers.h

@@ -22,8 +22,7 @@ namespace igl
     //   F  #F by 3 list of triangle indices into V
     //   N  #F by 3 list of per-face normals
     // Outputs:
-    //   odd  #F list of whether facet belongs to an odd iteration peel (otherwise
-    //     an even iteration peel)
+    //   I  #F list of which peel Iation a facet belongs 
     //   flip  #F list of whether a facet's orientation was flipped when facet
     //     "peeled" into its associated outer hull layer.
     // Returns number of peels
@@ -31,23 +30,23 @@ namespace igl
       typename DerivedV,
       typename DerivedF,
       typename DerivedN,
-      typename Derivedodd,
+      typename DerivedI,
       typename Derivedflip>
     IGL_INLINE size_t peel_outer_hull_layers(
       const Eigen::PlainObjectBase<DerivedV > & V,
       const Eigen::PlainObjectBase<DerivedF > & F,
       const Eigen::PlainObjectBase<DerivedN > & N,
-      Eigen::PlainObjectBase<Derivedodd > & odd,
+      Eigen::PlainObjectBase<DerivedI > & I,
       Eigen::PlainObjectBase<Derivedflip > & flip);
     template <
       typename DerivedV,
       typename DerivedF,
-      typename Derivedodd,
+      typename DerivedI,
       typename Derivedflip>
     IGL_INLINE size_t peel_outer_hull_layers(
       const Eigen::PlainObjectBase<DerivedV > & V,
       const Eigen::PlainObjectBase<DerivedF > & F,
-      Eigen::PlainObjectBase<Derivedodd > & odd,
+      Eigen::PlainObjectBase<DerivedI > & I,
       Eigen::PlainObjectBase<Derivedflip > & flip);
   }
 }

include/igl/mod.cpp

@@ -19,3 +19,13 @@ IGL_INLINE void igl::mod(
     *(B.data()+i) = (*(A.data()+i))%base;
   }
 }
+template <typename DerivedA>
+IGL_INLINE Eigen::PlainObjectBase<DerivedA> igl::mod(
+  const Eigen::PlainObjectBase<DerivedA> & A, const int base)
+{
+  Eigen::PlainObjectBase<DerivedA> B;
+  mod(A,base,B);
+  return B;
+}
+#ifdef IGL_STATIC_LIBRARY
+#endif

include/igl/mod.h

@@ -23,6 +23,9 @@ namespace igl
     const Eigen::PlainObjectBase<DerivedA> & A,
     const int base,
     Eigen::PlainObjectBase<DerivedB> & B);
+  template <typename DerivedA>
+  IGL_INLINE Eigen::PlainObjectBase<DerivedA> mod(
+    const Eigen::PlainObjectBase<DerivedA> & A, const int base);
 }
 #ifndef IGL_STATIC_LIBRARY
 #include "mod.cpp"

include/igl/point_mesh_squared_distance.cpp

@@ -23,9 +23,26 @@ IGL_INLINE void igl::point_mesh_squared_distance(
   Eigen::PlainObjectBase<DerivedC> & C)
 {
   using namespace std;
-  AABB<DerivedV,3> tree;
-  tree.init(V,Ele);
-  return tree.squared_distance(V,Ele,P,sqrD,I,C);
+  const size_t dim = P.cols();
+  assert((dim == 2 || dim == 3) && "P.cols() should be 2 or 3");
+  assert(P.cols() == V.cols() && "P.cols() should equal V.cols()");
+  switch(dim)
+  {
+    default:
+      // fall-through and pray
+    case 3:
+    {
+      AABB<DerivedV,3> tree;
+      tree.init(V,Ele);
+      return tree.squared_distance(V,Ele,P,sqrD,I,C);
+    }
+    case 2:
+    {
+      AABB<DerivedV,2> tree;
+      tree.init(V,Ele);
+      return tree.squared_distance(V,Ele,P,sqrD,I,C);
+    }
+  }
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/pseudonormal_test.cpp

@@ -62,3 +62,44 @@ IGL_INLINE void igl::pseudonormal_test(
   }
   s = (qc.dot(n) >= 0 ? 1. : -1.);
 }
+
+IGL_INLINE void igl::pseudonormal_test(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & E,
+  const Eigen::MatrixXd & EN,
+  const Eigen::MatrixXd & VN,
+  const Eigen::RowVector2d & q,
+  const int e,
+  const Eigen::RowVector2d & c,
+  double & s,
+  Eigen::RowVector2d & n)
+{
+  using namespace Eigen;
+  const auto & qc = q-c;
+  const double len = (V.row(E(e,1))-V.row(E(e,0))).norm();
+  // barycentric coordinates
+  RowVector2d b((c-V.row(E(e,1))).norm()/len,(c-V.row(E(e,0))).norm()/len);
+  // Determine which normal to use
+  const double epsilon = 1e-12;
+  const int type = (b.array()<=epsilon).cast<int>().sum();
+  switch(type)
+  {
+    case 1:
+      // Find vertex
+      for(int x = 0;x<2;x++)
+      {
+        if(b(x)>epsilon)
+        {
+          n = VN.row(E(e,x));
+          break;
+        }
+      }
+      break;
+    default:
+      assert(false && "all barycentric coords zero.");
+    case 0:
+      n = EN.row(e);
+      break;
+  }
+  s = (qc.dot(n) >= 0 ? 1. : -1.);
+}

include/igl/pseudonormal_test.h

@@ -41,6 +41,16 @@ namespace igl
     const Eigen::RowVector3d & c,
     double & s,
     Eigen::RowVector3d & n);
+  IGL_INLINE void pseudonormal_test(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & E,
+    const Eigen::MatrixXd & EN,
+    const Eigen::MatrixXd & VN,
+    const Eigen::RowVector2d & q,
+    const int i,
+    const Eigen::RowVector2d & c,
+    double & s,
+    Eigen::RowVector2d & n);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "pseudonormal_test.cpp"

include/igl/remove_duplicate_vertices.cpp

@@ -78,4 +78,5 @@ IGL_INLINE void igl::remove_duplicate_vertices(
 // Explicit instanciation
 template void igl::remove_duplicate_vertices<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::remove_duplicate_vertices<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+template void igl::remove_duplicate_vertices<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/signed_distance.cpp

@@ -26,22 +26,34 @@ IGL_INLINE void igl::signed_distance(
 {
   using namespace Eigen;
   using namespace std;
-  assert(V.cols() == 3 && "V should have 3d positions");
-  assert(P.cols() == 3 && "P should have 3d positions");
+  const int dim = V.cols();
+  assert((V.cols() == 3||V.cols() == 2) && "V should have 3d or 2d positions");
+  assert((P.cols() == 3||V.cols() == 2) && "P should have 3d or 2d positions");
+  assert(V.cols() == P.cols() && "V should have same dimension as P");
   // Only unsigned distance is supported for non-triangles
   if(sign_type != SIGNED_DISTANCE_TYPE_UNSIGNED)
   {
-    assert(F.cols() == 3 && "F should have triangles");
+    assert(F.cols() == dim && "F should have co-dimension 0 simplices");
   }
 
   // Prepare distance computation
-  AABB<MatrixXd,3> tree;
-  tree.init(V,F);
+  AABB<MatrixXd,3> tree3;
+  AABB<MatrixXd,2> tree2;
+  switch(dim)
+  {
+    default:
+    case 3:
+      tree3.init(V,F);
+      break;
+    case 2:
+      tree2.init(V,F);
+      break;
+  }
 
   Eigen::MatrixXd FN,VN,EN;
   Eigen::MatrixXi E;
   Eigen::VectorXi EMAP;
-  WindingNumberAABB<Eigen::Vector3d> hier;
+  WindingNumberAABB<Eigen::Vector3d> hier3;
   switch(sign_type)
   {
     default:
@@ -51,28 +63,72 @@ IGL_INLINE void igl::signed_distance(
       break;
     case SIGNED_DISTANCE_TYPE_DEFAULT:
     case SIGNED_DISTANCE_TYPE_WINDING_NUMBER:
-      hier.set_mesh(V,F);
-      hier.grow();
+      switch(dim)
+      {
+        default:
+        case 3:
+          hier3.set_mesh(V,F);
+          hier3.grow();
+          break;
+        case 2:
+          // no precomp, no hierarchy
+          break;
+      }
       break;
     case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
-      // "Signed Distance Computation Using the Angle Weighted Pseudonormal"
-      // [Bærentzen & Aanæs 2005]
-      per_face_normals(V,F,FN);
-      per_vertex_normals(V,F,PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE,FN,VN);
-      per_edge_normals(
-        V,F,PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM,FN,EN,E,EMAP);
-      N.resize(P.rows(),3);
+      switch(dim)
+      {
+        default:
+        case 3:
+          // "Signed Distance Computation Using the Angle Weighted Pseudonormal"
+          // [Bærentzen & Aanæs 2005]
+          per_face_normals(V,F,FN);
+          per_vertex_normals(V,F,PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE,FN,VN);
+          per_edge_normals(
+            V,F,PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM,FN,EN,E,EMAP);
+          break;
+        case 2:
+          FN.resize(F.rows(),2);
+          VN = MatrixXd::Zero(V.rows(),2);
+          for(int e = 0;e<F.rows();e++)
+          {
+            // rotate edge vector
+            FN(e,0) = -(V(F(e,1),1)-V(F(e,0),1));
+            FN(e,1) =  (V(F(e,1),0)-V(F(e,0),0));
+            FN.row(e).normalize();
+            // add to vertex normal
+            VN.row(F(e,1)) += FN.row(e);
+            VN.row(F(e,0)) += FN.row(e);
+          }
+          // normalize to average
+          VN.rowwise().normalize();
+          break;
+      }
+      N.resize(P.rows(),dim);
       break;
   }
 
   S.resize(P.rows(),1);
   I.resize(P.rows(),1);
-  C.resize(P.rows(),3);
+  C.resize(P.rows(),dim);
   for(int p = 0;p<P.rows();p++)
   {
-    const RowVector3d q = P.row(p);
+    RowVector3d q3;
+    RowVector2d q2;
+    switch(P.cols())
+    {
+      default:
+      case 3:
+        q3 = P.row(p);
+        break;
+      case 2:
+        q2 = P.row(p);
+        break;
+    }
     double s,sqrd;
-    RowVector3d c;
+    RowVectorXd c;
+    RowVector3d c3;
+    RowVector2d c2;
     int i=-1;
     switch(sign_type)
     {
@@ -80,23 +136,32 @@ IGL_INLINE void igl::signed_distance(
         assert(false && "Unknown SignedDistanceType");
       case SIGNED_DISTANCE_TYPE_UNSIGNED:
         s = 1.;
-        sqrd = tree.squared_distance(V,F,q,i,c);
+        sqrd = dim==3?
+          tree3.squared_distance(V,F,q3,i,c3):
+          tree2.squared_distance(V,F,q2,i,c2);
         break;
       case SIGNED_DISTANCE_TYPE_DEFAULT:
       case SIGNED_DISTANCE_TYPE_WINDING_NUMBER:
-        signed_distance_winding_number(tree,V,F,hier,q,s,sqrd,i,c);
+        dim==3 ? 
+          signed_distance_winding_number(tree3,V,F,hier3,q3,s,sqrd,i,c3):
+          signed_distance_winding_number(tree2,V,F,q2,s,sqrd,i,c2);
         break;
       case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
       {
-        Eigen::RowVector3d n;
-        signed_distance_pseudonormal(tree,V,F,FN,VN,EN,EMAP,q,s,sqrd,i,c,n);
+        RowVector3d n3;
+        RowVector2d n2;
+        dim==3 ?
+          signed_distance_pseudonormal(tree3,V,F,FN,VN,EN,EMAP,q3,s,sqrd,i,c3,n3):
+          signed_distance_pseudonormal(tree2,V,F,FN,VN,q2,s,sqrd,i,c2,n2);
+        Eigen::RowVectorXd n;
+        (dim==3 ? n = n3 : n = n3);
         N.row(p) = n;
         break;
       }
     }
     I(p) = i;
     S(p) = s*sqrt(sqrd);
-    C.row(p) = c;
+    C.row(p) = (dim==3 ? c=c3 : c=c2);
   }
 }
 
@@ -191,6 +256,25 @@ IGL_INLINE void igl::signed_distance_pseudonormal(
   pseudonormal_test(V,F,FN,VN,EN,EMAP,q,f,c,s,n);
 }
 
+IGL_INLINE void igl::signed_distance_pseudonormal(
+  const AABB<Eigen::MatrixXd,2> & tree,
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXd & FN,
+  const Eigen::MatrixXd & VN,
+  const Eigen::RowVector2d & q,
+  double & s,
+  double & sqrd,
+  int & f,
+  Eigen::RowVector2d & c,
+  Eigen::RowVector2d & n)
+{
+  using namespace Eigen;
+  using namespace std;
+  sqrd = tree.squared_distance(V,F,q,f,c);
+  pseudonormal_test(V,F,FN,VN,q,f,c,s,n);
+}
+
 IGL_INLINE double igl::signed_distance_winding_number(
   const AABB<Eigen::MatrixXd,3> & tree,
   const Eigen::MatrixXd & V,
@@ -210,12 +294,12 @@ IGL_INLINE void igl::signed_distance_winding_number(
   const AABB<Eigen::MatrixXd,3> & tree,
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
-  const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
-  const Eigen::RowVector3d & q,
+  const igl::WindingNumberAABB<Eigen::Matrix<double,3,1> > & hier,
+  const Eigen::Matrix<double,1,3> & q,
   double & s,
   double & sqrd,
   int & i,
-  Eigen::RowVector3d & c)
+  Eigen::Matrix<double,1,3> & c)
 {
   using namespace Eigen;
   using namespace std;
@@ -225,6 +309,26 @@ IGL_INLINE void igl::signed_distance_winding_number(
   s = 1.-2.*w;
 }
 
+IGL_INLINE void igl::signed_distance_winding_number(
+  const AABB<Eigen::MatrixXd,2> & tree,
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const Eigen::Matrix<double,1,2> & q,
+  double & s,
+  double & sqrd,
+  int & i,
+  Eigen::Matrix<double,1,2> & c)
+{
+  using namespace Eigen;
+  using namespace std;
+  using namespace igl;
+  sqrd = tree.squared_distance(V,F,q,i,c);
+  double w;
+  winding_number_2(V.data(), V.rows(), F.data(), F.rows(), q.data(), 1, &w);
+  s = 1.-2.*w;
+
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // This template is necessary for the others to compile with clang
 // http://stackoverflow.com/questions/27748442/is-clangs-c11-support-reliable

include/igl/signed_distance.h

@@ -77,6 +77,19 @@ namespace igl
   //   i  closest primitive
   //   c  closest point
   //   n  normal
+  IGL_INLINE void signed_distance_pseudonormal(
+    const Eigen::MatrixXd & P,
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const AABB<Eigen::MatrixXd,3> & tree,
+    const Eigen::MatrixXd & FN,
+    const Eigen::MatrixXd & VN,
+    const Eigen::MatrixXd & EN,
+    const Eigen::VectorXi & EMAP,
+    Eigen::VectorXd & S,
+    Eigen::VectorXi & I,
+    Eigen::MatrixXd & C,
+    Eigen::MatrixXd & N);
   IGL_INLINE void signed_distance_pseudonormal(
     const AABB<Eigen::MatrixXd,3> & tree,
     const Eigen::MatrixXd & V,
@@ -92,18 +105,17 @@ namespace igl
     Eigen::RowVector3d & c,
     Eigen::RowVector3d & n);
   IGL_INLINE void signed_distance_pseudonormal(
-    const Eigen::MatrixXd & P,
+    const AABB<Eigen::MatrixXd,2> & tree,
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
-    const AABB<Eigen::MatrixXd,3> & tree,
     const Eigen::MatrixXd & FN,
     const Eigen::MatrixXd & VN,
-    const Eigen::MatrixXd & EN,
-    const Eigen::VectorXi & EMAP,
-    Eigen::VectorXd & S,
-    Eigen::VectorXi & I,
-    Eigen::MatrixXd & C,
-    Eigen::MatrixXd & N);
+    const Eigen::RowVector2d & q,
+    double & s,
+    double & sqrd,
+    int & i,
+    Eigen::RowVector2d & c,
+    Eigen::RowVector2d & n);
 
   // Inputs:
   //   tree  AABB acceleration tree (see cgal/point_mesh_squared_distance.h)
@@ -124,12 +136,21 @@ namespace igl
     const AABB<Eigen::MatrixXd,3> & tree,
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
-    const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
-    const Eigen::RowVector3d & q,
+    const igl::WindingNumberAABB<Eigen::Matrix<double,3,1> > & hier,
+    const Eigen::Matrix<double,1,3> & q,
+    double & s,
+    double & sqrd,
+    int & i,
+    Eigen::Matrix<double,1,3> & c);
+  IGL_INLINE void signed_distance_winding_number(
+    const AABB<Eigen::MatrixXd,2> & tree,
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const Eigen::Matrix<double,1,2> & q,
     double & s,
     double & sqrd,
     int & i,
-    Eigen::RowVector3d & c);
+    Eigen::Matrix<double,1,2> & c);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/sort.cpp

@@ -231,6 +231,7 @@ template void igl::sort<long>(std::vector<long, std::allocator<long> > const&, b
 template void igl::sort<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::sort<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::sort<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::sort<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif
 
 #ifdef _WIN32

include/igl/tetgen/tetrahedralize.h

@@ -33,7 +33,7 @@ namespace igl
     // Outputs:
     //   TV  #V by 3 vertex position list
     //   TT  #T by 4 list of tet face indices
-    //   TF  #F by 3 list of trianlge face indices
+    //   TF  #F by 3 list of triangle face indices
     // Returns status:
     //   0 success
     //   1 tetgen threw exception

include/igl/triangle/triangulate.cpp

@@ -54,7 +54,7 @@ IGL_INLINE void igl::triangle::triangulate(
   using namespace Eigen;
 
   // Prepare the flags
-  string full_flags = flags + "pzBV";
+  string full_flags = flags + "pzB";
 
   // Prepare the input struct
   triangulateio in;

index.html

@@ -23,9 +23,7 @@
 <blockquote>
 <p>Get started with:</p>
 
-<pre><code class="bash">git clone https://github.com/libigl/libigl.git
-cd libigl
-git submodule update --init --recursive
+<pre><code class="bash">git clone --recursive https://github.com/libigl/libigl.git
 </code></pre>
 </blockquote>
 
@@ -154,6 +152,25 @@ Eigen::Matrix&lt;double, Eigen::Dynamic, 3&gt; A;
 them probably work just fine). This requires setting up unit testing, which is
 a major <em>todo</em> for our development.</p>
 
+<h2 id="gitsubmodules">Git Submodules</h2>
+
+<p>Libigl uses git submodules for its <em>optional</em> dependencies,
+in particular, those needed by the OpenGL viewer to run the examples in the
+<a href="tutorial/tutorial.html">tutorial</a>. Git submodules allow use to treat clones of
+other libraries as sub-directories within ours while separating our commits.
+Read the <a href="http://git-scm.com/docs/git-submodule">documentation</a> for a detailed
+explanation, but essentially our libigl repo stores a hash for each of its
+subrepos containing which version to update to. When a change is introduced in
+a dependencies repo we can incorporate that change by pulling in our sub-repo
+and updating (i.e. committing) that change to the hash.</p>
+
+<p>When pulling new changes to libigl it&#8217;s also a good idea to update changes to
+subrepos:</p>
+
+<pre><code class="bash">git pull
+git submodule update -- recursive
+</code></pre>
+
 <h2 id="howtocontribute">How to contribute</h2>
 
 <p>If you are interested in joining development, please fork the repository and
@@ -213,7 +230,8 @@ few labs/companies/institutions using libigl:</p>
 Jacobson</a> and <a href="http://www.inf.ethz.ch/personal/dpanozzo/">Daniele
 Panozzo</a>. Please <a href="&#109;&#97;&#105;&#108;&#x74;&#111;&#x3a;&#x61;&#x6c;&#101;&#x63;&#x6a;&#x61;&#x63;&#111;&#98;&#115;&#111;&#110;&#64;&#103;&#109;&#x61;&#105;&#108;&#x2e;&#x63;&#x6f;&#x6d;&#44;&#x64;&#x61;&#110;&#105;&#101;&#x6c;&#101;&#x2e;&#x70;&#x61;&#110;&#x6f;&#122;&#x7a;&#x6f;&#64;&#x67;&#x6d;&#x61;&#x69;&#x6c;&#x2e;&#99;&#x6f;&#x6d;">&#99;&#x6f;&#x6e;&#116;&#x61;&#99;&#x74;
 &#117;&#115;</a> if you have
-questions or comments. We are happy to get feedback!</p>
+questions or comments. For troubleshooting, please post an
+<a href="https://github.com/libigl/libigl/issues">issue</a> on github.</p>
 
 <p>If you&#8217;re using libigl in your projects, quickly <a href="&#x6d;&#x61;&#105;&#108;&#x74;&#x6f;&#58;&#97;&#x6c;&#101;&#99;&#106;&#x61;&#x63;&#111;&#x62;&#x73;&#x6f;&#110;&#64;&#103;&#x6d;&#x61;&#x69;&#x6c;&#46;&#99;&#x6f;&#x6d;&#x2c;&#100;&#x61;&#110;&#105;&#x65;&#108;&#x65;&#x2e;&#112;&#97;&#x6e;&#111;&#x7a;&#x7a;&#x6f;&#x40;&#103;&#109;&#x61;&#105;&#x6c;&#x2e;&#x63;&#111;&#109;">&#x64;&#x72;&#111;&#112; &#x75;&#x73; &#97;
 &#110;&#111;&#116;&#101;</a>. Tell us who you

matlab-to-eigen.html

@@ -255,7 +255,7 @@ IP.conservativeResize(stable_partition(
       <tr class=d0>
         <td><pre><code>B = A(R&lt;2,C&gt;1);</code></pre>
         <td><pre><code>B = igl::slice_mask(A,R.array()&lt;2,C.array()&gt;1);</code></pre>
-        <td>Where </td>
+        <td></td>
       </tr>
 
       <tr class=d1>
@@ -266,6 +266,12 @@ bool a = A.array().any();
         <td>Casts <code>Matrix::Scalar<code> to <code>bool</code>.</td>
       </tr>
 
+      <tr class=d0>
+        <td><pre><code>B = mod(A,2)</code></pre></td>
+        <td><pre><code>igl::mod(A,2,B)</code></pre></td>
+        <td></td>
+      </tr>
+
       <!-- Insert rows for each command pair -->
 
       <!-- Leave this here for copy and pasting

scripts/clone-and-build.sh

@@ -14,7 +14,7 @@
 #
 # and then add the line:
 #
-#     30 2 * * * /usr/local/igl/libigl/scripts/clone-and-build.sh
+#     30 2 * * * source /Users/ajx/.profile; /usr/local/igl/libigl/scripts/clone-and-build.sh
 #
 # replace the path above with the **full path** to this file.
 #