# Developability Approximation for Neural Implicits through Rank Minimization

Pratheba Selvaraju
University of Massachusetts, Amherst
<EMAIL_ADDRESS>
###### Abstract

Developability refers to the process of creating a surface without any tearing
or shearing from a two-dimensional plane. It finds practical applications in
the fabrication industry. An essential characteristic of a developable 3D
surface is its zero Gaussian curvature, which means that either one or both of
the principal curvatures are zero. This paper introduces a method for
reconstructing an approximate developable surface from a neural implicit
surface. The central idea of our method involves incorporating a
regularization term that operates on the second-order derivatives of the
neural implicits, effectively promoting zero Gaussian curvature. Implicit
surfaces offer the advantage of smoother deformation with infinite resolution,
overcoming the high polygonal constraints of state-of-the-art methods using
discrete representations. We draw inspiration from the properties of surface
curvature and employ rank minimization techniques derived from compressed
sensing. Experimental results on both developable and non-developable
surfaces, including those affected by noise, validate the generalizability of
our method.

Figure 1: Reconstructed developable approximation of implicit surface. We
propose Implicit developable surface approximation using rank minimization
heuristics as regularizer on second order differentiable of the objective
function. We are able to achieve approximate piecewise developable surfaces
with automatic crease emergence. The figure.1 shows (a) the ground truth
surface, (b) implicit (signed distance field - SDF) reconstruction without
regularizer, and (c) - (f) implicit reconstruction with regularizers detailed
in 4.

##  1 Introduction

Developable surfaces have practical applications in digital fabrication,
industrial design, architecture, and geometric surface abstraction. They are
used in the automotive industry for car body panels, in industrial design for
curved furniture elements, in computer graphics to simplify complex 3D shapes
to enhance rendering performance, and in architecture for building facades to
reduce material wastage. Research efforts are focused on developing computer-
aided algorithms to identify these developable patches on a surface, aiming to
minimize material wastage in industry applications during cutting and
reassembling processes. Several approaches exist for approximating developable
shapes, employing optimization techniques [35, 32, 3, 5], shape wrapping
techniques that utilize multiple developable patches [13, 15, 12] and methods
for synthesizing developable surfaces [6]. Existing methods for reconstructing
developable shapes often rely on fixed topology representations. In contrast,
implicit surfaces offer advantages such as smooth interpolation, deformation
capabilities, and the ability to handle topological changes naturally. In
recent years, there has been extensive research in neural implicit surface
representations for 3D reconstruction, from single shape reconstruction [4],
to data-driven methods [14, 39], and generative modeling by deformation [42].
These approaches achieve high-detail surface reconstruction, and
regularization techniques have been incorporated to encourage smoothness [1].
However, there is currently no implicit reconstruction method that
specifically promotes surface developability. Our paper presents a novel
approach where we introduce a regularization term into neural implicit
formulations to deform the surface into approximately piecewise developable
patches.

Our approach is driven by two key observations. Firstly, implicit surfaces
offer the advantage of providing access to gradients and higher-order
derivatives, allowing us to compute surface normals, curvature, and other
surface properties without additional computations. In our method, we utilize
the second-order derivatives of the implicit function with respect to input
coordinates, which are related to Gaussian curvature [31]. Secondly, surface
developability necessitates zero Gaussian curvature [21]. To achieve this, we
formulate the developability condition for implicits as a rank minimization
problem, inspired by the work of Sellan et al. [32], who applied the concept
to height fields derived from depth maps. We employ both gaussian curvature
minimization and implicit hessian rank minimization, combining them in an
objective that encourages both developability, resulting in the elimination of
Gaussian curvature and shape fitting of the implicit to the input point cloud.

Our implicit-driven approach has key advantages over current discrete
representation methods, which struggle with higher polygon counts. Discrete
optimization-based methods suffer from deformation-induced topology changes
leading to inaccuracies surface representation, our approach allows for easier
topological deformation while maintaining shape approximation. Additionally,
our method eliminates the need for specialized solvers by relying on a single
regularizer weight value to control the level of developability.

In summary, our main contribution is a novel approach introducing
regularization term and optimization procedure promoting developability in
neural implicit surface reconstruction. Qualitative and quantitative
evaluations in section 5 demonstrate its effectiveness for complex topologies,
robustness to noise, and superior shape preservation compared to existing
methods.

##  2 Related Works

Our approach leverages implicit functions and employs second-order
optimization techniques based on rank minimization to reconstruct an input
point cloud into a piecewise developable surface. In this section, we provide
a brief overview of prior research on developable surface reconstructions, and
applications of rank minimization approximations.

### Developable surface Approximation.

Existing developable surface reconstruction methods often use discrete
representations and rely on constrained optimization or patch-based wrapping.
For instance, Stein et al. [35] achieve piece-wise developable surfaces
through mesh vertex optimization. However, these methods assume noise-free
inputs and may require manual tuning, potentially getting trapped in local
minima due to non-convexity. Binninger and Verhoeven et al. [3] use Gaussian
image thinning for developability with automatic crease development. However,
increasing cone angles leads to significant shape deformation, deviating from
the ground truth shape. Gavriil et al. [19] propose a similar approach using
Gauss thinning. Sellen et al. [32] achieve developability with nuclear norm
minimization but focus on planar height fields. Computation time in these
methods scales with polygon count in the mesh. Patch wrapping methods
approximate a surface by fitting developable patches onto the input surface
geometry. Verhoeven et al. [38] use planar quad strips aligned to principal
curvatures for curved parts, while Rabinovich et al. [30] use orthogonal quad
meshes to optimize for developable surfaces with constrained overlapping quad
patches. Peternell et al. [22] fit a developable surface by estimating tangent
planes through curve approximation from data points. Spline-based methods
offer another discrete representation for smooth reconstruction. Tang et al.
[6] use cubic spline developable patches projected onto the surface and merge
them iteratively with proximity constraints. Leopoldseder et al. [34] employ
interpolation of tangent planes with right circular cones for appropriate
rulings. Gavriil et al. [19] propose B-spline surface optimization using Gauss
thinning for developability. Kenneth et al. [18] define a 3D polyline boundary
and generate a smooth discrete developable surface that interpolates this
boundary. Additionally, Solomon et al. [16] introduce a flexible structure for
developable surface representation.

On the other hand, our proposed method uses implicit functions and second-
order optimization based on rank minimization to achieve piecewise developable
surfaces from point cloud inputs.

### Rank minimization.

Stein et al. [27] proposes minimizing the $L_{1}$ norm of second derivatives
to reconstruct piecewise planar surfaces. Sellan et al. [32] utilizes nuclear
norm minimization for reconstructing piecewise developable surfaces from input
depth maps, but their method is restricted to heightfields. Liu et al. [11]
employ rank minimization using the $L_{1}$ norm to reconstruct surfaces in a
cube-like style while preserving the original shape’s content. Besides $L_{1}$
norm minimization, log-determinant minimization is utilized as rank
approximation for distance matrices [24], subspace clustering in [7]. $L_{0}$
minimization is employed for basis selection in [9] and for image editing in
[20]. Additionally, Oh et al. [37] uses partial sum nuclear norm minimization
for PCA.

##  3 Background

### Implicit surface representation.

In an implicit surface representation, a surface is defined as the zero level
set of the implicit function. The implicit function takes as input the
coordinates of a point $\mathbf{p}\in\mathcal{R}^{3}$ in space and returns a
scalar value $s=f(\mathbf{p})$ where $s\in\mathcal{R}$. Iso-levels of the
implicit function represent surfaces in $\mathcal{R}^{3}$. Points on the
surface are those where the implicit function evaluates to zero
$f(\mathbf{p})=0$. Points with negative scalar values $f(\mathbf{p})<0$
correspond to the shape interior, while points with positive values
$f(\mathbf{p})>0$ represent the shape exterior.

### Gaussian curvature and developability.

A developable surface is defined by its Gaussian curvature, which plays a key
role in its characterization. Gaussian curvature measures the curvature of a
surface at each point and is determined by the _principal curvatures_ [21]. In
the case of a developable surface, the Gaussian curvature is zero everywhere.
This means that either or both the principal curvatures, $K_{1}$ and $K_{2}$,
is zero. As a result, developable surfaces can be unrolled or flattened onto a
plane without any distortion or stretching. The principal curvatures of a
surface represent the extreme curvatures in different directions. They are
determined by the maximum and minimum values of the _normal curvature_ at a
given point on the surface. Normal curvature is the curvature of curves formed
by the intersection of normal planes and the surface at that point. The
direction with the maximum principal curvature ($K_{1}$) corresponds to the
maximum normal curvature, while the direction with the minimum principal
curvature ($K_{2}$) corresponds to the minimum normal curvature. For a deeper
understanding of the differential geometry of surfaces and its relation to
geometry processing applications, refer to the tutorial by Crane et al. [17].

### Gaussian curvature of implicits.

The Gaussian curvature, denoted as $K$, of a surface represented by implicit
function $f(\mathbf{p})$ can be defined by the first and second-order
derivatives of the implicit function with respect to the surface coordinates
$\mathbf{p}\in\mathcal{R}^{3}$. The first order derivative, known as the
gradient $\nabla f(\mathbf{p})=(\frac{\partial f(\mathbf{p})}{\partial
x},\frac{\partial f(\mathbf{p})}{\partial y},\frac{\partial
f(\mathbf{p})}{\partial z})$, yields a 3D vector pointing to the direction of
the surface normal. The second order derivative, referred to as the Hessian
$H_{f}(\mathbb{p})$ provides information about the rate of change of the
surface normal in different directions, which in turn determines the
curvature.

$\mathbf{H}_{f}(\mathbf{p})=\begin{bmatrix}\frac{\partial^{2}f(\mathbb{p})}{\partial^{2}x}&\frac{\partial^{2}f(\mathbb{p})}{\partial
x\partial y}&\frac{\partial^{2}f(\mathbb{p})}{\partial x\partial z}\\\
\frac{\partial^{2}f(\mathbb{p})}{\partial y\partial
x}&\frac{\partial^{2}f(\mathbb{p})}{\partial^{2}y}&\frac{\partial^{2}f(\mathbf{p})}{\partial
y\partial z}\\\ \frac{\partial^{2}f(\mathbf{p})}{\partial z\partial
x}&\frac{\partial^{2}f(\mathbf{p})}{\partial z\partial
y}&\frac{\partial^{2}f(\mathbf{p})}{\partial^{2}z}\\\ \end{bmatrix}$ (1)

To compute the Gaussian curvature at a specific point $\mathbf{p}$ we employ
the following procedure as outlined in the work by Goldman [31]:

$K(\mathbb{p})\\!\\!=\\!\\!-\frac{det(\hat{\mathbf{H}}_{f}(\mathbf{p}))}{\left|{\nabla
f(\mathbf{p})}\right|^{4}},\mbox{where}\,\,\hat{\mathbf{H}}_{f}(\mathbf{p})\\!\\!=\\!\\!\begin{bmatrix}\mathbf{H}_{f}(\mathbf{p})\\!&\\!\nabla{f(\mathbf{p})^{T}}\\\
\nabla f(\mathbf{p})\\!&\\!\mathbf{0}\end{bmatrix},$ (2)

For a smooth surface, zero Gaussian curvature at a point $\mathbf{p}$ entails
that $K(\mathbf{p})=0\Leftrightarrow K_{1}\cdot K_{2}=0\Leftrightarrow
K_{1}=0\text{ or }K_{2}=0$. Since the Gaussian surface is defined only on
points where the gradient $f(\mathbf{p})$ is non-zero, having zero Gaussian
curvature means that the determinant of the $4\times 4$ matrix in the
numerator of Eq.2 must be zero:
$\text{det}(\hat{\mathbf{H}}_{f}(\mathbf{p}))=0$. The equation for the
determinant of the matrix $\hat{\mathbf{H}}_{f}(\mathbf{p})$ can be expressed
as follows [31]:

$det(\hat{\mathbf{H}}_{f}(\mathbf{p}))=-\nabla
f(\mathbf{p})\cdot\text{Cof}(H_{f}(\mathbf{p}))\cdot\nabla f(\mathbf{p})^{T}$
(3)

By utilizing the properties of the cofactor matrix,
$\text{Cof}(\textbf{H}_{f}(\mathbf{p}))^{T}=det(\textbf{H}_{f}(\mathbf{p}))\textbf{H}_{f}(\mathbf{p})^{-1}$,
we can minimize the rank of $\mathbf{H}_{f}(\mathbf{p})$ than
$\hat{\mathbf{H}}_{f}(\mathbf{p})$.

### Norm minimization.

Equating the determinant of an $n\times n$ matrix $\mathbf{X}$ to zero is
essentially the same as ensuring that the matrix is not full rank, i.e.,
$rank(\mathbf{X})<n$. Considering that the absolute value of the determinant
of a square matrix is equivalent to the product of its singular values, i.e.,
$|det(A)|=\prod_{i=1}^{n}\sigma_{i}$, the rank minimization can be framed as a
problem of minimizing the $L_{0}$ norm of the singular values. To prevent
trivial solutions, the minimization process is subject to additional
constraints, often expressed as a linear system
$\mathbf{A}\mathbf{X}=\mathbf{B}$ [23, 32]:

$\min\limits_{\mathbf{A}\mathbf{X}=\mathbf{B}}\,rank(\mathbf{X})\Leftrightarrow\min\limits_{\mathbf{A}\mathbf{X}=\mathbf{B}}\,\lVert\sigma(\mathbf{X})\rVert_{0}$
(4)

,where the vector $\sigma(\mathbf{X})$ stores the singular values of the
matrix $\mathbf{X}$. The $L_{0}$ minimization problem is non-convex, non-
differentiable, and generally non-tractable (NP-hard) [26]. However, when
dealing with the hessian rank of the implicits with rank of 3, we can make use
of a smoother approximation of the $L_{0}$ minimization objective similar to
[8, 41]. Instead of employing the $L_{0}$ cardinality approximation, we
consider a relaxed alternative in the form of the nuclear norm (also known as
the $L_{1}$ norm) of the singular values which has been demonstrated as a
tight convex approximation to the rank function [25]. The nuclear norm
$\lVert\mathbf{X}\rVert_{*}$ of the matrix:

$\min\limits_{\mathbf{A}\mathbf{X}=\mathbf{B}}\,\lVert\sigma(\mathbf{X})\rVert_{1}\Leftrightarrow\min\limits_{\mathbf{A}\mathbf{X}=\mathbf{B}}\,\lVert\mathbf{X}\rVert_{*}$
(5)

### Alternative approximations.

There is a significant drawback in solving for nuclear norm minimization.
While minimizing the $L_{1}$ norm provides an approximation to the rank
minimization problem and leads to a low-rank matrix $X$, it also
simultaneously diminishes the high-variance information that includes
important details like the structure of the object. Thus we also explored the
use non-convex surrogate partial sum minimization method [36], minimizing the
sum of smaller singular values:

$\min\limits_{\mathbf{A}\mathbf{X}=\mathbf{B}}\,rank(\mathbf{X})=\sum_{i=r+1}^{n}\sigma_{i}(\mathbf{X})$
(6)

,where $\sigma_{i}(\mathbf{X})$ refers to the $i^{th}$ singular value,
arranged in decreasing order, and the parameter $r$ determines the number of
largest singular values to be excluded during the minimization process. To
mitigate the impact of large singular values, an alternative non-convex
surrogate involving the log determinant function was also employed [7]. Since
the hessian matrix $\mathbf{H}_{f}(\mathbf{p})$ of the implicit function may
not always be positive semi-definite, we express the rank minimization in the
form:

$\log(\det(\mathbf{X}^{T}\mathbf{X}+\mathbf{I}))=\sum_{i=1}^{n}\log(1+\sigma_{i}^{2})$
(7)

Some other alternative relaxation methods include Weighted Nuclear Norm
Minimization, Capped $L_{1}$ Norm, Schatten-p Norm, Truncated Rank
Minimization, which we could not cover in our experiments. We recommend [2]
for brief review of these approaches.

##  4 Method

### Overview.

Our approach takes a point cloud
$\mathcal{P}=\\{\mathbf{p}_{i},\mathbf{n}_{i}\\}_{i=1}^{N}$, as input, where
$\mathbf{p}_{i}\in\mathbb{R}^{3}$ represents the 3D position of a point,
$\mathbf{n}_{i}\in\mathbb{R}^{3}$ is its corresponding normal, and $N$ is the
total number of points. The output of our method is an implicit function
$f(\mathbf{p})$, which assigns a scalar value to input points
$\mathbf{p}\in\mathbb{R}^{3}$. The reconstructed surface is obtained by
extracting the zero iso-level ($s=0$) of the implicit function using the
marching cubes algorithm [40]. Our goal is to obtain an implicit function that
approximates the input point cloud and generates a surface that maximizes its
developability. This implies that the resulting surface points should ideally
possess negligible or zero Gaussian curvature while still retaining the
overall shape of the point cloud.

To achieve this objective, we formulate the problem as an optimization task
aimed at estimating the parameters $\boldsymbol{\theta}$ of a neural network
function $f(\mathbf{p};\boldsymbol{\theta})$ that represents the implicit
function (Fig.2). This optimization involves minimizing a loss function
consisting of two components: (a) a data fitting term
$L_{\text{data}}(\boldsymbol{\theta})$, which encourages the zero iso-surface
of the implicit function to closely match the input point cloud, and (b) a
regularizer term $L_{\text{*}}(\boldsymbol{\theta})$ that encourages surface
developability in the output. In the following sections, we elaborate on these
two terms, describe the network architecture, and outline our optimization
procedure.

Figure 2:  Our MLP architecture maps a point $\mathbf{p}\in\mathcal{R}^{3}$ to
an implicit function $f(\mathbf{p};\theta)$ with learned parameters $\theta$.
The parameters are optimized through a loss function comprising data term and
regularization term detailed in Section.4. Backpropagation (dotted arrows),
computes gradients and Hessians with respect to $\mathbf{p}$ coordinates for
fine-tuning $\theta$.

### Data term.

Single surface reconstruction methods involve fitting neural networks [4] to a
single input point cloud using data losses. The samples $\mathbf{p}_{j}$ used
for fitting are obtained by directly selecting input points from the point
cloud that have reference implicit values $s_{j}=0$ (”on-surface point
samples”), and also by perturbing points along the normals
$\mathbf{p}_{j}=\mathbf{p_{i}}+\epsilon\mathbf{n_{i}}$ resulting in
$s_{j}=\epsilon$ (”off-surface point samples”). In recent data-driven
approaches, surface reconstruction is performed by estimating the parameters
of the implicit function from a large dataset of different point clouds.
DeepSDF [14] is one such example, utilizing an auto-decoder architecture for
this purpose. In our method, we follow a similar architecture.

Our surface reconstruction network takes a point $\mathbf{p}$ as input and
maps it to an implicit function $f(\mathbf{p};\boldsymbol{\theta})$. During
training, the model parameters are optimized through backpropagation using
point samples $\\{\mathbf{p}_{j},s_{j}\\}_{j=1}^{K}$ taken around the input
point cloud [14]. The training process initially optimizes the parameters
based on the data fitting term. After achieving the fitting, the
developability regularizer is introduced for fine-tuning, resulting in a
surface close to being developable. The process of fitting an implicit
function to an input point cloud involves penalizing the differences between
estimated and reference implicit values at various sample points surrounding
the input point cloud. Specifically, given $K$ point samples
$\\{\mathbf{p}_{j}\\}_{j=1}^{K}$ with associated scalar signed distance values
$\\{s_{j}\\}_{j=1}^{K}$, the data term can be formulate as the $L_{1}$ loss or
clamped $L_{1}$ loss [14], which aims to make the parameter estimation more
sensitive to details near the surface as follows:

$L_{\text{data}}(\boldsymbol{\theta})=\sum_{j=1}^{K}|f(\mathbf{p}_{j};\boldsymbol{\theta})-\text{cl}(s_{j},\delta)|,$
(8)

, where $\text{cl}(\cdot,\delta)=\text{min}(\delta,\text{max}(-\delta,\cdot))$
and $\delta$ is a clamping parameter. We set $\delta=0.01$ in our experiments.

### Network architecture.

Our architecture is comprised of $8$ fully connected layers, each with 512
nodes, applying group normalization. To enable the computation of Hessians for
second-order differentiable regularization, we explored several second-order
activation functions [siLU[33], geLU[33], tanH, sine [39], and elu].

### Developability regularizer term.

Our regularizer term is motivated by rank minimization applied to the matrix
$\hat{\mathbf{H}}$ storing the gradients and Hessian of the implicit function
in Eq.2. We experimented with all rank minimization formulations except
$L_{0}$ as discussed in 3:

$\displaystyle
L_{\text{H}_{NN}}(\mbox{\boldmath$\theta$})=\sum_{i=1}^{N}\,\lVert\sigma(\mathbf{H}_{f}(\mathbf{p}_{i}))\rVert_{1}$
(9) $\displaystyle
L_{{\hat{\text{H}}_{det}}}(\mbox{\boldmath$\theta$})=\sum_{i=1}^{N}\det(\hat{\mathbf{H}}_{f}(\mathbf{p}_{i})),\forall\mbox{rank}(\hat{\mathbf{H}}_{f}(\mathbf{p}_{i}))=3$
(10) $\displaystyle
L_{\text{H}_{logdet}}(\mbox{\boldmath$\theta$})=\sum_{i=1}^{N}\log\det(\mathbf{H}_{f}^{\top}(\mathbf{p}_{i})\cdot\mathbf{H}_{f}(\mathbf{p}_{i})+\mathbf{I})$
(11) $\displaystyle
L_{\text{H}_{PNN}}(\mbox{\boldmath$\theta$})=\sum_{i=1}^{N}\sum_{o=r+1}^{o=3}\sigma_{o}(\mathbf{H}_{f}(\mathbf{p}_{i}))$
(12)

In our experiments, we set $r=1$ (i.e., we ignore the two largest singular
values of the Hessian). In our experiments, we minimized the determinants of
both $\hat{\text{H}}$ and H, and it was observed that $\hat{\text{H}}$
exhibited better developability.

Figure 3:  Visualization of Discrete Gaussian curvature of results of Table 1. (a) SDF Ground truth(230K$|$57K Vertices), (b) Binninger,Verhoeven et.al [3] on (a), (c) marching cube [40] reconstructed surface using $\hat{\text{H}}_{det}$ regularizer described in 4. _Note: Marching cube reconstruction lacks smooth edges and quality relies on voxel resolution (512 used here)._ Method |  Median K $\downarrow$ | Mean K $\downarrow$ | CD $\downarrow$
---|---|---|---
GT $\left|57K\right|$ | 0.004 | 0.012 | 0.0
[3] $\left|57K\right|$ | $3e^{-5}$ | 0.006 | 115.1
[3] $\left|230K\right|$ | $1e^{-5}$ | 0.003 | 102.1
$\text{SDF+H}_{det}\left|57K\right|$ | $1e^{-3}$ | 0.02 | 25.5
$\text{SDF+H}_{det}\left|950K\right|$ | $5e^{-5}$ | 0.003 | 25.6
$\text{SDF+}\hat{\text{H}}_{det}\left|3.8M\right|$ | $\boldsymbol{{1e^{-5}}}$ | 0.001 | 27.1
Table 1: Comparison of Discrete Gaussian curvature (K) and Chamfer Distance
(CD) for bunny 3. Ground-truth is the SDF surface reconstructed using marching
cubes. Chamfer distance evaluated as sum of square difference of the vertices.
Figure 4:  Histogram of Implicit Gaussian curvature of stanford bunny 6
reconstructed using (a) SDF shape approximation without any regularizer, in
comparison to the SDF Developable shape approximation using regularizer
variants 4, from (b) - (e).

### Minimization procedure.

Our optimization procedure aims to minimize an objective function combining
the data term and developability regularizer term:

$L(\mbox{\boldmath$\theta$})=L_{\text{data}}(\mbox{\boldmath$\theta$})+\lambda\cdot
L_{*}(\mbox{\boldmath$\theta$})$ (13)

where $L_{*}$ denotes one of the regularizer terms detailed above(4), and
$\lambda$ is set through hyperparameter tuning. Steps to minimizing the
objective involves (1) achieving good shape approximation to the point cloud,
(2) regularizing it for developability, and (3) computing its derivative with
respect to the input point and network parameters. We solve the first problem
by minimizing the data term $L_{data}(\mbox{\boldmath$\theta$})$ for the input
point cloud, yielding an approximate iso-surface. We then minimize the full
loss function $L(\mbox{\boldmath$\theta$})$ evaluated on the model parameters.
Computing the implicit’s Hessians needed in the developability term is
feasible due to the twice-differentiable feed-forward activation function used
in our network.

We use subgradient optimization for $L_{1}$ norm minimization. As for other
minimization terms, we leveraged pytorch’s autograd (_torch.autograd.grad_) to
calculate the gradients and their derivatives required for backpropagation. By
minimizing the complete loss function with backpropagation, we are able to
obtain piecewise developable surfaces with reduced Gaussian curvature and
automatic crease formation. The iterations were run until convergence.

Figure 5:  Ablation on regularizer weights, $\lambda$, and developability of
surface. Row (a), shows comparison for $\hat{\textbf{H}}_{det}$ regularized
surface, row (b), shows results for $\textbf{H}_{\small{PNN}}$ regularized
surface.

### Implementation details.

Prior to conducting the experiments, we standardized the inputs by fitting
point cloud to unit bounding box. For shape optimization, we set the learning
rate to $10^{-4}$ and adjust it to $10^{-5}$ during the fine-tuning stage for
developability. We employed the Adam minimizer [10] for optimization purposes.
Our implementation is in Pytorch. The source code and data will be made
available upon acceptance of the paper.

Method | Bunny | Horse | Dragon | Griffin
---|---|---|---|---
Med $\downarrow$ | Med $\downarrow$ | CD $\downarrow$ | Med $\downarrow$ | Med $\downarrow$ | CD $\downarrow$ | Med $\downarrow$ | Med $\downarrow$ | CD $\downarrow$ | Med $\downarrow$ | Med $\downarrow$ | CD $\downarrow$
$\text{K}_{min}$ | K |  | $\text{K}_{min}$ | K |  | $\text{K}_{min}$ | K |  | $\text{K}_{min}$ | K |
SDF | 1.6 | 10.2 | 14.4 | 2.0 | 19.9 | 1.3 | 2.9 | 37.4 | 2.9 | 1.6 | 18.5 | 2.5
SDF[3] | 0.2 | 0.8 | 430.7 | 0.2 | 1.27 | 170.5 | 0.3 | 1.8 | 860.2 | 0.2 | 0.9 | 174.4
SDF[35]($\left|3K\right|$) | 0.5 | 1.3 | 54.2 | 2.0 | 18.2 | 1.6 | - | - | - | - | - | -
SDF+$\text{H}_{NN}$ | 0.3 | 0.8 | 71.6 | 1.4 | 9.4 | 5.0 | 1.7 | 12.7 | 75.8 | 0.8 | 5.2 | 19.4
SDF+$\text{H}_{logdet}$ | 0.7 | 1.7 | 62.9 | 1.2 | 8.8 | 3.9 | 1.4 | 12.1 | 5.4 | 0.6 | 3.7 | 4.0
SDF+$\hat{\text{H}}_{det}$ | 0.08 | 0.5 | 110.4 | 0.2 | 1.31 | 18.6 | 0.5 | 3.5 | 51.0 | 0.2 | 0.8 | 11.5
SDF+$\text{H}_{PNN}$ | 0.03 | 0.2 | 134.7 | 0.3 | 2.1 | 10.2 | 0.4 | 4.5 | 63.2 | 0.2 | 1.7 | 15.6

Table 2:  Comparison of Implicit Curvature and Chamfer distance. Implicit
Curvature(K) is calculated according to Eq.2, and $\text{K}_{min}$ is
calculated using Eq.16, and their median value is taken (Med). Chamfer
distance (CD) uses sum of the squared distance metric for 500K points.
$\downarrow$ means lower values are better. For comparison, we employ SDF
reconstruction from the developable surface results of other methods (SDF[3],
SDF[35]). _Note: Only the implicit curvature values are measured from the SDF
reconstruction, Chamfer distance is measured from their discrete
representation by trying to match or lower the_ Med($\text{K}_{min}$)

.

Model |  Reg | Med $\downarrow$ | Mean $\downarrow$ | CD $\downarrow$
---|---|---|---|---
$\lambda=0$ | $\lambda$ | K | K |
| 0 | 14.5 | 2879.2 | 18.9
0.1 | 8.1 | 1227.9 | 20.7
0.5 | 3.9 | 1025.3 | 25.8
1 | 2.6 | 625.1 | 29.7
10 | 0.7 | 158.9 | 64.4
100 | 0.3 | 68.2 | 166.9
Table 3:  Quantitative result: Ablation study on regularizer weight,
$\lambda^{\prime}s$ effect on surface developability and chamfer distance.
Evaluation of implicit Gaussian curvature (K) for different $\lambda$ weights
on Albert Einstein by iczfirzis licensed under CC BY-SA.(Fig.5). Median(Med)
and mean of K as well as Chamfer Distance (CD) to ground-truth, are measured
using 250k surface points after ICP [29]. The best performance is highlighted
in bold.  Figure 6:  Comparison of developable surface reconstruction results
of two of our regularization methods (c) and (d) from 4 with the results of
(b) Binninger,Verhoeven et al. [3]. The insets shows that (b) fails to
preserve the structure of the shape while our methods does. Row 1 shows the
rotated image of back of the dragon, Row 3 shows drapes of back side of Lucy,
Row 4 the wings

##  5 Experiments

We now discuss our experiments for evaluating the effects of the proposed
developability regularization, the evaluation metrics, and finally show
results and comparisons.

### Competing variants.

We evaluate different variants based on the same network architecture
mentioned in Section 4. The evaluated variants include:

(a) Activation variants: _geLU_ , _siLU_ , _tanh_ , _elu_ used in the network
for both shape approximation and fine-tuning for developability (ref.
_supplement_)

(b) Regularization variants: $\text{H}_{\small{NN}}$, $\text{H}_{logdet}$,
$\text{H}_{PNN}$, $\hat{\text{H}}_{det}$, which perform fine-tuning of the
network parameters based on the total loss function (Eq.13). $\text{H}_{NN}$,
$\text{H}_{logdet}$ and $\text{H}_{PNN}$ minimizes the $L_{1}$ norm of the
singular values, the logarithmic determinant of squared hessian and the lowest
singular value of the implicit hessian $\mathbf{H}_{f}(\mathbf{p})$
respectively, while $\hat{\text{H}}_{det}$ minimizes the determinant of the
matrix $\hat{\mathbf{H}}_{f}(\mathbf{p})$ (Eq. 3). For fair comparison, these
variants were evaluated with the geLU activation function.

### Evaluation metrics.

We evaluate the regularizer variants using two metrics. Firstly, we employ the
Chamfer distance metric to assess the similarity between the reconstructed
surface and the ground-truth surface. This involves sampling 250k points on
both surfaces (using Poisson disk uniform sampling) followed by iterative
closest point (ICP) [29] for alignment. Then bidirectional Chamfer distance is
calculated between the point samples. Secondly, to measure surface
developability, we compute the median of the absolute minimum of the implicit
principal curvatures ($K_{\text{min}}$) at each point. Measuring the minimum
is crucial in cases where either one of $K_{1}$ or $K_{2}$ is close to 0 while
the other possesses a high value (such as creases). This allows us to
accurately assess the overall curvature, avoiding underestimation of
developability at points where the surface is nearly flat. To compute the
principal curvatures, and we follow Goldman [31] to obtain the implicit mean
curvature $M$:

$M(\mathbb{p})\\!\\!=\\!\\!\frac{\nabla{f(\mathbf{p})}\cdot\mathbf{H}_{f}(\mathbf{p})\cdot\nabla{f(\mathbf{p})^{T}}-\left|{\nabla
f(\mathbf{p})}\right|^{2}\cdot\text{Tr}(\mathbf{H}_{f}(\mathbf{p}))}{2{\left|{\nabla
f(\mathbf{p})}\right|^{3}}}$ (14)

Using Eq.2 and Eq.14, we get the principal curvatures values and their
absolute minimum as follows [31] :

$\displaystyle
K_{1},K_{2}\\!\\!=\\!\\!M(\mathbf{p})\pm\sqrt{M^{2}(\mathbf{p})-K(\mathbf{p})}$
(15) $\displaystyle
K_{\text{min}}\\!\\!=\\!\\!\text{min}(\left|K_{1}\right|,\left|K_{2}\right|)$
(16)

### Results.

Tables 1, 2 and 3 present our quantitative evaluation of the regularizer
variants. We find that all variants incorporating our developability term
significantly reduce Gaussian curvature on the reconstructed surface compared
to the ground truth. It is worth noting that while approaching developability,
the shape approximation of the reconstruction does not deviate much(Fig.6). It
shown in the Chamfer distance metric compared to other methods (1, 2).
Additionally, we observe that each variant offers a different level of
developability approximation (Fig.1), owing to the varied approaches and
relaxations used to minimize the hessian rank, but the $\text{H}_{NN}$ and
$\hat{\text{H}}_{det}$ variants provide consistent piecewise developable
patches with lower $K_{\text{min}}$ and lower chamfer distance. The same
conclusion can be observed qualitatively (1,3). $\hat{\text{H}}_{det}$ gives
piecewise planar developable patches, $\text{H}_{NN}$ gives planar and non-
planar developable patches, $\text{H}_{logdet}$ minimizes overall curvature,
but produces stronger crease lines, while $\text{H}_{NN}$ tries to minimize
both the principal curvatures but has adverse effect on shape approximation.

### Robustness.

We introduce a perturbation of $1\%$ to the positions of the input point cloud
while preserving the normal sign. In Figure 7, (b) shows the reconstruction of
noisy point cloud using SDF reconstruction without any regularizer, while (c)
and (d) depict the reconstructions using the $\text{H}_{NN}$ and
$\hat{\text{H}}_{det}$ regularizer variants for the same noisy input. Despite
the noise, the reconstructed surfaces remain close to the ground truth
surface. Note that, as the noise $\%$ increases, the surface SDF
reconstruction gets thicker.

Figure 7:  Comparison of developbale surface reconstruction results with noisy
input point cloud ($1\%$ noise) for and $\text{H}_{PNN}$ and
$\hat{\text{H}}_{det}$ regularizer variants.

### Ablation Results.

Our findings reveal that increasing the regularizer weight enhances
developability, as evident from the lower median and mean curvatures shown in
Table 3. However, this improvement comes with the trade-off of a greater
deviation from the shape approximation.

##  6 Conclusion

We have developed a method to approximate the developability of surfaces,
applicable to closed surfaces with varying levels of detail. Our approach
leverages the implicit representation of the surface and introduces a novel
regularizer term that acts on the implicit’s hessian and gradient, encouraging
the emergence of a piecewise developable surface with automatic crease
formation. Experimental results demonstrate the results of developability
achieved by our method while preserving the surface’s structural
characteristics better than alternative techniques.

Figure 8: Failed case using our method. (b) shows the reconstruction using
$\text{H}_{PNN}$ regularizer variant using larger regularizer weight
$\lambda$. Thin connectivity structure (dog ears) gets disconnected or
vanishes.

### Limitations and future work.

While our method shows promising results, it works only on single closed
surfaces and assumes correct normal signs in the provided point cloud. With
the current architecture the natural next step is to generalize this as data-
driven approach while extending it to handle open surfaces and noisy normals.
Our method is a global shape optimization method, thus does not adequately
preserve details in regions in a shape with different levels of detailing.
(Fig.6, Lucy model, where higher weights makes drapes with more developable,
but face details are not preserved). Moreover, there are potential failures
arising from marching cube reconstruction, which may not ensure connectivity
for shapes with thinner connectivity structures, as illustrated in Figure 8
with higher regularizer weights. Thus, automatically segmenting point clouds
and applying appropriate regularizer for individual segments is a promising
avenue for future research.

## References

  * Amos et al. [2020] Gropp Amos, Yariv Lior, Haim Niv, Atzmon Matan, and Lipman Yaron.  Implicit geometric regularization for learning shapes.  In _Proceedings of the 37th International Conference on Machine Learning_. JMLR.org, 2020.
  * April and John E. [2021] Sagan April and Mitchell John E.  Low-rank factorization for rank minimization with nonconvex regularizers.  _Computational Optimization and Applications_ , 79(2):273–300, 2021.
  * Binninger et al. [2021] Alexandre Binninger, Floor Verhoeven, Philipp Herholz, and Olga Sorkine-Hornung.  Developable approximation via gauss image thinning.  _Computer Graphics Forum (proceedings of SGP 2021)_ , 40(5):289–300, 2021.
  * Carr et al. [2001] J. C. Carr, R. K. Beatson, J. B. Cherrie, T. J. Mitchell, W. R. Fright, B. C. McCallum, and T. R. Evans.  Reconstruction and representation of 3d objects with radial basis functions.  In _Proceedings of the 28th Annual Conference on Computer Graphics and Interactive Techniques_ , page 67–76. Association for Computing Machinery, 2001.
  * Charlie and Kai [2004] Wang Charlie and Tang Kai.  Achieving developability of a polygonal surface by minimum deformation: A study of global and local optimization approaches.  _Visual Computer_ , 20:521–539, 2004.
  * Chengcheng et al. [2016] Tang Chengcheng, Bo Pengbo, Wallner Johannes, and Pottmann Helmut.  Interactive design of developable surfaces.  _ACM Trans. Graph._ , 35(2), 2016.
  * Chong et al. [2015] Peng Chong, Kang Zhao, Li Huiqing, and Cheng Qiang.  Subspace clustering using log-determinant rank approximation.  In _Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , page 925–934, New York, NY, USA, 2015\. Association for Computing Machinery.
  * Christos et al. [2018] Louizos Christos, Welling Max, and P. Kingma Diederik.  Learning sparse neural networks through l_0 regularization.  In _International Conference on Learning Representations_ , 2018.
  * David and Bhaskar [2004] Wipf David and Rao Bhaskar.  L0-norm minimization for basis selection.  In _Proceedings of the 17th International Conference on Neural Information Processing Systems_ , page 1513–1520, Cambridge, MA, USA, 2004. MIT Press.
  * Diederik P. and Jimmy [2015] Kingma Diederik P. and Ba Jimmy.  Adam: A method for stochastic optimization.  In _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , 2015.
  * Hsueh-Ti Derek and Alec [2019] Liu Hsueh-Ti Derek and Jacobson Alec.  Cubic stylization.  _ACM Trans. Graph._ , 38(6), 2019.
  * Idan et al. [2006] Shatz Idan, Tal Ayellet, and Leifman George.  Paper craft models from meshes.  _Vis. Comput._ , 22(9):825–834, 2006.
  * Ion et al. [2020] Alexandra Ion, Michael Rabinovich, Philipp Herholz, and Olga Sorkine-Hornung.  Shape approximation by developable wrapping.  _ACM Transactions on Graphics (proceedings of SIGGRAPH ASIA)_ , 39(6), 2020.
  * Jeong Joon et al. [2019] Park Jeong Joon, Florence Peter, Straub Julian, Newcombe Richard, and Lovegrove Steven.  Deepsdf: Learning continuous signed distance functions for shape representation.  In _The IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2019.
  * Jun and Hiromasa [2004] Mitani Jun and Suzuki Hiromasa.  Making papercraft toys from meshes using strip-based approximate unfolding.  _ACM Trans. Graph._ , 23(3):259–263, 2004.
  * Justin et al. [2012] Solomon Justin, Vouga Etienne, Wardetzky Max, and Grinspun Eitan.  Flexible Developable Surfaces.  _Computer Graphics Forum_ , 2012.
  * Keenan et al. [2013] Crane Keenan, de Goes Fernando, Desbrun Mathieu, and Schröder Peter.  Digital geometry processing with discrete exterior calculus.  In _ACM SIGGRAPH 2013 Courses_. Association for Computing Machinery, 2013.
  * Kenneth et al. [2007] Rose Kenneth, Sheffer Alla, Wither Jamie, Cani Marie-Paule, and Thibert Boris.  Developable surfaces from arbitrary sketched boundaries.  In _Proceedings of the Fifth Eurographics Symposium on Geometry Processing_ , page 163–172. Eurographics Association, 2007.
  * Konstantinos et al. [2019] Gavriil Konstantinos, Schiftner Alexander, and Pottmann Helmut.  Optimizing b-spline surfaces for developability and paneling architectural freeform surfaces.  _Computer-Aided Design_ , 111:29–43, 2019.
  * Li et al. [2011] Xu Li, Lu Cewu, Xu Yi, and Jia Jiaya.  Image smoothing via l0 gradient minimization.  _ACM Trans. Graph._ , 30(6):1–12, 2011.
  * Manfredo P. [1976] do Carmo Manfredo P.  _Differential Geometry of Curves and Surfaces_.  Prentice Hall, 1976.
  * Martin [2004] Peternell Martin.  Developable surface fitting to point clouds.  _Computer Aided Geometric Design_ , 21(8):785–803, 2004.
  * Maryam et al. [2001] Fazel Maryam, Hindi Haitham A., and Boyd Stephen P.  A rank minimization heuristic with application to minimum order system approximation.  In _Proceedings of the 2001 American Control Conference. (Cat. No.01CH37148)_ , pages 4734–4739 vol.6, 2001.
  * Maryam et al. [2003] Fazel Maryam, Hindi Haitham A., and Boyd Stephen P.  Log-det heuristic for matrix rank minimization with applications to hankel and euclidean distance matrices.  In _Proceedings of the 2003 American Control Conference, 2003._ , pages 2156–2162 vol.3, 2003.
  * Maryan [2002] Fazel Maryan.  _Matrix Rank Minimization with Applications_.  PhD thesis, Department of Electrical Engineering, Stanford University, 2002.
  * Natarajan [1995] B. K. Natarajan.  Sparse approximate solutions to linear systems.  _SIAM Journal on Computing_ , 24(2):227–234, 1995\.
  * Oded et al. [2018] Stein Oded, Grinspun Eitan, Wardetzky Max, and Jacobson Alec.  Natural boundary conditions for smoothing in geometry processing.  _ACM Trans. Graph._ , 37(2), 2018.
  * Paolo et al. [2008] Cignoni Paolo, Callieri Marco, Corsini Massimiliano, Dellepiane Matteo, Ganovelli Fabio, and Ranzuglia Guido.  MeshLab: an Open-Source Mesh Processing Tool.  In _Eurographics Italian Chapter Conference_. The Eurographics Association, 2008.
  * P.J. and Neil D. [1992] Besl P.J. and McKay Neil D.  A method for registration of 3-d shapes.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 14(2):239–256, 1992.
  * Rabinovich et al. [2018] Michael Rabinovich, Tim Hoffmann, and Olga Sorkine-Hornung.  Discrete geodesic nets for modeling developable surfaces.  _ACM Transactions on Graphics_ , 37(2), 2018.
  * Ron [2005] Goldman Ron.  Curvature formulas for implicit curves and surfaces.  22:632–658, 2005.
  * Silvia et al. [2020] Sellán Silvia, Aigerman Noam, and Jacobson Alec.  Developability of heightfields via rank minimization.  _ACM Transactions on Graphics_ , 2020.
  * Stefan et al. [2017] Elfwing Stefan, Uchibe Eiji, and Doya Kenji.  Sigmoid-weighted linear units for neural network function approximation in reinforcement learning.  _Neural networks : the official journal of the International Neural Network Society_ , 107:3–11, 2017.
  * Stefan and Helmut [1998] Leopoldseder Stefan and Pottmann Helmut.  Approximation of developable surfaces with cone spline surfaces.  _Computer-Aided Design_ , 30:571–582, 1998.
  * Stein et al. [2018] Oded Stein, Eitan Grinspun, and Keenan Crane.  Developability of triangle meshes.  _ACM Transaction on Graphics (TOG)_ , 37, 2018.
  * Tae-Hyun et al. [2016a] Oh Tae-Hyun, Tai Yu-Wing, Bazin Jean-Charles, Kim Hyeongwoo, and Kweon In So.  Partial sum minimization of singular values in robust pca: Algorithm and applications.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 38(4):744–758, 2016a.
  * Tae-Hyun et al. [2016b] Oh Tae-Hyun, Tai Yu-Wing, Bazin Jean-Charles, Kim Hyeongwoo, and Kweon In So.  Partial sum minimization of singular values in robust pca: Algorithm and applications.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 38(4):744–758, 2016b.
  * Verhoeven et al. [2022] Floor Verhoeven, Amir Vaxman, Tim Hoffmann, and Olga Sorkine-Hornung.  Dev2pq: Planar quadrilateral strip remeshing of developable surfaces.  _ACM Transactions on Graphics (TOG)_ , 41(3):29:1–18, 2022.
  * Vincent et al. [2020] Sitzmann Vincent, Martel Julien N.P., Bergman Alexander W., Lindell David B., and Wetzstein Gordon.  Implicit neural representations with periodic activation functions.  In _Proc. NeurIPS_ , 2020.
  * William E. and Harvey E. [1987] Lorensen William E. and Cline Harvey E.  Marching cubes: A high resolution 3d surface construction algorithm.  In _Proceedings of the 14th Annual Conference on Computer Graphics and Interactive Techniques_ , page 163–169. Association for Computing Machinery, 1987.
  * Yoshua et al. [2013] Bengio Yoshua, Léonard Nicholas, and C. Courville Aaron.  Estimating or propagating gradients through stochastic neurons for conditional computation.  _CoRR_ , abs/1308.3432, 2013.
  * Zhiqin and Hao [2019] Chen Zhiqin and Zhang Hao.  Learning implicit fields for generative shape modeling.  In _2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 5932–5941. IEEE Computer Society, 2019.

Supplementary Material

##  7 Implicit Curvature comparison

In Table.2, we compared the implicit curvature metric between the state-of-
the-art developable surface methods which operates on discrete representation
and our method. To facilitate this evaluation, we transformed the discrete
output to implicit representation by sampling 250k points from their surface
and fitted it using a SDF. To validate this evaluation approach, the following
table 4 showcases the calculated chamfer distances between the sampled points
from the discrete surfaces and those from the SDF-reconstructed surfaces.

Model |  Method |  CD $\downarrow$
---|---|---
($\text{n}_{samples}$ = 250k)
Bunny |  SDF[3] | 3.5
SDF[35] | 4.4
Horse |  SDF[3] | 1.3
SDF[35] | 2.5
Dragon |  SDF[3] | 3.8
Griffin |  SDF[3] | 2.2
Table 4:  Chamfer distance(CD) comparison of Sellan et al.[32] and
Binninger,Verhoeven et al.[3] with its corresponding SDF reconstruction. Grid
resolution of 512 is used for for marching cubes [40] reconstruction.

Note that Sellan et al. [32] required a reduced polygon count for compilation
and execution. Consequently, it was only tested on Bunny and Horse models,
achieved by reducing polygon numbers through meshlab [28]. Dragon and Griffin
models were excluded from this method, since reduction of polygons led to loss
of their surface details. The lower chamfer distance in Table 4 indicates that
the comparison of implicit curvature of our method against SOTA methods in
Table 2 demonstrates minimal deviation from the ground truth discrete
representation. The quantitative outcome is also illustrated in Figure 10.

Figure 9:  SDF reconstruction with different activation functions. (b) - (f)
shows the SDF reconstruction of stanford bunny using second order
differentiable activation functions under same constraints of epoch and
learning rate(Ref. 5).  Figure 10:  Comparison of discrete vs SDF
reconstructed developable surfaces of [3] and [32]. (a) shows the Ground
truth, (b) shows developable surface constructed using Binninger,Verhoeven et
al. [3] with maximum of 50k vertices, (d) shows developable surface
constructed using Sellen et al. [32] with maximum of $\approx$3k vertices, (c)
and (e) shows corresponding SDF reconstructed surface from (b) and (d) as
ground truth.

##  8 Activation comparison

We explored several second-order activation functions for SDF reconstruction.
The objective was to achieve both improved reconstruction accuracy in
comparison to the ground truth and faster convergence. During experimentation,
a learning rate of $1e^{-4}$ was used, and the training process was executed
until convergence. The $gelu$ activation function exhibited faster convergence
and better alignment with the ground truth, as assessed through the chamfer
distance metric. These results are depicted in Figure 9.

##  9 Developable surface results comparison

Figure.11 shows additional comparison of developable surface constructed by
Binninger,Verhoeven et al.[3] and the implicit reconstruction from our method
using all the variants described in Section.4. Figure.12 shows the developable
surface results varying the weight of the regularizer for $\text{H}_{PNN}$ and
$\hat{\text{H}}_{det}$ variant from our method.

Figure 11:  Comparison of developable surface reconstruction results of our
regularization methods (c) - (f) (Ref.4) with the results of (b)
Binninger,Verhoeven et al. [3] Figure 12:  Ablation on regularizer weights,
$\lambda$, and developability of surface. Row (a) and (c), shows comparison
for $\hat{\text{H}}_{det}$ regularized surface, row (b) and (d), shows results
for $\text{H}_{\small{PNN}}$ regularized surface. As the weight increases, the
surfaces gain more developability, which is a desired characteristic. However,
this improvement comes at the expense of preserving the original shape
structure.


# Learning Sensorimotor Primitives of Sequential Manipulation Tasks from
Visual Demonstrations

Junchi Liang, Bowen Wen, Kostas Bekris and Abdeslam Boularias This work was
supported by NSF awards IIS 1734492 and IIS 1846043. The authors are with the
Department of Computer Science, Rutgers University, NJ, USA. Emails:
{junchi.liang, bowen.wen, kostas.bekris<EMAIL_ADDRESS>
###### Abstract

This work aims to learn how to perform complex robot manipulation tasks that
are composed of several, consecutively executed low-level sub-tasks, given as
input a few visual demonstrations of the tasks performed by a person. The sub-
tasks consist of moving the robot’s end-effector until it reaches a sub-goal
region in the task space, performing an action, and triggering the next sub-
task when a pre-condition is met. Most prior work in this domain has been
concerned with learning only low-level tasks, such as hitting a ball or
reaching an object and grasping it. This paper describes a new neural network-
based framework for learning simultaneously low-level policies as well as
high-level policies, such as deciding which object to pick next or where to
place it relative to other objects in the scene. A key feature of the proposed
approach is that the policies are learned directly from raw videos of task
demonstrations, without any manual annotation or post-processing of the data.
Empirical results on object manipulation tasks with a robotic arm show that
the proposed network can efficiently learn from real visual demonstrations to
perform the tasks, and outperforms popular imitation learning algorithms.

##  I Introduction

Complex manipulation tasks are performed by combining low-level sensorimotor
primitives, such as grasping, pushing and simple arm movements, with high-
level reasoning skills, such as deciding which object to grasp next and where
to place it. While low-level sensorimotor primitives have been extensively
studied in robotics, learning how to perform high-level task planning is
relatively less explored.

High-level reasoning consists of appropriately chaining low-level skills, such
as picking and placing. It determines when the goal of a low-level skill has
been reached, and the pre-conditions for switching to the next skill are
satisfied. This work proposes a unified framework for learning both low and
high level skills in an end-to-end manner from visual demonstrations of tasks
performed by people. The focus is on tasks that require manipulating several
objects in a sequence. Examples include stacking objects to form a structure,
as in Fig. 1, removing lug nuts from a tire to replace it, and dipping a brush
into a bucket before pressing it on a surface for painting. These tasks are
considered in the experimental section of this work. For all of these tasks,
the pre-conditions of low-level skills depend on the types of objects as well
as their spatial poses relative to each other, in addition to the history of
executed actions. To support the networks responsible for the control
policies, this work uses a separate vision neural network to recognize the
objects and to track their 6D poses both over the demonstration videos as well
as during execution. The output of the vision network is the semantic category
of each object and its 6D pose relative to other objects. This output along
with the history of executed actions is passed to a high-level reasoning
neural network, which selects a pair of two objects that an intermediate level
policy needs to focus its attention on.

Figure 1: System overview: (top left) Video demonstrations of sequential
manipulation tasks performed by a person. (bottom left) The manipulated
objects are tracked in the input video to automatically identify sub-tasks.
(center) Given the tracking, policy networks are trained to perform high-level
reasoning and compute low-level controls. (right) The output of the policy
networks is forwarded to a robot, which manipulates the same objects as in the
demonstrations without other manual annotation. It is important to note that
the objects are always randomly placed on the table at the beginning of each
demonstration and each test scenario.

The first object is referred to as the tool, and to the second one as the
target. In a stacking task, the tool is the object grasped by the robot and
the target is a table or any other object on top of which the tool will be
placed. In the painting task, the tool is the brush and the target is the
paint bucket or the canvas. If no object is grasped, then the tool is the
robot’s end-effector and the target is the next object that needs to be
grasped or manipulated. An intermediate-level network receives the pair of
objects indicated by the high-level reasoning network, their 6D poses relative
to each other, and a history of executed actions. The intermediate-level
network returns a sub-goal state, defined as a way-point in $SE(3)$. Finally,
a low-level neural network generates the end-effector’s motion to reach the
way-point. The policy neural networks are summarized in Fig. 2.

While the proposed formulation is not exhaustive, it allows to cast a large
range of manipulation tasks, and to use the same network to learn them. The
proposed architecture requires only raw RGB-D videos, without the need to
segment them into sub-tasks, or even to indicate the number of sub-tasks. The
efficacy of the method is demonstrated in extensive experiments using real
objects in visual demonstrations, as well as both simulation and a real robot
for execution.

##  II Related Work

Most of the existing techniques in imitation learning in robotics are related
to learning basic low-level sensorimotor primitives, such as grasping, pushing
and simple arm movements [1, 2]. The problem of learning spatial preconditions
of manipulation skills has been addressed in some prior works [3, 4]. Random
forests were used [3] to classify configurations of pairs of objects, obtained
from real images, into a family of low level skills. However, the method
presented in [3] considers only static images where the objects are in the
vicinity of each other [3], in contrast to the proposed model, which
continuously predicts low-level skills while the objects are being manipulated
and moved by the robot. Moreover, it does not consider complex tasks that are
composed of several low-level motor primitives [3].

A closely related line of work models each sub-task as a funnel in the state
space, wherein the input and output regions of each sub-task are modeled as a
multi-modal Gaussian mixture [4, 5], and learned from explanatory data through
an elaborate clustering process. Explicit segmentation and clustering have
also been used [6]. Compared to these methods, the proposed approach is simple
to reproduce and uses significantly less hyper-parameters since it does not
involve any clustering process. Our approach trains an LSTM to select and
remember pertinent past actions. The proposed approach also aims for data-
efficiency through an attention mechanism provided by the high-level network.
Hierarchical imitation learning with high and low level policies is
investigated in recent work [7, 8]. These methods require ground-truth
labeling of each sub-task to train the high-level policy, while the proposed
method is unsupervised.

Skill chaining was considered in other domains, such as 2D robot navigation
[9]. Long-horizon manipulation tasks have also been solved by using symbolic
representations via Task and Motion Planning (TAMP) [10, 11, 12, 13].
Nevertheless, all the variables of the reward function in these works are
assumed to be known and fully observable, in contrast to the proposed
approach. A finite-state machine that supports the specification of reward
functions was presented and used to accelerate reinforcement learning of
structured policies [14]. In contrast to the proposed method, the structure of
the reward machine was assumed to be known. A similar idea has also been
adopted in other efforts [15, 16].

While 6D poses and labels of objects are provided from a vision module [17] in
the proposed approach, other recent works have shown that complex tasks can be
completed by learning directly from pixels [18, 19, 20, 21, 22, 23, 24]. This
objective is typically accomplished by using compositional policy structures
that are learned by imitation [18, 19], or that are manually specified [20,
21]. Some of these methods have been used for simulated control tasks [25, 26,
27]. These promising end-to-end techniques still require orders of magnitude
more training trajectories compared to methods like the one proposed here,
which separates the object tracking and policy learning problems.

Figure 2: Left: robot performing a stacking task. Right: policy networks. A
high-level policy network uses the latest observation, current progress
towards the current way-point and memory to select a pair of tool/target
objects and their relative poses, which are then passed along to the
intermediate policy for generating a next way-point. The low-level policy is
responsible of generating the actual motion of the end-effector toward the
way-point.

##  III Problem Formulation and Architecture

This approach employs a hierarchical neural network for learning to perform
tasks that require consecutive manipulation of multiple objects. The
assumption is that each scene contains at most $n$ objects from a predefined
set $\mathcal{O}=\\{o^{1},o^{2},\dots,o^{n}\\}$. The robot’s end-effector is
included as a special object in $\mathcal{O}$. The robot receives as inputs at
each time-step $t$ sensory data as an observation $z_{t}=(e_{t},\langle
l_{t}^{1},\dots,l_{t}^{m}\rangle,p_{t})$, where
$e_{t}\in\mathbb{R}^{3}\times\mathbb{SO}(3)$ is the 6D pose of the end-
effector in the world frame, $m$ is the maximum number of objects present in
the scene, $l_{t}^{i}$ is the semantic label of object $o^{i}$, and $p_{t}$ is
a $7n\times 7n$ matrix that contains the 6D poses of all objects relative to
each other, i.e., $p_{t}[o^{i},o^{j}]$ is a $7$-dim. vector that represents
$o^{i}$’s orientation and translation in the frame of object $o^{j}$. The
objects have known geometric models and have fixed frames of reference defined
by their centers and 3 principal axes. The objects are detected and tracked
using the technique presented in Section V-B. The maximum number of objects
$n$ is fixed a priori.

The system returns at each time-step $t$ an action
$a_{t}\in\mathbb{R}^{3}\times\mathbb{SO}(3)$, i.e., a desired change in the
pose $e_{t}$ of the robotic end-effector. An individual low-level sub-task is
identified by a tool denoted by $o_{t}^{+}$ and a target denoted by
$o_{t}^{*}$, along with a way-point $w_{t}$. The tool is the object being
grasped by the robot at time $t$, the target is the object to manipulate using
the grasped tool and the predicted way-point is the desired pose of the tool
in the target’s frame at the end of the sub-task. The way-point $w_{t}$ is a
function of time as it changes based on the current pose of the tool relative
to the target. Several way-points are often necessary to perform even simple
tasks. For instance, in painting, a brush is the tool $o_{t}^{+}$ and a paint
bucket is the target $o_{t}^{*}$. To load a brush with paint, several way-
points in the bucket’s frame need to be predicted. The first way-point can be
when the brush touches the paint, while the second way-point is slightly above
the paint. The tool $o_{t}^{+}$ and target assignment $o_{t}^{*}$ are also
functions of time $t$, and change as the system switches from one sub-task to
the next, based on the current observation $z_{t}$ and on what has been
accomplished so far. For instance, after loading the brush, the robot switches
to the next sub-task wherein the brush is still the tool object, but the
painting canvas or surface becomes the new target object.

In the proposed model, observations are limited to $6D$ poses of objects and
their semantic labels. These observations are often insufficient by themselves
for determining the current stage of the task, for deciding to terminate the
current sub-task and for selecting the next sub-task. For instance, in the
painting example, the vision module does not provide information regarding the
current status of the brush. Therefore, the robot needs to remember whether it
has already dipped the brush in the paint. Since it is not practical to keep
the entire sequence of past actions in memory, the approach uses a Long Short-
Term Memory (LSTM) to compress the history $h_{t}$ of the actions that the
robot has performed so far, and use it as an input to the system along with
observation $z_{t}$. The LSTM is trained along with the other parameters of
the neural network.

The following describes the three levels of the hierarchical network
architecture as depicted in Figure 2.

A high-level policy, denoted by $\pi_{h}$, returns a probability distribution
over pairs $(o^{+}_{t},o^{*}_{t})$ of objects, wherein $o^{+}_{t}$ is the
predicted tool at time $t$, and $o^{*}_{t}$ is the predicted target at time
$t$. The high-level policy takes as inputs 6D poses of objects and their
semantic labels, along with the pose of the robot’s end-effector.
Additionally, the high-level policy $\pi_{h}$ receives as inputs a history
$h_{t}=(o^{+}_{k},o^{*}_{k},w_{k})_{k=0,\dots,t-1}$ of pairs of tools and
targets $(o^{+}_{k},o^{*}_{k})$ and way-points $w_{k}$ at different times $k$
in the past, compressed into an LSTM unit, as well as a progress vector
$g_{t}=p_{t}[o^{+}_{t-1},o^{*}_{t-1}]-w_{t-1}$ that indicates how far is the
tool $o^{+}_{t-1}$ from the previous desired way-point $w_{t-1}$ with respect
to the target $o^{*}_{t-1}$.

An intermediate-level policy, denoted by $\pi_{m}$, receives as inputs the
current tool $o^{+}_{t}$ and target $o^{*}_{t}$, the pose
$p_{t}[o^{+}_{t},o^{*}_{t}]$ of $o^{+}_{t}$ relative to $o^{*}_{t}$, in
addition to history $h_{t}$ and progress vector $g_{t}$. Both tool $o^{+}_{t}$
and target $o^{*}_{t}$ are predicted by the high-level policy $\pi_{h}$, as
explained above. The intermediate-level policy returns a way-point
$w_{t}\in\mathbb{R}^{3}\times\mathbb{SO}(3)$, expressed in the coordinates
system of the target object $o^{*}_{t}$.

A low-level policy, denoted by $\pi_{l}$, receives as inputs the current pose
$p_{t}[o^{+}_{t},o^{*}_{t}]$ of the current tool $o^{+}_{t}$ relative to the
current target $o^{*}_{t}$, in addition to the way-point $w_{t}$ predicted by
the intermediate policy, and returns a Gaussian distribution on action
$a_{t}\in\mathbb{R}^{3}\times\mathbb{SO}(3)$ that corresponds to a desired
change in the pose $e_{t}$ of the robotic end-effector.

##  IV Learning Approach

In the proposed framework, an RGB-D camera is used to record a human
performing an object manipulation task multiple times with varying initial
placements of the objects. The pose estimation and tracking technique,
explained in Section V-B, is then used to extract several trajectories of the
form $\tau=(z_{1},z_{2},\dots,z_{H})$, wherein $z_{t}=(e_{t},\langle
l_{t}^{1},\dots,l_{t}^{n}\rangle,p_{t})$ is the observed $6D$ poses of all
objects at time $t$, including the end-effector’s pose $e_{t}$. The goal of
the learning process is to learn parameters of the three policy neural
networks $\pi_{h},\pi_{m}$ and $\pi_{l}$ that maximize the likelihood of the
data $\tau$ and the inferred way-points, tools and targets, so that the system
can generalize to novel placements of the objects that did not occur in the
demonstrations. The likelihood is given by:

$\displaystyle
P\big{(}(z_{t},w_{t},o^{+}_{t},o^{*}_{t})_{t=0:H}\big{)}=\Pi_{t=1}^{H-1}P_{A,t}P_{B,t}\text{,
with}$
$\begin{cases}P_{A,t}&\overset{\Delta}{=}P(w_{t},o^{+}_{t},o^{*}_{t}|z_{t})\\\
&=\pi_{h}(o^{+}_{t},o^{*}_{t}|z_{t},h_{t},g_{t})\pi_{m}(w_{t}|h_{t},g_{t},p_{t}[o^{+}_{t},o^{*}_{t}]),\\\
P_{B,t}&\overset{\Delta}{=}P(z_{t+1}|z_{t},w_{t},o^{+}_{t},o^{*}_{t})\\\
&=\pi_{l}(e_{t+1}-e_{t}|w_{t},o^{+}_{t},o^{*}_{t},p_{t}[o^{+}_{t},o^{*}_{t}]),\end{cases}\vspace{-.05in}$

wherein $h_{t}=(w_{t},o^{+}_{t},o^{*}_{t})_{i=0}^{t-1}$ is the history and
$g_{t}=p_{t}[o^{+}_{t-1},o^{*}_{t-1}]-w_{t-1}$ is the progress vector.

The principal challenge here lies in the fact that the sequence
$(w_{t},o^{+}_{t},o^{*}_{t})_{t=0}^{H}$ of way-points, tools, and targets is
unknown, since the proposed approach uses as inputs only 6D poses of objects
at different time-steps and does not require any sort of manual annotation of
the data.

To address this problem, an iterative learning process performed in three
steps is proposed. First, the low-level policy is initialized by training on
basic reaching tasks. The intermediate and high-level policies are initialized
with prior distributions that simply encourage time continuity and proximity
of way-points to target objects. Then, an expectation-maximization (EM)
algorithm is devised to infer the most likely sequence
$(w_{t},o^{+}_{t},o^{*}_{t})_{t=0}^{H}$ of way-points, tools and targets in
the demonstration data $(z_{t})_{t=0}^{H}$. Finally, the three policy networks
are trained by maximizing the likelihood of the demonstration data
$(z_{t})_{t=0}^{H}$ and the pseudo ground-truth data
$(w_{t},o^{+}_{t},o^{*}_{t})_{t=0}^{H}$ obtained from the EM algorithm. This
process is repeated until the inferred pseudo ground-truth data
$(w_{t},o^{+}_{t},o^{*}_{t})_{t=0}^{H}$ become constant across iterations.

###  IV-A Prior Initialization

This section first explains how the low-level policy $\pi_{l}$ is initialized.
The most basic low-level skill is moving the end-effector between two points
in $\mathbb{R}^{3}\times\mathbb{SO}(3)$ that are relatively close to each
other. We therefore initialize the low-level policy by training the policy
network, using gradient-ascent, to maximize the likelihood of straight-line
movements between consecutive poses $e_{t+1}$ and $e_{t}$ of the end-effector
while aiming at way-points
$\hat{w}_{t}\overset{\Delta}{=}p_{t+1}[o^{+}_{t},o^{*}_{t}]$. Therefore, the
objective of the initialization process is given as
$\max_{\theta_{l}}\sum_{t=1}^{H}\pi_{l}(e_{t+1}-e_{t}|o^{+}_{t},o^{*}_{t},p_{t}[o^{+}_{t},o^{*}_{t}],\hat{w}_{t})$,
wherein each $\hat{w}_{t}$ is expressed in the frame of the target
$o^{*}_{t}$, and $\theta_{l}$ are the parameters of the neural network
$\pi_{l}$. Both $o^{+}_{t}$ and $o^{*}_{t}$ are also chosen randomly in this
initialization phase. The goal is to learn simple reaching skills, which will
be refined and adapted in the learning steps to produce more complex motions,
such as rotations.

The intermediate policy $\pi_{m}$ is responsible for selecting way-point
$w_{t}$ given history $h_{t}$. It is initialized by constructing a discrete
probability distribution over points
$(\hat{w}_{t},\hat{w}_{t+1},\hat{w}_{t+2},\dots,\hat{w}_{H})$, defined as
$\hat{w}_{t}\overset{\Delta}{=}p_{t+1}[o^{+}_{t},o^{*}_{t}]$. Poses $p_{t+1}$
used as way-points $\hat{w}_{t}$ are obtained directly from demonstration data
$(z_{t})_{t=0}^{H}$. Specifically, we set
$\pi_{m}(\hat{w}_{k}|h_{t},o^{+}_{t},o^{*}_{t},p_{t}[o^{+}_{t},o^{*}_{t}])=0$
for $k<t$,
$\pi_{m}(\hat{w}_{k}|h_{t},o^{+}_{t},o^{*}_{t},p_{t}[o^{+}_{t},o^{*}_{t}])\propto\exp{(-\alpha\|\hat{w}_{k}\|_{2})}$
for $k=t$, and
$\pi_{m}(\hat{w}_{k}|h_{t},o^{+}_{t},o^{*}_{t},p_{t}[o^{+}_{t},o^{*}_{t}])\propto\exp{(-\alpha\|\hat{w}_{k}\|_{2})}\frac{1-\beta}{H-t}$
for $k>t$, where $\alpha$ and $\beta$ are predefined fixed hyper-parameters,
and $\hat{w}_{k}$ is expressed in the coordinates system of the target
$o^{*}_{t}$. This distribution encourages way-points to be close to the target
at time $t$. This distribution is constructed for each candidate target
$o^{*}_{t}\in\mathcal{O}$ at each time-step $t$, except for the robot’s end-
effector, which cannot be a target.

High-level policy $\pi_{h}$ is responsible for selecting tools and targets
$(o^{+}_{t},o^{*}_{t})$ as a function of context. It is initialized by setting
the tool as the object with the most motion relative to others:
$o^{+}_{t}=\arg\max_{o^{i}\in\mathcal{O}\setminus\\{o^{e}\\}}\sum_{o^{j}\in\mathcal{O}}\|p_{t+1}[o^{i},o^{j}]-p_{t}[o^{i},o^{j}]\|$,
excluding the end-effector (or human hand) $o^{e}$. If all the objects besides
the end-effector are stationary relative to each other, then no object is
being used, and the end-effector is selected as the tool. Once the tool
$o^{+}_{t}$ is fixed, the prior distribution on the target $o^{+}_{t}$ is set
as: $\pi_{h}(o^{+}_{t},o^{*}_{t}|h_{t},p_{t},g_{t})=0$ if $o^{+}_{t}=o^{e}$
(the end-effector cannot be a target),
$\pi_{h}(o^{+}_{t},o^{*}_{t}|h_{t},p_{t},g_{t})=\gamma$ if
$o^{+}_{t}=o^{+}_{t-1}$, and
$\pi_{h}(o^{+}_{t},o^{*}_{t}|h_{t},p_{t},g_{t})=\frac{1-\gamma}{n-2}$ if
$o^{+}_{t}\neq o^{+}_{t-1}$, where $o^{+}_{t-1}$ is obtained from history
$h_{t}$, $n$ is the number of objects and $\gamma$ is a fixed hyper-parameter,
set to a value close to $1$ to ensure that switching between targets does not
occur frequently in a given trajectory.

###  IV-B Pseudo Ground-Truth Inference

After initializing $\pi_{h},\pi_{m}$ and $\pi_{l}$ as in Section IV-A, the
next step consists of inferring from the demonstrations $(z_{t})_{t=1}^{H}$ a
sequence $(o^{+}_{t},o^{*}_{t},w_{t})_{t=1}^{H}$ of tools, targets and way-
points that has the highest joint probability
$P\big{(}(z_{t},w_{t},o^{+}_{t},o^{*}_{t})_{t=0:H}\big{)}$ (Algorithm 1, lines
2-15). This problem is solved by using the Viterbi technique. In a forward
pass (lines 2-12), the method computes the probability of the most likely
sequence up to time $t-1$ that results in a choice
$(o^{+}_{t},o^{*}_{t},w_{t})$ at time $t$. The log of this probability,
denoted by $F_{t}[(o^{+}_{t},o^{*}_{t},w_{t})]$, is computed by taking the
product of three probabilities: (i) $\pi_{h}$: the probability of switching
from $o^{+}_{t-1}$ and $o^{*}_{t-1}$ as tools and targets to $o^{+}_{t}$ and
$o^{*}_{t}$, given the progress vector $g_{t}$ and the object poses relative
to each other provided by the matrix $p_{t}$ (which is obtained from
observation $z_{t}$); (ii) $\pi_{m}$: the probability of selecting as a way-
point a future pose $p_{k}[o^{+}_{t},o^{*}_{t}]$ (denoted as $w_{k}$, $k\geq
t$) for the tool relative to the target in the demonstration trajectory; this
probability is also conditioned on choices made at the previous time step
$t-1$; (iii) the likelihood of the observed movement of the objects at time
$t$ in the demonstration, given the choice $(o^{+}_{t},o^{*}_{t},w_{k})$ and
the relative poses of the objects with respect to each other (given by matrix
$p_{t}$). For each candidate $(o^{+}_{t},o^{*}_{t},w_{k})$ at time $t$, we
keep in $R_{t}$ the trace of the candidate at time $t-1$ that maximizes their
joint probability. The backward pass (lines 13-15) finds the most likely
sequence $(o^{+}_{t},o^{*}_{t},w_{t})_{t=1}^{H}$ by starting from the end of
the demonstration and following the trace of that sequence in $R_{t}$. The
last step is to train $\pi_{m}$, $\pi_{l}$ and $\pi_{h}$ using the most likely
sequence $(o^{+}_{t},o^{*}_{t},w_{t})_{t=1}^{H}$ as a pseudo ground-truth for
the tools, targets and way-points.

Input: A set of $n$ objects $\mathcal{O}=\\{o^{1},o^{2},\dots,o^{n}\\}$; one
or several demonstration trajectories $\\{z_{t}\\}_{t=1}^{H}$, wherein
$z_{t}=(e_{t},\langle l_{t}^{1},\dots,l_{t}^{n}\rangle,p_{t})$, $e_{t}$ is the
end-effector’s pose at time $t$, $p_{t}[o^{i},o^{j}]$ is the 6D pose of
$o^{i}\in\mathcal{O}$ relative to $o^{j}\in\mathcal{O}$,
$\forall(o^{i},o^{j})\in\mathcal{O}\times\mathcal{O}$;

Output: High-level, intermediate-level, and low-level policies $\pi_{h}$,
$\pi_{m}$ and $\pi_{l}$;

1 Initialize $\pi_{h}$, $\pi_{m}$ and $\pi_{l}$ (Section IV-A);
$F_{0}[:]\leftarrow-\infty$;

2 for  _$t=1;t\leq H;t\leftarrow t+1$ _ do

3 foreach
_$(o^{+}_{t},o^{*}_{t},k)\in\mathcal{O}\times\mathcal{O}\times\\{t,{t+1},\dots,H\\}$_
do

4 $x_{t}\overset{\Delta}{=}(o^{+}_{t},o^{*}_{t},k)$;
$w_{k}\overset{\Delta}{=}p_{k}[o^{+}_{t},o^{*}_{t}]$;

5 $\Delta
p_{t}\overset{\Delta}{=}p_{t+1}[o^{+}_{t},o^{*}_{t}]-p_{t}[o^{+}_{t},o^{*}_{t}]$;

6 foreach
_$(o^{+}_{t-1},o^{*}_{t-1},k^{\prime})\in\mathcal{O}\times\mathcal{O}\times\\{t-1,\dots,H\\}$_
do

7 $x_{t-1}\overset{\Delta}{=}(o^{+}_{t-1},o^{*}_{t-1},k^{\prime})$;
$w_{k^{\prime}}\overset{\Delta}{=}p_{k^{\prime}}[o^{+}_{t-1},o^{*}_{t-1}]$;

8 $h_{t}\leftarrow(o^{+}_{t-1},o^{*}_{t-1},w_{k^{\prime}})$;

9 $g_{t}\leftarrow p_{t}[o^{+}_{t},o^{*}_{t}]-w_{t}$;

10
$Q[x_{t-1},x_{t}]\leftarrow\log\big{(}\pi_{h}(o^{+}_{t},o^{*}_{t}|h_{t},p_{t},g_{t})\big{)}+\log\big{(}\pi_{m}(w_{k}|h_{t},o^{+}_{t},o^{*}_{t},p_{t}[o^{+}_{t},o^{*}_{t}])\big{)}+\log\big{(}\pi_{l}(\Delta
p_{t}[o^{+}_{t},o^{*}_{t}]|o^{+}_{t},o^{*}_{t},p_{t}[o^{+}_{t},o^{*}_{t}],w_{k})\big{)}+F_{t-1}[x_{t-1}]$;

11 $F_{t}[x_{t}]\leftarrow\max_{x_{t-1}}Q[x_{t-1},x_{t}]$;

12 $R_{t}[x_{t}]\leftarrow\arg\max_{x_{t-1}}Q[x_{t-1},x_{t}]$;

13

14

/* Construct the most likely sequence */

15 $(o^{+}_{H},o^{*}_{H},k)\leftarrow\operatorname*{argmax}_{x}R_{t}[x]$;
$w_{H}\leftarrow p_{k}[o^{+}_{H},o^{*}_{H}]$;

16 for  _$t=H-1;t >0;t\leftarrow t-1$_ do

17 $(o^{+}_{t},o^{*}_{t},k)\leftarrow\operatorname*{argmax}_{x}R_{t+1}[x]$;
$w_{t}\leftarrow p_{k}[o^{+}_{t},o^{*}_{t}]$;

18

19Train the policy networks $\pi_{h}$, $\pi_{m}$ and $\pi_{l}$ with
$(o^{+}_{t},o^{*}_{t},w_{t})_{t=1}^{H}$ and $\\{z_{t}\\}_{t=1}^{H}$;

Optional: Go to 1 and repeat with updated policies $\pi_{h}$, $\pi_{m}$,
$\pi_{l}$;

Algorithm 1 Learning Policies from Visual Demonstrations

###  IV-C Training the Policy Networks

To train $\pi_{m}$, $\pi_{l}$ and $\pi_{h}$ using the pseudo ground-truth
$(o^{+}_{t},o^{*}_{t},w_{t})_{t=1}^{H}$, obtained as explained in the previous
section, we apply the stochastic gradient-descent technique to simultaneously
optimize the parameters of the networks by minimizing a loss function $L$
defined as follows. $L$ is defined as the sum of multiple terms. The first two
are $L_{o^{+}}=\sum_{t}CE(\hat{o^{+}}_{t},o^{+}_{t})$ and
$L_{o^{*}}=\sum_{t}CE(\hat{o^{*}}_{t},o^{*}_{t})$ where $CE$ is the cross
entropy, and $(\hat{o^{+}}_{t},\hat{o^{*}}_{t})$ is the current prediction of
$\pi_{h}$. The third term is $L_{w}=\sum_{t}MSE(\hat{w}_{t},w_{t})$ where
$\hat{w}_{t}$ is the output of $\pi_{m}$ and $MSE$ is the mean square error.
The next term is
$L_{action}=-log[\pi_{l}(a_{t}|p_{t}[o^{+}_{t},o^{*}_{t}],w_{t},o^{+}_{t},o^{*}_{t})]$,
which corresponds to the log-likelihood of the low-level actions in the
demonstrations. To further facilitate the training, two auxiliary losses are
introduced. The first one is to encourage consistency within each sub-task. As
the role of the memory in the architecture is to indicate the sequence sub-
tasks that have been already performed, it should not change before
$(o^{+}_{t},o^{*}_{t},w_{t})$ changes. If we denote the LSTM’s output as
$M(h_{t})$ at time-step $t$, the consistency loss is defined as
$L_{mem}=\sum_{t}\|M(h_{t})-M(h_{t-1})\|_{2}\times
I_{o^{+}_{t-1}=o^{+}_{t-1},o^{*}_{t-1}=o^{*}_{t},w_{t-1}=w_{t}}$ where $I$ is
the indicator function. The last loss term, $L_{ret}$, is used to ensure that
memory $M(h_{t})$ retains sufficient information from previous steps. Thus, we
train an additional layer $(ret_{o^{*}},ret_{o^{+}},ret_{w})$ directly after
$M(h_{t})$ to retrieve at step $t$ the target, tool and way-point of step
$t-1$. $L_{ret}$ is defined as
$L_{ret}=CE(ret_{o^{*}}(M(h_{t})),o^{*}_{t-1})+CE(ret_{o^{+}}(M(h_{t})),o^{+}_{t-1})+MSE(ret_{w}(M(h_{t})),w_{t-1}).$
As a result, the complete proposed architecture is trained with the loss
$L=L_{o^{+}}+L_{o^{*}}+L_{w}+L_{action}+L_{mem}+L_{ret}$.

##  V Experimental Results

###  V-A Data collection

We used an Intel RealSense 415 camera to record several demonstrations of a
human subject performing three tasks. The first task consists in inserting a
paint brush into a bucket, then moving it to a painting surface and painting a
virtual straight line on the surface. The poses of the brush, bucket and
painting surface are all tracked in real-time using the technique explained in
V-B. The second task consists in picking up various blocks and stacking them
on top of each other to form a predefined desired pattern. The third task is
similar to the second one, with the only difference being the desired stacking
pattern. Additionally, we use the PyBullet physics engine to simulate a Kuka
robotic arm and collect data regarding a fourth task. The fourth task consists
in moving a wrench that is attached to the end-effector to four precise
locations on a wheel, sequentially, rotating the wrench at each location to
remove the lug-nuts, then moving the wrench to the wheel’s center before
finally pulling it.

Figure 3: Tasks considered in the experiments: painting (left), stacking
(middle), and tire removal (right).

Figure 4: Average success rate (left) and episode length (right) as a function
of the number of training trajectories

###  V-B Object Pose Parsing from Demonstration Video

In each demonstration video, 6D poses $\xi\in SE(3)$ and semantic labels of
all relevant objects are estimated in real-time and used to create
observations $(z_{t})_{t=1}^{H}$ as explained in Section III. Concretely, a
scene-level multi-object pose estimator [28] is leveraged to compute globally
the relevant objects’ 6D poses in the first frame. It starts with a pose
sampling process and performs Integer Linear Programming to ensure physical
consistency by checking collisions between any two objects as well as
collisions between any object and the table. Next, the poses computed from the
first frame are used to initialize the se(3)-TrackNet [17], which returns a 6D
pose for each object in every frame of the video. The 6D tracker requires
access to the objects’ CAD models. For the painting task, the brush and the
bucket are 3D-scanned as in [29], while for all the other objects, models are
obtained from CAD designs derived from geometric primitives. During inference
on the demonstration videos, the tracker operates in $90Hz$, resulting in an
average processing time of $13.3s$ for a 1 min demonstration video. The entire
video parsing process is fully automated, and did not require any human input
beyond providing CAD models of the objects offline to the se(3)-TrackNet in
order to learn to track them.

###  V-C Training and architecture details

The high-level, intermediate-level, and low-level policies are all neural
networks. In the high-level policy, the progress vector $g_{t}$ is embedded by
a fully connected layer with $64$ units followed by a ReLU layer. An LSTM
layer with $32$ units is used to encode history $h_{t}$. Observation $z_{t}$
is concatenated with the LSTM units and the embedded progress vector and fed
as an input to two hidden layers with $64$ units followed by a ReLU layer.
From the last hidden layer, target and tool objects are predicted by a fully
connected layer and a softmax. The intermediate-level policy network consists
of two hidden layers with $64$ units. The low-level policy concatenates into a
vector four inputs: 6D pose of the target in the frame of the tool object,
way-point, and the semantic labels of the target and tool objects. After two
hidden layers of $64$ units and ReLU, the low-level policy outputs a Gaussian
action distribution. The number of training iterations for all tasks is
$20,000$, the batch size is $2,048$ steps, the learning rate is $1e-4$, and
the optimizer is Adam. The hyper-parameters used in Section IV-A are set as
$\gamma=0.95$, $\alpha=100$, and $\beta=0.95$.

The proposed method is compared to three other techniques. Generative
Adversarial Imitation Learning (GAIL) [30], Learning Task Decomposition via
Temporal Alignment for Control (TACO) [8], and Behavioral Cloning [31] where
we train the policy network of [30] directly to maximize the likelihood of the
demonstrations without learning rewards. We also compare to a supervised
variant of our proposed technique where we manually provide ground-truth tool
and target objects and way-points for each frame. The supervised variant
provides an upper bound on the performance of our unsupervised algorithm. Note
also that TACO [8] requires providing manual sketches of the sub-tasks,
whereas our algorithm is fully unsupervised.

###  V-D Evaluation

Except for the tire removal, all tasks are evaluated using real demonstrations
and a real Kuka LBR robot equipped with a Robotiq hand. The policies learned
from real demonstrations are also tested extensively in the PyBullet simulator
before testing them on the real robot.

A painting task is successfully accomplished if the brush is moved into a
specific region of a radius of $3cm$ inside the paint bucket, the brush is
then moved to a plane that is $10cm$ on top of the painting surface, and
finally the brush draws a virtual straight line of $5cm$ at least on that
plane. A tire removal task is successfully accomplished if the robot removes
all bolts by rotating its end-effector on top of each bolt (with a toleance of
$5mm$) with at least $30^{\circ}$ counter-clockwise, and then moves to the
center of the wheel. A stacking task is successful if the centers of all the
objects in their final configuration are within $0.5cm$ of the corresponding
desired locations.

| Painting | Stacking 1 | Stacking 2
---|---|---|---
Proposed (unsupervised) | $5/5$ | $5/5$ | $5/5$
Behavioral Cloning [31] | $0/5$ | $0/5$ | $0/5$
GAIL [30] | $0/5$ | $0/5$ | $0/5$
TACO [8] | $1/5$ | $1/5$ | $2/5$
TABLE I: Success rates on the real Kuka robot

Figure 4 shows the success rates of the compared methods for the four tasks,
as well as the length of the generated trajectories while solving these tasks
in simulation, as a function of the number of demonstration trajectories
collected as explained in Section V-A. The results are averaged over $5$
independent runs, each run contains $50$ test episodes that start with random
layouts of the objects. Table I shows the success rates of the compared
methods on the real Kuka robot, using the same demonstration trajectories that
were used to generate Figure 4 ($70$ trajectories for painting and $40$ for
each of the remaining tasks). These results show clearly that the proposed
approach outperforms the compared alternatives in terms of success rates and
solves the four tasks with a smaller number of actions. The performance of our
unsupervised approach is also close to that of the supervised variant. The
proposed approach outperforms TACO despite the fact that TACO requires a form
of supervision in its training. We also note that both our approach and TACO
outperform GAIL and the behavioral cloning techniques, which clearly indicates
the data-efficiency of compositional and hierarchical methods. Videos and
supplementary material can be found at https://tinyurl.com/2zrp2rzm.

##  VI Conclusion

We presented a unified neural-network framework for training robots to perform
complex manipulation tasks that are composed of several sub-tasks. The
proposed framework employs the principal of attention by training a high-level
policy network to select a pair of tool and target objects dynamically,
depending on the context. The proposed method outperformed alternative
techniques for imitation learning, without requiring any supervision beyond
recorded demonstration videos. While the current video parsing module requires
the objects’ CAD models beforehand, it is possible in future work to leverage
model-free 6D pose trackers [32] for learning from demonstration involving
novel unknown objects. We will also explore other applications of the proposed
framework, such as real-world assembly tasks.

## References

  * [1] S. Calinon, _Robot Programming by Demonstration_ , 1st ed. USA: CRC Press, Inc., 2009.
  * [2] T. Osa, J. Pajarinen, and G. Neumann, _An Algorithmic Perspective on Imitation Learning_. Hanover, MA, USA: Now Publishers Inc., 2018.
  * [3] O. Kroemer and G. S. Sukhatme, “Learning spatial preconditions of manipulation skills using random forests,” in _Proceedings of the IEEE-RAS International Conference on Humanoid Robotics_ , 2016. [Online]. Available: http://robotics.usc.edu/publications/954/
  * [4] A. S. Wang and O. Kroemer, “Learning robust manipulation strategies with multimodal state transition models and recovery heuristics,” in _Proceedings of (ICRA) International Conference on Robotics and Automation_ , May 2019, pp. 1309 – 1315.
  * [5] Z. Wang, C. R. Garrett, L. P. Kaelbling, and T. Lozano-Pérez, “Learning compositional models of robot skills for task and motion planning,” _Int. J. Robotics Res._ , vol. 40, no. 6-7, 2021. [Online]. Available: https://doi.org/10.1177/02783649211004615
  * [6] Z. Su, O. Kroemer, G. E. Loeb, G. S. Sukhatme, and S. Schaal, “Learning to switch between sensorimotor primitives using multimodal haptic signals,” in _Proceedings of International Conference on Simulation of Adaptive Behavior (SAB ’16): From Animals to Animats 14_ , August 2016, pp. 170 – 182.
  * [7] H. Le, N. Jiang, A. Agarwal, M. Dudík, Y. Yue, and H. Daumé, “Hierarchical imitation and reinforcement learning,” in _International conference on machine learning_. PMLR, 2018, pp. 2917–2926.
  * [8] K. Shiarlis, M. Wulfmeier, S. Salter, S. Whiteson, and I. Posner, “Taco: Learning task decomposition via temporal alignment for control,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 4654–4663.
  * [9] G. Konidaris and A. Barto, “Skill discovery in continuous reinforcement learning domains using skill chaining,” in _Advances in Neural Information Processing Systems_ , Y. Bengio, D. Schuurmans, J. Lafferty, C. Williams, and A. Culotta, Eds., vol. 22. Curran Associates, Inc., 2009.
  * [10] M. Toussaint, K. R. Allen, K. A. Smith, and J. B. Tenenbaum, “Differentiable physics and stable modes for tool-use and manipulation planning,” in _Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, IJCAI 2019, Macao, China, August 10-16, 2019_ , S. Kraus, Ed. ijcai.org, 2019, pp. 6231–6235. [Online]. Available: https://doi.org/10.24963/ijcai.2019/869
  * [11] L. P. Kaelbling, “Learning to achieve goals,” in _IN PROC. OF IJCAI-93_. Morgan Kaufmann, 1993, pp. 1094–1098.
  * [12] L. P. Kaelbling and T. Lozano-Pérez, “Hierarchical task and motion planning in the now,” in _Proceedings of the 1st AAAI Conference on Bridging the Gap Between Task and Motion Planning_ , ser. AAAIWS’10-01. AAAI Press, 2010, pp. 33–42.
  * [13] C. R. Garrett, R. Chitnis, R. Holladay, B. Kim, T. Silver, L. P. Kaelbling, and T. Lozano-Pérez, “Integrated task and motion planning,” 2020.
  * [14] R. T. Icarte, T. Klassen, R. Valenzano, and S. McIlraith, “Using reward machines for high-level task specification and decomposition in reinforcement learning,” ser. Proceedings of Machine Learning Research, J. Dy and A. Krause, Eds., vol. 80. StockholmsmÃ€ssan, Stockholm Sweden: PMLR, 10–15 Jul 2018, pp. 2107–2116.
  * [15] R. Toro Icarte, E. Waldie, T. Klassen, R. Valenzano, M. Castro, and S. McIlraith, “Learning reward machines for partially observable reinforcement learning,” _Advances in Neural Information Processing Systems_ , vol. 32, pp. 15 523–15 534, 2019.
  * [16] A. Camacho, R. T. Icarte, T. Q. Klassen, R. A. Valenzano, and S. A. McIlraith, “Ltl and beyond: Formal languages for reward function specification in reinforcement learning.”
  * [17] B. Wen, C. Mitash, B. Ren, and K. E. Bekris, “se (3)-tracknet: Data-driven 6d pose tracking by calibrating image residuals in synthetic domains,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, pp. 10 367–10 373.
  * [18] D. Kalashnikov, A. Irpan, P. Pastor, J. Ibarz, A. Herzog, E. Jang, D. Quillen, E. Holly, M. Kalakrishnan, V. Vanhoucke, and S. Levine, “Scalable deep reinforcement learning for vision-based robotic manipulation,” ser. Proceedings of Machine Learning Research, A. Billard, A. Dragan, J. Peters, and J. Morimoto, Eds., vol. 87. PMLR, 29–31 Oct 2018, pp. 651–673.
  * [19] R. Fox, R. Shin, S. Krishnan, K. Goldberg, D. Song, and I. Stoica, “Parametrized hierarchical procedures for neural programming,” in _International Conference on Learning Representations_ , 2018. [Online]. Available: https://openreview.net/forum?id=rJl63fZRb
  * [20] D. Xu, S. Nair, Y. Zhu, J. Gao, A. Garg, L. Fei-Fei, and S. Savarese, “Neural task programming: Learning to generalize across hierarchical tasks.” _CoRR_ , vol. abs/1710.01813, 2017. [Online]. Available: http://dblp.uni-trier.de/db/journals/corr/corr1710.html#abs-1710-01813
  * [21] D. Huang, S. Nair, D. Xu, Y. Zhu, A. Garg, L. Fei-Fei, S. Savarese, and J. C. Niebles, “Neural task graphs: Generalizing to unseen tasks from a single video demonstration,” _CoRR_ , vol. abs/1807.03480, 2018. [Online]. Available: http://arxiv.org/abs/1807.03480
  * [22] S. Nair, M. Babaeizadeh, C. Finn, S. Levine, and V. Kumar, “TRASS: time reversal as self-supervision,” in _2020 IEEE International Conference on Robotics and Automation, ICRA 2020, Paris, France, May 31 - August 31, 2020_. IEEE, 2020, pp. 115–121. [Online]. Available: https://doi.org/10.1109/ICRA40945.2020.9196862
  * [23] M. Andrychowicz, F. Wolski, A. Ray, J. Schneider, R. Fong, P. Welinder, B. McGrew, J. Tobin, P. Abbeel, and W. Zaremba, “Hindsight experience replay,” in _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , ser. NIPS’17. Curran Associates Inc., 2017, pp. 5055–5065.
  * [24] S. Nair and C. Finn, “Hierarchical foresight: Self-supervised learning of long-horizon tasks via visual subgoal generation,” _CoRR_ , vol. abs/1909.05829, 2019. [Online]. Available: http://arxiv.org/abs/1909.05829
  * [25] P.-L. Bacon, J. Harb, and D. Precup, “The option-critic architecture,” ser. AAAI’17. AAAI Press, 2017, pp. 1726–1734.
  * [26] O. Nachum, S. Gu, H. Lee, and S. Levine, “Data-efficient hierarchical reinforcement learning,” in _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , ser. NIPS’18, 2018, pp. 3307–3317.
  * [27] B. Eysenbach, A. Gupta, J. Ibarz, and S. Levine, “Diversity is all you need: Learning skills without a reward function,” in _International Conference on Learning Representations_ , 2019. [Online]. Available: https://openreview.net/forum?id=SJx63jRqFm
  * [28] C. Mitash, B. Wen, K. Bekris, and A. Boularias, “Scene-level pose estimation for multiple instances of densely packed objects,” in _Conference on Robot Learning_. PMLR, 2020, pp. 1133–1145.
  * [29] A. S. Morgan, B. Wen, J. Liang, A. Boularias, A. M. Dollar, and K. Bekris, “Vision-driven compliant manipulation for reliable, high-precision assembly tasks,” _RSS_ , 2021.
  * [30] J. Ho and S. Ermon, “Generative adversarial imitation learning,” _Advances in neural information processing systems_ , vol. 29, pp. 4565–4573, 2016.
  * [31] D. A. Pomerleau, “Alvinn: An autonomous land vehicle in a neural network,” Carnegie-Mellon University Pittsburgh PA, Tech. Rep., 1989.
  * [32] B. Wen and K. Bekris, “Bundletrack: 6d pose tracking for novel objects without instance or category-level 3d models,” _IROS_ , 2021.