MainReferences SuppSupplementary References

# Towards 3D Vision with Low-Cost Single-Photon Cameras

Fangzhou Mu†, Carter Sifferman†, Sacha Jungerman, Yiquan Li‡, Mark Han‡,
Michael Gleicher, Mohit Gupta, Yin Li
† co-first author, ‡ equal contribution
University of Wisconsin-Madison

###### Abstract

We present a method for reconstructing 3D shape of arbitrary Lambertian
objects based on measurements by miniature, energy-efficient, low-cost single-
photon cameras. These cameras, operating as time resolved image sensors,
illuminate the scene with a very fast pulse of diffuse light and record the
shape of that pulse as it returns back from the scene at a high temporal
resolution. We propose to model this image formation process, account for its
non-idealities, and adapt neural rendering to reconstruct 3D geometry from a
set of spatially distributed sensors with known poses. We show that our
approach can successfully recover complex 3D shapes from simulated data. We
further demonstrate 3D object reconstruction from real-world captures,
utilizing measurements from a commodity proximity sensor. Our work draws a
connection between image-based modeling and active range scanning, and offers
a step towards 3D vision with single-photon cameras. Our project webpage is at
https://cpsiff.github.io/towards_3d_vision/.

Figure 1: We demonstrate that measurements from spatially distributed low-cost
single-photon proximity sensors (left) can be used to reconstruct 3D shape of
real world objects (right). Our method combines a differentiable image
formation model and neural rendering to recover 3D geometry based on
measurements (transient histograms) from sensors with known poses. This is
done by minimizing the difference between the observed and rendered sensor
measurements. For clarity, a subset of sensor poses and measurements are
shown.

##  1 Introduction

Reconstructing 3D shape of real objects remains a central problem in vision,
solutions to which have evolved into two parallel branches. Image-based
modeling szeliski2022cv leverages a plethora of visual cues from multiple
photographs (_e.g_., stereo, motion, shading), leading to problems including
multi-view stereo seitz2006mvs, photometric stereo ackermann2015ps and the
more recent neural radiance fields (NeRF) Mildenhall2020nerf. Conversely,
active range scanning jarvis1983range combines an active light source with an
imaging sensor, giving rise to imaging techniques such as structured light
geng2011sl, and time-of-flight hansard2012tof. Conventional wisdom suggests
that range scanning yields more precise 3D geometry than image-based modeling
at the cost of using specialized, expensive hardware.

An emerging approach for range scanning is direct time-of-flight imaging with
active _single-photon cameras_ , a form of time-resolved image sensor. This
approach couples a pico-to-nanosecond detector with a fast coherent light
source, illuminates the scene with a very short pulse of light, and measures
the intensity of the light over time as it reflects back from the scene. The
resulting incident wavefront is recorded and quantized, forming a transient
histogram. A special case of this approach is single-photon LiDAR, in which
the light source (laser) is highly focused, the detector finds the peak in the
histogram, and the sensor reports a single distance value per detector pixel.
When using a diffuse light source, these time-resolved sensors capture visual
information beyond the distance measurements extracted in LiDAR. Transient
histograms in this case record distributions of times-of-flight, encoding the
product of scene geometry and reflectance over each imaged scene patch
Jungerman_ECCV_22.

The entirety of information in the transient histogram has been previously
utilized in applications such as fluorescence lifetime imaging lagarto2020real
and non-line-of-sight imaging velten2012recovering. Recently, time-resolved
measurements from such sensors have been explored for line-of-sight 3D
reconstruction and novel view synthesis Malik2023transient. Unfortunately,
these systems, due to their prototype-like nature, exist in a bulky benchtop
form factor, require significant power, and are costly (over $10,000 USD).
They are thus not suited for many applications including autonomous drones,
wearable computing, and augmented reality, in which low-cost, energy
efficient, miniature sensors are required.

Low-cost single-photon cameras have recently become available in the form of
active single photon avalanche diodes (SPADs). They include one or more SPADs
paired with an eye-safe diffuse light source (_e.g_., an infrared VCSEL
laser), are very small ($<$20 mm3), inexpensive ($<$$5 USD), and power
efficient ($<$10 milliwatts per measurement) TMF8820, VL6180X. These are sold
as proximity sensors and some can be configured to report transient
histograms. Compared to laboratory-grade systems, however, they lack precise
optics, calibration, and timing characteristics and have an order of magnitude
lower spatial and temporal resolution. Still, these sensors have proven
successful for material classification Becker2023plastic, human pose
recognition Ruget2022pixels2pose, and simple shape recovery (_i.e_., a planar
surface) Jungerman_ECCV_22,Sifferman2023unlocking.

In this work, we address the problem of reconstructing 3D shape of arbitrary
Lambertian objects from a set of spatially distributed low-cost single-photon
cameras with known poses. We present an approach that combines a neural signed
distance field surface representation, a differentiable transient formation
model for practical active single-photon cameras, and an optimization scheme
following the analysis-by-synthesis pipeline. Fig. 1 illustrates our sensor,
imaging setting, and approach. We show that our approach can successfully
recover complex 3D shapes from simulated data. We further demonstrate 3D
object reconstruction from real-world captures, utilizing measurements from a
low-cost, off-the-shelf proximity sensor.

Our approach draws a connection between image-based modeling and active range
scanning while avoiding some of their pitfalls. By using an active light
source and SPADs, our approach illuminates and images an object from multiple
angles, resembling the imaging setup of multiview photometric stereo. Our
model also echoes the design of NeRF. The key difference is that our model
considers an input of low spatial resolution (_i.e_., hundreds of pixels per
scene), and leverages rich temporal information encoded in the transient
histograms. Unlike image-based modeling, our approach can operate under low
light conditions, and reconstruct objects without texture. Unlike most range
scanning methods, our approach performs favorably in the presence of strong
ambient flux and is robust to mildly specular surfaces (see top of the spray
bottle in Fig. 1).

Practical Applications. Our approach presents a compact, energy-efficient, and
low-cost solution for 3D sensing when relative sensor poses are known. We
envision applying our method to settings where relative sensor poses are
fixed, and 3D geometry can be estimated on a per-frame basis by combining
measurements from a distributed set of sensors. This includes applications
like wearable hand tracking devrio2022discoband,hu2020fingertrak, and
collision avoidance and mapping for drones, mobile robots, and robot
manipulators. This work provides a first step towards enabling such
applications.

Scope and limitation. Similarly to other methods for active range scanning,
the proposed method fails on highly specular objects and requires hundreds of
views to reconstruct a complex 3D shape. Our imaging system is tailored for
low-cost sensors with limited range (4-6 meters) and temporal resolution (over
100 picoseconds), and thus it is unproven at imaging larger space or finer
details. While methods exist for fast NeRF muller2022instant, our approach is
not optimized for speed and takes a few hours to reconstruct a single object.
Instead, we focus on demonstrating the feasibility of reconstructing real-
world 3D objects using commodity sensors.

##  2 Related Work

Time-resolved imaging. Time-resolved sensors have long been used in
applications such as non-line-of-sight (NLOS) imaging
xin2019theory,Faccio2020non, where a scene is recovered from around the
corner, fluorescence lifetime imaging lee2023caspi, a microscopy technique for
characterizing a biological sample, and to measure distance using direct time-
of-flight heide2018sub,gupta2019photon. Recent techniques Malik2023transient
used high-end time-resolved sensors to directly recover scene geometry.
However, these methods rely on hardware prototypes which are often
prohibitively expensive and not accessible to consumers. For example,
Transient NeRF Malik2023transient performs view synthesis and 3D
reconstruction using a lab prototype with no ambient light and $2$-$5$ views
of $512\times 512$ transient histograms, each with 1500 bins. In contrast, in
this work, we use commodity hardware with ambient light and $128$-$240$ total
transients, each comprised of $128$ bins.

3D Imaging with low-cost SPADs. With cheap time-resolved proximity sensors
becoming commonplace, recent works have investigated their use for 3D
reconstruction. Callenberg _et al_. Callenberg2021CheapSPAD demonstrate that,
with some additional hardware, high-resolution depth imaging from a single
viewpoint is possible. A low-cost SPAD has been used to augment an RGB SLAM
system tofslam. Other works utilize supervised machine learning to recover
geometric information from single low-cost SPAD measurements, such as 3D human
pose Ruget2022pixels2pose or high-resolution depth images
Jungerman_ECCV_22,deltar. Jungerman _et al_. Jungerman_ECCV_22 use
differentiable rendering to recover two degrees of freedom of a planar surface
from a single low-cost SPAD transient histogram. Sifferman _et al_.
Sifferman2023unlocking extend this method to fully recover a planar surface
using a low-cost SPAD with multiple detector pixels.

Neural implicit representations. Neural representations, as popularized by
NeRF Mildenhall2020nerf, enable novel view synthesis and 3D reconstruction by
representing the scene as a neural network. While the original NeRF
representation encoded view-dependent volumetric effects, alternative
encodings have been proposed to better model geometry and reconstruct
surfaces. NeuS Wang2021neus represents the scene as a level set, allowing for
better modeling of surfaces at the expense of not being able to represent
volumetric effects. Many works extend these ideas to work with different
sensing modalities and external supervision, such as depth queues from
structure-from-motion deng2022dsnerf, RGB images plus continuous-wave time-of-
flight sensors attal2021torf, only depth information
OrtiziSDF2022,liu2023multimodal, or more recently using only transients
Huang2023nfl,Malik2023transient. In this work, we perform 3D reconstruction
using only transient histograms as captured using commodity hardware by
adapting the implicit surface representation introduced by NeuS to render
transients.

Multiview photometric stereo. Our approach illuminates and images an object
from multiple views in order to reconstruct its 3D shape, having conceptual
similarity to image-based modeling methods. This concept has previously been
explored as multiview photometric stereo (MVPS)
hernandez2008mvps,cheng2021multi: reconstructing 3D geometry given distributed
views of the scene and various lighting conditions. More recent approaches
resort to deep learning models zhao2023mvpsnet and consider NeRF-based
representations kaya2023multi,yang2022psnerf.

##  3 3D Vision with Single-Photon Cameras

Figure 2: Method Overview: The scene is modeled as a neural implicit surface
in the form of an SDF. To render a transient, we approximate Eq. 8 by sampling
rays within each pixel’s FoV, and subsequently points on those rays. This
idealized transient waveform is then convolved with the sensor’s laser impulse
response to model the transient histogram formation. Finally, we optimize the
scene representation by minimizing a loss between the rendered transients and
the observations.

We propose to reconstruct 3D scene geometry using a sparse set of measurements
from low-cost time-resolved SPAD sensors. Our approach assumes a distributed
set of $N$ SPAD sensors with known pose, each comprising a single-pixel
detector co-located with a diffuse laser. The time of flight, $t$, of
returning photons, received over the detector FoV 111We assume that the
illumination and detector FoV coincide and refer to both as the sensor FoV in
the rest of the paper. is recorded and binned into _transient histograms_
Jungerman_ECCV_22,Sifferman2023unlocking, $\\{\mathbf{h}_{j}\\}_{j=1}^{N}$:
time-resolved measurements that encode rich information about the scene
geometry and reflectance. Our goal is to recover scene geometry from these
measurements.

Traditional depth ranging. Our setup differs from the conventional use of
SPADs for depth ranging heide2018sub,gupta2019photon,Malik2023transient, which
assumes a collimated laser and focused detector. There, the depths of
individual scene points can be independently estimated from separate
transients, and a dense scan of thousands of scene points is needed to obtain
full scene geometry. These systems require high-quality beam lasers and
detectors with high temporal resolution. In contrast, our imaging
configuration is indicative of the real-world use cases of low-cost SPAD
sensors with diffuse lasers, wide FoVs, and low-resolution detectors. In our
case, the transient encodes rich information beyond its peak location.

Low-cost SPAD sensing. Our setup extends recent systems
Jungerman_ECCV_22,Sifferman2023unlocking built around low-cost SPAD sensors
for 3D sensing in two ways. First, our system consists of multiple posed
sensors as opposed to a single sensor. Second, it is capable of capturing and
reconstructing complex, non-parametric scenes as opposed to the parametric
geometry of a single plane. With these upgrades, our system is reminiscent of
multi-view stereo and NeRF-like systems for conventional RGB cameras, and is a
substantial step towards practical 3D vision using commodity proximity
sensors.

###  3.1 Transient Formation Model

We model the transient formation process in two steps. We first derive the
_transient waveform_ of a scene by modeling the interaction of light with
scene geometry. The next step uses this waveform as the input to a sensor
model that accounts for non-linearities that occur when capturing _transient
histograms_ with real hardware.

Transient waveform. As in Jungerman_ECCV_22,Sifferman2023unlocking, we ignore
high-order light paths and only model direct reflection from the object
surface. The transient waveform $\tau_{\mathbf{s}}(t)$ of a sensor
$\mathbf{s}$ located at the origin can be written in path integral form
veach1998robust as

$\begin{split}\tau_{\mathbf{s}}(t)&=\int_{\mathcal{S}}\rho(\mathbf{x})\frac{G(\mathbf{x},\boldsymbol{\omega})}{\|\mathbf{x}\|^{4}}\delta\left(\|\mathbf{x}\|_{2}-\frac{ct}{2}\right)\mathrm{d}A(\mathbf{x}),\\\
G(\mathbf{x},\boldsymbol{\omega})&=f_{r}(\mathbf{x},\mathbf{n}_{\mathbf{x}},\boldsymbol{\omega})V(\mathbf{x})\langle\boldsymbol{\omega},\mathbf{n}_{\mathbf{x}}\rangle^{2},\end{split}$
(1)

where $\mathrm{d}A(\mathbf{x})$ is an infinitesimal area around point
$\mathbf{x}$ on the object surface $\mathcal{S}$, $\rho$ is the albedo,
$\delta$ is the Dirac delta function, $c$ the speed of light, and
$\boldsymbol{\omega}$ is a ray direction from $\mathbf{x}$ to $\mathbf{s}$.
Time-independent effects such as $V$, the visibility function governed by
scene geometry and sensor FoV, $f_{r}$ the bidirectional reflectance
distribution function (BRDF), and foreshortening effects where
$\mathbf{n}_{\mathbf{x}}$ is the surface normal at $\mathbf{x}$, are modeled
as $G(\mathbf{x},\boldsymbol{\omega})$.

Sensor model. The sensor model accounts for laser and detector characteristics
when converting a transient waveform into a transient histogram
hernandez2017computational. Our sensor model considers the laser pulse, laser
power, detector quantum efficiency, ambient photon flux, internal detector
noise, pile-up effect, and time jitter. In practice, the _laser pulse_ is not
a perfect impulse and, despite bandpass filters, the measured transient also
captures some constant ambient light. To model this, we convolve $\tau(t)$
222We drop the subscript $\mathbf{s}$ hereafter for clarity. with the laser’s
impulse response $g(t)$, scaled by $\phi^{scale}$ which absorbs laser power
and quantum efficiency of the detector, and then offset its intensity by
$\phi^{bkgd}$ which encapsulates ambient photon flux and internal detector
noise:

$\tilde{\tau}(t)=\phi^{scale}(\tau*g)(t)+\phi^{bkgd}.$ (2)

$\tilde{\tau}(t)$ is subsequently discretized into a histogram of Poisson
rates $\mathbf{r}=[r_{1},...,r_{B}]$ with $B$ bins. The probability $q_{i}$ of
at least one photon falling inside the $i^{\mathrm{th}}$ bin is given by
coates1968correction

$q_{i}=1-\exp(-r_{i}).$ (3)

In practice, SPADs aggregate photon counts over $C$ laser cycles, with only
the first incident photon being detected in each cycle. This results in _pile-
up_ , a nonlinear distortion of transients, leaving photons arriving at a
later timestamp less likely to be detected pediredla2018pileup. Specifically,
the probability $p_{i}$ of detecting a photon in the $i^{\mathrm{th}}$ bin in
a cycle is given by pediredla2018pileup

$p_{i}=q_{i}\Pi_{k=1}^{i-1}(1-q_{k}).$ (4)

The photon counts $[h_{1},...,h_{B}]$ in a transient histogram
$\tilde{\mathbf{h}}$ follow a multinomial distribution:

$[h_{1},...,h_{B+1}]\sim\mathrm{Multinomial}(C,(p_{1},...,p_{B+1})),$ (5)

where $p_{B+1}=1-\sum_{i=1}^{B}p_{i}$, and $h_{B+1}$ counts the number of
cycles without detected photons. $\tilde{\mathbf{h}}$ is subsequently
convolved with a discretized _time jitter_ kernel $\mathbf{s}$ to account for
the temporal uncertainty of photon detection events, yielding the final
histogram $\mathbf{h}$ measured by a SPAD detector:

$\mathbf{h}[b]=\Sigma_{k}\mathbf{h}[k]\mathbf{s}[b-k].$ (6)

We now present a reconstruction algorithm based on the analysis-by-synthesis
principle to recover scene geometry from the transients. The input to this
system is a distributed set of posed transients
$\\{\mathbf{h}_{j}\\}_{j=1}^{N}$ captured of the scene.

###  3.2 Neural Scene Reconstruction from SPADs

The reconstruction problem that we attempt to solve is extremely challenging.
Unlike in conventional depth ranging where the depth of a scene point can be
directly determined from the histogram via peak finding, a transient from our
system represents the superposition of light reflected from numerous scene
points, as illuminated by a diffuse laser, and is further contaminated by non-
idealities of the detector (_e.g_., pile-up). The direct inversion of the
signal is thus a highly ill-posed problem. Further, we cannot adapt methods
from the NLOS imaging literature o2018lct,lindell2019fk,xin2019theory as they
only support dense 2D scans, whereas our system uses a distributed, sparse and
unstructured set of measurements.

To overcome these challenges, we resort to an analysis-by-synthesis approach
based on differentiable rendering. Our approach allows flexible positioning of
sensors and accurate modeling of histogram formation, thereby enabling high-
quality reconstruction of scene geometry. We now describe our reconstruction
algorithm in detail.

Neural scene representation. Following NeuS Wang2021neus, we represent the
scene geometry as a signed distance function (SDF), parameterized as a multi-
layer perceptron (MLP) $f_{\theta}:\mathbb{R}^{3}\rightarrow\mathbb{R}$.
$f_{\theta}$ maps the position-encoded (PE) $xyz$-coordinates of a point
$\mathbf{x}$ to its signed distance $d$:

$d=f_{\theta}(\mathrm{PE}(\mathbf{x})).$ (7)

Compared to Jungerman_ECCV_22,Sifferman2023unlocking, this neural SDF allows
our method to represent scene geometry beyond simple parametric shapes as the
level set
$\mathcal{S}=\\{\mathbf{x}\in\mathbb{R}^{3}|f_{\theta}(\mathbf{x})=0\\}$.

Transient volume rendering. The key idea behind our analysis-by-synthesis
approach is to render $f_{\theta}$ into transients and compare them with those
captured by our system. To adapt Equation 1 for the rendering of $f_{\theta}$,
we first rewrite it in angular integral form as

$\begin{split}\tau(t)=\int_{\Omega}&\frac{\rho}{\pi}\frac{V(\mathbf{x})\langle-\boldsymbol{\omega},\mathbf{n}_{\mathbf{x}}\rangle}{\|\mathbf{x}\|^{2}}\delta\left(\|\mathbf{x}\|_{2}-\frac{ct}{2}\right)\mathrm{d}\boldsymbol{\omega},\end{split}$
(8)

where $\boldsymbol{\omega}$ are ray directions in the sensor FoV $\Omega$, and
$\mathbf{x}$ the point where $\boldsymbol{\omega}$ intersects with the object
surface $\mathcal{S}$ ($\infty$ if no intersection). For simplicity, we assume
a learned spatially uniform albedo $\rho$ and Lambertian BRDF $f_{r}=1/\pi$
333This assumption may be relaxed to allow more expressive BRDF models such as
the Phong reflection model bui1975phong..

Inspired by NeRF Mildenhall2020nerf and NeuS Wang2021neus, we approximate
Equation 8 via volume rendering to resolve surface discontinuities, enabling
the optimization of $\theta$ via gradient descent:

$\hat{\tau}(t)=\int_{\Omega}\frac{\rho}{\pi}\frac{T^{2}(t)\sigma(\mathbf{p}(\boldsymbol{\omega},t))\langle-\boldsymbol{\omega},\mathbf{n}_{\mathbf{p}}\rangle}{\|\mathbf{p}(\boldsymbol{\omega},t)\|^{2}}\mathrm{d}\boldsymbol{\omega}.$
(9)

Here,
$\mathbf{p}(\boldsymbol{\omega},t)=\nicefrac{{ct}}{{2}}\,\boldsymbol{\omega}$
are points along $\boldsymbol{\omega}$, the volume density $\sigma$ is a
function of $f_{\theta}$ as in NeuS Wang2021neus, and the transmittance $T$ is
given by

$T(t)=\exp{\left(-\int_{0}^{t}\sigma(\mathbf{p}(u))\mathrm{d}u\right)}.$ (10)

In practice, we discretize $\hat{\tau}(t)$ over the transient bin intervals
$\\{[t_{i},t_{i+1})\\}_{i=1}^{B}$ and work with the histogram
$\hat{\boldsymbol{\tau}}=[\hat{\tau}_{1},...\hat{\tau}_{B}]$, where

$\hat{\tau}_{i}=\int_{\Omega}\frac{\rho}{\pi}\int_{t_{i}}^{t_{i+1}}\frac{T^{2}(t)\sigma(\mathbf{p}(\boldsymbol{\omega},t))\langle-\boldsymbol{\omega},\mathbf{n}_{\mathbf{p}}\rangle}{\|\mathbf{p}(\boldsymbol{\omega},t)\|^{2}}\mathrm{d}t\,\mathrm{d}\boldsymbol{\omega}.$
(11)

We estimate the intractable Equation 11 via Monte Carlo sampling of
$\boldsymbol{\omega}$ and subsequently of $\mathbf{p}(\boldsymbol{\omega},t)$.

Bilevel importance sampling. Similar to Mildenhall2020nerf,Wang2021neus, the
sampling of $\mathbf{p}(\boldsymbol{\omega},t)$ is weighted by a probability
density function (PDF) over the equally sized bin intervals. This PDF is
proportional to the per-bin weights $w_{i}$ given by

$w_{i}=\exp{\left(-\Sigma_{j=1}^{i-1}\sigma_{j}\Delta\right)}(1-\exp{(-\sigma_{i}\Delta})),$
(12)

where $\sigma_{i}$ is evaluated at the mid-point of the $i^{\mathrm{th}}$ bin,
and $\Delta$ is the bin size in distance.

We extend this idea to the importance sampling of $\boldsymbol{\omega}$.
Specifically, the sampling PDF over a uniform partitioning of FoV $\Omega$ is
proportional to the cumulative weights $w^{(k)}$ over rays
$\boldsymbol{\omega}^{(k)}$ drawn from each partition $k$:

$w^{(k)}=\Sigma_{i=1}^{B}w^{(k)}_{i}.$ (13)

Intuitively, this allows us to point more rays at high-density regions
occupied by the object surface.

Differentiable sensor modeling. Modeling sensor behavior is particularly
important for our analysis-by-synthesis approach. This is because the
synthesis targets $\boldsymbol{\tau}$ are not determined by the scene geometry
alone but reflect the complex interplay of geometry with sensor non-idealities
including pulse shape, pile-up and time jitter. To this end, we cascade
$\hat{\boldsymbol{\tau}}$ to a _differentiable sensor model_ $\Gamma$ to
simulate the transformation applied by the sensor to raw waveforms.

Specifically, $\Gamma$ closely follows the sensor model in Section 3.1;
Equations 2-4 are differentiable and applied sequentially on
$\hat{\boldsymbol{\tau}}$, yielding per-bin photon detection probabilities
$\hat{\mathbf{p}}=[\hat{p}_{1},...,\hat{p}_{B}]$. Instead of sampling photon
counts using Equation 5, we directly convolve $\hat{\mathbf{p}}$ with the
jitter kernel as in Equation 6. This allows us to sidestep the non-
differentiable sampling step while producing an unbiased estimate of the
transient $\hat{\mathbf{h}}=\Gamma(\hat{\boldsymbol{\tau}})$ for loss
evaluation.

Loss functions. The optimization of $\theta$ is driven by three loss terms: a
histogram reconstruction loss $\mathcal{L}_{hist}$ that minimizes the L1
distance between $\hat{\mathbf{h}}$ and $\mathbf{h}$, an Eikonal loss
gropp2020eikonal $\mathcal{L}_{Eikonal}$ that regularizes the SDF, and a total
variation regularizer mu2022physics $\mathcal{L}_{TV}$ that penalizes floaters
in empty space. The combined loss function $\mathcal{L}$ is thus given by

$\mathcal{L}=\mathcal{L}_{hist}+\lambda_{Eikonal}\mathcal{L}_{Eikonal}+\lambda_{TV}\mathcal{L}_{TV},$
(14)

where $\lambda_{Eikonal}$ and $\lambda_{TV}$ are the respective loss weights.

Comparison to NLOS-NeuS fujimura2023nlos-neus. NLOS-NeuS is a concurrent work
for scene reconstruction using SPAD sensors. Despite similarities in the
imaging model and reconstruction algorithm, the two works differ substantially
in application (direct line-of-sight _vs_. non-line-of-sight), laser
characteristics (eye-safe diffuse laser _vs_. high-energy beam laser), sensor
quality and cost (commodity, low-cost _vs_. laboratory-grade, high-cost),
transient modeling (with _vs_. without sensor non-idealities) and scan pattern
(distributed and sparse _vs_. dense and structured 2D grid). We believe the
two works complement each other and together unveil an exciting avenue toward
single-photon 3D vision systems.

##  4 Experiments

We demonstrate the effectiveness of our method for 3D geometric reconstruction
of various objects in simulation, and in the real world with a low-cost SPAD.
We provide qualitative and quantitative results for both settings. See the
supplement for implementation details (_e.g_. learning rates).

|  Chamfer Distance (mm) $\downarrow$
---|---
Method | Armadillo | Bear | Bunny | Digit | Einstein | Skull | Soap | Sphere
Reprojection (Peak) | 54.29 | 40.36 | 34.95 | 55.85 | 43.25 | 48.90 | 51.71 | 51.07
Reprojection (Threshold) | 65.43 | 60.72 | 54.05 | 60.64 | 61.31 | 65.14 | 68.74 | 63.16
Space Carving | 34.78 | 24.53 | 22.29 | 45.44 | 26.60 | 25.49 | 21.44 | 25.47
Ours | 3.93 | 5.95 | 3.84 | 3.27 | 3.51 | 3.22 | 3.23 | 3.77
Table 1: Quantitative results on simulated data. Our method more accurately
recovers 3D shapes than baselines across $8$ objects. Figure 3: Qualitative
results on simulated data. Our method reconstructs dense and detailed 3D
shapes. Space carving provides only hulls of a target shape, and is prone to
carving away extra space when thin structures are present. Reprojection yields
sparse points.

Baselines. We compare our method to two baselines: reprojection
heide2018sub,gupta2019photon and space carving
kutulakos2000space,tsai2017geometry.

_Reprojection_ , also known as back-projection, reconstructs a scene as a
point cloud and is the _de facto_ standard for depth ranging. We compare to
two forms of reprojection. The peak method finds the distance $d$
corresponding to the histogram bin with the highest intensity. For a sensor at
position $\mathbf{s}$ with an outwards pointing optical axis $\mathbf{u}$, a
point is placed in the scene at position $\mathbf{s}+d\mathbf{u}$. The
threshold method works in the same way but finds $d$ by locating the lowest-
index bin with intensity above a threshold $t_{p}$. If no bin passes the
threshold, no point is projected. We do not apply surface reconstruction to
the generated point clouds, as it is often unreliable for poorly imaged
scenes. See the supplement for further discussion.

_Space carving_ reconstructs a scene as a voxel grid. Like thresholded
reprojection, it finds the distance $d$ corresponding to the lowest-index bin
with intensity above a threshold $t_{s}$. All voxels in the sensor’s FoV and
nearer than $d$ are marked empty, along with voxels outside the FoV. Voxels in
the FoV and further than $d$ are marked as occupied. The carved scene is the
union of the occupied set for all sensors.

To ensure strong baselines for real-world experiments, we perform a brute-
force search over $t_{p}$ and $t_{s}$ and choose values that minimize Chamfer
distance over the entire real-world dataset. Space carving voxel size was set
to $1.0\mathrm{cm}$.

Evaluation protocol. Following NeuS Wang2021neus, we evaluate all methods
using Chamfer distance. We report standard (two-way) Chamfer on simulated
data. For real-world captures, we report Chamfer in both directions to
evaluate the quality of reconstruction. Prior to Chamfer calculation, we
convert ground-truth meshes and reconstructions from our method to point
clouds by drawing $5$ million points uniformly at random on the mesh surface.
For space carving, occupied voxels are converted to points if they touch
unoccupied space, excluding the edge of the grid.

|  Chamfer Distance (mm) $\downarrow$ Rec $\rightarrow$ GT / GT $\rightarrow$
Rec
---|---
Method |  Big Box∗ |  Block∗ |  Pyramid∗ |  Toy Container† |  Cereal Box†
Reprojection (Peak) | 77.4 | 24.9/ | 52.5 | 51.8 | 12.8/ | 39.0 | 94.7 | 17.5/ | 77.1 |  71.0 | 24.9/ | 46.0 |  49.3 | 17.3/ | 31.9
Reprojection (Threshold) | 67.5 | 14.8/ | 52.7 | 52.3 | 8.5/ | 43.8 | 75.4 | 5.9/ | 69.5 | 52.4 | 8.9/ | 43.5 | 51.6 | 19.2/ | 32.4
Space Carving | 67.9 | 35.1/ | 32.8 | 69.2 | 33.4/ | 35.8 | 80.1 | 39.5/ | 40.6 | 98.9 | 52.8/ | 46.0 | 44.1 | 24.4/ | 19.6
Ours | 12.5 | 6.1/ | 6.4 | 9.8 | 5.6/ | 4.2 | 18.4 | 9.0/ | 9.3 | 11.5 | 5.8/ | 5.6 | 16.3 | 8.3/ | 8.0
Table 2: Quantitative results on real-world captures. Our method more
accurately reconstructs real-world objects with homogeneous (*) and rich (†)
texture. Reprojection yields sparse and unevenly distributed points, harming
one-way Chamfer from GT to reconstruction.

###  4.1 Simulated Experiments

Experiment setup. We simulate transients for eight scenes of varying
complexity using the image formation model in Section 3.1. The objects are
centered on the ground plane ($z=0$) with the largest dimension $\approx
0.3\mathrm{m}$. Sensors with a conical FoV are uniformly distributed on a
hemisphere at the origin with a radius of $0.5\mathrm{m}$, and are all pointed
at the origin. In our simulation, $N=256$, $B=256$, $\Delta=5\mathrm{mm}$,
$\mathrm{FoV}=30^{\circ}$, $\phi^{scale}=1$, $\phi^{bkgd}=0.001$, $C=5000$ and
$\rho=0.8$. The laser pulse, $g$, has a full-width-at-half-maximum (FWHM) of
$50\mathrm{ps}$, and $\mathbf{s}$ is a tabulated PDF obtained from experiments
hernandez2017computational. The sensor parameters are deliberately chosen to
reflect the characteristics of low-cost sensors. See the supplement for a
sensitivity study on how these parameters affect the quality of
reconstruction.

Results. Table 1 summarizes the quantitative results of all methods. Our
method achieves an average Chamfer distance of $<5\mathrm{mm}$, an order of
magnitude lower than all baselines. A key reason is that the baselines only
use depth information from a single histogram bin, whereas our method makes
effective use of the entire waveform, which contains rich geometry cues about
a large scene patch.

We provide visualizations of our results in Figure 3. Our method recovers
global scene structure as well as local geometry details. In contrast,
reprojection yields sparse point clouds without sufficient coverage of the
scene. While space carving produces dense reconstructions, the occupancy grid
only represents an envelope of the scene, leaving it difficult to recognize
the precise shape of an object.

Figure 4: To capture real-world data from a wide set of viewpoints, we mount
the TMF8820 proximity sensor to a robot arm. Forward kinematics of the robot
are used to gather sensor pose. Figure 5: Qualitative results on real-world
captures. Our method again attains the highest reconstruction quality. Poses
in column two are subsampled by a factor of two for clarity. See supplement
for additional qualitative results.

###  4.2 Hardware Prototype

We use the SPAD-based AMS TMF8820 proximity sensor TMF8820, which retails for
$10 USD. We connect the sensor to a microcontroller via I2C and use the AMS-
provided driver to extract transient histograms.

The sensor contains a total of 216 SPADs, which are pooled onboard the sensor
into $3\times 3$ zones, each of which images a different FoV. The sensor
captures one transient histogram for each zone. We pool histograms from all
zones, which is equivalent to capturing one wide-FoV histogram per-measurement
as SPADs do not suffer from readout noise zappa2007principles. In doing so, we
avoid inter-histogram interference previously observed by
Sifferman2023unlocking and avoid the need to model individual fields-of-view
of the sensor, which we empirically observed to have soft and poorly specified
boundaries. We slightly modify our method to accommodate the AMS TMF8820
sensor used in real-world experiments.

Laser Impulse. The laser impulse response of the TMF8820 is not Gaussian and
varies slightly between measurements. Fortunately, the sensor captures the
shape of its laser impulse for each measurement in a “reference histogram”. We
record this histogram for each measurement and incorporate it into our forward
model by cross-correlating the idealized scene response with this recorded
reference histogram. We observe that the bin size $\Delta_{r}$ of the
reference histogram is smaller than the bin size $\Delta$ of the transient
histograms captured by the sensor. To account for this, we scale the reference
histogram in the temporal dimension by a factor $\Delta_{r}/\Delta$ before
cross-correlation. Further, we find that it is necessary to temporally shift
the reference histogram by a fixed amount $\phi^{delay}$ before correlation.

To calibrate the parameters $\Delta$, $\Delta_{r}$, and $\phi^{delay}$, we
perform the one-off intrinsic calibration procedure separately introduced by
Sifferman2023unlocking. The TMF8820 sensor is pointed at a planar surface from
a range of known distances and angles-of-incidence. A differentiable render-
and-compare method is used to optimize for the unknown sensor intrinsic
parameters given known planar geometry.

Pile-up Correction. While our forward model assumes that the target transients
exhibit nonlinear distortion due to pile-up, the TMF8820 sensor performs pile-
up correction on-sensor, and it cannot be disabled. To accommodate this, we
incorporate the differentiable Coates’ correction coates1968correction as a
final step in the forward model.

Other Sensors. Our design can be easily applied to any SPAD with a co-located
diffuse illumination source, including other low-cost sensors VL53L8CH, high-
end setups Jungerman_ECCV_22, Nishimura2020disambiguating, and things in
between. The only modification necessary is the calibration of sensor
intrinsics (_e.g_., using Sifferman2023unlocking).

###  4.3 Real-world Experiments

Experiment Setup. We capture a real-world tabletop dataset of eight
objects444See supplement for full results over the entire dataset. of varying
geometry and texture. To capture many posed views of the target object, we
attach the sensor to a Universal Robots UR5 robot arm. We program the arm to
automatically move to a set of poses and record sensor measurements at each
pose. To obtain sensor poses, we use the forward kinematics of the robot,
which are accurate to $\pm 0.5$mm pollak2020measurement. Each object is
captured from between 128 and 240 viewpoints. Five of the objects are simple
geometric primitives, for which we manually generate ground-truth meshes based
on the dimensions of the target object and measurements of its position from
the robot’s forward kinematics. Meshes are trimmed to an axis-aligned bounding
box 16cm larger than the target object in each dimension before the Chamfer
distance calculation.

Results. As seen in Table 2, our method outperforms all baselines by a wide
margin as measured by two-way Chamfer distance. While reprojection is at times
competitive in one-way distance from reconstruction to ground truth, it
performs poorly in the opposite direction due to the sparse and unevenly
distributed point cloud generated, as visualized in Figure 5. While space
carving outperforms reprojection on simulated data under a highly structured
sensor pose distribution (_i.e_. all sensors are facing the center of the
object), it yields poor results on real-world scenes, in which we vary sensor
orientation by $\pm 10^{\circ}$ to emulate real-world capture conditions and
increase coverage. By contrast, our method benefits from the more varied
sensor poses as is shown in Figure 6.

Figure 6: Space carving performs more poorly when sensor poses include some
variation in target point, while our system takes advantage of the increased
view diversity and coverage. Figure 7: Because our reconstruction method
assumes a Lambertian surface, it fails to reconstruct highly specular scenes,
such as a glossy white bust (left) or mirror-finish kettle (right).

Further, our method is surprisingly robust to violation of assumptions made
about surface reflectance; it successfully reconstructs non-Lambertian objects
with rich texture despite assuming Lambertian BRDF with a spatially uniform
albedo. These include both simple shapes (Toy Container in Figure 5 and Cereal
Box) and challenging objects with complex geometry (Spray Bottle in Figure 5).
We hypothesize that the overlapping FoVs of distributed sensors help constrain
the optimization of our model and encourage a plausible reconstruction that
best explains all transients. Our strong results on simulated and real-world
data validate our modeling approach and demonstrate a single-photon 3D vision
system for real-world scene reconstruction.

##  5 Discussion and Future Work

Despite assuming spatial uniform reflectance and albedo, our method is robust
to rich textures (Fig. 5) and compares favorably to baseline methods for
reconstructing challenging scenes with high specularities (Fig. 7). Future
work will investigate recovering spatially varying reflectance or incremental
learning of geometry as new measurements become available OrtiziSDF2022,
Sucar2021imap, enabling applications like real-time mapping and SLAM. Our
method may be particularly relevant in applications such as robotics and
wearable computing, where the small size, low power requirements, and robust
hardware of proximity sensors are very valuable.

Acknowledgements: This work was supported by Los Alamos National Lab and the
Department of Energy, Wisconsin Alumni Research Foundation, NSF CAREER Award
1943149, NSF grant CNS-2107060, and a SONY Faculty Innovation Award. Gleicher
and Gupta hold concurrent appointments at Amazon and Cruise, respectively.
This work is not associated with either company.

ieeenat_fullname main

Supplementary Material

In this supplement, we describe (1) The implementation details of our
reconstruction pipeline (Section A); (2) a sensitivity analysis to understand
the robustness of our method w.r.t. imaging parameters (Section B); (3) A
discussion of the difficulty of point cloud surface reconstruction (Section
C); (4) additional qualitative results on simulated data and real world
captures (Section D); and (5) further discussion of our work (Section E). We
also include three short videos, which provide animated $360^{\circ}$ views of
qualitative reconstruction results. The videos can also be viewed on the
project web page: https://cpsiff.github.io/towards_3d_vision/.

For sections, figures and equations, we use numbers (_e.g_., Sec. 1) to refer
to the main paper and capital letters (_e.g_., Sec. A) to refer to this
supplement.

##  A Implementation Details

For both simulated and real world experiments, we use an $8$-layer MLP with
$256$ hidden units as our SDF, $f_{\theta}$, and initialize it as a sphere,
centered at the origin with radius $0.3$m, using geometric initialization
atzmon2020sal_. For each transient, we sample $256$ rays $\boldsymbol{\omega}$
over $\Omega$ and sample $256$ points per ray. We set $\lambda_{Eikonal}$ to
$0.1$ across all experiments and set $\lambda_{TV}$ to $0$ and $0.01$
respectively for the simulated and real-world experiments. We train
$f_{\theta}$ for $300\mathrm{K}$ steps using Adam kingma2014adam_ with a mini-
batch size $2$, a learning rate $0.0005$, and cosine decay. The learned SDFs
are converted to meshes using Marching Cubes lorensen1998marching_.

##  B Sensitivity Analysis

Figure A: Sensitivity analysis of our method compared to baselines across a
range of imaging parameters. In almost every case, our method outperforms
baseline methods on Chamfer distance. Missing datapoint in (d) indicates that
our method failed to converge. Illumination power (e) is unit-less as it also
absorbs factors like quantum efficiency and does not map directly to any real
world parameter.

We perform extensive experiments to understand the robustness of our method in
comparison to baselines under varying sensor parameters in simulation. All
experiments are based on the Bunny scene and the parameters are varied one at
a time while other parameters remain fixed at the base condition (as described
in Section 4 of the main paper). To ensure strong baselines for every sensor
configuration, we calibrate the thresholds $t_{p}$ and $t_{s}$ for the
projection (threshold) and space carving baselines respectively per
reconstruction. We perform a brute force search over possible thresholds and
report the best Chamfer achieved. As this amounts to calibrating on the test
set, the numbers reported represent the best possible performance of the
baseline methods on the given data. The results of this sensitivity analysis
are presented in Figure A. In what follows we discuss some of the main
findings.

Sensor Placement. We study two key parameters that control sensor placement:
the number of views and the minimum elevation angle at which the sensors are
placed. Our method consistently outperforms all baselines in Chamfer distance
by an order of magnitude across a broad range of parameter choices. In
particular, our method readily supports as few as $128$ views above a
considerably large elevation angle of $30^{\circ}$ without harming
reconstruction quality. This robust gain in performance confirms that our
method takes advantage of broad-band signal in transients not exploited by the
baseline methods.

Temporal Resolution. Our system takes advantage of the temporal information in
transient histograms, and therefore benefits when that information is present
at a high resolution. Because of this, our method outperforms baselines by a
very wide margin at a small bin size, but the margin vanishes as bins become
wider than $2\mathrm{cm}$ (equivalently $66\mathrm{ps}$), because decomposing
the temporal signal becomes impractical beyond this limit. Fortunately,
today’s commodity SPADs operate at a smaller bin size ($\sim 40\mathrm{ps}$).
Baseline methods show no performance gain at small bin sizes, as they do not
take advantage of the temporal resolution.

Angular Resolution. Our system resolves spatial resolution from wide-FoV
sensors by taking advantage of the time dimension. In this regime, the optimal
sensor field-of-view size is not obvious: a smaller FoV means more highly
constrained geometry, as each histogram images a smaller region, but too small
of a field-of-view means a lack of coverage and under-constrained geometry. We
find that an angular resolution in the $30^{\circ}$ to $60^{\circ}$ range is
optimal for reconstructing 3D geometry with our method on the bunny scene.
Reprojection based methods benefit more from a smaller field-of-view, while
space carving performs best with a wider field-of-view so that space is
sufficiently carved away. In every case, our method outperforms baselines by a
wide margin.

Signal-to-noise ratio (SNR). We consider three parameters that jointly impact
SNR: illumination power, ambient flux, and number of illumination cycles. Our
method again outperforms all baselines by a significant margin across all test
conditions. Notably, the baselines fail or perform considerably worse under
high ambient flux, as signal photons are blocked by background photons due to
pile-up. By contrast, our method is robust against a broad range of ambient
flux levels, as we model the effects of ambient flux directly.

Figure B: Off-the-shelf algorithms for surface reconstruction (Poisson Surface
Reconstruction kazhdan2006poisson_ and The Ball-Pivoting Algorithm
Bernardini1999ball_) do not perform well on point clouds generated by
reprojection.

##  C Point Cloud Surface Reconstruction

We do not apply surface reconstruction to the point clouds computed by
reprojection because off-the-shelf reconstruction techniques do not perform
reliably on the generated point clouds, as shown in Figure B. Additionally,
calculating Chamfer distance to the computed point cloud ensures that we are
capturing the efficacy of reprojection rather than a given surface
reconstruction method. We include one-way reconstruction to ground truth
Chamfer distance in Table 2 to provide a metric which does not penalize the
sparsity of point clouds produced by reprojection.

##  D Additional Qualitative Results

More reconstruction results. We present additional qualitative results on
simulated data (Figure D) and real-world captures (Figure C). These results
were omitted from the main paper due to lack of space.

Surface normal visualization. Moving beyond the 3D shapes, we further examine
the surface normal of our reconstructed 3D objects. The surface normal of a
point $\mathbf{x}$ on the reconstructed mesh is estimated as

$\tilde{\mathbf{n}}_{\mathbf{x}}=\frac{\nabla_{\mathbf{x}}(f_{\theta}(\mathrm{PE(\mathbf{x})}))}{\|\nabla_{\mathbf{x}}(f_{\theta}(\mathrm{PE(\mathbf{x})}))\|},$
(A)

where $f_{\theta}$ is the learned SDF and PE denotes the positional encoding
function. The error $e_{\mathbf{x}}$ w.r.t. the ground-truth normal
$\mathbf{n}_{\mathbf{x}}$ is given by

$e_{\mathbf{x}}=|\langle\mathbf{n}_{\mathbf{x}},\tilde{\mathbf{n}}_{\mathbf{x}}\rangle|.$
(B)

Surface normal results. We provide visualizations of surface normals for
simulated data in Figure E. Our method can successfully recover smoothly
varying normals. Error typically occurs at edges and depth discontinuities
with fast-changing normals. We hypothesize that sensors with higher temporal
and spatial resolution are needed for more accurate surface normal
reconstruction.

##  E Further Discussion

Beyond Lambertian objects. Our method assumes a spatially uniform Lambertian
BRDF, but in practice can effectively reconstruct objects with spatially
varying albedo and slightly glossy appearance (_e.g_. the spray bottle). In
theory, our method can easily be adapted to incorporate a parametric lighting
model. Recovery of the parameters of such a model are likely possible because,
by sharing information among many observations, the BRDF is effectively
sampled at many incident and exitant angles. An intriguing direction for
future work is investigating which BRDF parameterizations can be recovered
with our imaging setup, and the effect of the reflectance model on
reconstruction quality. We suspect that a non-parametric NeRF-like BRDF would
not be suitable as it does not sufficiently constrain the optimization. A
parametric lighting model, _e.g_. Phong bui1975phong_ or Oren-Nayar
oren1994oren-nayar_ may appropriately constrain the optimization while
allowing the model to learn a more accurate scene representation.

Runtime efficiency. Our method takes on the order of hours to reconstruct a
scene, making it unsuitable for real-time applications in its current state.
Future work should investigate ways to speed up forward rendering and model
training. Improved importance sampling would likely yield modest improvements
in convergence time. Another option is to render only summary statistics of
the histogram (_e.g_. mean, peak locations or widths) rather than the entire
histogram, which would likely be faster to render at the expense of yielding a
lower-quality reconstruction.

Sensor pose. In this work, we used an industrial robot arm to gather posed
sensor measurements. We chose this modality as it is guaranteed to provide
highly accurate sensor poses, and allows control over precise sensor
placement. For applications like wearable computing and distributed sensing
for robotics, camera poses might be pre-calibrated and remain fixed relative
to each other during operation. Alternatively, the low-cost single-photon
camera could be combined with a sensor-based localization system (_e.g_., an
IMU based yi2007imu_ or a camera based mur2015orb_ system) to recover camera
pose, a setup which is standard in related works Mildenhall2020nerf_,
OrtiziSDF2022_. Such a capture setup would allow capture of more organic and
large scale scenes, which more closely mimic the potential use cases of the
sensor (_e.g_. on mobile robots and drones).

Comparison to other 3D imaging modalities. Our work provides a low-cost 3D
imaging system using single-photon cameras. We provide detailed comparisons
between our method and baseline methods, but do not compare our
reconstructions to those gathered from other 3D modalities, such as continuous
wave time-of-flight attal2021torf or LiDAR Huang2023nfl. Future work should
provide a comparison to these other modalities to provide insights into the
niche (in terms of accuracy, size, power, etc.) filled by each.

Comparison to a NeRF with a miniature RGB camera. Using multiple calibrated
RGB cameras, NeRF-like approaches can leverage photometric cues (_e.g_.
correspondence and shading) to recover 3D shapes from a distributed set of
cameras. _Conceptually_ , NeRF-like approaches fall short under non-ideal
lighting (_e.g_. low-light) or insufficient correspondence (_e.g_. textureless
objects), while our method with active lighting and depth remains effective.
With sufficient light and distinct texture patterns, NeRF-like approaches on
RGB images will yield a higher quality reconstruction than our method, due to
significantly increased data rates. Without compression, a VGA-resolution
image with 8-bit color channels contains 7Mbits of information, while an image
(_i.e_., transient histogram) from our sensor contains 3Kbits of information.
With $2400\times$ the information per-view, it is not surprising that NeRF
with RGB cameras could outperform our method given the same number of input
views. _Practically_ , existing low-cost RGB cameras are larger and less power
efficient than the SPAD sensors that we utilize. While miniature RGB cameras
do exist (_e.g_., those used for endoscopy), they are $>$20$\times$ the cost
of a SPAD.

Commodity sensors. One challenge for future work is a lack of hardware support
for measurement and use of transient histograms. Very few low-cost sensors
allow access to transient histograms, and those that do often perform pre-
processing that is proprietary or undocumented. Additionally, most sensors are
equipped with very low-bandwidth I2C interfaces, limiting their effective FPS.
We hope that manufacturers will see value in users having access to transient
histogram data and support the use of this data with documentation, low-level
access, and high-bandwidth interfaces in the future.

Ethical concerns. Our work presents a new method for imaging 3D objects with
low-cost single-photon cameras. We do not anticipate major ethical concerns.

Figure C: Additional qualitative results on real-world captures. Our method
achieves the highest reconstruction quality. Poses in column two are
subsampled by a factor of two for clarity. Figure D: Additional qualitative
results on simulated data. Our method achieves the highest reconstruction
quality. Space carving captures an envelope of the shape, and may carve away
occupied areas in concave shapes (_e.g_. Armadillo). Reprojection gives a
sparse reconstruction of convex shapes, (_e.g_. skull, soap, sphere), the
scale of which may be distorted due to biases introduced by the wide field-of-
view of the sensor. Figure E: Visualizations of surface normals for simulated
data. Our method correctly estimates surface normals in flat regions. Error
mainly occurs at edges and depth discontinuities. We hypothesize that sensors
with higher temporal and spatial resolution are needed to detect rapid changes
in surface normals.

ieeenat_fullname main


11institutetext:  Université publique, Marseille, France 22institutetext: Aix-
Marseille Univ., CNRS, LIS, Marseille, France

# Reducing block-sequential automata networks to _smaller_ parallel networks
with isomorphic limit dynamics

Pacôme Perrotin  1122 Sylvain Sené  1122

# Turning block-sequential automata networks into smaller parallel networks
with isomorphic limit dynamics

Pacôme Perrotin  1122 Sylvain Sené  1122

###### Abstract

We state an algorithm that, given an automata network and a block-sequential
update schedule, produces an automata network of the same size or smaller with
the same limit dynamics under the parallel update schedule. Then, we focus on
the family of automata cycles which share a unique path of automata, called
tangential cycles, and show that a restriction of our algorithm allows to
reduce any instance of these networks under a block-sequential update schedule
into a smaller parallel network of the family and to characterize the number
of reductions operated while conserving their limit dynamics. We also show
that any tangential cycles reduced by our main algorithm are transformed into
a network whose size is that of the largest cycle of the initial network. We
end by showing that the restricted algorithm allows the direct
characterization of block-sequential double cycles as parallel ones.

##  1 Introduction

Automata networks are classically used to model gene regulatory networks [9,
16] [10, 2, 4]. In these applications the dynamics of automata networks help
to understand how the biological systems might evolve. As such, there is
motivation in improving our computation and characterization of automata
networks dynamics. This problem is a difficult one to approach considering the
vast diversity of network structures, local functions and update schedules
that are studied. Rather than considering the problem in general, we look for
families or properties which allow for simpler dynamics that we might be able
to characterize [7, 8].

We are interested in studying the limit dynamics of automata networks, that
is, the limit cycles and fixed points that they adopt over time, notably since
these asymptotic behaviors of the underlying dynamical systems may correspond
to real biological phenomenologies such as the genetic expression patterns of
cellular types, tissues, or paces. More precisely, we are less interested in
the possible configurations themselves than in the information that is being
transfered and computed in networks over time. As such, given families of
networks, one of our objectives is to count the fixed points and limit cycles
they possess.

In this paper, we provide an algorithm that, given an automata network and a
block-sequential update schedule, produces an automata network of the same
size or smaller with isomorphic limit dynamics under the parallel update
schedule. After definitions in Section 2, this algorithm is detailed in
Section 3. In Section 4, the feasibility of the algorithm on _Tangential
Cycles_ (TC) is studied, a TC being a set of cycles that intersect on a shared
path of automata.

_Why focusing on TCs?_ Cycles are fundamental retroactive patterns that are
necessary to observe complex dynamics [14]. They are present in many
biological regulation networks [17] and are perfectly understood in isolation
[6, 12]. In theory, cycles generate an exponential amount of limit cycles,
which is incoherent with the observed behavior of biological systems [9]. The
only way to reduce the amount of limit cycles is to constrain the degrees of
freedom induced by isolated cycles, which can only be done by intersecting
cycles from the purely structural standpoint. This leads us naturally to TCs,
as a simple intersection case. Double cycles (intersections of two isolated
cycles) in particular are the largest family of intersecting cycles for which
a complete characterization exists [12, 5]; the present paper generalizes this
result to block-sequential update schedules. Moreover, from the biological
standpoint, double cycles are also observed in biological regulation networks,
in which they seem to serve as inhibitors of their limit behavior [3].

##  2 Definitions

Let $\Sigma$ be a finite alphabet. We denote by $\Sigma^{n}$ the set of all
words of size $n$ over the alphabet $\Sigma$, such that for all $1\leq i\leq
n$ and $x\in\Sigma^{n}$, $x_{i}$ is the $i$th letter of that word. An
_automata network (AN)_ is a function $F:\Sigma^{n}\to\Sigma^{n}$, where $n$
is the size of the network. A configuration of $F$ is a word over
$\Sigma^{n}$. The global function $F$ can be divided into functions that are
local to each automaton: $\forall k,f_{k}:\Sigma^{n}\to\Sigma$, and the global
function can be redefined as the parallel application of every local function:
$\forall 1\leq i\leq n,F(x)_{i}=f_{i}(x)$. For convenience, the set of
automata $\\{1,\ldots,n\\}$ is denoted by $S$, and will sometimes be
considered as a set of letters rather than numbers. For questions of
complexity, we consider that _local functions are always encoded as circuits._

For $(i,j)$ any pair of automata, $i$ is said to _influence_ $j$ if and only
if there exists a configuration $x\in\Sigma^{n}$ in which there exists a state
change of $i$ that changes the state of $f_{j}(x)$. More formally, $i$
influences $j$ if and only if there exists $x,x^{\prime}\in\Sigma^{n}$ such
that $\forall k,x_{k}=x^{\prime}_{k}\Leftrightarrow k\neq i$ and $f_{j}(x)\neq
f_{j}(x^{\prime})$.

It is common to represent an automata network $F$ as the digraph with its
automata as nodes so that $(i,j)$ is an edge if and only if $i$ influences
$j$. This digraph is called the _interaction digraph_ and is denoted by
$G_{I}(F)=(S,E)$, with $E$ the set of edges. The automata network described in
Example 1 is illustrated as an interaction digraph in Figure 1.

###### Example 1

Let $F:\mathbb{B}^{3}\to\mathbb{B}^{3}$ be an AN with local functions

$\displaystyle f_{a}(x)$ $\displaystyle=\neg x_{b}\vee x_{c}$ $\displaystyle
f_{b}(x)$ $\displaystyle=x_{a}$ $\displaystyle f_{c}(x)$ $\displaystyle=\neg
x_{b}$

acb Figure 1: Interaction digraph of the AN detailed in Example 1.

An _update schedule_ is an infinite sequence of non-empty subsets of $S$,
called blocks. Such a sequence describes in which order the local functions
are to be applied to update the network, and there are uncountably infinitely
many of them. A _periodic update schedule_ is an infinite periodic sequence of
non-empty subsets of $S$, which we directly define by its period. The
application of an update schedule on a configuration of a network is the
parallel application of the local functions of the subsets in the sequence,
each subset being applied one after the other.

For example, the sequence $\pi=(S)$ is the parallel update schedule. It is
periodic, and its application on a configuration is undistinguishable from the
application of $F$. The sequence $(\\{1\\},\ldots,\\{n\\})$ is also a periodic
update schedule, and implies the application of every local function in order,
one at a time.

Formally, the application of a periodic update schedule $\Delta$ to a
configuration $x\in\Sigma^{n}$ is denoted by the function $F_{\Delta}$, and is
defined as the composition of the applications of the local functions in the
order specified by $\Delta$. For any subset $X\subseteq S$, updating $X$ into
$x$ is denoted by $F_{X}(x)$ and is defined as

$\forall i\in S,\ F_{X}(x)_{i}=\begin{array}[]{\\{}{ll}.f_{i}(x)&\text{if
}i\in X\\\ x_{i}&\text{otherwise}\end{array}\text{.}$

Example 2 provides an example of the execution of the network detailed in
Example 1 under some non-trivial update schedule.

###### Example 2

Let $\Delta=(\\{b,c\\},\\{a\\},\\{a,b\\})$ be a periodic update schedule, and
let $x=000$ be an initial configuration. For $F$ the AN detailed in Example 1,
we have that:

$\displaystyle F_{\Delta}(000)$ $\displaystyle=(F_{\\{a,b\\}}\circ
F_{\\{a\\}}\circ F_{\\{b,c\\}})(000)$ $\displaystyle=(F_{\\{a,b\\}}\circ
F_{\\{a\\}})(001)$ $\displaystyle=F_{\\{a,b\\}}(101)=111.$

A _block-sequential update schedule_ is a periodic update schedule where all
the subsets in a period form a partition of $S$; that is, every automaton is
updated exactly once in the sequence. If every subset in the sequence is of
cardinality $1$, the update schedule is said to be sequential. For any AN with
automata $S$, both the parallel update schedule and the $|S|!$ different
sequential update schedules are block-sequential. Block-sequential update
schedules are _fair_ update schedules, in the sense that applying it updates
each automaton the same amount of times.

The application of a block-sequential update schedule on an AN can be
otherwise represented as an update digraph, introduced in [15, 1]. For $F$ an
AN and $\Delta$ a block-sequential update schedule, the _update digraph_ of
$F_{\Delta}$, denoted by $G_{U}(F_{\Delta})$, is an annotation of the
network’s interaction digraph, where any edge $(u,v)$ is annotated with $<$ if
$u$ is updated strictly before $v$ in $\Delta$, and with $\geqslant$
otherwise. An update digraph of the AN detailed in Example 1 is illustrated in
Figure 2.

acb$<$$<$$\geqslant$$\geqslant$ Figure 2: Update digraph of the AN detailed in
Example 1, for $\Delta=(\\{a\\},\\{b\\},\\{c\\})$.

Given an automata network $F$ and a periodic update schedule $\Delta$, we
define the _dynamics_ of $F_{\Delta}$ as the digraph with all configurations
$x\in\Sigma^{n}$ as nodes, so that $(x,y)$ is an edge of the dynamics if and
only if $F_{\Delta}(x)=y$. We call _limit cycle of length $k$_ any sequence of
unique configurations $(x_{1},x_{2},\ldots,x_{k})$ such that
$F_{\Delta}(x_{i})=x_{i+1}$ for all $1\leq i<k$, and
$F_{\Delta}(x_{k})=x_{1}$. A limit cycle of length one is called a _fixed
point_. The _limit dynamics_ of $F_{\Delta}$ is the subgraph which contains
only the limit cycles and the fixed points of the dynamics. The limit dynamics
of the network defined in Example 1 are emphasized in Figure 3.

000001100010101011110111

000001100010101011110111

Figure 3: Two dynamics of the AN $F$ detailed in Example 1. On the left, the
dynamics of $F$ under the parallel update schedule. On the right, the dynamics
of $F$ under the update schedule $\Delta=(\\{a\\},\\{b\\},\\{c\\})$. The limit
dynamics are depicted with bold arrows.

Since the dynamics of a network is a graph that is exponential in size
relative to the number of its automaton, naively computing the limit dynamics
of a family of network is a computationally expensive process.

##  3 The algorithm

In this section, we look at an algorithm that can turn any automata network
$F$ with a block-sequential update schedule $\Delta$ into another automata
network $F^{\prime}$, such that the limit dynamics of $F_{\Delta}$ stays
isomorphic to the limit dynamics of $F^{\prime}$ under the parallel update
schedule $\pi$. Furthermore, the size of $F^{\prime}$ will always be the size
of $F$, _or less_.

This algorithm is built from two parts: first, we parallelize the network
thanks to a known algorithm in the folklore of automata networks. Second, we
remove automata from the networks based on redundancies created in the first
step.

First, let us state the usual algorithm that, given an automata network $F$
and a block-sequential update schedule $\Delta$, provides a new automata
network $F^{\prime}$ defined on the same set of automata, such that
$F_{\Delta}$ and $F^{\prime}_{\pi}$ have the same exact dynamics.

Algorithm 1 Parallelization algorithm of $F_{\Delta}$

Input

$F$local functions of a network over $S$, encoded as circuits

$\Delta$block-sequential update schedule over $S$

Output

$F$local functions of a parallel network over $S$, encoded as circuits

for $(u,v)$ such that $u$ precedes $v$ in $\Delta$ do

apply the substitution $x_{u}\mapsto\theta_{u}$ in $f_{v}$ $\triangleright$
$\theta$ is a temporary symbol

let $X\leftarrow S$

while $|X|>0$ do

let $s\in X$ such that $f_{s}$ contains no $\theta$ symbol

$X\leftarrow X\setminus\\{s\\}$

for $s^{\prime}\in X$ do

if $f_{s^{\prime}}$ contains $\theta_{s}$ then

apply the substitution $\theta_{s}\mapsto f_{s}$ in $f_{s^{\prime}}$

return $F$

Algorithm 1 proceeds with two waves of subsitutions. First, for every $<$-edge
$(u,<,v)$, the influencing automaton $u$ is replaced in the local function of
$v$ by a token symbol $\theta_{u}$. All of these token symbols are then
replaced by the corresponding local functions (in this case, $f_{u}$) in the
correct order: that is, no function is ever used in a substitution if it
contains a token character. This way, even if the network contains a complex
tree of $<$-edges, the substitutions will be applied in the correct order. It
holds that this algorithm always returns, and runs in polynomial time.

###### Property 1

Algorithm 1 always returns, and does so in polynomial time.

###### Proof

Let us denote by _$ <$-graph_ the subset of the graph $G_{U}(F_{\Delta})$
where only the $<$-edges have been preserved. The $<$-graph of $F_{\Delta}$ is
always a tree (or multiple disconnected trees): if this wasn’t the case, there
would be a cycle of $<$-edges in $G_{U}(F_{\Delta})$, which would mean a cycle
of automata that are all updated strictly before their out-neighbor, which is
impossible.

Algorithm 1 will place a $\theta$ symbol for every edge in the $<$-graph. In
the second loop, the selected $s$ is always a leaf of one of the trees
contained in the $<$-graph. The applied substitution removes that leaf from
the $<$-graph. By the structure of a tree, all the $<$-edges will be removed
and the algorithm terminates.

To see that this algorithm can be performed in polynomial time, consider that
all of the local functions are encoded as circuits. As such, it is enough to
prepare a copy of each local function into one large circuit, on which every
substitution will be performed. Any substitution $x_{u}\mapsto\theta_{u}$ is
performed by renaming the corresponding input gate. Any substitution
$\theta_{s}\mapsto f_{s}(x)$ is performed by replacing the input gate which
corresponds to $\theta_{s}$ by a connection to the output gate of the circuit
that computes the local function $f_{s}$. These substitutions are performed
for every $<$-edge in the update digraph of $F_{\Delta}$, which can be done by
doing one substitution for every pair in the partial order provided by
$\Delta$, which is never more than $n^{2}$. The resulting circuit is then
duplicated for every automaton in the output network, which leads to a total
size of no more than $k^{2}$, for $k$ the size of the input. ∎

###### Remark 1

This algorithm is not polynomial if the local functions are encoded as
formulæ, which is a detail often overlooked in the literature where this
parallelization algorithm is always assumed to be polynomial.

###### Theorem 3.1

For any $F_{\Delta}$ Algorithm 1 returns a network $F^{\prime}$ such that the
dynamics of $F_{\Delta}$ is equal to that of $F^{\prime}_{\pi}$.

###### Proof

Let us consider some configuration $x\in\Sigma^{n}$, and let us compute its
image $x^{\prime}$ in both systems. Let us consider the initial block $X_{0}$
in $\Delta$. For any automaton in $X_{0}$, its local function is untouched in
$F^{\prime}$, and thus $F_{\Delta}(x)|_{X_{0}}=F^{\prime}(x)|_{X_{0}}$.
Suppose that $F_{\Delta}(x)|_{X_{0}\cup\ldots\cup
X_{k}}=F^{\prime}(x)|_{X_{0}\cup\ldots\cup X_{k}}$ for some $k$, let us prove
that is true when including the next block $X_{k+1}$.

Let $v\in X_{k+1}$. By the nature of updates in $\Delta$, $f_{v}$ will be
updated using the values in $F_{\Delta}(x)$ for any $x_{u}$ such that $u\in
X_{0}\cup\ldots\cup X_{k}$, and in $x$ otherwise. In $F^{\prime}$, in the
local function $f^{\prime}_{v}$ and for any $u\in X_{0}\cup\ldots\cup X_{k}$
that influences $v$, a substitution has replaced $x_{u}$ by
$f^{\prime}_{u}(x)$, which implies that the value of $v$ will be updated using
a value of $u$ in $F^{\prime}(x)$. Pulling this together, we obtain that
$f_{v}(x)=f^{\prime}_{v}(x)$ and $F_{\Delta}(x)|_{X_{0}\cup\ldots\cup
X_{k+1}}=F^{\prime}(x)|_{X_{0}\cup\ldots\cup X_{k+1}}$, and the recurrence
yields $F_{\Delta}(x)=F^{\prime}(x)$ for any $x$. ∎

Algorithm 2 Parallelization algorithm of $F_{\Delta}$, with a possible
reduction in size

Input

$F$local functions of a network over $S$, encoded as circuits

$\Delta$update schedule over $S$

Output

$F^{\prime}$local functions of a parallel network over a subset of $S$,

encoded as circuits

let $F^{\prime}\leftarrow$ apply Algorithm 1 to $F_{\Delta}$

let $G_{I}(F^{\prime})\leftarrow$ the interaction digraph of $F^{\prime}$

for $(u,v)\in S^{2}$ do

if $\forall x\in\Sigma^{n},f_{u}(x)=g(f_{v}(x))$ then $\triangleright$ for
some $g:\Sigma\to\Sigma$

for $(u,w)\in E(G_{I}(F^{\prime}))$ do

apply the substitution $x_{u}\mapsto g(x_{v})$ in $f_{w}$

while $\exists u\in S$ such that $u$ has no accessible neighbor in
$G^{\prime}_{I}$ do

$S\leftarrow S\setminus\\{u\\}$ $\triangleright$ $u$ is removed from the
network

return $F^{\prime}$

Algorithm 2 is our contribution to this process, and removes automata that are
not necessary for the limit dynamics of the network. It proceeds in two steps:
first, the algorithm identifies pairs of automata with equivalent local
functions, up to some function. In other terms, if one automaton $u$ can be
computed as a function $g$ of the local function of another automaton $v$,
then $u$ is not necessary and all references to $x_{u}$ in the network can be
replaced by $g(x_{v})$ for an identical result. Of course, this only works
under the hypothesis that $u$ and $v$ are updated synchronously, which is the
case after the application of Algorithm 1. Second, the algorithm iteratively
removes any automaton that has no influence in the network, that is, that has
no accessible neighbor in the interaction graph of the network. These automata
are not part of cycles and do not lead to cycles, and as such have no impact
on the attractors.

Algorithm 2 is non-deterministic, and when the local functions of any pair of
automata $(u,v)$ are shown to be equivalent up to some reversible function
$g:\Sigma\to\Sigma$, either automata could replace the influence of the other
without preference. As such, more than one result network is possible, but all
are equivalent in their limit dynamics, as will be shown later.

While it is clear that Algorithm 2 always terminates, its complexity is out of
the deterministic polynomial range, as applying it implies solving the coNP-
complete decision problem of testing if two Boolean formulæ are equal for all
possible pairs of automata and for every possible function
$g:\Sigma\to\Sigma$. As such, a polynomial implementation of this algorithm
would (at least) imply P $=$ NP. This drastic conclusion is softened when
looking at restricted classes of networks where redundancies can be easily
pointed out, which is the case for the rest of the paper.

###### Theorem 3.2

For any $F_{\Delta}$, Algorithm 2 returns a network $F^{\prime}$ such that the
limit dynamics of $F_{\Delta}$ and $F^{\prime}_{\pi}$ are isomorphic.

###### Proof

By Theorem 3.1, the network $F^{\prime}$ returned by the application of
Algorithm 1 to $F_{\Delta}$ has identical dynamics to $F_{\Delta}$.

Algorithm 2 operates two kinds of modifications.

The first operation is replacing the influence of any automaton $u$ by another
automaton $v$ if they are found to have equivalent local function up to some
$g:\Sigma\to\Sigma$, that is, $f_{u}=g\circ f_{v}$. For any configuration $x$,
the value of $f_{u}(x)$ and $g(f_{v}(x))$ are always equal. Thus, substituting
the variable $x_{u}$ by $g(x_{v})$ in the local functions of every out-
neighbor of $u$ will lead to an identical limit behavior. After this
substitution, the automaton $u$ does not have any influence over the network.
Morever, all its previous out-neighbors in $G_{I}(F^{\prime})$ are now the
out-neighbors of $v$.

The second operation is iteratively removing automata that do not influence
any automaton. Let $u$ be such an deleted automaton. Consider a limit cycle
$(x^{1},x^{2},\ldots,$ $x^{k})$ in $G$. By definition of a limit cycle,
$G(x^{i})=x^{i+1}$ for any $i$, $G(x^{k})=x^{1}$, and $x^{i}=x^{j}\Rightarrow
i=j$. Consider the component $x_{u}^{i}$ for some $i$. Since $u$ does not
influence any automaton, $x^{i+1}$ is a function of
$x^{i}|_{S\setminus\\{u\\}}$. As the entire sequence is aperiodic, the
sequence of the subconfigurations $x^{i}|_{S\setminus\\{u\\}}$ is also
aperiodic, and the attractor is preserved in $F^{\prime}$. ∎

##  4 Reductions in size of tangential cycles

In this section, we characterize the reduction in size that our algorithm
provides on a specific family of networks. We call _tangential cycles_ (TC)
any AN composed of any number of cycles $\\{C_{1},C_{2},\ldots,C_{k}\\}$ such
that a unique path of automata, called the _tangent_ , is shared by all of the
cycles. The first automaton of the tangent is the only automaton with more
than one in-neighbor, and is called the _central automaton_. A TC is
represented as part of Figure 4, which contains three cycles and a tangent of
length $0$ (only one node is shared between the cycles).

###  4.1 Reducing block-sequential TCs

The reduction in size provided by Algorithm 2 can be quite large on TCs, as
even TCs updated in parallel can be reduced in size by merging the different
cycles as much as possible. As such, the reduction power of this algorithm is
greater than just removing the redundancies inherent to the block-sequential
to parallel update translation. Indeed, Figure 4 provides an example of a
parallel TC, the size of which is greatly reduced by the application of
Algorithm 2. But, by this process, the final result of Algorithm 2 is no
longer a TC.

As explained above, TCs are studied as the next simplest cases of complex ANs
that make biological sense, after automata cycles. Both isolated cycles and
double cycles are examples of TCs. To show that the study of TCs under block-
sequential update schedules can be directly reduced to the study of TCs under
the parallel update schedule, we provide an algorithm that transforms any TC
under a block-sequential update schedule into a TC under the parallel update
schedule, such that their limit dynamics are isomorphic, and the local
functions of their central automaton equivalent. This is done by simply
stopping the process of Algorithm 2 earlier to preserve the TC shape of the
network.

Algorithm 3 Parallelization algorithm of a TC $F$ under the block-sequential
update schedule $\Delta$, with a possible reduction in size

Input

$F$local functions of a network over $S$, encoded as circuits

$\Delta$update schedule over $S$

Output

$F^{\prime}$local functions of a parallel network over a subset of $S$,

encoded as circuits

let $F^{\prime}\leftarrow$ apply Algorithm 1 to $F_{\Delta}$

let $G_{I}(F^{\prime})\leftarrow$ the interaction digraph of $F^{\prime}$

for $(u,v)\in S^{2}$, such that either $u$ or $v$ has more than one in-
neighbor in $G_{I}(F^{\prime})$ do

if $\forall x\in\Sigma^{n},f_{u}(x)=g(f_{v}(x))$ then $\triangleright$ for
some $g:\Sigma\to\Sigma$

for $(u,w)\in E(G_{I}(F^{\prime}))$ do

apply the substitution $x_{u}\mapsto g(x_{v})$ in $f_{w}$

for $u\in S$ do

if $u$ has no accessible neighbors in $G_{I}(F^{\prime})$ then

$S\leftarrow S\setminus\\{u\\}$ $\triangleright$ $u$ is removed from the
network

return $F^{\prime}$

abcdeh$<$$\geqslant$$\geqslant$$\geqslant$$\geqslant$$\geqslant$$<$$\geqslant$

$\displaystyle f_{a}(x)$ $\displaystyle=x_{a}\vee x_{d}\vee\neg x_{h}$
$\displaystyle f_{b}(x)$ $\displaystyle=\neg x_{a}$ $\displaystyle f_{c}(x)$
$\displaystyle=x_{b}$ $\displaystyle f_{d}(x)$ $\displaystyle=x_{c}$
$\displaystyle f_{e}(x)$ $\displaystyle=x_{a}$ $\displaystyle f_{h}(x)$
$\displaystyle=x_{e}$

Initial network, $\Delta=(\\{h\\},\\{a,c,d,e\\},\\{b\\})$

N / A

$\displaystyle f_{a}(x)$ $\displaystyle=x_{a}\vee x_{d}\vee\neg\theta_{h}$
$\displaystyle f_{b}(x)$ $\displaystyle=\neg\theta_{a}$

$\ldots$

Algorithm 1, first loop (c, d, e and h are unchanged)

abcdeh

$\displaystyle f^{\prime}_{a}(x)$ $\displaystyle=x_{a}\vee x_{d}\vee\neg
x_{e}$ $\displaystyle f^{\prime}_{b}(x)$ $\displaystyle=\neg(x_{a}\vee
x_{d}\vee\neg x_{e})$ $\displaystyle f^{\prime}_{c}(x)$ $\displaystyle=x_{b}$
$\displaystyle f^{\prime}_{d}(x)$ $\displaystyle=x_{c}$ $\displaystyle
f^{\prime}_{e}(x)$ $\displaystyle=x_{a}$ $\displaystyle f^{\prime}_{h}(x)$
$\displaystyle=x_{e}$

Algorithm 1, second loop

adce

$\displaystyle f^{\prime}_{a}(x)$ $\displaystyle=x_{a}\vee x_{d}\vee\neg
x_{e}$ $\displaystyle f^{\prime}_{c}(x)$ $\displaystyle=\neg x_{a}$
$\displaystyle f^{\prime}_{d}(x)$ $\displaystyle=x_{c}$ $\displaystyle
f^{\prime}_{e}(x)$ $\displaystyle=x_{a}$

Algorithm 2 after the deletion of $h$ and the merge of $b$ into $a$. This is
where Algorithm 3 ends

adc

$\displaystyle f^{\prime}_{a}(x)$ $\displaystyle=x_{a}\vee x_{d}\vee\neg
x_{c}$ $\displaystyle f^{\prime}_{c}(x)$ $\displaystyle=x_{a}$ $\displaystyle
f^{\prime}_{d}(x)$ $\displaystyle=\neg x_{c}$

End of Algorithm 2, after the merge of $e$ into $c$

Figure 4: Application of Algorithm 1, 2 and 3 on an example network. Different
steps of the algorithm are represented and separated using horizontal lines.
At each step, the interaction graph or update graph and the local functions
are the result of the operations indicated on the right. As the initial
network is a TC, the fourth step represents the result returned by Algorithm
3, which is a TC of smaller size. The fifth step represents the result
returned by Algorithm 2, which is not a TC. Dashed lines in the interaction
digraph connect automata the local function of which are equivalent up to a
negation. Only the first graph is represented as an update digraph, as all the
other networks are updated in parallel.

The only difference between Algorithms 2 and 3 is that the latter restricts
the reductions it operates. If two local functions are found to be equivalent
up to some function $g$, Algorithm 3 removes a node if and only if these local
functions are duplicates of the previous local function of the central
automaton of the network. Removing duplications of any function that is part
of a cycle would merge two cycles and the network would no longer be
tangential cycles, in a way that is harder to count the reductions for. Since
Algorithm 3 is a variation of Algorithm 2 that only does less reductions,
Theorem 3.2 still applies in its case. An application of Algorithms 2 and 3 is
illustrated in Figure 4 and the difference between the algorithms is
highlighted.

###### Theorem 4.1

Let $F$ be a TC and $\Delta$ a block-sequential update schedule. The amount of
reductions in size that Algorithm 3 operates on $F_{\Delta}$ is the number of
$<$-edges in the update digraph of $F_{\Delta}$, and the result is a TC.

###### Proof

Algorithm 1 operates a substitution for every $<$-edge in the update digraph
of $F_{\Delta}$. In this proof, we will show that each of the possible
transformations implies the removal of exactly one node from the network.

For any such edge $(u,<,v)$, there are two cases. Either $u$ is the central
automaton, or not. In any case, $u\neq v$ since the contrary would imply that
an automaton is updated strictly before itself.

If we suppose that $u$ is the central automaton, this means that $f_{v}$ is a
local function that only depends on $x_{u}$. It can thus be written
$f_{v}(x)=g(x_{u})$ for some $g:\Sigma\to\Sigma$. After the application of
Algorithm 1, we thus obtain $f_{v}(x)=g(f_{u}(x))$, which implies the removal
of either $u$ or $v$ (but at this point, not both) by Algorithm 3.

If we suppose that $u$ is not the central automaton, this means that $f_{v}$
is an arbitrary formula which contains $x_{u}$, and $f_{u}$ is a function of
the form $f_{u}(x)=g(x_{w})$ for some $g$ and some $w\in S$. Note that $w\neq
u$ by the hypothesis that $F$ is a TC, either $w$ is the previous automaton in
the path, or it is the central automaton $v$. As such, applying Algorithm 1
substitutes any mention of $x_{u}$ in $f_{v}$ by $g(x_{w})$. Previously, $u$
only had one accessible neighbor, as it was part of a path connecting to the
central automaton. This leaves $x_{u}$ without any accessible neighbors in the
interaction digraph of $F$, which means that it is removed by Algorithm 3. If
the removed edge is part of a cycle, this means that this cycle will be
reduced in size. If the edge is part of the tangent, this means that the
tangent will be reduced in size.

We thus obtain that the number of reductions is at least the number of
$<$-edges in the update digraph of the network. Suppose now that some extra
automaton $u$ is removed on top of any $<$-edge related reduction. First
observe that if $u$ has no accessible neighbor, it must have had none from
before the application of Algorithm 1, since in none of the two cases are
external automaton disconnected from each other. Now suppose that $f_{u}$ is
equivalent to some $f_{v}$ up to some $g$. Neither $u$ nor $v$ can be the
central automaton, as any duplication of that function is handled in the first
case. This proves that the number of reductions is exactly the number of
$<$-edges in the update digraph of $(F,\Delta)$.

Let us now show that the result of Algorithm 3 is a TC. If the initial network
had a central automata, there still exists a unique central automata at the
end of the algorithm, even if the original central automata was removed in a
chosen reduction. Paths that exit the central automata in the previous network
still exit the central automata in the result, in the same number, and still
share some tangent. The paths can be smaller in size, as well as the tangent,
but they still end in the central automata.∎

If Algorithm 2 cannot be polynomial in the worst case under the hypothesis
that $\text{P}\neq\text{NP}$, Algorithm 3 can be simplified to the following
rule: taking a TC with a block-sequential update schedule, we obtain the
equivalent parallel TC by reducing each cycle by the number of $<$-edges that
its update digraph contains. This process is quadratic, since we only need to
check the possible $<$-edges defined by the partial order defined by $\Delta$,
which are no more than $n^{2}$.

###  4.2 Reducing parallel Boolean TCs further

Applying Algorithm 2 to its full extent to a Boolean TC (That is, a TC defined
over the Boolean alphabet) may result in a larger reduction in size. As any
automaton that is not the central one has a unary function as its local
function, any pair of non-central local functions is equivalent up to some
$g:\Sigma\to\Sigma$ if they are influenced by the same automaton. For example,
if the central automaton influences three other automata that represent the
start of three chains, these three automata can be merged into one. Continuing
this zipping process yields a final network only as large as the longest cycle
of the initial TC.

This process is not straightforward for non-Boolean TCs, as the local
functions along the chains can be non-reversible using modular arithmetics,
for example. Optimizing these networks is still possible, but requires a more
complex set of substitutions to do so. It has been proven in general using
modules and output functions [13]. The following theorem corresponds to the
Boolean case, proven with more classical means. An example of its application
is illustrated in the two last steps of Figure 4.

###### Theorem 4.2

Let $F$ be a Boolean TC. Applying Algorithm 2 to $F_{\pi}$ generates a network
$F^{\prime}$ whose size is that of the largest cycle in $F$.

###### Proof

Starting from the initial TC $F$, all of the automata directly influenced by
the automata at the end of the tangent $u$ (but that are not $u$) have local
functions $f_{v}(x)=g(x_{u}),f_{w}(x)=h(x_{u}),\ldots$ for
$g,h,\ldots:\Sigma=\\{0,1\\}\to\Sigma$. All these functions $g,h,\ldots$ are
not constant functions, since the automata that they represent are influenced
by an automaton by hypothesis. Thus, they can only be defined as the identity
or the negation of $x_{u}$. As a consequence, all but one of these automata
will be removed by the algorithm as they are all equivalent up to some $g$.

This same argument can be repeated by taking all the automata influenced by
the only automaton resulting from the previous iteration, excluding the
central automaton. At each step, all of the automata at the same distance from
the central automaton are merged. Hence, at the end of this process, whatever
the choices made for merging automata along the iterative process, the
resulting AN will be compoesed of $k$ automata, with $k$ the length of the
largest cycle of $F$. ∎

##  5 An application: disjunctive double cycles

As an application of this algorithm and as an example to the algorithm’s
capacities to reduce the size of the provided network, we turn to the family
of disjunctive double cycles. Notice that the result still holds for
conjunctive double cycles since conjuctive and disjunctive cycles have
isomorphic dynamics [12, 11].

In disjunctive automata networks, an edge $(u,v)$ is signed positively if the
$x_{u}$ appears as a positive variable in $f_{v}$. An edge $(u,v)$ is signed
negatively if $x_{u}$ appears as a negative variable in $f_{v}$. A cycle is
said to be positive if it contains an even number of negative edges, and
negative otherwise.

A _disjunctive double cycle_ is an automata network with an interaction
digraph that is composed of two automata cycles that intersect in one
automaton. The local function of this central automaton is a disjunctive
clause. This family of networks is very simple to define, and is a simple and
intuitive next step after the family of Boolean automata cycles, which are
composed of a single cycle.

Both families have been characterized under the parallel update schedule [12,
7]; that is to say, given basic parameters concerning the size of the cycles,
their sign, and any integer $k$, an explicit formula (defined as a
polynomially computable function) has been given among other to count the
number of limit cycles of size $k$ of such networks under the parallel update
schedule. In this section, we extend this characterization to the block-
sequential equivalents by showing how applying our algorithm reduces the
network to a smaller instance of the same family of networks.

Furthermore, as Boolean automata cycles and disjunctive double cycles are TCs,
our method can be simplified to the simple following rules: given a TC $F$, a
block-sequential update schedule $\Delta$, count the number of $<$-edges in
the update digraph $G_{U}(F_{\Delta})$; for every cycle, substract to its size
the number of such edges it contains, while keeping its sign; the final
network, under the parallel update schedule, and the initial network under
$\Delta$ have isomorphic limit dynamics. This is a simple application of
Theorem 4.1, and of the rule of thumb deduced from Algorithm 3.

We denote by $DC(s,s^{\prime},a,b)$ the disjunctive double cycles with cycle
sizes $a,b$ and signs $s,s^{\prime}$.

###### Theorem 5.1

Let $D=DC(s,s^{\prime},a,b)$ be disjunctive double cycles, $\Delta$ a block-
sequential update schedule. For $A$ ($B$ respectively) the number of $<$-edges
on the cycle of size $a$ ($b$ respectively) in $G_{U}(F_{\Delta})$, the limit
dynamics of $D_{\Delta}$ is isomorphic to that of $D^{\prime}_{\pi}$, where
$D^{\prime}=DC(s,s^{\prime},a-A,b-B)$.

###### Proof

This is a straightforward application of Theorem 4.1.∎

##  6 Conclusion

In this paper we provide a novel algorithm which allows the reduction in size
of automata networks, in particular by passing the network from a block-
sequential to a parallel update schedule, while keeping isomorphic limit
dynamics. While this algorithm is too computationally expensive for the
general case, we study the specific family of intersection of automata cycles,
on which this algorithm is easily applied. This study allows the discovery
that all block-sequential tangential cycles have isomorphic limit dynamics to
parallel tangential cycles. Finally, we apply this fact to Boolean automata
double cycles to characterize their behavior under block-sequential update
schedules.

It seems now clear to us that the difference between the parallel update
schedule and block-sequential update schedules is that the latter changes the
timing of the information along sections of the network. In particular,
structures such as tangential cycles can be directly translated into an
equivalent parallel network with shorter cycles. We are interested in seeing
what effects this translation could have in a more general set of families of
networks, and if there exists other families in which block-sequential update
schedules lead to equivalent parallel networks which are still part of the
family.

As a perspective, we would like to characterize more redundancies that can be
removed from networks to help with the computation of their dynamics. For
example, we are currently interested in more complex compositions of automata
cycles, and have already found equivalences that show that many networks are
equivalent in their limit dynamics where complex parts of automata networks
can be moved alongside cycles without affecting the network’s limit dynamics.

Isolated paths are also a strong candidate for size reduction. Isolated paths
are paths that lead from cycles to other cycles but can only be crossed once.
Our current algorithms conserve such paths, despite it being possible to
reduce them completely without changing the limit dynamics of the network in
many cases, for example when an isolated path is the only way to go from one
part to another. We have to be careful when multiple isolated paths exit from
and join onto the same parts, as the synchronicity of the information in the
entire network must be preserved.

Acknowledgments. This work has been partially funded by ANR-18-CE40-0002 FANs
project (PP & SS), ECOS-Sud CE19E02 SyDySy project (PP & SS), and STIC-AmSud
22-STIC-02 CAMA project (SS).

## References

  * [1] J. Aracena, E. Goles, A. Moreira, and L. Salinas.  On the robustness of update schedules in Boolean networks.  Biosystems, 97:1–8, 2009.
  * [2] M.I. Davidich and S. Bornholdt.  Boolean network model predicts cell cycle sequence of fission yeast.  PLoS One, 3:e1672, 2008.
  * [3] J. Demongeot, A. Elena, M. Noual, S. Sené, and F. Thuderoz.  ”Immunetworks”, intersecting circuits and dynamics.  J. Theor. Biol., 280:18–33, 2011.
  * [4] J. Demongeot, E. Goles, M. Morvan, M. Noual, and S. Sené.  Attraction basins as gauges of robustness against boundary conditions in biological complex systems.  PLoS One, 5:e11793, 2010.
  * [5] J. Demongeot, T. Melliti, M. Noual, D. Regnault, and S. Sené.  Automata and Complexity: Essays presented to Eric Goles on the occasion of his 70th birthday, volume 42 of Emergence, Complexity, Computation, chapter On Boolean automata isolated cycles and tangential double-cycles dynamics, pages 145–178.  Springer, 2022.
  * [6] J. Demongeot, M. Noual, and S. Sené.  On the number of attractors of positive and negative Boolean automata circuits.  In Proceedings of WAINA’10, pages 782–789. IEEE Press, 2010.
  * [7] J. Demongeot, M. Noual, and S. Sené.  Combinatorics of Boolean automata circuits dynamics.  Discr. Appl. Math., 160:398–415, 2012.
  * [8] Z. Gao, X. Chen, and T. Başar.  Stability structures of conjunctive boolean networks.  Automatica, 89:8–20, 2018.
  * [9] S. A. Kauffman.  Metabolic stability and epigenesis in randomly constructed genetic nets.  J. Theor. Biol., 22:437–467, 1969.
  * [10] L. Mendoza and E. R. Alvarez-Buylla.  Dynamics of the genetic regulatory network for _Arabidopsis thaliana_ flower morphogenesis.  J. Theor. Biol., 193:307–319, 1998.
  * [11] M. Noual.  Dynamics of circuits and intersecting circuits.  In Proceedings of LATA’12, volume 7183 of LNCS, pages 433–444. Springer, 2011.
  * [12] M. Noual.  Updating Automata Networks.  PhD thesis, École Normale Supérieure de Lyon, 2012.
  * [13] P. Perrotin.  Simulation entre modèles de calcul naturel et modularité des réseaux d’automate.  PhD thesis, Aix-Marseile Université, 2021.
  * [14] F. Robert.  Itérations sur des ensembles finis et automates cellulaires contractants.  Linear Algebra and its Applications, 29:393–412, 1980.
  * [15] L. Salinas.  Estudio de modelos discretos: estructura y dinámica.  PhD thesis, Universidad de Chile, Santiago, Chile, 2008.
  * [16] R. Thomas.  Boolean formalization of genetic control circuits.  J. Theor. Biol., 42:563–585, 1973.
  * [17] R. Thomas.  On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations.  In Numerical methods in the study of critical phenomena, pages 180–193. Springer, 1981.