# Shared Control for Bimanual Telesurgery with Optimized Robotic Partner

Ziwei Wang, Yanpei Huang, Xiaoxiao Cheng,
Pakorn Uttayopas, Etienne Burdet This research was partly supported by the
European Commission grants H2020 PH-CODING (FETOPEN 829186), NIMA (FETOPEN
899626) and UK EPSRC grant FAIR-SPACE EP/R026092/1. (Corresponding author:
Yanpei Huang and Xiaoxiao Cheng)The authors are with the Department of
Bioengineering, Imperial College of Science, Technology and Medicine, SW7 2AZ,
United Kingdom (e-mail: {ziwei.wang, yanpei.huang,
xiaoxiao.cheng,pakorn.uttayopas18, e.burdet}@imperial.ac.uk).

###### Abstract

Traditional telesurgery relies on the surgeon’s full control of robot on the
patient’s side, which tends to increase surgeon fatigue and may reduce the
efficiency of the operation. This paper introduces a Robotic Partner (RP) to
facilitate intuitive bimanual telesurgery, aiming at reducing the surgeon
workload and enhancing surgeon-assisted capability. An interval type-2
polynomial fuzzy-model-based learning algorithm is employed to extract expert
domain knowledge from surgeons and reflect environmental interaction
information. Based on this, a bimanual shared control is developed to interact
with the other robot teleoperated by the surgeon, understanding their control
and providing assistance. As prior information of environment model is not
required, it reduces reliance on force sensors in control design. Experimental
results on the DaVinci Surgical System show that the RP could assist peg-
transfer tasks and reduce the surgeon’s workload by 51% in force-sensor-free
scenarios.

###### Index Terms:

Shared control, force estimate, polynomial fuzzy model, bimanual telesurgery,
DaVinci surgical system.

##  I Introduction

Robotic technology and its integration with surgical procedures have largely
improved the operation intuitiveness for the surgeons [1]. Most of the current
surgical robotic systems are conducted in master-slave mode, in which the
surgeon could comfortably sit at a master console, remotely controlling the
slave robotic arm and holding surgical tools through master interfaces and
high-definition camera view [2, 3]. These teleoperation systems tend to
simplify the surgeon operation to some extent. However, the surgeon’s
performance would suffer from fatigue during long time operation. Also, the
surgeon’s performance is limited by delayed visual and tactile information
from the interaction environment due to the communication channel. Therefore,
sharing the surgical robot with a degree of autonomy could facilitate the
telesurgery, easing both physical and mental workloads for surgeons at the
master console.

To this end, the slave surgical robot needs to perceive interaction forces in
terms of telepresence. Nevertheless, it is challenging to obtain interactive
forces in contact with patient tissue directly due to the mounting position of
the force sensors conflicting with the surgical tool. Some studies [4, 5]
integrated the environment information into the control design by force
sensors or model calculation to enable the robot to conduct force interaction
with environment automatically. These methods require accurate model
structures of remote interactive dynamics or force sensors on the end-effector
of the robot. It is revealed in [6] that accurate dynamics model is difficult
to build in practical teleoperation subject to uncertainties, which may limit
the merits of the above methods in telesurgery.

Another critical issue to realise the target is how the surgical robot can
predict the human intention and give assistance, thus reducing the workload.
The current methods mostly focus on training the robot to carry out repetitive
tasks [7, 8], which lack the control flexibility to reflect on different
scenarios. To assist the surgeon better, the surgical robot is required to
adapt to the cooperated surgeon, react to the interacted environment and even
improve the performance of the surgeon through sensory augmentation [9].

Based on this idea, we propose an intelligent Robotic Partner (RP) in bimanual
shared control scheme to solve the two issues mentioned above. The RP is
trained by interval type-2 polynomial fuzzy-model-based learning (IT2 PFMBL)
algorithm according to human-human demonstration and robot-environment
interaction, which enables the learning speed with few-shot data through
virtue of the fitting ability of polynomials and the nonlinear expression
ability of fuzzy models. The RP could then estimate interaction force without
the exact physical model nor force sensors. Moreover, the proposed RP provides
intuitive assistance to perform a cooperation task (e.g. peg-transfer) with
human operator.

##  II Method

IT2 PFMBL algorithm is designed to estimate the human motions (i.e., movement
trajectories and gesture actions) and interactive force with environment
simultaneously based on the human-robot and robot-environment interaction
data.

Let collected input-output pair data be
$[x_{i}(k),v_{i}(k),y_{i}(k)]_{k=1}^{N}$, where $i\in\\{h_{mt},h_{ga},e\\}$
stands for human movement trajectories, gesture actions or interactive forces
with environment, respectively.
$x_{i}(k),v_{i}(k),y_{i}(k)\in\mathbb{R}^{n_{i}}$ are position, velocity and
the associated output in Cartesian space at the $k$th sampling. $N$ is the
number of sampling pairs. Fuzzy C-Means clustering [10] and subtractive
clustering [11] are utilised to obtain the number of fuzzy sets $a_{i}$, the
centre of each fuzzy set ${\mathcal{B}}_{\nu}^{i}$ and the corresponding
membership function
${\mu}_{\mathcal{B}_{\nu}^{i}}\mathcal{A}_{i}^{\nu}(z_{i}(k))$, where
$\nu\in\mathbb{N}_{+}{[1,a]}$,
$z_{i}(k)=[x_{i}(k)\;v_{i}(k)\;v_{i}(k+1)]\in\mathbb{R}^{3n_{i}}$ and
$\mathcal{A}_{i}^{\nu}(z_{i}(k))$ is the premise variable. The membership
function
$0\leq\underline{\mu}_{\mathcal{B}_{\nu}^{i}}\mathcal{A}_{i}^{\nu}(z_{i}(k))\leq{\mu}_{\mathcal{B}_{\nu}^{i}}\mathcal{A}_{i}^{\nu}(z_{i}(k))\leq\bar{\mu}_{\mathcal{B}_{\nu}^{i}}\mathcal{A}_{i}^{\nu}(z_{i}(k))\leq
1$ describes the degree to which the input data matches the fuzzy set.

With the help of the obtained fuzzy sets, the inner relationship between the
input and output can be represented by the following polynomial discrete-time
fuzzy logic system

$\displaystyle\text{Rule}\;l:$
$\displaystyle\text{IF}\;\mathcal{A}_{i}^{1}(z_{i}(k))\;\text{is}\;{\mathcal{B}}_{1}^{il}\;\text{AND},...,\text{AND}\;\mathcal{A}_{i}^{a}(z_{i}(k))\;\text{is}\;{\mathcal{B}}_{a}^{il}$
(1)
$\displaystyle\text{THEN}\;y_{i}(k)=M_{i}^{l}(z_{i}(k))\frac{v_{i}(k+1)-v_{i}(k)}{\Delta
T}$
$\displaystyle+C_{i}^{l}(z_{i}(k))v_{i}(k)+K_{i}^{l}(z_{i}(k))x_{i}(k)+f_{i}^{l}(z_{i}(k)),$

where ${\Delta T}$ is the sampling time; ${\mathcal{B}}_{\nu}^{il}$ is a
polynomial fuzzy set of Rule $l$ corresponding to the premise variable
$\mathcal{A}_{i}^{\nu}(z_{i}(k))$, $l\in\mathbb{N}_{+}{[1,p]}$;
$M_{i}^{l}(z_{i}(k))\in\mathbb{R}^{n_{i}\times n_{i}}$,
$C_{i}^{l}(z_{i}(k))\in\mathbb{R}^{n_{i}\times n_{i}}$,
$K_{i}^{l}(z_{i}(k))\in\mathbb{R}^{n_{i}\times n_{i}}$ and
$f_{i}^{l}(z_{i}(k))\in\mathbb{R}^{n_{i}}$ are the polynomial system matrices,
where each element is a polynomial in $z_{i}(k)$. For example, in the case of
$i=e$, the first three system matrices donate the equivalent inertia, damping
and stiffness ones while the last the bounded exogenous force concerned with
the position of the contact point. Robust linear least-square approach is
employed in each Rule to identify the polynomial coefficients of the
associated system matrices.

Considering the measurement noise in $z_{i}(k)$ and the system uncertainty,
the firing strength of the Rule $l$ is of the following interval set

$W_{i}^{l}(z_{i}(k))=[\underline{\omega}_{i}^{l}(z_{i}(k)),\bar{\omega}_{i}^{l}(z_{i}(k))],l\in\mathbb{N}_{+}{[1,p]},$
(2)

where
$\underline{\omega}_{i}^{l}(z_{i}(k))=\prod_{\nu=1}^{a}\underline{\mu}_{\mathcal{B}_{\nu}^{il}}\mathcal{A}_{i}^{\nu}(z_{i}(k))$,
$\bar{\omega}_{i}^{l}(z_{i}(k))=\prod_{\nu=1}^{a}\bar{\mu}_{\mathcal{B}_{\nu}^{il}}\mathcal{A}_{i}^{\nu}(z_{i}(k))$,
in which
$\underline{\mu}_{\mathcal{B}_{\nu}^{il}}\mathcal{A}_{i}^{\nu}(z_{i}(k))$ and
$\bar{\mu}_{\mathcal{B}_{\nu}^{il}}\mathcal{A}_{i}^{\nu}(z_{i}(k))$ denote the
lower and upper membership function such that
$0\leq\underline{\mu}_{\mathcal{B}_{\nu}^{il}}\mathcal{A}_{i}^{\nu}(z_{i}(k))\leq\bar{\mu}_{\mathcal{B}_{\nu}^{il}}\mathcal{A}_{i}^{\nu}(z_{i}(k))\leq
1$; $\underline{\omega}_{i}^{l}(z_{i}(k))$ and
$\bar{\omega}_{i}^{l}(z_{i}(k))$ are the lower and upper grade of membership,
respectively, satisfying
$0\leq\underline{\omega}_{i}^{l}(z_{i}(k))\leq\bar{\omega}_{i}^{l}(z_{i}(k))$.
Hence, the data-driven IT2 PFMBL system is described by

$\displaystyle y_{i}(k)=$
$\displaystyle\sum_{l=1}^{p}\tilde{\omega}_{i}^{l}(z_{i}(k))\Bigg{(}M_{i}^{l}(z_{i}(k))\frac{v_{i}(k+1)-v_{i}(k)}{\Delta
T}$ (3)
$\displaystyle+C_{i}^{l}(z_{i}(k))v_{i}(k)+K_{i}^{l}(z_{i}(k))x_{i}(k)+f_{i}^{l}(z_{i}(k))\Bigg{)},$

and

$\tilde{\omega}_{i}^{l}(z_{i}(k))=\underline{\omega}_{i}^{l}(z(k))\underline{b}_{i}^{l}(z_{i}(k))+\bar{\omega}_{l}(z_{i}(k))\bar{b}_{i}^{l}(z_{i}(k)),$
(4)

with $\underline{b}_{i}^{l}(z_{i}(k))$, $\bar{b}_{i}^{l}(z_{i}(k))\in[0,1]$
such that $\underline{b}_{i}^{l}(z_{i}(k))+\bar{b}_{i}^{l}(z_{i}(k))=1$.

Remark 1. A unified data-driven learning paradigm is developed without the
prior information (i.e., deterministic model structure) in [12, 13]. For
$i\in\\{h_{mt},h_{ga},e\\}$, $y_{i}$ represents the human motion trajectories
in Cartesian space, gripper angles and interactive forces with the
environment, leading to different system matrices for different tasks.

##  III Experimental validation and results

In this section, we have verified the proposed method by the following
experiments, which provides implementation foundation for the RP.

###  III-A Prediction of robot-environment interaction force

Figure 1: H-Man planar robot-environment interaction platform with robotic
finger.

We firstly apply the proposed IT2 PFMBL method to predict the interaction
force in a planar robot-environment interaction platform using modified H-Man
system [14], which is shown in Fig. 1. The environment is selected as hard
silicone having similar mechanical characteristics with the organ tissues. The
robotic finger equipped with an ATI 6-axis force/torque sensor will press the
silicon at three different levels of depth: shallow, medium and deep, and scan
the whole environment surface to precept the overall material properties.
Interaction force is collected during 50 interactive trials for each type of
pressing, totalling 150 trials for each sampling point. The proposed IT2 PFMBL
method is advantageous at producing high-prediction accuracy with a small
training set by introducing a type-2 fuzzy rule. We randomly allocate
available data into 10/90 proportions of training and testing sets. Fig. 2a
shows that the root-mean-square-error (RMSE) of the IT2 PFMBL method decreases
rapidly with few-shot training data.

Figure 2: (a) RMSE of interactive force estimation using IT2 PFMBL algorithm.
(b) Performance comparison in terms of SRMSE and SAME (The smaller these two
indicators are, the better). Figure 3: Experimental setup. (a) Master
interfaces: dual Omega 7 devices with visual feedback. (b) Slave robot:
surgical robotic manipulators. (c) Task: peg transfer.

To validate the merits of the IT2 PFMBL method, we compare its performances
with other three common interaction force models, namely linear Kelvin-Voigt
(LKV) model [15], T-S FMB [16] and PFMB model. Due to the improved non-linear
representation of material stiffness properties, Fig. 2b indicates the force
estimation of IT2 PFMBL surpasses the above models on the sum of RMSE (SRMSE)
and mean absolute error (SMAE) in testing trials.

###  III-B Human-Robot shared control using IT2 PFMBL

To further assess the performance of the IT2 PFMBL method, we implemented it
in a surgical robotic teleoperation system to train the RP to cooperate with a
human operator. The setup of the surgical robotic system is shown in Fig. 3.
The master console (Fig. 3a) consists of two hand haptic interfaces (Force
Dimension Inc., Switzerland) and a monitor providing visual feedback of the
patient side. Two patient-side manipulators of Da Vinci robotic system
(Intuitive Surgical, USA) are used as slave robotic arms as shown in Fig. 3b.

The end-effectors of the robotic arms are connected to surgical grippers to
perform a peg-transfer task (Fig. 3c), which is a standard laparoscopic
training task from Fundamentals of Laparoscopic Surgery (FLS) [17]. The left
gripper is controlled by a human operator and the right gripper by a RP.
During the task, the human operator navigates the left gripper to the initial
peg with an orange block, grasps the block and moves it to conduct a handover
to the RP controlled gripper. After the handover, the RP moves the block to
the target peg (Fig. 3c) (see supplementary video 1). The manipulation data of
two operators were used as the input to the IT2 PFMBL method, resulting in a
final interaction environment model including the motion trajectories, and
gripper movements. After training, the RP could imitate the human partner and
cooperate intelligently with another human operator.

Figure 4: Experimental results from (a-b) human-human training to (c-d) human-
robot testing. Motion learning from human operator and human partner:(a)
Cartesian trajectories in human-human telesurgery (x, y and z denote the right
hand movement along the associated direction; the dashed lines are the
estimate with respect to the associated direction using IT2 PFMBL). (b)
Gripper grasping $\theta$: the dashed line stands for the estimate of gripper
movement. $\theta_{\text{LH}}$ and $\theta_{\text{RH}}$ represent the angle at
which the left and right hand grippers open. Motion learning from human
operator and RP: (c) Cartesian trajectories in human-robot telesurgery (x, y
and z denote the RP movement along the associated direction). (d) Assistive
gripper grasping: $\theta_{\text{RP}}$ is the angle of RP gripper.

Fig. 4 a&b show the comparison of motion trajectories between human-human
manipulation and our training results from IT2 PFMBL method. It indicates that
the derived models could effectively represent the human partner’s behavior
with minimal errors. We implemented these models in the surgical robotic
system to validate their performance. The results (Fig. 4c,d) show that the
trained RP could successfully complete the transfer task with closely
cooperating with a human partner (see supplementary video 2). The task
completion of human operator with RP is within 27 seconds.

##  IV Discussion and Conclusion

An innovative Robotic Partner concept is proposed in shared control for
bimanual telesurgery to reduce the workload of the surgeon. A novel interval
type-2 polynomial fuzzy model-based method is designed to realize precise
estimation of contact force and few-shot training for intelligent surgeon
assistance. The experimental results illustrate that the proposed IT2 PFMBL
method achieves superior nonlinear representation capability and fitting
performance in terms of SRMSE and SAME. Compared to T-S FMB method, PFMBL
method can better characterise the nonlinearities of interaction force because
its model parameters draw on the ability to fit polynomials. Besides, IT2
fuzzy sets are effective to capture the noise and system uncertainty by using
the upper and lower bound of membership functions, thereby demonstrating
satisfactory estimation robustness. The IT2 PFMBL is therefore able to
estimate force accurately and provide haptic feedback to telesurgery system,
avoiding installing force sensors on delicate surgical tools.

Furthermore, the experiment results showed that the RP could reduce the
intensity of the surgeon’s manipulation effectively with the implemented IT2
PFMBL method, by assisting and adapting to surgeon behaviour autonomously.
Trained by trajectories and gripper movements of two human operators, the
intuitive human-robot cooperation is demonstrated in a peg-transfer task with
self-determining movements from the RP. With the assistance of the RP, the
surgeon could teleoperate the other surgical robot in the peg-transfer task
alone avoiding the communication errors of two surgeons and saving the
manpower.

## References

  * [1] J. Reynoso, A. Meyer, J. Unnirevi, and D. Oleynikov, _Robotics for minimally invasive surgery (MIS) and natural orifice transluminal endoscopic surgery (NOTES)_. Elsevier Ltd, 2012, pp. 210–223.
  * [2] Y. Huang, W. Lai, L. Cao, J. Liu, E. Burdet, and S. J. Phee, “A three-limb teleoperated robotic system with foot control for flexible endoscopic surgery,” _Annals of Biomedical Engineering_ , 2021.
  * [3] S. Treratanakulchai and F. Baena, “A passive decoupling mechanism for misalignment compensation in master-slave teleoperation,” _IEEE Transactions on Medical Robotics and Bionics_ , vol. 3, pp. 285–288, 2021.
  * [4] C. Yang, G. Peng, L. Cheng, J. Na, and Z. Li, “Force sensorless admittance control for teleoperation of uncertain robot manipulator using neural networks,” _IEEE Transactions on Systems, Man, and Cybernetics: Systems_ , vol. 51, no. 5, pp. 3282–3292, 2021.
  * [5] J.-Z. Xu, M.-F. Ge, Z.-W. Liu, W.-Y. Zhang, and W. Wei, “Force-reflecting hierarchical approach for human-aided teleoperation of nrs with event-triggered local communication,” _IEEE Transactions on Industrial Electronics_ , pp. 1–1, 2021.
  * [6] Z. Wang, H.-K. Lam, B. Xiao, Z. Chen, B. Liang, and T. Zhang, “Event-triggered prescribed-time fuzzy control for space teleoperation systems subject to multiple constraints and uncertainties,” _IEEE Transactions on Fuzzy Systems_ , pp. 1–1, 2020.
  * [7] G. A. Fontanelli, G.-Z. Yang, and B. Siciliano, “A comparison of assistive methods for suturing in MIRS,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2018, pp. 4389–4395.
  * [8] H. Ishida, M. M. Marinho, K. Harada, J. Gao, and M. Mitsuishi, “Virtual-fixtures for robotic-assisted bi-manual cutting using vector-field inequalities,” in _2020 IEEE/SICE International Symposium on System Integration (SII)_ , 2020, pp. 395–400.
  * [9] A. Takagi, G. Ganesh, T. Yoshioka, M. Kawato, and E. Burdet, “Physically interacting individuals estimate the partner’s goal to enhance their movements,” _Nature Human Behaviour_ , vol. 1, no. 3, pp. 1–6, 2017.
  * [10] J. C. Bezdek, _Pattern Recognition with Fuzzy Objective Function Algorithms_. Springer Science & Business Media, 2013.
  * [11] S. L. Chiu, “Fuzzy model identification based on cluster estimation,” _Journal of Intelligent & Fuzzy Systems_, vol. 2, no. 3, pp. 267–278, 1994\.
  * [12] A. Shahdi and S. Sirouspour, “Adaptive/robust control for time-delay teleoperation,” _IEEE Transactions on Robotics_ , vol. 25, no. 1, pp. 196–205, 2009.
  * [13] Z. Lu, P. Huang, and Z. Liu, “High-gain nonlinear observer-based impedance control for deformable object cooperative teleoperation with nonlinear contact model,” _International Journal of Robust and Nonlinear Control_ , vol. 30, no. 4, pp. 1329–1350, 2020.
  * [14] D. Campolo, P. Tommasino, K. Gamage, J. Klein, C. M. Hughes, and L. Masia, “H-man: A planar, h-shape cabled differential robotic manipulandum for experiments on human motor control,” _Journal of Neuroscience Methods_ , vol. 235, pp. 285–297, 2014.
  * [15] H. Li and K. Kawashima, “Achieving stable tracking in wave-variable-based teleoperation,” _IEEE/ASME Transactions on Mechatronics_ , vol. 19, no. 5, pp. 1574–1582, 2014.
  * [16] Z. Wang, H.-K. Lam, Z. Chen, B. Liang, and T. Zhang, “Event-triggered interval type-2 fuzzy control for uncertain space teleoperation systems with state constraints,” in _2020 IEEE International Conference on Fuzzy Systems (FUZZ-IEEE)_ , 2020, pp. 1–8.
  * [17] H.-C. Hur, D. Arden, L. E. Dodge, B. Zheng, and H. A. Ricciotti, “Fundamentals of laparoscopic surgery: A surgical skills assessment tool in gynecology,” _JSLS Journal of the Society of Laparoscopic & Robotic Surgeons_, vol. 15, no. 1, pp. 21–26, 2011.


Figure 1. We propose a novel method to correctly differentiate a neural SDF
renderer by reparameterizing the pixel integral. Direct application of
automatic differentiation to the renderer fails because of discontinuities
like silhouette boundaries. In this work we show that, by carefully designing
a discontinuity-aware warp function $\mathcal{V}(u;\theta)$ to reparameterize
the pixel domain, we can remove these discontinuities, and the reparameterized
integral is amenable to automatic differentiation. We demonstrate the benefits
of our method on inverse rendering problems. Starting from a multiview dataset
of real photos (a), our reparameterized renderer (c) can optimize a neural SDF
that closely matches the input data, and generalizes to novel views. Our
renderer matches or outperforms prior SDF renderers (Yariv et al., 2020) (b),
while doing away with their need for additional geometric supervision in the
form of per-view masks, which can be unreliable for real-world data. We show
additional surface reconstructions obtained with our inverse renderer in (d).

# Differentiable Rendering of Neural SDFs through Reparameterization

Sai Praveen Bangaru  MIT CSAILUSA<EMAIL_ADDRESS>,  Michaël Gharbi  Adobe
ResearchUSA<EMAIL_ADDRESS>,  Tzu-Mao Li  UC San DiegoUSA<EMAIL_ADDRESS>,
Fujun Luan  Adobe ResearchUSA<EMAIL_ADDRESS>,  Kalyan Sunkavalli  Adobe
ResearchUSA<EMAIL_ADDRESS>,  Miloš Hašan  Adobe ResearchUSA
<EMAIL_ADDRESS>,  Sai Bi  Adobe ResearchUSA<EMAIL_ADDRESS>,  Zexiang Xu
Adobe ResearchUSA<EMAIL_ADDRESS>,  Gilbert Bernstein  MIT CSAIL & UC
BerkeleyUSA<EMAIL_ADDRESS>and  Frédo Durand  MIT CSAILUSA<EMAIL_ADDRESS>
###### Abstract.

We present a method to automatically compute _correct_ gradients with respect
to geometric scene parameters in neural SDF renderers. Recent physically-based
differentiable rendering techniques for meshes have used _edge-sampling_ to
handle discontinuities, particularly at object silhouettes, but SDFs do not
have a simple parametric form amenable to sampling. Instead, our approach
builds on _area-sampling_ techniques and develops a continuous warping
function for SDFs to account for these discontinuities. Our method leverages
the distance to surface encoded in an SDF and uses quadrature on sphere tracer
points to compute this warping function. We further show that this can be done
by subsampling the points to make the method tractable for neural SDFs. Our
differentiable renderer can be used to optimize neural shapes from multi-view
images and produces comparable 3D reconstructions to recent SDF-based inverse
rendering methods, without the need for 2D segmentation masks to guide the
geometry optimization and no volumetric approximations to the geometry.

††submissionid: 322

##  1\. Introduction

Differentiable rendering algorithms have become crucial tools in solving
challenging inverse problems (Zhao et al., 2020), thanks to their ability to
compute the derivatives of images with respect to arbitrary scene parameters.
Naive differentiation of rendering algorithms does not handle discontinuities
caused by visibility changes and object boundaries correctly. Previous work
has observed that the discontinuities can be handled by properly handling the
Dirac delta signals, and derived algorithms for explicit geometry
representations like triangle meshes (Li et al., 2018; Zhang et al., 2020).

On the other hand, implicit representations like signed distance fields (SDFs)
are appealing since they do not require the initialized geometry to have the
right topology. Recent work has demonstrated the use of SDFs—usually
parameterized using multi-layer perceptron networks—for the task of
reconstructing shape and surface reflectance from images. However, these
methods either require additional geometric supervision such as segmentation
masks (Yariv et al., 2020; Zhang et al., 2021b) or make approximations to the
geometry using volumetric models (Yariv et al., 2021; Oechsle et al., 2021)
that limit their applicability.

In this paper, we derive an algorithm to automatically compute _correct_
gradients with respect to geometric scene parameters in neural SDF renderers.
Previous methods that rely on silhouette sampling are not directly applicable
to SDFs since direct sampling of boundaries of implicit functions is
challenging. Instead, we build on the _reparameterization_ approaches (Loubet
et al., 2019; Bangaru et al., 2020), which removes discontinuities through
reparameterization while preserving the integral values. These methods do not
require explicit sampling along discontinuities. Previous reparameterization
methods focused on triangle meshes, and require new derivation for
reparameterizing SDF rendering.

Figure 2. Overview. Without proper case a SDF rendering pipeline is
discontinuous, which means there are points $u$ where the rendering function
$L(u;\theta)$ is not differentiable in $\theta$, highlighted in red (a). Our
method uses intermediate points from a sphere tracer (b) applied to an SDF
$f$, to compute a warp function $\mathcal{V}$ (c). Using this warp, we
reparameterize the integration domain to avoid discontinuities (d), which
allows us to compute correct gradients of the rendering equation. The key to
achieving this is to design the warp $\mathcal{V}$ so it is continuous in $u$
everywhere, and satisfies some consistency criterion on the geometric
boundaries.

Specifically, we construct a silhouette-aware reparameterization similar to
that of Loubet et al. (2019), but following the equivalent unbiased warp
definition that Bangaru et al. (2020) used to produce correct gradients for
triangle-meshes. We leverage the fact that SDFs naturally encode the distance
to the surface, and develop a practical algorithm that uses a quadrature on
sphere tracing (Hart, 1996) samples to construct a reparameterization that
removes the discontinuities. We further show that this can be computed using
only a subset of sphere tracing samples, reducing computational burden of the
backward pass for bulky neural SDFs111We only became aware of the concurrent
work ”Differentiable Signed Distance Function Rendering” [Vicini et al. 2022]
a week before submission. They use a warping reparameterization similar to
ours that also builds on Bangaru et al. (2020), but focus on voxelized SDFs.
We show that reparameterization works with neural SDFs. Our method also uses a
more efficient ‘top-k’ weighting scheme..

Our algorithm produces correct geometry gradients for SDFs. It does away with
the segmentation masks and depth guidance required by previous techniques
(Yariv et al., 2020), without making a volumetric approximation to the
geometry (Yariv et al., 2021; Oechsle et al., 2021). We show that our
differentiable renderer can be used to optimize neural shapes from multi-view
images, with no additional information beyond the RGB data and the
corresponding camera parameters. Our focus is on occlusion discontinuities, so
the rest of the paper assumes a differentiable shading model.

##  2\. Related Work

We focus on work that recover the latent 3D scene from images through
differentiable rendering. We categorize them by the type of scene
representation.

##### Meshes.

To account for discontinuities, earlier work focused on approximating the
derivatives of mesh rendering by smoothing the geometry (de La Gorce et al.,
2011; Loper and Black, 2014; Rhodin et al., 2015; Kato et al., 2018; Liu et
al., 2019). Alternatively, some work derived correct analytical derivatives
under simplified assumptions (Arvo, 1994; Zhou et al., 2021). Li et al. (2018)
noticed that the differentiation of discontinuities caused by the visibility
and geometric boundaries lead to Dirac delta signals, and can be integrated by
the pixel antialiasing integral or the rendering equation. They proposed an
_edge sampling_ algorithm to explicitly sample the Dirac delta on triangle
mesh silhouettes. Importance sampling the silhouettes can be difficult,
therefore Loubet et al. and Bangaru et al. (2019; 2020) later proposed to
convert the silhouette integral into an area integral. Loubet et al.
formulated the conversion using a reparametrization, and derived an
approximated reparametrization to remove discontinuities. Bangaru et al. built
on Loubet et al.’s work and derived an unbiased estimator by showing the
equivalence between the reparametrization and divergence theorem. On the other
hand, Zhang et al. (2020) showed that directly sampling silhouette in path-
space (Veach, 1998) can also be done efficiently. Directly sampling the
silhouette for SDFs is difficult. Our work extends the reparametrization
approach to handle SDFs, including approximate SDFs defined by neural
networks.

##### Level sets and signed distance fields.

A level set defines a surface using the roots of a 3D implicit function. A
signed distance field is a specific kind of level set where the implicit
function defines the distance of a 3D point to the surfaces, where the sign is
negative when the point is inside the object. SDFs can be represented using
polynomials (Blinn, 1982), voxels (Izadi et al., 2011), or neural networks
(Park et al., 2019). Differentiable rendering for SDFs has been discussed in
computer vision and used for 3D surface reconstruction (Niemeyer et al., 2020;
Jiang et al., 2020; Yariv et al., 2020; Zhang et al., 2021b; Kellnhofer et
al., 2021), but currently methods all ignore the discontinuities when
differentiating and require 2D object masks to converge. An alternative way to
render the signed distance field is to convert it to another format such as a
thin participating medium (Oechsle et al., 2021; Yariv et al., 2021; Wang et
al., 2021), a mesh (Remelli et al., 2020), or a point cloud (Cole et al.,
2021). These methods all introduce approximation. Instead, we focus on
deriving accurate gradients without approximation.

##### Volumes.

A scene can also be represented as participating media instead of solid
surfaces. Gkioulekas et al. (2013) pioneered the use of differentiable volume
rendering for inverse problems. Zhang et al. (2019; 2021c) tackled
discontinuities at volumetric boundaries. Recently, there has been a surging
interest in using volumetric representations—parameterized either as
discretized grids or neural networks—for view synthesis (Lombardi et al.,
2019; Mildenhall et al., 2020; Liu et al., 2020; Xie et al., 2022). These
volumetric representations allow for a trivially differentiable rendering
model and can achieve high-quality novel view synthesis and appearance
acquisition (Bi et al., 2020b, a). However, it is still a challenge to extract
high-quality surface geometry from these methods, and while the trade-offs
between surface and volume representations is an interesting research topic,
we focus on surface representations.

##### Light transport.

In addition to handling discontinuities, recent work also studies the
reduction of variance and memory consumption for Monte Carlo rendering
(Zeltner et al., 2021; Zhang et al., 2021a; Nimier-David et al., 2020; Vicini
et al., 2021). Earlier rendering work used derivatives for forward rendering
(Ward and Heckbert, 1992; Ramamoorthi et al., 2007; Li et al., 2015; Luan et
al., 2020). Our work is largely orthogonal to these.

##  3\. Method

Our method computes the correct gradient of a rendering function (i.e. the
pixel integral of the radiance function on the camera image plane) with
respect to geometric parameters, in the presence of primary visibility
discontinuities, for scenes where the geometry is represented by a signed
distance field $f$, parameterized by $\theta$ (e.g., the weights of neural
network). Our approach builds on Bangaru et al. (2020). We show how to extend
their warp function to SDFs in order to reparameterize an intractable boundary
integral. We summarize the necessary background in § 3.1. We then derive a
warp function for SDFs that is continuous and boundary consistent (§ 3.2) as
an integral along camera rays, and show how to compute it via quadrature using
sphere tracer points (§ 3.3). In Section 3.4, we finally give an unbiased
approximation for this warp that is tractable for use with neural SDFs, and we
show the criteria of unbiasedness and a sketch of proof. Section 3.5 provides
details on how to use our approach to solve inverse rendering problems.

###  3.1. Background: boundary-aware warping

Without loss of generality, assume a box pixel filter, so that
$\mathcal{U}\subset\mathbb{R}^{2}$ is the image plane region corresponding to
the pixel of interest. Let $L(u;\theta)$ denote the radiance along the ray
from $u\in\mathcal{U}$, a point on the image plane, and denote
$\theta\in\mathbb{R}^{N}$ the vector of geometric scene parameters (e.g.
neural network weights). In matrix expressions below, we will assume vector
quantities ($u$, $x$, $\theta$) to be row vectors, and gradients with respect
to $\theta$ to be column vectors.

We aim to compute the gradient of the rendering integral $I$ with respect to
parameters $\theta$:

(1)
$\partial_{\theta}I=\frac{\partial}{\partial\theta}\int_{\mathcal{U}}L(u;\theta)du.$

Primary visibility discontinuities make the radiance function non-
differentiable along occlusion boundaries (Fig. 3). Denoting
$\mathcal{U}_{\text{sil}}(\theta)\subset\mathcal{U}$ the set of object
silhouettes, for a point $u_{\text{sil}}\in\mathcal{U}_{\text{sil}}$, the
radiance $L(u_{\text{sil}};\theta)$ is discontinuous in $\theta$. This makes
naive automatic differentiation methods applied to the Monte Carlo sampling of
$I$ produce incorrect gradients since they ignore the Dirac delta that arises
from the differentiation.

Li et al.(2018; 2019) and Zhang et al. (2019) showed that Eq. (1) can be split
into two terms: an interior integral, for contributions away from the
discontinuities; and a boundary integral, along the discontinuities:

(2)
$\partial_{\theta}I=\int_{\mathcal{U}}\frac{\partial}{\partial\theta}L(u;\theta)du+I_{\text{sil}}.$

The second integral $I_{\text{sil}}$ is harder to compute because sampling the
boundary is generally technically difficult. This is particularly true for
SDFs whose surface boundaries admit no easy parametric form. We will not cover
boundary sampling in detail, since we will not use it; instead, we will use a
result from Bangaru et al. (2020), who showed, using the divergence theorem,
that this boundary term can be turned into an integral over the interior
$\mathcal{U}\setminus\mathcal{U}_{\text{sil}}(\theta)$, which is easier to
sample:

(3)
$I_{\text{sil}}=\int_{\mathcal{U}\setminus\mathcal{U}_{\text{sil}}(\theta)}\nabla_{u}\cdot\left(L(u;\theta)\mathcal{V}(u;\theta)\right)du.$

Here $\nabla_{u}\cdot$ is the divergence operator, and
$\mathcal{V}(u;\theta)\in\mathbb{R}^{N\times 2}$ is a warping function
required to satisfy two properties:

  1. (1)

continuity: $\mathcal{V}(\cdot;\theta)$ is continuous on $\mathcal{U}$, and

  2. (2)

boundary consistency: $\mathcal{V}$ agrees with the derivative of the
discontinuity points when $u$ approaches the discontinuity. That is,
$\lim_{u\to
u_{\text{sil}}}\mathcal{V}(u;\theta)=\partial_{\theta}u_{\text{sil}}$ for
$u_{\text{sil}}\in\mathcal{U}_{\text{sil}}(\theta)$.

Bangaru et al. further show that the area integral is equivalent to applying
the change of variable (Loubet et al., 2019) $u\mapsto
T(u,\theta)=u+(\theta-\theta_{0})\mathcal{V}(u;\theta)$ in Eq. (1), where
$\theta_{0}$ is the point at which the derivative is computed but
$\partial_{\theta}\theta_{0}=0$. Applying the reparameterization we obtain:

(4)
$\displaystyle\partial_{\theta}I=\int_{\mathcal{U}}\partial_{\theta}\left[L(T(u,\theta),\theta)\left|\det\left(\partial_{u}T(u,\theta)\right)\right|\right]\mathrm{d}u.$

Expanding $T$ and using Eq. (3), one can show that Eq. (4) indeed computes
$\partial_{\theta}I$. Intuitively, the reparameterization $T$ moves each point
on the boundary locally at the velocity of their derivatives, essentially
removing the discontinuities, while the determinant term accounts for the
change of measure.

The main goal of this paper is to derive a suitable form for
$\mathcal{V}(u;\theta)$ for SDFs, that can be tractably computed, so that we
can evaluate Eq. (4) using Monte Carlo estimation.

##### Rendering.

To render an SDF $f$ and compute $L(u;\theta)$, we need to find the closest
intersection point $x(u,t)\in\mathbb{R}^{3}$ such that $f(x;\theta)=0$, where
$t$ is the distance along the primary ray associated with pixel location $u$.
To find the intersection distance, we use sphere tracing (Hart, 1996), which
applies a fixed-point iteration to generate a sequence of points
$x_{n}\in\mathcal{T}(u)$, such that $\lim_{n\rightarrow\infty}x_{n}=x$.

###  3.2. Continuous boundary-consistent warp for SDFs

In this section, we construct an idealized warp function
$\mathcal{V}^{\text{int}}$ that satisfies the continuity and boundary-
consistency conditions of Section 3.1. First, we derive the boundary gradient
$\partial_{\theta}u_{\text{sil}}$ with which the warp should agree at
silhouette points (§ 3.2.1). We then smoothly extrapolate this gradient using
a weighted integral along the primary ray passing through $u$, to obtain our
warp function (§ 3.2.2). We show necessary and sufficient conditions on the
weighs to make the warp continuous and boundary-consistent (§ 3.2.3).

Figure 3.  As geometric parameters $\theta$ vary, visibility creates
discontinuities in the rendering function $L(u;\theta)$ which traditional
automatic differentiation cannot handle (left). These discontinuities
correspond to continuous changes in our SDF representation $f(x;\theta)$
(right). We compute the derivative $\partial_{\theta}u_{sil}$ of a pixel-space
silhouette point w.r.t. the geometry parameters $\theta$ by computing the
derivative $G(x;\theta)$ of the corresponding 3D scene point $x$, and
projecting it onto the screen space $\mathcal{U}$ through the inverse
Jacobian.

####  3.2.1. Boundary consistency for implicit functions

The boundary consistency condition in Section 3.1, requires that, at a
discontinuity point $u_{\text{sil}}$ the warp agrees with
$\partial_{\theta}u_{\text{sil}}$. The derivation proposed by Bangaru et al.
(2020) does not apply directly to implicit surfaces, so we derive this
boundary derivative using the implicit function theorem. Specifically, the
derivative of a scene point $x\in\mathbb{R}^{3}$ on the surface, i.e.,
$f(x;\theta)=0$, w.r.t. parameters $\theta\in\mathbb{R}^{N}$ is given by:

(5)
$G(x;\theta)\vcentcolon=\partial_{\theta}x=-\frac{\partial_{\theta}f\partial_{x}f^{T}}{\|\partial_{x}f\|^{2}}\in\mathbb{R}^{N\times
3}.$

The above directly follows from the implicit function theorem applied to
$f(x;\theta)=0$. To get the derivative in pixel coordinates
$\partial_{\theta}u=\partial_{\theta}x\cdot\partial_{x}u\in\mathbb{R}^{N\times
2}$, we need to project this derivative by the Jacobian
$\partial_{x}u\in\mathbb{R}^{3\times 2}$, which for a perspective camera can
be easily derived by hand. For more generality, we can obtain this Jacobian as
the pseudo-inverse $\dagger$ of the forward Jacobian:

(6) $\partial_{x}u=(\partial_{u}x(u,t))^{\dagger}.$

Taken together, the derivative at a silhouette point $u_{\text{sil}}$, with
corresponding 3D position $x_{\text{sil}}=x(u_{\text{sil}},t_{\text{sil}})$,
is then:

(7) $\partial_{\theta}u_{\text{sil}}=G(x_{\text{sil}};\theta)\,\partial_{x}u.$

Figure 3 illustrates the geometric configuration.

####  3.2.2. Extending to a smooth warp $\mathcal{V}^{\text{int}}(u;\theta)$
by integration along the ray

Now that we have an expression for the warp at silhouette points, we extend it
to all points, by smoothing this term in a consistent manner. Our method takes
advantage of the fact that our implicit SDF $f(x;\theta)$ is continuous in 3D
space and achieves the smoothing by convolving _along the ray_ (Fig. 4 (b)).
This avoids casting expensive additional rays which are needed by Bangaru et
al. (2020), and also propagates gradients to points in free space near the
boundary points. This can have a stabilizing effect on the optimization of
neural SDFs, as noted by Oechsle et al. (2021) and Wang et al. (2021); note
that, while they adapt a volumetric rendering model to achieve better
convergence, we do so while computing correct boundary gradients for a
surface-based representation.

Our proposed warp function smoothly extends Eq. (7) to non-boundary points as
follows:

(8)
$\mathcal{V}^{\text{int}}(u;\theta)=\frac{\int^{t_{0}}_{t=0}w\left(x(u,t)\right)\,G(x;\theta)\,\partial_{x}u\,dt}{\int^{t_{0}}_{t=0}w\left(x(u,t)\right)\,dt},$

with $t_{0}$ the distance to the closest intersection, $t_{0}=\infty$ when the
ray does not intersect.

####  3.2.3. Choice of weights

In order to satisfy the boundary consistency criteria, the weights need to
asymptotically satisfy the limit:

(9) $\lim_{u\to
u_{\text{sil}}}\frac{w\left(x(u,t)\right)}{\int^{t_{0}}_{t^{\prime}=0}w\left(x(u,t^{\prime})\right)\,dt^{\prime}}=\delta(t-t_{\text{sil}}),$

where $\delta$ is the Dirac delta operator.

From Eq. (9), we see that our weights have to depend on some notion of
distance to the silhouette. For an implicit function $f$ that is at least
$C_{1}$ continuous, the following constraints implicitly characterize the
silhouette points (Gargallo et al., 2007):

(10) $\displaystyle f(x(u,t);\theta)=0,$
$\displaystyle\partial_{x}f(x(u,t);\theta)^{T}\partial_{t}x(u,t)=0.$

The first condition requires the point to be on the surface, and the second
condition requires the SDF gradient to be perpendicular to ray direction
(Hertzmann, 1999). We can use these equations to build a silhouette
characteristic function $\mathcal{S}(x)$, which takes value $0$ whenever $x$
is a silhouette point, and is continuous everywhere. Specifically, we define:

(11)
$\mathcal{S}(x)=\left|f(x;\theta)\right|+\lambda_{d}\,\left|\partial_{x}f(x;\theta)^{T}\partial_{t}x\right|,$

where $\lambda_{d}>0$. This characteristic function is similar to the boundary
test function used by Bangaru et al. (2020) for meshes. However, unlike their
boundary test, $\mathcal{S}(x)$ is defined everywhere in the SDF’s 3D domain,
not just the surface points. This allows us to use these weights for our
integral along any ray.

Our final harmonic weights are given by:

(12) $w(x)=\mathcal{S}(x)^{-\gamma},\gamma>2.$

For $\gamma>2$, our weights satisfy the limit in Eq. (9). Intuitively, this is
because the $w(x)\to\delta(x-x_{\text{sil}})$ as $x\to x_{\text{sil}}$. See
our supplementary material for a discussion of correctness, and derivation of
$\gamma>2$.

Figure 4. (a) Warped-area sampling uses additional random rays around the
primal ray and averages $\partial_{\theta}{x(u;t)}$ using boundary-aware
harmonic weights. (b) Our method instead takes a weighted average along the
ray, repurposing the SDF $f$ into weights $S(x)^{-\gamma}$. Figure 5. Weight
visualization. A contour plot of a sample 2D SDF (first row). We use an
orthographic camera for illustration, so camera rays are parallel to the
horizontal axis. We show our three weighting schemes in unnormalized (second
row) and normalized (third row) form. Our proposed harmonic weights (a) for
$\gamma=4.0$, $\lambda_{d}=1e-1$ are well approximated by a trapezoidal
quadrature on the sphere tracer points (b). The blank regions with no weight
can be excluded from the computation, which leads to our proposed our proposed
top-$k$ subset weights (c), for $k=8$. This reduces both the compute and
memory burden of the backward pass. We visualize the weight in a _symlog_
plot, values are linear in $[0,10^{1}]$ and $[0,10^{-3}]$ for the unnormalized
and normalized weights, respectively.

Fig. 5(a) shows our weight distribution along the ray for all $u$ in a 1D
example sphere tracer.

###  3.3. Estimating the warp through its quadrature $\mathcal{V}^{q}$

We now have a clear form for our warp function that can be used to
reparameterize and differentiate the rendering function. Unfortunately, the
asymptotical sharpness of our weights required to obtain a valid warp, also
makes the integral (8) very difficult to sample. For $u$ close to the
silhouette $\mathcal{U}_{\text{sil}}$, the weights become very concentrated
near the surface boundary, presenting a tricky integrand if we were to
uniformly sample along the ray.

Careful importance sampling of areas near the boundary could remedy this, but
there is unfortunately no straightforward way to implement this: the weight
distribution depends heavily on the configuration of silhouettes near $u$,
dictated by the SDF.

Our approach foregoes stochastic sampling altogether. We construct a
trapezoidal quadrature on the series of intermediate points
$x_{n}\in\mathcal{T}(u)$ generated by the sphere tracer, shown in Fig. 5(b).
This quadrature estimator for the warp is given by:

(13)
$\displaystyle\mathcal{V}^{q}(u;\theta)=\frac{\sum_{x_{i}\in\mathcal{T}(u)}w^{q}\left(x_{i}\right)\,G(x;\theta)\,\partial_{x}u}{\sum_{x_{i}\in\mathcal{T}(u)}w^{q}\left(x_{i}\right)},$
$\displaystyle\text{where }w^{q}(x_{i})=w(x_{i})\frac{(t_{i-1}-t_{i+1})}{2},$

and $t_{i}$ is the distance along the ray to sphere tracer point $x_{i}$.
Assuming the underlying SDF $f(x;\theta)$ is $C_{1}$ continuous, the
intermediate points of the sphere tracer are continuous at all
$u\not\in\mathcal{U}_{\text{sil}}$. By composition of continuous functions,
$\mathcal{V}^{q}(\cdot;\theta)$ is also continuous.

Our quadrature warp $\mathcal{V}^{q}$ satisfies the continuity and boundary
consistency condition (§3.1). Since we apply trapezoidal quadrature,
$\mathcal{V}^{q}(u;\theta)$ is in general a biased estimator of integral
$\mathcal{V}^{\text{int}}(u;\theta)$. However, the two terms are equal in the
limit as $u$ approaches the silhouette, i.e., for
$u_{\text{sil}}\in\mathcal{U}_{\text{sil}}$, $\lim_{u\to
u_{\text{sil}}}\mathcal{V}^{q}(u;\theta)=\lim_{u\to
u_{\text{sil}}}\mathcal{V}^{\text{int}}(u;\theta)$, and since the right-hand
side is boundary consistent, so is our quadrature warp $\mathcal{V}^{q}$. See
supplemental for a sketch proof of correctness.

###  3.4. Top-$k$ subset weighting $\bar{w}_{k}$ to reduce memory use

For complex SDFs such as a neural network, our quadrature warp
$\mathcal{V}^{\text{q}}$ has the caveat that it requires back-propagating
through _every_ sphere tracer point. Previous work like IDR (Yariv et al.,
2020) do not have this issue since their (biased) gradient is only computed at
the intersection point, and they exclude other points from the gradient
computation. Our approach, on the other hand, uses a weighted sum, so we
cannot discard intermediate points.

However, as shown in Fig. 5(b), the vast majority of sphere tracer points have
negligible weight, and most of the mass is concentrated close to the
silhouette. We exploit this by only using the subset of points with the
highest weight in our warp estimation. That is, instead of using all of
$T(u)$, we can instead use a top-k subset $T_{k}(u)$. Selecting the top-$k$
weights requires adjusting them to ensure that they remain continuous. For a
subset size of $k$, our weights are222Note that even though we consider the
top $k$ weights, only $k-1$ weights actually have a non-zero contribution to
the $\theta$-derivative.:

(14) $\bar{w}_{k}(x)=\begin{cases}w^{\text{q}}(x)-\min\limits_{x\in
T_{k}(u)}w^{\text{q}}(x),&\text{if }x\in T_{k}(u)\\\
0&\text{otherwise.}\end{cases}$

The weights $\bar{w}_{k}(x)$ still produce a continuous warp field (see
supplemental for a sketch of proof). Intuitively, even though the set of
points change as a function of $u$, whenever this change occurs, the points
that swap in or out of the set always have weight 0.

###  3.5. Inverse Rendering Details

In this section, we briefly discuss some details that make our inverse
rendering pipeline tractable.

##### Implementation

Our method requires 3 nested derivative passes to (i) compute normals
$\partial_{\mathbf{x}}f$, (ii) compute Jacobian of the transformation
$\partial_{u}T$ and (iii) to compute derivatives of the full pipeline
$\partial_{\theta}\left[L(T(u,\theta))\left|\det\left(\partial_{u}T(u,\theta)\right)\right|\right]$.
We use the Python JAX automatic differentiation system (Bradbury et al.,
2018), which supports nested forward+backward differentiation. We use forward-
mode for (i) and (ii), and reverse-mode for (iii).

##### Network architecture

For our inverse rendering results, we use the network architecture shown in
Fig. 11 of Yariv et al. (2020). Since our method is slightly more memory-
intensive (even with top-$k$ subset weights), we reduce the width of the SDF
network to 256 channels per layer. In this architecture, the shading network
predicts the final radiance based on the position, viewing direction and a
geometric feature vector. However, in contrast to NeRF (Mildenhall et al.,
2020)-like methods, the shading network is only evaluated at _surface_ points.
We use 6-levels of positional encoding on the input position $x$ to allow the
network to reconstruct fine geometry.

##### Pixel sampling

Similar to Yariv et al. (2020) and other neural methods, we sample a subset of
pixels for each iteration since it can be computationally prohibitive to trace
the entire image when using a deep neural representation. However, unlike
Yariv et al. (2020), which works with a single ray at the center of the pixel,
our approach must integrate the spatially-varying warp $\mathcal{V}$ over each
pixel. We achieve this by Monte-Carlo sampling within each pixel. Appendix A
discusses how we incorporate pixel filters.

##### Multi-level optimization

Since we only use a subset of pixels, the likelihood of sampling a pixel with
silhouette gradient is fairly low. For unbiased derivatives, only pixels that
are partially covered by a surface have a non-zero boundary contribution. This
is in contrast to approximate derivatives (e.g., (Liu et al., 2019), (Yariv et
al., 2020)) that have a wider spatial footprint. To alleviate this issue, we
use a multi-scale pyramid of the target image throughout our optimization to
help with gradient stability.

##### Initialization

We use the geometric network initialization (Atzmon and Lipman, 2020) which
approximately produces a spherical SDF. We also initializes the weights of the
positional encoding layer to 0 (Yariv et al., 2020). We found this subtle
modification implicitly enforces a coarse-to-fine mechanism that yields
significantly better generalization to novel views.

##### Eikonal constraint

We represent our SDF $f$ using a neural network, which does not necessarily
satisfy the distance property. We adopt the Eikonal regularization loss (Gropp
et al., 2020) to explicitly enforce this. In-spite of the additional loss, $f$
is still an approximation of an SDF, and therefore we pad our weights with a
small $\epsilon$ in order to avoid infinities.

##  4\. Results

Figure 6. Gradient quality. We compare the image gradients computed naïvely
without reparameterization (Yariv et al., 2020) and with our method against
the “ground truth” gradient computed with finite differences for three scenes.
Our method properly handles boundary discontinuities both due to object edges
(in purple insets) and self-occlusions (in green insets).

###  4.1. Ground truth gradient comparions

We first evaluate the correctness of our gradient by visualizing gradients on
three different scenes (illustrated in Fig. 6). For Torus—a analytical torus
model textured with a diffuse Perlin noise albedo—we visualize the gradients
w.r.t the outer radius (distance from the center to the center of the ring).
(Santa and Kitty) are 3D models that we represent as neural SDFs. We take the
parameters of the neural SDF from an intermediate iteration during an inverse
rendering optimization, and visualize the gradient w.r.t the bias parameter of
the last layer output (i.e. the level set perturbation). We also compute the
gradient without reparameterization; this is similar to the gradient used in
previous SDF-based inverse rendering methods (Yariv et al., 2020). Note that
the interior gradient is largely unaffected by reparameterization, with the
gradient at the silhouettes being the largest benefit of our method,
especially at self-occlusions. In the next subsection, we show that this
boundary gradient is critical and without it, the inverse rendering diverges.

santa |  |  |  |  |  |
---|---|---|---|---|---|---
kitty |  |  |  |  |  |
duck |  |  |  |  |  |
pony |  |  |  |  |  |
dragon |  |  |  |  |  |
| Input image | IDR Radiance | IDR Depth | Ours Radiance | Ours Depth | Ours 3D geometry
Figure 7. Neural SDF reconstruction. We compare with IDR (Yariv et al., 2020)
on three synthetic scenes (top three rows) and two real captured scene (bottom
rows). IDR requires 2D mask supervision, without which it completely diverges.
Thanks to our accurate gradient computation, our reconstructions are on par
with IDR, without requiring any additional supervision beyond the input images
and camera poses. In fact, on the real scenes, our reconstructions (without
masks) outperform IDR _with_ masks (see the head of the pony, or the legs,
tails and wings of the dragon) because of errors in automatic 2D segmentation.

###  4.2. Comparisons against IDR

We compare our reconstructions against the SDF-based inverse rendering method
of IDR (Yariv et al., 2020). IDR does not correctly account for the boundary
term of gradient of the rendering integral and requires additional
supervision, in the form of accurate 2D segmentation masks. We implement IDR
in our pipeline to ensure that the only difference is our reparameterization.
We use the same network architecture for both methods (See Sec. 3.5 for
details), and report results after roughly 25,000 network updates. Note that
our method uses more samples ($2$ in the interior + $4$ on each pixel
boundary) since we use a Monte-Carlo approach to estimate the warp. IDR only
requires one sample, fixed at the center of the pixel.

Figure 7 shows that, on three synthetic scenes (Santa, Kitty and Duck), our
method _without_ any 2D masks supervision obtains comparable depth and RGB
reconstruction as IDR _with_ (perfect) mask supervision. We also show
reconstructions of a captured real scene (Pony from Bi et al. (2020b)). Here,
we provide IDR with 2D masks derived from a COLMAP reconstruction, which has
errors. As a result, our reconstruction outperforms IDR on this scene.

We also tried to compare with IDR without mask supervision. In most cases, IDR
without masks diverges completely because of the lack of gradients from the
silhouette. This is similar to the observation made by Oechsle et al. (2021).

###  4.3. Ablation study: Subset size

Our top-k weighting scheme reduces the memory footprint of our optimization,
but this comes at a cost. The smaller $k$, the sharper the weight landscape.
This can cause high variance that can impede the optimization of fine details.
We explore this through an ablation study on the Santa dataset, varying $k$
shown below. We use 36 views for this study, and report results after $20,000$
network updates. Details are resolved for $k\geq 14$.

|  |  |
---|---|---|---
$k=4$ | $k=8$ | $k=15$ | $k=22$

##  5\. Conclusion

We have presented a novel method to correctly differentiate neural SDFs.
Unlike prior work that relies on accurate masks or biased approximations of
the boundary gradients, we reparameterize the pixel coordinates with
discontinuity-aware transformation such that geometric derivatives exist
everywhere and can be correctly computed via standard automatic
differentiation. We have validated the correctness of our approach by
comparing our gradients with finite difference ground truth, and demonstrated
superior optimization convergence comparing state-of-the-art neural SDF
baseline.

While we have focused on primary visibility in this work, our formulation can
be extended to global light transport. In particular, we expect to be able to
model light rays and jointly optimize for geometry (represented as an SDF) as
well as surface reflectance (instead of the radiance we are currently
reconstructing) and illumination. Modeling full global light transport
(interreflections) with neural SDFs will require extensions or approximations
to be computationally tractable. Finally, inverse rendering under unknown,
natural illumination is ill-posed and it would be interesting to explore
geometry, material and illumination priors that can be combined with our
differentiable rendering formulation.

###### Acknowledgements.

This work was partially completed during an internship at Adobe Research and
subsequently funded by the Toyota Research Institute and the National Science
Foundation (NSF 2105806).

## References

  * (1)
  * Arvo (1994) James Arvo. 1994\.  The Irradiance Jacobian for Partially Occluded Polyhedral Sources. In _SIGGRAPH_. ACM Press/Addison-Wesley Publishing Co., 343–350.
  * Atzmon and Lipman (2020) Matan Atzmon and Yaron Lipman. 2020.  SAL: Sign Agnostic Learning of Shapes From Raw Data. In _IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_.
  * Bangaru et al. (2020) Sai Bangaru, Tzu-Mao Li, and Frédo Durand. 2020.  Unbiased Warped-Area Sampling for Differentiable Rendering.  _ACM Trans. Graph._ 39, 6 (2020), 245:1–245:18.
  * Bi et al. (2020a) Sai Bi, Zexiang Xu, Pratul P. Srinivasan, Ben Mildenhall, Kalyan Sunkavalli, Milos Hasan, Yannick Hold-Geoffroy, David J. Kriegman, and Ravi Ramamoorthi. 2020a.  Neural Reflectance Fields for Appearance Acquisition.  _CoRR_ abs/2008.03824 (2020).
  * Bi et al. (2020b) Sai Bi, Zexiang Xu, Kalyan Sunkavalli, Milos Hasan, Yannick Hold-Geoffroy, David Kriegman, and Ravi Ramamoorthi. 2020b.  Deep Reflectance Volumes: Relightable Reconstructions from Multi-View Photometric Images. In _Proceedings of the European Conference on Computer Vision (ECCV)_.
  * Blinn (1982) James F Blinn. 1982\.  A generalization of algebraic surface drawing.  _ACM Trans. Graph._ 1, 3 (1982), 235–256.
  * Bradbury et al. (2018) James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang. 2018\.  _JAX: composable transformations of Python+NumPy programs_.  http://github.com/google/jax
  * Cole et al. (2021) Forrester Cole, Kyle Genova, Avneesh Sud, Daniel Vlasic, and Zhoutong Zhang. 2021.  Differentiable surface rendering via non-differentiable sampling. In _International Conference on Computer Vision_. 6088–6097.
  * de La Gorce et al. (2011) Martin de La Gorce, David J Fleet, and Nikos Paragios. 2011\.  Model-based 3D hand pose estimation from monocular video.  _IEEE Trans. Pattern Anal. Mach. Intell._ 33, 9 (2011), 1793–1805.
  * Gargallo et al. (2007) Pau Gargallo, Emmanuel Prados, and Peter Sturm. 2007\.  Minimizing the reprojection error in surface reconstruction from images. In _International Conference on Computer Vision_. 1–8.
  * Gkioulekas et al. (2013) Ioannis Gkioulekas, Shuang Zhao, Kavita Bala, Todd Zickler, and Anat Levin. 2013.  Inverse Volume Rendering with Material Dictionaries.  _ACM Trans. Graph._ 32, 6 (nov 2013), 162:1–162:13.
  * Gropp et al. (2020) Amos Gropp, Lior Yariv, Niv Haim, Matan Atzmon, and Yaron Lipman. 2020.  Implicit Geometric Regularization for Learning Shapes.  In _Proceedings of Machine Learning and Systems 2020_. 3569–3579.
  * Hart (1996) John C Hart. 1996\.  Sphere tracing: A geometric method for the antialiased ray tracing of implicit surfaces.  _The Visual Computer_ 12, 10 (1996), 527–545.
  * Hertzmann (1999) Aaron Hertzmann. 1999\.  Introduction to 3D Non-Photorealistic Rendering: Silhouettes and Outlines.  In _SIGGRAPH Course Notes. Course on Non-Photorelistic Rendering_ , Stuart Green (Ed.). ACM Press/ACM SIGGRAPH, New York.
  * Izadi et al. (2011) Shahram Izadi, David Kim, Otmar Hilliges, David Molyneaux, Richard Newcombe, Pushmeet Kohli, Jamie Shotton, Steve Hodges, Dustin Freeman, Andrew Davison, et al. 2011\.  KinectFusion: real-time 3D reconstruction and interaction using a moving depth camera. In _symposium on User interface software and technology_. 559–568.
  * Jiang et al. (2020) Yue Jiang, Dantong Ji, Zhizhong Han, and Matthias Zwicker. 2020\.  SDFDiff: Differentiable rendering of signed distance fields for 3D shape optimization. In _Computer Vision and Pattern Recognition_. 1251–1261.
  * Kato et al. (2018) Hiroharu Kato, Yoshitaka Ushiku, and Tatsuya Harada. 2018\.  Neural 3D Mesh Renderer. In _Computer Vision and Pattern Recognition_. IEEE, 3907–3916.
  * Kellnhofer et al. (2021) Petr Kellnhofer, Lars C Jebe, Andrew Jones, Ryan Spicer, Kari Pulli, and Gordon Wetzstein. 2021\.  Neural lumigraph rendering. In _Computer Vision and Pattern Recognition_. 4287–4297.
  * Li (2019) Tzu-Mao Li. 2019\.  _Differentiable Visual Computing_.  Ph.D. Dissertation. Massachusetts Institute of Technology.  Advisor(s) Durand, Frédo.
  * Li et al. (2018) Tzu-Mao Li, Miika Aittala, Frédo Durand, and Jaakko Lehtinen. 2018.  Differentiable Monte Carlo Ray Tracing through Edge Sampling.  _ACM Trans. Graph. (Proc. SIGGRAPH Asia)_ 37, 6 (2018), 222:1–222:11.
  * Li et al. (2015) Tzu-Mao Li, Jaakko Lehtinen, Ravi Ramamoorthi, Wenzel Jakob, and Frédo Durand. 2015\.  Anisotropic Gaussian Mutations for Metropolis Light Transport through Hessian-Hamiltonian Dynamics.  _ACM Trans. Graph. (Proc. SIGGRAPH Asia)_ 34, 6 (2015), 209:1–209:13.
  * Liu et al. (2020) Lingjie Liu, Jiatao Gu, Kyaw Zaw Lin, Tat-Seng Chua, and Christian Theobalt. 2020.  Neural Sparse Voxel Fields.  _NeurIPS_ (2020).
  * Liu et al. (2019) Shichen Liu, Tianye Li, Weikai Chen, and Hao Li. 2019\.  Soft Rasterizer: A Differentiable Renderer for Image-based 3D Reasoning.  _International Conference on Computer Vision_ (2019).
  * Lombardi et al. (2019) Stephen Lombardi, Tomas Simon, Jason Saragih, Gabriel Schwartz, Andreas Lehrmann, and Yaser Sheikh. 2019\.  Neural Volumes: Learning Dynamic Renderable Volumes from Images.  _ACM Trans. Graph._ 38, 4, Article 65 (July 2019), 65:1–65:14 pages.
  * Loper and Black (2014) Matthew M. Loper and Michael J. Black. 2014.  OpenDR: An Approximate Differentiable Renderer. In _European Conference on Computer Vision_ , Vol. 8695. ACM, 154–169.
  * Loubet et al. (2019) Guillaume Loubet, Nicolas Holzschuch, and Wenzel Jakob. 2019\.  Reparameterizing discontinuous integrands for differentiable rendering.  _ACM Trans. Graph. (Proc. SIGGRAPH Asia)_ 38, 6 (2019), 228\.
  * Luan et al. (2020) Fujun Luan, Shuang Zhao, Kavita Bala, and Ioannis Gkioulekas. 2020\.  Langevin Monte Carlo Rendering with Gradient-based Adaptation.  _ACM Trans. Graph. (Proc. SIGGRAPH)_ 39, 4 (2020).
  * Mildenhall et al. (2020) Ben Mildenhall, Pratul P. Srinivasan, Matthew Tancik, Jonathan T. Barron, Ravi Ramamoorthi, and Ren Ng. 2020\.  NeRF: Representing Scenes as Neural Radiance Fields for View Synthesis. In _European Conference on Computer Vision_.
  * Niemeyer et al. (2020) Michael Niemeyer, Lars Mescheder, Michael Oechsle, and Andreas Geiger. 2020.  Differentiable volumetric rendering: Learning implicit 3D representations without 3D supervision. In _Computer Vision and Pattern Recognition_. 3504–3515.
  * Nimier-David et al. (2020) Merlin Nimier-David, Sébastien Speierer, Benoît Ruiz, and Wenzel Jakob. 2020.  Radiative Backpropagation: An Adjoint Method for Lightning-Fast Differentiable Rendering.  _ACM Trans. Graph. (Proc. SIGGRAPH)_ 39, 4 (July 2020).  https://doi.org/10.1145/3386569.3392406
  * Oechsle et al. (2021) Michael Oechsle, Songyou Peng, and Andreas Geiger. 2021\.  UNISURF: Unifying Neural Implicit Surfaces and Radiance Fields for Multi-View Reconstruction. In _International Conference on Computer Vision (ICCV)_.
  * Park et al. (2019) Jeong Joon Park, Peter Florence, Julian Straub, Richard Newcombe, and Steven Lovegrove. 2019\.  DeepSDF: Learning continuous signed distance functions for shape representation. In _Computer Vision and Pattern Recognition_. 165–174.
  * Ramamoorthi et al. (2007) Ravi Ramamoorthi, Dhruv Mahajan, and Peter Belhumeur. 2007\.  A First-order Analysis of Lighting, Shading, and Shadows.  _ACM Trans. Graph._ 26, 1 (2007), 2.
  * Remelli et al. (2020) Edoardo Remelli, Artem Lukoianov, Stephan Richter, Benoît Guillard, Timur Bagautdinov, Pierre Baque, and Pascal Fua. 2020.  MeshSDF: Differentiable iso-surface extraction. In _Advances in Neural Information Processing Systems_ , Vol. 33. 22468–22478.
  * Rhodin et al. (2015) Helge Rhodin, Nadia Robertini, Christian Richardt, Hans-Peter Seidel, and Christian Theobalt. 2015\.  A Versatile Scene Model with Differentiable Visibility Applied to Generative Pose Estimation. In _International Conference on Computer Vision_. IEEE, 765–773.
  * Veach (1998) Eric Veach. 1998\.  _Robust Monte Carlo Methods for Light Transport Simulation_.  Ph.D. Dissertation. Stanford University.  Advisor(s) Guibas, Leonidas J.
  * Vicini et al. (2021) Delio Vicini, Sébastien Speierer, and Wenzel Jakob. 2021.  Path replay backpropagation: differentiating light paths using constant memory and linear time.  _ACM Trans. Graph. (Proc. SIGGRAPH)_ 40, 4 (2021).
  * Wang et al. (2021) Peng Wang, Lingjie Liu, Yuan Liu, Christian Theobalt, Taku Komura, and Wenping Wang. 2021\.  NeuS: Learning Neural Implicit Surfaces by Volume Rendering for Multi-view Reconstruction. In _Advances in Neural Information Processing Systems_.
  * Ward and Heckbert (1992) Greg Ward and Paul Heckbert. 1992.  Irradiance Gradients. In _Eurographics Workshop on Rendering_. Eurographics Association, 85–98.
  * Xie et al. (2022) Yiheng Xie, Towaki Takikawa, Shunsuke Saito, Or Litany, Shiqin Yan, Numair Khan, Federico Tombari, James Tompkin, Vincent Sitzmann, and Srinath Sridhar. 2022\.  Neural Fields in Visual Computing and Beyond.  _Comput. Graph. Forum. State of the Art Report_ (2022).
  * Yariv et al. (2021) Lior Yariv, Jiatao Gu, Yoni Kasten, and Yaron Lipman. 2021\.  Volume rendering of neural implicit surfaces. In _Advances in Neural Information Processing Systems_.
  * Yariv et al. (2020) Lior Yariv, Yoni Kasten, Dror Moran, Meirav Galun, Matan Atzmon, Basri Ronen, and Yaron Lipman. 2020.  Multiview Neural Surface Reconstruction by Disentangling Geometry and Appearance.  _Advances in Neural Information Processing Systems_ 33 (2020).
  * Zeltner et al. (2021) Tizian Zeltner, Sébastien Speierer, Iliyan Georgiev, and Wenzel Jakob. 2021.  Monte Carlo estimators for differential light transport.  _ACM Trans. Graph. (Proc. SIGGRAPH)_ 40, 4 (2021).
  * Zhang et al. (2021a) Cheng Zhang, Zhao Dong, Michael Doggett, and Shuang Zhao. 2021a.  Antithetic Sampling for Monte Carlo Differentiable Rendering.  _ACM Trans. Graph. (Proc. SIGGRAPH)_ 40, 4 (2021).
  * Zhang et al. (2020) Cheng Zhang, Bailey Miller, Kai Yan, Ioannis Gkioulekas, and Shuang Zhao. 2020.  Path-Space Differentiable Rendering.  _ACM Trans. Graph. (Proc. SIGGRAPH)_ 39, 6 (2020), 143:1–143:19.
  * Zhang et al. (2019) Cheng Zhang, Lifan Wu, Changxi Zheng, Ioannis Gkioulekas, Ravi Ramamoorthi, and Shuang Zhao. 2019\.  A differential theory of radiative transfer.  _ACM Trans. Graph. (Proc. SIGGRAPH Asia)_ 38, 6 (2019), 227\.
  * Zhang et al. (2021c) Cheng Zhang, Zihan Yu, and Shuang Zhao. 2021c.  Path-Space Differentiable Rendering of Participating Media.  _ACM Trans. Graph. (Proc. SIGGRAPH)_ 40, 4 (2021), 76:1–76:15.
  * Zhang et al. (2021b) Kai Zhang, Fujun Luan, Qianqian Wang, Kavita Bala, and Noah Snavely. 2021b.  PhySG: Inverse rendering with spherical gaussians for physics-based material editing and relighting. In _Computer Vision and Pattern Recognition_. 5453–5462.
  * Zhao et al. (2020) Shuang Zhao, Wenzel Jakob, and Tzu-Mao Li. 2020.  Physics-Based Differentiable Rendering: From Theory to Implementation. In _ACM SIGGRAPH 2020 Courses_. Article 14, 30 pages.
  * Zhou et al. (2021) Yang Zhou, Lifan Wu, Ravi Ramamoorthi, and Ling-Qi Yan. 2021\.  Vectorization for Fast, Analytic, and Differentiable Visibility.  _ACM Trans. Graph._ 40, 3 (2021).

##  Appendix A Pixel boundary sampling

Figure 8.  On the left, existing methods (Bangaru et al., 2020; Loubet et al.,
2019) use an unbounded Gaussian filter to avoid the need to handle the
boundary of the pixel filter $\mathcal{U}_{b}$, but this comes at the cost of
increased variance in the interior due to the derivative of the Gaussian
weights. On the right, our method uses a box filter and explicitly removes the
discrepancy in the warp field $\mathcal{V}$ through a boundary integral over
$\mathcal{U}_{b}$.

In Eqn. 3, we use a box which implies the pixel domain $\mathcal{U}$ is
bounded. An implication of this is that we must also consider the boundaries
of the pixel filter support (denoted by
$\mathcal{U}_{b}\subset\mathtt{R}^{2}$) as discontinuities in
$\mathcal{U}_{\text{sil}}$. Previously, to avoid this additional complexity,
Bangaru et al. (2020) used a Gaussian filter that has infinite support. We
have found that this introduces extra variance due to the variation in the
pixel filters in the divergence.

We instead keep the box filter as well as exclude the pixel boundary from the
area integral $\mathcal{U}_{\text{sil}}$. This means that Eqn. 3 is no longer
valid since the product $(L\mathcal{V})$ does not vanish smoothly at the pixel
filter boundary $\mathcal{U}_{\text{b}}$. We must instead rewrite the integral
domain as an unbounded space $\mathcal{U}_{\infty}$. We can further split the
unbounded integral into two parts, one inside the pixel filter domain
$\mathcal{U}$ and one outside (we omit parentheses here for brevity)

(15)
$I_{\text{sil}}=\int_{\mathcal{U}\setminus\mathcal{U}_{\text{sil}}}\nabla\cdot\left(L\mathcal{V}\right)+\int_{(\mathcal{U}_{\infty}\setminus\mathcal{U})\setminus\mathcal{U}_{\text{b}}}\nabla\cdot\left(L\mathcal{V}\right).$

We can then use the divergence theorem on the second area integral to turn it
into a boundary integral over $\mathcal{U}_{\text{b}}$

(16)
$I_{\text{sil}}=\int_{\mathcal{U}\setminus\mathcal{U}_{\text{sil}}}\nabla\cdot\left(L\mathcal{V}\right)-\oint_{\mathcal{U}_{\text{b}}}L\left(\mathcal{V}\cdot\mathbf{n}_{b}\right),$

where $\mathbf{n}_{b}$ is the outward pointing normal of the pixel filter
boundary, and the negative sign comes from the fact that we consider regions
outside of the pixel filter instead of inside. Unlike silhouette boundaries in
$\mathcal{U}_{\text{sil}}$, $\mathcal{U}_{b}$ is easy to sample since it only
contains axis-aligned line segments of equal length. Fig. 8 illustrates the
difference between using a smooth unbounded filter and using a box filter with
pixel boundary sampling.