

\section{Spherical Shape Analogies}
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/analogies.jpg}
  \vspace{-5pt}
  \caption{We generate an output shape $B'$ that relates to the input $B$ in the same way as how the surface normal of a given primitive $A'$ relates to the surface normal of a sphere $A$. } 
  \label{fig:analogies}
  \vspace{-10pt}
\end{figure}
Our main idea is to use surface normals to capture the style of 3D objects: if two shapes share a similar normal ``profile'', we consider them to exhibit the same geometric style. 
Centered around this observation, as discussed in \refsec{intro}, we propose an analogy-based stylization method to translate the relationship between the normals of $A', A$ to create a stylized output shape $B'$ (see \reffig{analogies}). 
Throughout the paper, we use green color to denote the target style shape $A'$, gray color to denote the input shape $B$, and blue color to denote the output stylized shape $B'$.

\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/analogies_gauss.jpg}
  \vspace{-5pt}
  \caption{Our algorithm defines the analogous relation based on the surface normals. We first map the normals of the style shape $N_{A'}$ to a unit sphere to obtain $\widetilde{N}_{A'}$ (top row), transfer the relationship between $N_A$ and $\widetilde{N}_{A'}$ to the input shape to obtain the target normals $T$ (middle row), then optimize the input shape $B$ so that the actual output normals are aligned with the target normal $T$ (bottom row). } 
  \label{fig:analogies_gauss}
  \vspace{-10pt}
\end{figure}

Our algorithm consists of three simple steps, described in \reffig{analogies_gauss}: (1) we map the surface normal of $A'$ to a unit sphere $A$ in order to compute target normals on a sphere $\widetilde{N}_{A'}$, (2) we construct analogous target normals $T$ that relate to $N_{B}$ the same way $\widetilde{N}_{A'}$ relate to $N_{A}$, (3) we take $B, T$ as inputs and generate the stylized shape $B'$ whose normals approximate $T$ via optimization. 

\subsection{Generating $\widetilde{N}_{A'}$}\label{sec:generateNS}
Depending on the provided style shape $A'$ or user preferences, we consider three ways to get a set of target normals on a sphere $\widetilde{N}_{A'}$.

\textit{1. Closest normals.}\ 
The simplest case is when the style shape $A'$ is a simple convex shape with only few distinct face normals (e.g., icosahedron). We compute $\widetilde{N}_{A'}$ simply via snapping the normals of the sphere $N_{A}$ to the nearest normal in the style shape $N_{A'}$.

\textit{2. Spherical parameterization.}\ 
For a generic genus-0 shape (e.g., smooth or concave), we compute its parameterization to a sphere using, for example, conformalized mean curvature flow \cite{KazhdanSB12}. Then $\widetilde{N}_{A'}$ can be computed from the spherical parameterization.
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/meanCurvatureFlow.jpg}
  \vspace{-5pt}
  \caption{Given a style shape $A'$, we run mean curvature flow \cite{KazhdanSB12} to map the normals of style shape $N_{A'}$ to a sphere as $\widetilde{N}_{A'}$.} 
  \label{fig:meanCurvatureFlow}
  \vspace{-10pt}
\end{figure}

\textit{3. Normal Capture.}\ 
If one desires more control, one can manually specify $\widetilde{N}_{A'}$, in the spirit of how \emph{MatCap} (material capture \cite{SloanMGG01}) is used in rendering. We can then skip the first step in \reffig{analogies_gauss} and move on to the second step using the user-provided $\widetilde{N}_{A'}$.

\subsection{Generating $T$}\label{sec:generateT}
\begin{wrapfigure}[6]{r}{1.33in}
	\raggedleft
  \vspace{-12pt}
 
	\hspace*{-0.7\columnsep}
	\includegraphics[width=1.48in, trim={6mm 0mm -1mm 0mm}]{figures/gaussmap.jpg}
	\label{fig:gaussmap}
\end{wrapfigure}
Generating target normals $T$ on the input shape $B$ using analogy requires the correspondences between $A$, $B$. We compute the map using the \emph{Gauss map}, leveraging the fact that our $A$ 
is always a unit sphere (see the inset, where we use colors to visualize the correspondences).
\update{Specifically, the unit normal vector of each element (e.g., vertex or face) on the input shape $B$ can be equivalently interpreted as a point on the unit sphere $A$. Thus, we can easily map signals from $A$ back to $B$.}
Once the correspondences are obtained via the normals of input shape $N_{B}$, we can trivially compute $T$ by ``pasting'' $\widetilde{N}_{A'}$ on top of $B$.

\subsection{Generating $B'$}
After obtaining a set of target normals $T = \{ {\vecFont{t}}_k \}$ for each vertex $k$, our goal is to obtain a deformed output shape $B'$ whose surface normals approximate $T$. 
Let ${\matFont{V}}$ be a matrix of vertex locations with size $|{\matFont{V}}|$-by-3 and ${\matFont{F}}$ be the face list with size $|{\matFont{F}}|$-by-3 of the input shape $B$. Our output shape $B'$ is a deformed version of the input shape and we use $\widetilde{{\matFont{V}}}$ to denote the $|{\matFont{V}}|$-by-3 matrix of the deformed vertex locations. We formulate the normal-driven deformation as an energy optimization in the following form:
\begin{align}
  \min_{\widetilde{{\matFont{V}}}} \sum_{k \in {\matFont{V}}} E_R({\vecFont{v}}_k,{\vecFont{u}}_k) + \lambda a_k \| \hat{n}_k(\widetilde{{\matFont{V}}}) - {\vecFont{t}}_k \|_2^2,
\end{align} 
where $E_R$ denotes a regularization energy to preserve the details of the input mesh, and the second part measures the squared distance from the output unit surface normal $\hat{{\vecFont{n}}}_k(\widetilde{{\matFont{V}}})$ to the target output normal ${\vecFont{t}}_k$ at vertex $k$. We use $a_k$ to denote the Voronoi area of the vertex $k$, $\lambda$ is a weighting parameter to control the balance between the two terms, and ${\vecFont{v}}_k$ (${\vecFont{u}}_k$) is the input (output) location of vertex $k$.
In \reffig{lambda}, we can observe that using a small $\lambda$, the method preserves the input shape $B$. Using a bigger $\lambda$, the method favors in deforming the shape more into the style of $A'$.
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/lambda_more.jpg}
  \vspace{-5pt}
  \caption{\update{The $\lambda$ parameter in \refequ{arapNormalDriven} controls the balance between preserving the original shape and satisfying the desired style. We show different stylization results with increasing $\lambda$. \copyright Spiral Light Bulb (top) by benglish under CC BY-SA.}} 
  \label{fig:lambda}
  \vspace{-10pt}
\end{figure} 
The choice of $E_R$ depends on the user's intent. One can apply different regularizations to obtain different results.
For the purposes of this exposition, we introduce our optimization based on {\textsc{arap}}\xspace regularization in \refsec{normalDriven_arap}. We discuss how to extend to other regularizations in \refsec{regularizations}. 

\subsection{Normal-Driven Optimization with {\textsc{arap}}\xspace}\label{sec:normalDriven_arap}
We use ${\vecFont{e}}_{ij} \coloneqq {\vecFont{v}}_j - {\vecFont{v}}_i \in \mathbb{R}^3$ to denote the edge vector between vertices $i, j$ on the original mesh, and ${{\vecFont{e}}'}_{ij} \coloneqq \u_j - \u_i$ for the edge vectors on the deformed mesh. We can write down the energy that uses {\textsc{arap}}\xspace regularization as
\begin{align}\label{equ:arapOutputNormal}
  \min_{\widetilde{{\matFont{V}}},{\matFont{R}}}\ \sum_{k \in  V} \underbrace{\sum_{i,j \in \mathcal{N}_k} w_{ij} \| {\matFont{R}}_k {\vecFont{e}}_{ij} -  {{\vecFont{e}}'}_{ij} \|^2_2}_{E_{\textsc{arap}}\xspace} + \lambda a_k \| \hat{n}_k(\widetilde{{\matFont{V}}}) - {\vecFont{t}}_k \|_2^2,
\end{align}

\begin{wrapfigure}[6]{r}{0.9in}
	\raggedleft

	\vspace{-14pt} 
	\hspace*{-0.9\columnsep}
	\includegraphics[width=1.25in, trim={5mm 0mm 1mm 0mm}]{figures/spokes_rims.jpg} 
\end{wrapfigure}
We use $\mathcal{N}_k$ to denote the edge vectors of the \emph{spokes and rims} at vertex $k$ (see the inset) \cite{ChaoPSS10}, ${\matFont{R}}_k \in \text{SO(3)}$ to denote a 3-by-3 rotation matrix defined on $k$, and $w_{ij}$ is the cotangent weight of edge $i,j$ \cite{PinkallP93}. 
However, this energy is difficult to optimize because the term $\hat{n}_k(\widetilde{{\matFont{V}}})$ is non-linear in $\widetilde{{\matFont{V}}}$. 

We adapt the observation made in \cite{LiuCubic19} that the space of unit vectors can be captured by rotations. Thus, we can perform a change of variables by replacing $\hat{n}_k(\widetilde{{\matFont{V}}})$ with the \emph{rotated} unit normal of the input mesh ${\matFont{R}}_k \hat{{\vecFont{n}}}_k$ as
\begin{align} \label{equ:arapNormalDriven}
  \boxed{
  \min_{\widetilde{{\matFont{V}}},{\matFont{R}}}\ \sum_{k \in V} \sum_{i,j \in \mathcal{N}_k} w_{ij} \| {\matFont{R}}_k {\vecFont{e}}_{ij} -  {{\vecFont{e}}'}_{ij} \|^2_2  + \lambda a_k \| {\matFont{R}}_k \hat{{\vecFont{n}}}_k  - {\vecFont{t}}_k \|_2^2,}
\end{align}
where $\hat{{\vecFont{n}}}_k$ is the $k$th unit vertex normal of the input mesh computed via area-weighted average of face normals, which is constant throughout the optimization. This ${\matFont{R}}_k \hat{{\vecFont{n}}}_k$ can be perceived as an approximation of the area-weighted vertex normals of the output mesh $\hat{n}_k(\widetilde{{\matFont{V}}})$. In \reffig{RN_difference}, we visualize the difference between the output normals $\hat{n}_k(\widetilde{{\matFont{V}}})$ and the rotated input normals ${\matFont{R}}_k \hat{{\vecFont{n}}}_k$. We can notice that ${\matFont{R}}_k \hat{{\vecFont{n}}}_k$ is a decent approximation of the output vertex normals computed via area-weighted average. 
\update{We can observe that error tend to concentrate on high-curvature regions because discrete vertex normals are less accurate along on those regions and the {\textsc{arap}}\xspace regularization encourages smooth deformation.}
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/RN_difference.jpg}
  \vspace{-5pt}
  \caption{Empirically, we show that rotated input normals ${\matFont{R}}_k \hat{{\vecFont{n}}}_k$ is a good approximation of the area-weighted output vertex normals $\hat{n}_k(\widetilde{{\matFont{V}}})$. We can observe that the error mostly occurs on the high-curvature regions (right). \copyright Proto Paradigm under CC BY.} 
  \label{fig:RN_difference}
  \vspace{-10pt}
\end{figure} 
This change of variables allows us to solve for ${\matFont{R}}_k$s in parallel and make this energy quadratic in $\widetilde{{\matFont{V}}}$. In addition, the fact that ${\matFont{R}}_k$ is shared across the {\textsc{arap}}\xspace term and the normal term enables us to jointly consider both the regularization and the normal terms when obtaining the deformed vertex locations $\widetilde{{\matFont{V}}}$.

We minimize this energy via the local/global strategy \cite{SorkineA07}, where the local step involves solving a set of small \emph{Orthogonal Procrustes} problems and the global step amounts to a linear solve. For the sake of reproducibility, we reiterate the local-global steps for our energy in \refapp{localStep}, \refappnum{globalStep}. 
Non-linear methods, such as Newton’s method, could be applied to our scenario. It is however far slower than the local-global optimization since a single iteration of the Newton's method could be more expensive than 100 iterations of the local-global iterations (see \cite{LiuBK17}). Thus, it is less suitable for our interactive applications. 
Further accelerating our solver using other optimization methods (e.g., \cite{KovalskyGL16, PengDZGQL18, ZhuBK18}) should be possible, but is left as future work.

\section{Introduction}\label{sec:intro}



Analogies of the form $A:A'\dblcolon B:B'$ is a reasoning process that conveys \emph{$A$ is to $A'$ as $B$ is to $B'$}. This formulation has become a core technique for creating artistic 2D digital content, such as image analogies \cite{HertzmannJOCS01} in Photoshop \cite{photoshop} for image stylization and the Lit Sphere \cite{SloanMGG01} (a.k.a. MatCap) in ZBrush \cite{zbrush} for non-photorealistic renderings.
However, leveraging analogies to stylize 3D geometry is still at a preliminary stage because defining the analogy relationship on surface meshes requires dealing with irregular discretizations, curved metrics, and different topologies. 


In this paper, we introduce a step towards a more general 3D shape analogies, named \emph{spherical shape analogies}. We consider a specific case where $A$ is a unit sphere. This restriction enables us to operate on an input mesh $B$ with arbitrary topologies, boundaries, and geometric complexity. 
While not fully general, because $A$ is restricted to be a sphere, we demonstrate that this formulation can immediately achieve different geometric styles within a single framework. In \reffig{teaser}, we show that by providing different target style shapes $A'$ to the algorithm, we can turn the input shape $B$ into different styles.
In addition to stylization, our method can encompass many existing applications, such as developable surface approximation and PolyCube deformation. 

One key observation in our spherical shape analogies is that the surface normal is an effective instrument to capture geometric styles. Thus, we define the analogy relationship based on normals: we optimize a stylized shape $B'$ such that the relationship between the surface normals of $B$ and $B'$ is the same as the relationship between the surface normals of $A$ and $A'$

We realize this by casting it as a simple and effective normal-driven shape optimization problem which aims at deforming the input shape towards a set of desired normals. However, such an optimization problem is often difficult due to the nonlinearity of unit normals. We draw inspiration from previous works and apply a change of variables to accelerate the computation: instead of directly optimizing the vertex positions, we optimize a set of rotations that rotate the normals of the input mesh to the set of desired normals. 
Our simple formulation with the change of variables results in a generic stylization algorithm that runs at interactive rates.   





\section{Extensions \& Analysis}\label{sec:analysis} 
In this section, we introduce its extensions to different regularizations and how to handle cases where target normals $T(B')$ are a function of output geometry.   

\subsection{Different Regularizations}\label{sec:regularizations}
In addition to $E_{\textsc{arap}}\xspace$, the normal-driven optimization supports different regularization energies for different modeling intents. 
One could use {\textsc{arap}}\xspace when the goal is to produce a smooth deformation that preserves surface details.
If one wants to produce a non-smooth deformation (e.g., sharp creases) while preserving local rigidity, one could instead use a \emph{face-only} {\textsc{arap}}\xspace energy $E_{\textsc{farap}}\xspace$ \update{discussed in \cite{ZhaoG16a, LeviG15}} which consists of only the membrane term.
If one is interested in preserving the textures and allowing local scaling, one could use an \emph{as-conformal-as-possible} energy $E_{\textsc{acap}}\xspace$ \cite{BouazizDSWP12}.

\begin{wrapfigure}[4]{r}{0.53in}
	\raggedleft 
	\vspace{-16pt}

	\hspace*{-0.9\columnsep}
	\includegraphics[width=0.85in, trim={3mm 0mm 1mm 0mm}]{figures/face_only.jpg} 
	\label{fig:stretchElement} 
\end{wrapfigure} 
\textit{Face-only {\textsc{arap}}\xspace.}\ 
The core idea is to remove the bending term from {\textsc{arap}}\xspace and only measure the membrane term \cite{TerzopoulosPBF87}, so that two adjacent triangles can bend freely. We achieve this by applying the idea from \cite{ZhaoG16a, LeviG15} which only measures the {\textsc{arap}}\xspace energy over the three edge vectors of a face $f_k$ (see the inset), instead of the spokes and rims $\mathcal{N}_k$. Precisely, we can write this ``face-only'' {\textsc{arap}}\xspace regularization $E_{\textsc{farap}}\xspace$ as
\begin{align}\label{equ:farap}
  E_{\textsc{farap}}\xspace(\widetilde{{\matFont{V}}}, {\matFont{R}}) =  \sum_{k \in F} \sum_{i,j \in f_k} w_{ij} \| {\matFont{R}}_k {\vecFont{e}}_{ij} -  {{\vecFont{e}}'}_{ij} \|^2_2,
\end{align}

\textit{As-conformal-as-possible.}\ 
If the goal is to create novel geometric details, it is crucial to allow non-rigid deformations. However, an arbitrary deformation may lead to undesirable behaviors, such as badly shaped triangles. Thus constraining the angle preservation, a.k.a. conformality, will be a suitable regularization.
Specifically, we use the {\textsc{acap}}\xspace energy $E_{\textsc{acap}}\xspace$ in \cite{BouazizDSWP12} as our regularization 
\begin{align}\label{equ:acap}
  E_{\textsc{acap}}\xspace(\widetilde{{\matFont{V}}}, {\matFont{R}}, {\vecFont{s}}) =  \sum_{k \in F} \sum_{i,j \in \mathcal{N}_k} w_{ij} \| s_k {\matFont{R}}_k {\vecFont{e}}_{ij} -  {{\vecFont{e}}'}_{ij} \|^2_2,
\end{align}
where $s_k$ is a scalar representing the scaling of local patch. One can compute the optimal $s_k$ analytically via the method by Sch{\"o}nemann et al. \cite{schonemann1970fitting} (see \refapp{regularizations}).

Deploying these regularizations $E_{\textsc{farap}}\xspace, E_{\textsc{acap}}\xspace$ requires only a few changes in the optimization steps. Deploying $E_{\textsc{farap}}\xspace$ only involves changing the incidence matrix. Deploying $E_{\textsc{acap}}\xspace$ only requires adding one more line of code in the local step to solve an \emph{isotropic orthogonal Procrustes} problem \cite{schonemann1970fitting}. We detail such changes in \refapp{regularizations}. In \reffig{regularizations}, we apply the same deformation to a sheet but with different regularizations. We can perceive that different regularizations favor drastically different solutions.  

Our framework allows one to easily plug-and-play different regularizations. Specifically, we use $E_{\textsc{arap}}\xspace$ for applications that favor smooth deformation (e.g., \reffig{primStyles}), $E_{\textsc{farap}}\xspace$ for creating sharp creases (\reffig{generalized_PolyCube}, \reffignum{developable}), and $E_{\textsc{acap}}\xspace$ when one wants to manipulate geometric details such as in \reffig{synthesis}.
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/regularizations.jpg}
  \vspace{-5pt}
  \caption{Different regularizations favor different behaviors. Given a sheet (gray), we pull up the center part (central blue dots) and shrink the boundary (blue dots on the boundary), then we minimize each regularization energy to determine the unconstrained vertices. We can observe that $E_{\textsc{arap}}\xspace$ favors rigid and smooth interpolation, $E_{\textsc{farap}}\xspace$ favors sharp bending between triangles, and $E_{\textsc{acap}}\xspace$ favors to preserve angles while allowing local scaling. } 
  \label{fig:regularizations}
  \vspace{-10pt}
\end{figure}
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/update_t.jpg}
  \vspace{-5pt}
  \caption{Setting the target normal $T$ as a constant or treating it as a function of the output mesh $T(B')$ leads to different local minima. In many cases (left pair), both options lead to similar looking results. But setting $T$ as a constant may result in an undesirable local minimum in some cases (right), such as the ears of the bunny.} 
  \label{fig:update_t}
  \vspace{-10pt}
\end{figure} 

\subsection{Dynamic Target Normals}
Our method converges to a local minimum. 
Empirically, we observe that treating the target normal $T$ as a constant throughout the optimization may work fine perceptually in many cases (see the left pair in \reffig{update_t}). However, constant $T$ may lead to an undesirable local minimum due to a sub-optimal assignment of $T$ (see the right pair in \reffig{update_t}).
Inspired by Projective Dynamics \cite{BouazizMLKP14}, a simple solution to avoid such local minima is to treat $T$ as a function of $B'$ ($N_{B'}$ specifically), and update $T$ at every iteration.
We summarize the pseudo code in \refalg{normalDriven}. If $T$ is a constant throughout the optimization, one can simply skip the optional step at line 9.
\IncMargin{1em}
\begin{algorithm}[t]
    \SetKwInOut{Input}{Input}
    \SetKwInOut{Output}{Output}
    \caption{Normal-driven optimization}
    \label{alg:normalDriven}
    \Indentp{-1em}
        \Input{$\,$A triangle mesh ${\matFont{V}}, {\matFont{F}}$ and a weight $\lambda$}
        \Output{$\,$Deformed vertex positions $\widetilde{{\matFont{V}}}$}
        \BlankLine
    \Indentp{1em}
        compute $\widetilde{N}_{A'}$ \hfill \text{\color{iglGreen} // step 1, Sec. 3.1}\\
        compute $T$ from $\hat{n}({\matFont{V}})$ \hfill \text{\color{iglGreen} // step 2, Sec. 3.2}\\
        $\hat{{\vecFont{n}}} \leftarrow \hat{n}({\matFont{V}})$ \hfill \text{\color{iglGreen} // compute input surface normals}\\
        ${\matFont{Q}}, {\matFont{K}} \leftarrow precompute({\matFont{V}},{\matFont{F}})$ \hfill \text{\color{iglGreen} // see App. B}\\
        $\widetilde{{\matFont{V}}} \leftarrow {\matFont{V}}$\\
        \While{\textit{not converge}}{
          ${\matFont{R}} \leftarrow local\text{-}step\,(\widetilde{{\matFont{V}}}, \hat{{\vecFont{n}}}, T, \lambda)$  \hfill \text{\color{iglGreen} // App. A}\\
          $\widetilde{{\matFont{V}}} \leftarrow global\text{-}step\,({\matFont{R}}, {\matFont{Q}},{\matFont{K}})$ \hfill \text{\color{iglGreen} // App. B}\\
          compute $T$ from $\hat{n}(\widetilde{{\matFont{V}}})$ \hfill \text{\color{iglGreen} // (optional) for dynamic $T$}\\
        }
\end{algorithm} 
\DecMargin{1em}

\begin{wrapfigure}[7]{r}{1.33in}
	\raggedleft
 
  \vspace{-12pt}
 
	\hspace*{-0.7\columnsep}
	\includegraphics[width=1.55in, trim={6mm 0mm -1mm 0mm}]{figures/convergence.jpg}
	\label{fig:convergence}
\end{wrapfigure}
In terms of convergence, in the case where $T$ is constant, the convergence behaves the same as the original {\textsc{arap}}\xspace \cite{SorkineA07}, where the energy decreases monotonically. In the case where $T$ is dependent to $B'$, we do not guarantee a monotonic decrease in energy, but the optimization still converges in our experiments. In the inset, we visualize the convergence plot for examples in \reffig{update_t}. 

\begin{wrapfigure}[7]{r}{1.33in}
	\raggedleft
  \vspace{-12pt}
 
	\hspace*{-0.7\columnsep}
	\includegraphics[width=1.55in, trim={6mm 0mm -1mm 0mm}]{figures/runtime.jpg}
	\label{fig:runtime}
\end{wrapfigure}
We implement our algorithm in C++ with Eigen \cite{eigenweb} and evaluate our method on a MacBook Pro with an Intel i5 2.3GHz processor.
Our method runs 24 iterations per second for a mesh with around 20k vertices. We report a complete picture of our runtime in the inset. The local step will be the computation bottleneck for meshes with less than 20k vertices, but further acceleration can be achieved via the method by Zhang et al. \cite{ZhangFast2021}. 
Typically, within the first 10 iterations, our method can achieve a visually similar result compared to the converged solution. This property enables us to build an interactive tool for users to play with different style shapes $A'$ or artistic controls.


\section{Applications}
The major benefit of our analogy-based stylization method is that one can plug-and-play different style shapes to obtain different results. \update{When one provides} convex primitives with few distinct face normals, we can simply use the method discussed in \refsec{generateNS} to turn an input shape into the style of the primitive (see \reffig{primStyles}, \reffignum{primStyles_bd}). In \reffig{quan_normal}, we also quantitatively show that our method can effectively reduce the difference between mesh normals and the normals of a primitive. 
If the provided style shape is smooth or non-convex, where the simple closest normal may fail to capture the style, one could use a spherical parameterization described in \refsec{generateNS} to achieve the stylization.
If desiring more user controls, one could ``paint'' the desired surface normals on a unit sphere (see \refsec{generateNS}), and then transfer the style of the painted normals directly to the input (see \reffig{manual_gauss})

\subsection{PolyCube Deformation}
If one is interested in PolyCube maps \cite{TariniHCM04}, we can adapt normal driven editing to create PolyCube maps,following the observation in \cite{ZhaoLLZXG17}. Specifically, we need to use a cube as a style shape and move the pre-computation step in \refalg{normalDriven} to the optimization loop. 
Moving the pre-computation in the loop would no longer preserve the original details, which is desirable for creating PolyCube shapes. This modification may also lead to badly shaped triangle. 
When these faces appear, a quick solution is to move the vertex towards the 1-ring average by a small amount to improve triangle quality.
For the sake of comparison, we use the same PolyCube segmentation as in \cite{FuBL16} and show that we can achieve comparable results in \reffig{PolyCube_comparison}. 
We can further generalize the PolyCube map to other polygonal boxes by specifying non-cube normals (see \reffig{generalized_PolyCube}).

\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/generalized_PolyCube.jpg}
  \vspace{-5pt}
  \caption{\update{By using} different sets of normals, we can generalize the PolyCube method (left) to create polygonal boxed maps.} 
  \label{fig:generalized_PolyCube}
 
  \vspace{-10pt} 
\end{figure} 
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/control_developable.jpg}
  \vspace{-5pt}
  \caption{Stein et al. \cite{SteinGC18} control the patches on the developable surfaces via remeshing the input. We, instead, can control the amount of creases (middle, right) by tuning a single parameter (see \refapp{developable}).}
 
  \label{fig:control_developable}
 
  \vspace{-10pt} 
\end{figure}

\subsection{Developable Surface Approximation}
So far we have only considered an explicit shape or a set of painted normals as our style shape. Here we further extend our method to support an energy that describes a certain style. In particular, we consider the target normal $T$ is computed via an optimization
\begin{align}
  T = \argmin_T f(B'),
\end{align}
and, similar to the case where $T$ is dependent to $B'$, we update $T$ at every iteration in the local/global solve. 

\begin{wrapfigure}[7]{r}{1.33in}
	\raggedleft
    \vspace{-14pt}
	\hspace*{-0.7\columnsep}
	\includegraphics[width=1.48in, trim={6mm 0mm -1mm 0mm}]{figures/K_hist}
	\label{fig:K_hist}
\end{wrapfigure}
We evaluate this extension via setting $f$ to be the discrete developability energy proposed in \cite{SteinGC18}, with details provided in \refapp{developable}. Compared to the original method, our approach contains a regularization term in addition to the developable energy, thus our optimization requires no remeshing and results in the faster optimization (see \reffig{developable_comaprison}).  
In \reffig{control_developable}, we further show that our framework enables one to control the number of creases in the results.
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/developable.jpg}
  \vspace{-5pt}
  \caption{Our normal driven editing can be used to create many piece-wise developable surfaces (blue). Our method requires no remeshing and is fast enough for interactive modeling. \copyright cerberus333 (third) under CC BY-NC.} 
  \label{fig:developable}
  \vspace{-10pt}
\end{figure}
With our framework one can interactively create a variety of piece-wise developable shapes (see \reffig{developable}).
In \reffig{K_developable}, we evaluate our results by visualizing the discrete Gaussian curvature before and after running our developable flow. We can observe that the Gaussian curvature concentrates along the creases and results in a piece-wise developable surface. In the inset, we quantitatively demonstrate that our method effectively increases the developability of the mesh in \reffig{K_developable}.
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/K_developable.jpg}
  \vspace{-5pt}
  \caption{We use our normal driven editing to deform the input shape (gray) into a piece-wise developable approximation (blue). In the bottom row, we visualize the Gaussian curvature concentrates on the creases after the deformation, leading to a piece-wise developable shape. \copyright Oliver Laric under CC BY-NC-SA.} 
  \label{fig:K_developable}
  \vspace{-10pt}
\end{figure}

\section{Related Work}
Our work shares similar motivations to computer-assisted image stylization pioneered by Haeberli \cite{Haeberli90}. But since our outputs are stylized 3D geometries, we focus the discussion on geometric stylization and geometric deformation methods. 

\subsubsection*{Analogy-based Geometric Stylization}
Many generative models have been proposed for creating stylized 3D objects, such as collage art \cite{gal20073d, Theobalt2007Collage}, manga style \cite{ShenLHC12}, cubic style \cite{LiuCubic19}, and neuronal homunculus \cite{ReinertRS12}. However, these methods are tailor-made for only a specific style.

\emph{Analogy} $A:A'\dblcolon B:B'$ is a powerful idea to achieve different stylization results within a single framework. 
This idea has inspired several design tools for images \cite{HertzmannJOCS01}, non-photorealistic renderings \cite{SloanMGG01,Shechtman16}, and curves \cite{HertzmannOCS02}. 
\update{Beyond 2D data, the idea of analogy has also been used for transferring 3D geometric details from one shape to another. 
We omit the discussion on methods that are not based on analogies, such as mesh cloning \cite{ZhouHWTDGS06, TakayamaSSIBS11} and geometric learning \cite{LiuKCAJ20, HertzHGC20, WangAKCS20, chen2021decor, li_cvpr21}, and focus on analogy-based techniques.}
Ma et al. \cite{MaHSKW14} propose a method for 3D style transfer based on patch-based assembly. However, their method cannot handle free-form deformations and requires the source and the exemplar shape to share a similar structure in order to compute high-quality correspondences. 
\update{Bhat et al. \cite{BhatIT04} propose a voxel-based texture synthesis method for transferring geometric details encoded in the volumetric grid. Berkiten et al. \cite{BerkitenHSMLR17} use metric learning for details represented as displacement maps. These methods are designed for high-frequency details (e.g., wrinkles on the surface).} In contrast, our spherical shape analogies focuses on larger scale free-form deformations. 
Albeit limited --- in our analogies $A$ is restricted to the unit sphere --- our method enables a first step in this exciting direction.

\subsubsection*{Surface Normals in Shape Deformation}
A key insight of our spherical shape analogies is to leverage surface normals to capture geometric styles.
The surface normal is a fundamental geometric quantity and is ubiquitous in geometry processing.
A representative example is in the \emph{PolyCube} deformation \cite{TariniHCM04} where the goal is to optimize surface normals to be axis-aligned.
Gregson et al. \cite{GregsonSZ11} and Zhao et al. \cite{ZhaoLLZXG17} use the closest rotation from the surface normal to an axis-aligned direction to drive the PolyCube deformation. Huang et al. \cite{HuangJSTBD14} and Fu et al. \cite{FuBL16} propose to minimize energies defined on normals to create PolyCube shapes.
In architectural geometry design, surface normals are a main ingredient to characterize polygon meshes with planar faces. The methods proposed by Deng et al. \cite{DengPW11} and Poranne et al. \cite{PoranneOG13} utilize normals to formulate a \emph{distance-from-plane} constraint to encourage planarity. Tang et al. \cite{TangSGWP14} use the dot product between a face normal and its adjacent edge vectors to determine whether the vertices of a polygon are coplanar. 
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/silvia.jpg}
  \vspace{-5pt}
  \caption{Sell{\'{a}}n et al. \cite{Sellan2020} propose a technique to make 2D heightfields developable (purple). In contrast, our method can create developable approximations for surface meshes in 3D (blue). } 
  \label{fig:silvia}
  \vspace{-10pt}
\end{figure}
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/developable_comaprison.jpg}
  \vspace{-5pt}
  \caption{Compared to the method proposed by Stein et al. \cite{SteinGC18} for creating developable approximations (left), our method can create visually comparable results (right) with significant speed-ups.} 
  \label{fig:developable_comaprison}
  \vspace{-10pt}
\end{figure}
Characterizing whether a mesh can be flattened to 2D without stretching or shearing, a.k.a. \emph{developability}, also relies on surface normals. Stein et al. \cite{SteinGC18} characterize discrete developability based on the 1-ring face normals, and propose an algorithm to compute piecewise developable surfaces. Sell{\'{a}}n et al. \cite{Sellan2020} reformulate the developable energy into a convex semidefinite program for finding piecewise developable heightfields.
In addition to these examples, deforming shapes into the cubic style \cite{LiuCubic19, Fumero0R20}, constructing shape abstractions \cite{alexa2021polycover}, surface parameterization \cite{ZhaoSLZYLWGG20}, \update{and interactive mesh editing \cite{YuZXSBGS04, SorkineCLARS04}} are all related to surface normals. Many more examples can be found in the design of geometric filters, such as the Guided filter \cite{ZhangDZBL15}, the Shock filter \cite{PradaK15}, the Bilateral normal filter \cite{ZhengFAT11}, and the Total Variation mesh denoising \cite{ZhangWZD15}.

Our method can be adapted to these normal-based deformations. Compared to the PolyCube method \cite{FuBL16}, we achieve comparable quality (see \reffig{PolyCube_comparison}), but we can further generalize to polytopes (see \reffig{generalized_PolyCube}).
Compared to \cite{LiuCubic19} in cubic stylization (see \reffig{cubic_comparison}), we can achieve similar performance, but we can further generalize to many styles other than the cubic style (see \reffig{teaser}).
In developable surface approximation, in contrast to the method by Sell{\'{a}}n et al. \cite{Sellan2020}, our method can be applied to surface triangle meshes (see \reffig{silvia}) and is significantly faster than the method by Stein et al. \cite{SteinGC18} (see \reffig{developable_comaprison}).

\subsubsection*{Shape Deformation}
Our geometric stylization method can also be perceived as \update{a type of shape deformation method}. We share technical similarities with methods that deform a shape while addressing given modeling constraints.  
A common choice is to minimize the \emph{as-rigid-as-possible} ({\textsc{arap}}\xspace) energy \cite{SorkineA07,IgarashiMH05, ChaoPSS10} while satisfying the constraints. This {\textsc{arap}}\xspace energy measures the rigidity of local surface patches and favors detail-preserving smooth deformations.
In the case where locally rigid deformations are too constrained, the conformal energy \cite{CranePS11, VaxmanMW15} which preserves angles is commonly used.  In contrast to {\textsc{arap}}\xspace, the conformal energy often triggers larger deformations as it allows both local uniform scaling and rigid transformations. 
In addition to mesh deformations, similar energies have also been used for parameterization \cite{LiuZXGG08}, shape optimization \cite{BouazizDSWP12}, and simulating mass-spring systems \cite{LiuBOK13}.
The {\textsc{arap}}\xspace and conformal energies are also commonly used as regularization terms in mesh optimization problems, such as reconstruction \cite{IzadiF14}, surface registration \cite{HuangAWG08, YoshiyasuMYK14}, PolyCube construction \cite{HuangJSTBD14}, and surface stylization \cite{LiuCubic19}. 
Their popularity comes from the property that they favor smooth deformations and are amenable to fast optimizations.
For the same reasons, we also use these as our regularization energies for interactive modeling tasks (see \reffig{regularizations}). 




























\section{Limitations \& Future Work}\label{sec:limitations}
Our method draws inspiration from Projective Dynamics \cite{BouazizMLKP14} to handle the case where target normals $T$ are a function of output shape $B'$ (e.g., \reffig{update_t}, \reffignum{developable}). Although being fast and suitable for our intended interactive applications, it often struggles to converge to a highly accurate solution. 
Extending our optimization to, for example, Newton's method would be desirable for applications that desire highly accurate solutions.

Our approach is restricted to a sphere as our reference shape $A$, and uses the Gauss map to determine the correspondences between $A$ and the input $B$. As the Gauss map purely relies on surface normals to determine the map, the resulting map is ignorant to area distortion. 
This characteristic is beneficial to handle input shapes $B$ that are very different (e.g., different genus) from a sphere because in these cases it is challenging to obtain a map with low area distortion. 
However, the price we have to pay is that we cannot support structured 
\begin{wrapfigure}[8]{r}{0.62in}
	\raggedleft 


	\vspace{-6pt}
	\hspace*{-0.7\columnsep}
	\includegraphics[width=0.75in, trim={4mm 0mm 1mm 0mm}]{figures/printedHead_cut} 
\end{wrapfigure} 
and high-frequency patterns (e.g., geometric texture synthesis). Thus, if one is interested in stylizing shapes with detailed textures, we suggest to first synthesize target normals on the surface directly \cite{WieLKT09} then perform the normal-driven optimization (\refsec{normalDriven_arap}). 
In \reffig{synthesis} we demonstrate this alternative by unbaking an existing normal map for manufacturing purposes (see the inset) and synthesizing normal textures from an image.
\begin{figure}
  \centering
  \includegraphics[width=3.33in]{figures/synthesis.jpg}
  \vspace{-5pt}
  \caption{If one is interested in creating high-frequency geometric textures, we recommend to compute target normals via texture synthesis and then optimize the geometry via the normal-driven optimization. We demonstrate an example of ``unbaking'' normal maps (left) and an example of geometric texture synthesis (right).}
   
  \label{fig:synthesis}
  \vspace{-10pt}
\end{figure} 

Our method currently supports manifold triangle meshes. 
Extending to non-manifold meshes, polygon meshes, volumetric meshes, and point clouds could be beneficial to handle real-world geometric data. 
Not every shape or normal capture sphere is valid to serve as the style shape of our algorithm. Discovering the validity of a style shape is important to understand the behavior of these novel modeling methods.
Removing the assumption about the source shape being a sphere could lead to a more general analogy-based shape editing. 
Based on the observation that surface normals are a promising geometric quantity to capture the style of a shape. Developing a better categorization of styles based on normals or exploring learning-based techniques on normals (instead of vertices) could lead to novel stylization methods.

\section*{Acknowledgements}
  Our research is funded in part by NSERC Discovery (RGPIN2017–05235, RGPAS–2017–507938), New Frontiers of Research Fund (NFRFE–201), the Ontario Early Research Award program, the Canada Research Chairs Program, the Fields Centre for Quantitative Analysis and Modelling and gifts by Adobe Systems, Autodesk and MESH Inc. 
  %
  We thank Sheldon Andrews, Abhishek Madan, Silvia Sell{\'{a}}n, Oded Stein, Li-Yi Wei for helps on experiments. We thank members of Dynamic Graphics Project at the University of Toronto; Sarah Kushner, Abhishek Madan, Silvia Sell{\'{a}}n, Let{\'{i}}cia Mattos Da Silva, Towaki Takikawa for proofreading; John Hancock for the IT support. We thank all the artists for sharing a rich variety of 3D models.
  %
\section{Local Step with $E_{\textsc{arap}}\xspace$} \label{app:localStep}
Given a fixed $\widetilde{{\matFont{V}}}$, we obtain the optimal rotation for each vertex $k$ by solving the following minimization problem
\begin{align*}
  {\matFont{R}}_k = \argmin_{{\matFont{R}}_k \in\text{SO(3)}} \sum_{i,j \in \mathcal{N}_k} w_{ij} \| {\matFont{R}}_k {\vecFont{e}}_{ij} -  {{\vecFont{e}}'}_{ij} \|^2_2 + \lambda a_k\| {\matFont{R}}_k \hat{{\vecFont{n}}}_k - {\vecFont{t}}_k \|_2^2
\end{align*}
The above optimization is an instance of the \emph{orthogonal Procrustes} which finds the best rotation matrix ${\matFont{R}}_k$ to map a set of vectors (${\vecFont{e}}_{ij}, \hat{{\vecFont{n}}}_k$) to another set of vectors (${{\vecFont{e}}'}_{ij}, {\vecFont{t}}_k$). We can re-write it into a more compact expression as:
\begin{align}\label{equ:localStep}
  &{\matFont{R}}_k^\star = \argmax_{{\matFont{R}}_k \in \text{SO}(3)} \ \Tr({\matFont{R}}_k {\matFont{X}}_k) \\
  &{\matFont{X}}_k = 
  \begin{bmatrix}
    {\matFont{E}}_k & \hat{{\vecFont{n}}}_k
  \end{bmatrix}
  \begin{bmatrix}
    {\matFont{W}}_k & \vspace{3.5pt}\\
    & \lambda a_k
  \end{bmatrix}
  \begin{bmatrix}
    {{\matFont{E}}'}_k^\top  \vspace{2pt}\\
    {\vecFont{t}}_k^\top
  \end{bmatrix}.
\end{align}
where ${\matFont{W}}_k$ is a $|\mathcal{N}_k|$-by-$|\mathcal{N}_k|$ diagonal matrix of the cotangent weights $w_{ij}$, ${\matFont{E}}_k$ and ${{\matFont{E}}'}_k$ are $3$-by-$|\mathcal{N}_k|$ matrices concatenating the edge vectors of the face one-ring at the rest and deformed states, respectively. 
One can then derive the optimal ${\matFont{R}}_k$ from the SVD of ${\matFont{X}}_k = \mathcal{U}_k \sum_k \mathcal{V}_k^\top$
\begin{align}\label{equ:optRotation}
  {\matFont{R}}_k = \mathcal{V}_k \mathcal{U}_k^\top,
\end{align}
up to changing the sign of the column of $\mathcal{U}_k$ so that $\det({\matFont{R}}_k ) > 0$.

\section{Global Step with $E_{\textsc{arap}}\xspace$} \label{app:globalStep}
The global step updates the deformed vertex positions $\widetilde{{\matFont{V}}}$ from a fixed set of rotations ${\matFont{R}}$ obtained via the local step. This boils down to solving the following problem 
\begin{align*}
  \widetilde{{\matFont{V}}}^\star = \argmin_{\widetilde{{\matFont{V}}}} \sum_{k \in V} &\sum_{i,j \in \mathcal{N}_k} w_{ij} \| {\matFont{R}}_k {\vecFont{e}}_{ij} -  {{\vecFont{e}}'}_{ij} \|^2_2
\end{align*}
We can expand this energy as
\begin{align*}
  &\sum_{k \in V} \sum_{i,j \in \mathcal{N}_k} w_{ij} \| {\matFont{R}}_k {\vecFont{e}}_{ij} -  {{\vecFont{e}}'}_{ij} \|^2_2 \\
  &\quad = \sum_{k \in V} \sum_{i,j \in f_k} w_{ij} {{\vecFont{e}}'}_{ij}^\top {{\vecFont{e}}'}_{ij} - 2 w_{ij} {{\vecFont{e}}'}_{ij}^\top {\matFont{R}}_k {\vecFont{e}}_{ij} + \text{constant} 
\end{align*}
It is often convenient to express the summation in terms of matrices. 
We introduce a directed incidence matrix ${\matFont{A}}_k$ with size $|V|$-by-$|\mathcal{N}_k|$ to  represent the edge vectors in $\mathcal{N}_k$ as ${\matFont{V}}^\top {\matFont{A}}_k$, and we use ${\matFont{M}}_k$ to represent a $|\mathcal{N}_k|$-by-$|\mathcal{N}_k|$ diagonal matrix of the weights $w_{ij}$.  
Then we can re-write the energy in terms of matrices as
\begin{align}\label{equ:arapGlobalStep}
  &\sum_{k \in V} \Tr({\matFont{M}}_k {\matFont{A}}_k^\top \widetilde{{\matFont{V}}} \widetilde{{\matFont{V}}}^\top {\matFont{A}}_k) - 2 \Tr({\matFont{M}}_k {\matFont{A}}_k^\top \widetilde{{\matFont{V}}} {\matFont{R}}_k {\matFont{V}}^\top {\matFont{A}}_k ) \nonumber\\
  &= \sum_{k \in V} \Tr(\widetilde{{\matFont{V}}}^\top {\matFont{A}}_k{\matFont{M}}_k {\matFont{A}}_k^\top \widetilde{{\matFont{V}}}) - 2 \Tr({\matFont{R}}_k{\matFont{V}}^\top{\matFont{A}}_k{\matFont{M}}_k {\matFont{A}}_k^\top \widetilde{{\matFont{V}}}) \nonumber\\
  &= \Tr\Big(\widetilde{{\matFont{V}}}^\top \Big(\textstyle\sum_k {\matFont{A}}_k{\matFont{M}}_k {\matFont{A}}_k^\top \Big) \widetilde{{\matFont{V}}} \Big) -2 \Tr\Big(\Big(\textstyle\sum_k {\matFont{R}}_k{\matFont{V}}^\top{\matFont{A}}_k{\matFont{M}}_k {\matFont{A}}_k^\top\Big) \widetilde{{\matFont{V}}}\Big) \nonumber\\
  &= \Tr(\widetilde{{\matFont{V}}}^\top {\matFont{Q}} \widetilde{{\matFont{V}}}) -2 \Tr({\matFont{R}} {\matFont{K}} \widetilde{{\matFont{V}}}),
\end{align}

where ${\matFont{R}} = \{ {\matFont{R}}_k \}$ is the concatenation of all the rotations, ${\matFont{Q}}$ is a $|V|$-by-$|V|$ symmetric matrix, and ${\matFont{K}}$ is a $|9V|$-by-$|3V|$ matrix stacking the constant terms which can be computed during the precomputation. We can then find the optimal $\widetilde{{\matFont{V}}}$ by solving a linear system
\begin{align*}
  {\matFont{Q}} \widetilde{{\matFont{V}}} = {\matFont{K}}^\top {\matFont{R}}^\top
\end{align*}
As we know from \cite{SorkineA07}, ${\matFont{Q}}$ is the cotangent Laplacian \cite{PinkallP93}. We can pre-factorize ${\matFont{Q}}$ to speed up runtime performance. With these pieces in hand, we can minimize our energy \refequ{arapNormalDriven} by iteratively performing the local and the global steps (see \refalg{normalDriven}).

\section{Generalize to $E_{\textsc{farap}}\xspace$ and $E_{\textsc{acap}}\xspace$}\label{app:regularizations}
Changing the regularization from $E_{\textsc{arap}}\xspace$ to the membrane-only regularization $E_{\textsc{farap}}\xspace$ (\refequ{farap}) requires to re-define ${\matFont{R}}$ on each face and change the set of edge vectors to the three edge vectors of a triangle. These changes would lead us to replace the ${\matFont{E}}_k, {{\matFont{E}}'}_k$ in the local step \refequ{localStep} to the three edge vectors of a face, and $a_k$ to the face area. In the global step, one only needs to update the incidence matrix ${\matFont{A}}_k$ in \refequ{arapGlobalStep} to a $|V|$-by-$|f_k|$ matrix containing the three edge vectors information. 

Deploying the as-conformal-as-possible regularization $E_{\textsc{acap}}\xspace$ (\refequ{acap}) changes the local step to solve an instance of the isotropic orthogonal Procrustes problem, where an analytical solution has been derived in \cite{schonemann1970fitting}. In short, one can obtain the optimal rotation the same way as \refequ{optRotation}, and compute the optimal scaling $s_k$ analytically as
\begin{align*}
  s_k = \frac{\Tr({\matFont{W}}_k {{\matFont{E}}'}^\top_k {\matFont{R}}_k {\matFont{E}}_k) + \lambda a_k \hat{{\vecFont{n}}}_k^\top {\vecFont{t}}_k}{\Tr({\matFont{W}}_k {\matFont{E}}^\top_k {\matFont{E}}_k) + \lambda a_k \hat{{\vecFont{n}}}_k^\top {\vecFont{t}}_k}.
\end{align*}
When assembling the matrices for the global step, using $E_{\textsc{acap}}\xspace$ would require replacing ${\matFont{R}}_k$ with $s_k{\matFont{R}}_k$.

\section{Projective Dynamics for Dynamic Target Normals}\label{app:PD}
We draw inspiration from \emph{projective dynamics} \cite{BouazizMLKP14} to handle cases where the target normal $T$ is a function of output geometry $B'$. Let us first define 
\begin{align*} 
  T = \argmin E_N(\widetilde{{\matFont{V}}})
\end{align*}
as a minimizer of an energy $E_N$ defined on the output shape. In our cases, $E_N$ could be the distance to the closet normals or the developable energy \cite{SteinGC18}. With this definition, we re-write \refequ{arapNormalDriven} as
\begin{align*} 
  &\min_{\widetilde{{\matFont{V}}},{\matFont{R}}}\ \sum_{k \in V} \sum_{i,j \in\mathcal{N}_k} w_{ij} \| {\matFont{R}}_k {\vecFont{e}}_{ij} -  {{\vecFont{e}}'}_{ij} \|^2_2  + \lambda a_k \| {\matFont{R}}_k \hat{{\vecFont{n}}}_k  - {\vecFont{t}}_k \|_2^2,\\
  &\text{ subject to}\ T = \argmin E_N(\widetilde{{\matFont{V}}})
\end{align*}
This reformulation allows us to directly deploy the projective dynamics solver by first projecting $T = \{ {\vecFont{t}}_k\}$ to the ``constraint'' $E_N$, fixing ${\vecFont{t}}_k$, and solving the original problem as \refequ{arapNormalDriven} via the local/global solver to get $\widetilde{{\matFont{V}}}$ at the next iteration. We then iterate this procedure (see \refalg{normalDriven}) until convergence. 
This expression enables us to plug-and-play different $E_N$ for different modeling objectives.

\section{Normal Driven Developable Surfaces} \label{app:developable}
Our normal-driven editing can be used to create developable surfaces by specifying a set of target normals that are developable.
Stein et al. \cite{SteinGC18} propose a characterization of discrete developability based on face normals of a vertex one-ring.
In short, if all the one-ring face normals correspond to a common plane or two planes, then this local one-ring is piecewise developable. 

With this characterization, we can easily get a set of ``developable'' face normals by (1) visiting all the one-ring faces of a vertex, (2) performing a small principle component analysis on the face normals for each one-ring, and (3) projecting the normals to one or two common planes by zeroing out the components correspond to the smallest eigenvalues. 
By using a different threshold to decide whether to zero out the smallest or the smallest two components, we can control the amount of creases in the developable approximation (see \reffig{control_developable}).
As each face will receive three (possibly) different developable normals from the previous procedure, we simply average them to get the target face normals. 
We perform this developable normal computation at each iteration in parallel, which corresponds to the Line 9 of \refalg{normalDriven}.