rowsep=1pt SMReferences

#  OpEnCam: Lensless Optical Encryption Camera

Salman S. Khan, Xiang Yu, Kaushik Mitra, Manmohan Chandraker, Francesco
Pittaluga S. S. Khan is with Rice University, Houston, TX, USA - 77005X. Yu is
with Amazon, USAK. Mitra is with the Indian Institute of Technology Madras,
Chennai, India - 600036 M. Chandraker and F. Pittaluga are with the NEC Labs
America, San Jose, CA, USA - 95110This paper was produced by the IEEE
Publication Technology Group. They are in Piscataway, NJ.

###### Abstract

Lensless cameras multiplex the incoming light before it is recorded by the
sensor. This ability to multiplex the incoming light has led to the
development of ultra-thin, high-speed, and single-shot 3D imagers. Recently,
there have been various attempts at demonstrating another useful aspect of
lensless cameras - their ability to preserve the privacy of a scene by
capturing encrypted measurements. However, existing lensless camera designs
suffer numerous inherent privacy vulnerabilities. To demonstrate this, we
develop the first comprehensive attack model for encryption cameras, and
propose OpEnCam– a novel lensless optical encryption camera design that
overcomes these vulnerabilities. OpEnCam encrypts the incoming light before
capturing it using the modulating ability of optical masks. Recovery of the
original scene from an OpEnCam measurement is possible only if one has access
to the camera’s encryption key, defined by the unique optical elements of each
camera. Our OpEnCam design introduces two major improvements over existing
lensless camera designs - (a) the use of two co-axially located optical masks,
one stuck to the sensor and the other a few millimeters above the sensor and
(b) the design of mask patterns, which are derived heuristically from signal
processing ideas. We show, through experiments, that OpEnCam is robust against
a range of attack types while still maintaining the imaging capabilities of
existing lensless cameras. We validate the efficacy of OpEnCam using simulated
and real data. Finally, we built and tested a prototype in the lab for proof-
of-concept.

###### Index Terms:

Lensless imaging, visual privacy, inverse problems.

##  I Introduction

Figure 1: Optical Encryption. Optical encryption camera encodes the scene
optically prior to capture. Decryption is possible computationally using a key
that is unique to the hardware of each encryption camera. Existing optical
encryption designs are susceptible to powerful learning-based attacks even
without the key. Proposed OpEnCam design prohibits such learning-based attacks
effectively.

Lensless imaging is a computational imaging modality that replaces the lens of
a conventional camera with a thin optical mask and computation. The addition
of the mask, which no longer has the focusing ability of the lens, leads to
multiplexing of the incoming light prior to capture. This multiplexing of
light has been exploited to develop ultra-thin cameras for 2D [1, 2] and 3D
imaging [2, 3] high-speed imaging [4] and hyper-spectral imaging [5]. Another
interesting aspect of lensless cameras that has been relatively less explored
is their ability to preserve privacy in the optical domain itself before
capture. This is made possible by the ability of lensless cameras to perform
computation in the optical space, with operations defined by their unique mask
patterns and arrangements.

With camera-based technologies becoming increasingly integrated into every
aspect of our lives, the potential for leakage of sensitive visual information
is ever-increasing. This is leading to significant privacy and security
concerns from consumers, citizens, and governments. The standard approach to
address these concerns is to encrypt sensitive image data at the sensor level,
after image capture, via software or specialized hardware. However, most
cameras lack such specialized hardware. Further, both software and hardware-
based solutions may still be susceptible to data sniffing attacks that gain
access to sensitive data before encryption. Thus, protecting the privacy of a
scene through optical processing before capture is currently the need of the
hour.

While optical image encryption has been around for many decades [6], most
existing solutions require coherent illumination for encoding and decoding,
making them unsuitable for deployment in unconstrained real-world
environments. Further, incoherent illumination solutions that have been
proposed are susceptible to ciphertext-only attacks[7] or have yet to be
realized with real hardware due to the experimental nature of the optics[8].

A few works like [9, 10, 11] have used existing lensless cameras [1] to show
optical encryption-based privacy-preserving applications. However, these works
are fundamentally limited by the inherent design flaws of the existing
lensless cameras, which were designed for imaging applications with no privacy
considerations. In this work, we look at the design of existing lensless
cameras and make two significant modifications to it that allows us to improve
their privacy-preserving ability prior to capture. Our novel lensless camera
design is called OpEnCam – Optical Encryption Camera. In OpEnCam, we first
introduce a novel double mask design using two co-axially placed masks – one
at the top of the bare sensor placed flush against it and the other a few
millimeters above it. Such a double mask design makes lensless cameras
significantly less susceptible to various attacks without compromising their
imaging ability. Secondly, we introduce a novel design strategy based on
signal processing heuristics for the two masks used in our system that further
enhances the privacy of our system. Finally, we show through experiments that
our proposed OpEnCam design still has the imaging ability of existing lensless
cameras while preserving privacy.

To summarize, we make the following contributions:

  1. 1.

A novel framework for generating lensless computational cameras that encrypt
image data before image capture, using two co-axially placed masks - a scaling
mask placed flush against the bare sensor, and a multiplexing mask placed a
few millimeters from the sensor.

  2. 2.

Novel heuristically designed random key-generators for the two masks,
providing a distinct mask pattern for each camera. We show that these key-
generators are robust against various forms of ciphertext-only attacks and
known plaintext attacks.

  3. 3.

A novel transformer-based keyed decryption approach.

  4. 4.

We validate our results both in simulation and real data collected using an
inexpensive proof-of-concept prototype built in our lab.

##  II Background

###  II-A Optical Cryptography

A cryptographic system is judged to be reliably secure only if it can
successfully survive a rigorous evaluation via cryptanalysis. The goal of
cryptanalysis is to identify potential weaknesses in encryption mechanisms
that would allow someone to decode the encrypted data, even without possessing
the confidential key. Technological progress in the realm of optical
cryptography has spurred the introduction of various specialized cryptanalysis
approaches, such as chosen plaintext attacks (CPA) [12, 13], known plaintext
attacks (KPA) [14, 15], and ciphertext-only attacks (COA) [16, 17, 18].

In a CPA scenario, the attacker can choose arbitrary plaintexts to be
encrypted and then examine the resulting ciphertexts. The primary objective is
to use this information to discover a weakness in the encryption algorithm or
even to deduce the secret key used for encryption. In contrast, in a KPA
scenario, the attacker does not have the luxury of choosing the plaintexts
that are encrypted. Instead, the attacker must work with pre-existing pairs of
plaintexts and ciphertexts. Despite this limitation, a successful KPA could
potentially reveal enough about the encryption scheme or key to decrypt other
ciphertexts encrypted with the same key or to weaken the overall security of
the cryptosystem. Finally, in a COA, the attacker has access solely to the
ciphertext—that is, the encrypted data—without any accompanying plaintext or
other additional information. The goal of the attacker is to deduce either the
plaintext, the encryption key, or details about the encryption algorithm,
based solely on the available ciphertexts. Ciphertext-only attacks are
considered one of the most challenging forms of cryptanalysis because the
attacker has minimal information to work with. Unlike chosen CPA or KPA, where
the attacker has more control or information, COA provides the least amount of
leverage for the attacker. If a cryptographic algorithm can resist a
ciphertext-only attack, it’s generally considered to be quite secure.

One of the earliest works on optical encryption is [6]. The authors place
independent white uniformly distributed phase masks in the input and Fourier
planes of a 4f correlator to encode images. However, this approach only works
for coherent illumination, and is also susceptible to autocorrelation-based
COAs and impulse-based CPAs [19]. [19] improves on [6] by adding a third
independent white uniformly distributed phase mask in the image plane of the
4f correlator to inhibit impulse-based CPAs [19]. In [20], the authors place a
random white uniformly distributed phase mask in the aperture plane of a
camera to produce an incoherent imaging system for optical image encryption
and demonstrate that such a system is susceptible to autocorrelation-based
COAs. In [8], the authors simulate an experimental 3D printed optical fiber
bundle [21] at the imaging plane of a camera to produce pixel shuffling.
However, such optical fiber bundles have only been shown for parallel fibers,
which are not suitable for pixel shuffling [21]. Our framework employs
standard optical masks that are cost-effective and easily mass-produced, and
we construct a prototype encryption camera. Works like [10, 7, 11, 9] use
single optical masks for encryption. Although these methods are practical,
they are not robust to various decryption attacks. In comparison, our solution
is not only practical but also robust to these attacks. It should also be
noted that none of the above existing works provide a complete analysis of the
different forms of attacks that need to be considered for an optical
encryption system, especially the powerful learning-based attacks that exploit
hidden priors in data.

###  II-B Lensless Imaging

In a conventional single-mask lensless camera like [1, 2], the lens of the
camera is replaced by a thin optical mask placed at some distance from the
sensor that modulates the incoming light. For a scene $X$, the measurement
recorded by a sufficiently large sensor $Y$ is given by:

$Y=P*X+N,$ (1)

where $*$ is the full-size convolutional operator (no cropping due to finite
sensor size), $P$ is the point spread function (PSF), and $N$ is additive
noise. The PSF is the response of the camera to a point source, and it depends
on the mask pattern. The optical mask can be implemented using an amplitude
mask that attenuates the incoming light [1] or a phase mask that modulates
based on diffraction [2]. Due to the large PSF size and limited angular
response of the pixels[2, 3], the convolution follows a zero-padded boundary
condition. Existing works like [2, 1, 3] use a similar model for thin 2D and
3D lensless imaging. These lensless imagers were developed for imaging
applications, and their ability to perform optical encryption has not been
fully explored.

FlatCam [1] is a lensless camera that places a separable coded amplitude mask
above a bare sensor array to enable a thin and flat form-factor imaging
device, which can simulate a conventional camera by reconstructing
conventional images from coded measurements. DiffuserCam [3] and PhlatCam [2]
replace the coded amplitude mask from FlatCam [1] with a coded phase mask for
improved light efficiency and reconstruction quality. Spectral DiffuserCam[5]
exploits the multiplexing ability of lensless imagers to do hyper-spectral
imaging. OpEnCam is similar to these methods in that a coded optical mask is
used to capture coded measurements. However, OpEnCam design goes a step
further - our aim is not only to enable high-quality reconstruction of
conventional images but also to prevent decryption attacks. To this end, we
employ a second optical scaling mask positioned flush against the bare sensor
array and propose novel mask designs for the multiplexing and scaling masks.

###  II-C Image Reconstruction

Image reconstruction is a core problem in computational imaging, and plays a
key role in lensless imaging. The problem of scene reconstruction involves the
computational recovery of the latent scene from a lensless capture.
Conventional approaches for lensless image reconstruction are mostly based on
regularized least squares [22, 1, 3, 4]. More recently, deep-learning-based
reconstruction methods have proposed [23, 24, 25, 26]. We employ both
regularized least squares and deep-learning-based reconstruction methods.

##  III Attack Model

Since optical elements behave as linear operators [27], optical encryption
cameras will always be susceptible to chosen-plaintext attacks in which an
attacker uses a large set of ciphertext-plaintext pairs to estimate the
encryption function. Naturally, this limits the practical use of optical
encryption cameras to settings where attackers lack physical access to the
camera location. Accordingly, in this paper, we assume that attackers do not
have physical access to the camera locations, but can access a camera’s
stream. In other words, attackers have access to ciphertexts (measurements)
from a camera, but not their corresponding plaintext (raw images), except for
three special cases which we detail below. This is a reasonable assumption for
many settings, such as security cameras in homes and offices.

In this paper, we consider four types of attacks: ciphertext-only attacks
(COA) and three special cases of known plaintext attacks, namely, impulse
known plaintext attack (I-KPA), uniform known plaintext attack (U-KPA), and
uniform-impulse known plaintext attack (UI-KPA). The reason we consider the
known impulse/uniform plaintext settings is that bright impulse-like
illumination sources and uniform backgrounds often appear naturally in scenes
and may be recognizable in the corresponding encrypted sensor measurement.
Hence, attackers may gain access to an approximate impulse or uniform response
of the sensor, without having physical access to the sensor.

Let $C:\mathbb{R}^{W_{1}\times H_{1}\times 3}\times\mathbb{R}^{W_{2}\times
H_{2}\times M}\rightarrow\mathbb{R}^{W_{3}\times H_{3}\times 3}$ denote an
optical encryption camera that maps a plaintext (scene)
$X\in\mathbb{R}^{W_{1}\times H_{1}\times 3}$ and key
$K\in\mathbb{R}^{W_{2}\times H_{2}\times M}$ to the ciphertext (measurement)
$Y=C(K,X)$. The goal of an attacker is to recover the plaintext $X$ from
ciphertext $Y$ with knowledge of function $C$, but partial or no knowledge of
key $K$.

For each of the four attack types, we train a Dense Prediction Transformer
(DPT) [28] $D$ to learn the inverse mapping from $B=\\{Y,\hat{X},\hat{K}\\}$
for arbitrary key $K$ to the corresponding plaintext $X$, where $Y$ denotes a
single cipherext, $\hat{X}$ an estimate of the plaintext, and $\hat{K}$ an
estimate of the key. For the COA, no estimates $\hat{X}$ and $\hat{K}$ are
available, so $B$ reduces to $Y$. For the KPAs, the methods for estimating
$\hat{X}$ and $\hat{K}$ depend on the design of the targeted optical
encryption camera (See section V for additional details.)

To train $D$, we generate a random key for each ciphertext sample in each
batch and use a loss consisting of a combination of an L1 pixel loss and an L2
perceptual loss (as in [29, 30]) over the outputs of layers _relu1_1_ ,
_relu2_2_ , and _relu3_3_ of VGG16 [31] pre-trained for image classification
on the ImageNet [32] dataset. Concretely, the loss is given by

$\mathcal{L}_{D}=w_{1}||D(B)-X||_{1}+w_{2}\sum_{i=1}^{3}||\phi_{i}(D(B))-\phi_{i}(X)||_{2}^{2},$
(2)

where $w_{1}=0.5$, $w_{2}=1.2$, and $\phi_{1}:\mathbb{R}^{H\times W\times
3}\rightarrow\mathbb{R}^{\frac{H}{2}\times\frac{W}{2}\times 64}$,
$\phi_{2}:\mathbb{R}^{H\times W\times
3}\rightarrow\mathbb{R}^{\frac{H}{4}\times\frac{W}{4}\times 128}$, and
$\phi_{3}:\mathbb{R}^{H\times W\times
3}\rightarrow\mathbb{R}^{\frac{H}{8}\times\frac{W}{8}\times 256}$ denote the
layers _relu1_1_ , _relu2_2_ , and _relu3_3_ , respectively, of the pre-
trained VGG16 network.

Figure 2: OpEnCam Framework. Our Optical Encryption Cameras encrypt image data
before image capture using optical masks. Recovery of decrypted image data is
possible only if you have access to the camera’s encryption key, which is
defined by the unique optical elements of each camera. The optical mask
patterns are generated by a novel mask generator that is robust against
various ciphertext-only attacks.

##  IV OpEnCam

###  IV-A Imaging Model (Encryption Process)

We propose a novel double-mask sensor design with two co-axially placed masks
- an optical scaling mask $S\in\mathbb{R}^{W_{3}\times H_{3}\times 3}$
positioned flush against the sensor array and an optical multiplexing mask
positioned a few millimeters above the scaling mask. The scaling mask is
implemented via an amplitude mask. For the multiplexing mask, an amplitude or
phase mask can be used, as described in the previous section. Let
$P\in\mathbb{R}^{W_{2}\times H_{2}\times 3}$ denote the point-spread-function
(PSF) produced by the multiplexing mask. Then, $P$ and $S$ constitute the
“encryption key” of the system, and the encrypted measurements
$Y\in\mathbb{R}^{W_{3}\times H_{3}\times 3}$ captured by the sensor are given
by

$Y=S\cdot(P\ast X)+N,$ (3)

where $X\in\mathbb{R}^{W_{1}\times H_{1}\times 3}$ denotes the scene and
$N\in\mathbb{R}^{W_{3}\times H_{3}\times 3}$ the sensor noise. It should be
noted again that $\ast$ represents a full-sized convolution with zero-padded
boundary conditions similar to the model used in existing thin lensless
cameras.

Compared to single-mask lensless designs [20, 2, 10, 11, 1], our novel double-
mask design has a larger key space, which makes decryption attacks more
challenging, and also inhibits impulse and uniform known plain text attacks.
For the single mask design, the encryption key (i.e., PSF $P$) is equal to the
impulse response of the camera, so the IKPA completely compromises the system.
This is not the case for our double-mask design, as the scaling mask makes the
system shift-variant, so the impulse response only reveals the “effective“
encryption key for a single scene pixel.

Imaging a uniform scene. We refer to the response of the camera to a uniform
scene as Uniform Scene Response (USR). Because of the structure of OpEnCam
forward imaging model, it is natural to wonder if a USR, due to scenes like a
plain wall, can reveal the scaling mask $S$ up to a scale. For a uniform scene
to reveal the scaling mask $S$, $P*X$ must also be uniform. However, due to
full-sized, linear convolution with zero-padded boundary conditions, $P*X$ is
never uniform, even for a plain wall scene. Therefore, it is not possible to
reveal $S$ from a USR.

###  IV-B Mask Design (Key Generation)

For both single and double-mask optical encryption cameras, the design of the
masks is critical for good performance, i.e., for enabling high-quality keyed
decryption while also thwarting decryption attacks. Consider the design with a
single multiplexing mask. The mask should have the following desirable
properties:

  1. 1.

PSF should contain directional filters for all angles to enable high-quality
lensless reconstruction [2].

  2. 2.

The autocorrelation of the PSF should not be an impulse. The PSFs for which
the autocorrelation is an impulse, are susceptible to autocorrelation-based
COAs [20]. These attacks exploit the fact that taking the autocorrelation of
encoded measurements eliminates the mask component if the autocorrelation of
the PSF is an impulse. This reduces the COA problem to recovering an image
from its autocorrelation.

  3. 3.

The PSF produced by the multiplexing mask should not be binary. Binary PSFs
may be revealed to an attacker when a point source appears in the scene.
Although our double-mask design inherently provides protection against such
attacks, it is still not desirable for the PSF to be compromised, as it
reduces the key size significantly. Consider a 1D case of binary PSF $p$, and
positive scaling mask $s$. The measurement $y(n)$ recorded at the sensor due
to a single point source is given by,

$y(n)=\begin{cases}s(n),&\text{if }n\in\mathcal{J}\\\
0,&\text{otherwise},\end{cases}$ (4)

where $\mathcal{J}$ is the set of all pixel locations for which $p(n)$ is 1.
Simple thresholding of $y$ i.e. $y>0$, reveals $p$ and reduces the size of the
key.

Keeping the above criteria in mind, we propose a novel design for the
multiplexing mask. The PSF due to OpEnCam multiplexing mask is given by,

$P=\alpha P_{colr}+(1-\alpha)P_{cont},$ (5)

where $P_{cont}$ is a binary contour PSF obtained from Perlin noise[33]
proposed in [2], $P_{colr}$ is a colored noise with certain roll-off. Perlin
contours are known for high-quality lensless imaging as shown in [2, 23] but
have impulse-like autocorrelations. On the hand, $P_{colr}$ is smoother and
has non-impulse autocorrelation with varying roll-off. We fix the Perlin
feature size, randomize the permutation vector for Perlin noise, and
correspondingly the contour PSF $P_{cont}$. The length of the permutation
vector is the same as the height or width of the PSF. To generate $P_{colr}$,
we generate colored noise with Power Spectral Density (PSD) given by,

$H_{\beta}(u,v)=\frac{1}{(u^{2}+v^{2})^{(\frac{\beta}{2})}}.$ (6)

The corresponding colored noise $P_{colr}(\beta)$ is given by,

$P_{colr}(\beta)=\mathcal{F}^{-1}(\sqrt{H_{\beta}}e^{j\theta}).$ (7)

Here $\mathcal{F}^{-1}$ is the inverse FFT, $\theta$ is random white noise
sampled from $\mathcal{U}(0,1)$ and $\beta$ is a scalar hyper-parameter that
controls the roll-off of the colored noise. We randomly sample $\beta$ from
$\mathcal{U}(1,10)$. Finally, the linear combination co-efficient $\alpha$ is
chosen uniformly at random. Apart from the above criteria, it should also be
noted that PSF must be non-negative as we can not subtract light.

For the scaling mask $S$, we use colored noise again without the Perlin
contours. For a given pair of $(S,P)$, the color of the noise ($\beta$) is the
same for both $S$ and $P_{colr}$ to avoid attacks due to filtering of the
either component. It should be noted that $S$ has to be positive as we do not
want to throw information from any of the sensor pixels. Some samples of the
generated mask patterns are shown in fig. 3.

Figure 3: Generated Mask Samples. Each column shows a scaling mask and
multiplexing mask pair generated from our key generator. Each of these is
randomly generated. Moreover, the underlying Perlin contours in the
multiplexing mask PSF are also randomized. Figure 4: Autocorrelation Analysis.
It is not possible to reconstruct the underlying plaintext/scene from the
measurement autocorrelation for the OpEnCam design $P_{O}$. While the
reconstruction is possible for $P_{W}$ (Perlin Contours + White noise).

###  IV-C Effect of Colored Noise PSF - Autocorrelation Analysis

In the above subsection, we claimed that autocorrelation of the multiplexing
mask should not be an impulse and, as a result, settled on a combination of
Perlin contour and colored noise for the multiplexing PSF of OpEnCam. In this
experiment, we validate this claim by comparing the performance of two
different multiplexing mask designs under autocorrelation-based ciphertext-
only attack. The mask designs being compared are (a) $P_{O}=\alpha
P_{colr}+(1-\alpha)P_{cont}$, (b)$P_{W}=\alpha P_{white}+(1-\alpha)P_{cont}$,
where $P_{white}$ is white noise PSF with impulse autocorrelation. $P_{O}$ is
the OpEnCam design while $P_{W}$ is a corresponding version with white noise
used instead of colored noise for PSF. To perform this experiment, we first
simulate measurements using each of the mask generators and perform
autocorrelation of the simulated measurements. We then feed this
autocorrelation to a U-Net, which then learns to estimate the underlying scene
or plaintext from the autocorrelation. Given, that the autocorrelation of
$P_{W}$ is close to an impulse (because autocorrelations of $P_{white}$ and
$P_{cont}$ are impulse), U-Net finds it easy to learn the mapping from the
measurement autocorrelation to the underlying scene for $P_{W}$, while it
fails to do so for $P_{O}$. We perform this experiment on MNIST dataset. Fig.
4 shows the visual results for the experiment. As can be seen, it is possible
to reconstruct the underlying scene from the autocorrelation of the
measurement for $P_{W}$, while this is not the case for $P_{O}$ i.e. our
OpEnCam design.

Figure 5: Keyed Decryption. The first five columns show the keyed-decryption
performance for existing single mask systems. Keyed decryption using PhlatCam
(Boominathan et al.[2]) is currently state of the art for lensless imaging.
OpEnCam keyed decryption performance matches that of existing mask designs for
lensless imaging and optical encryption despite the addition of a scaling mask
to enhance privacy.

Encryption Design PSNR $\uparrow$ SSIM $\uparrow$ LPIPS $\downarrow$ PhlatCam
[2] 23.10 0.81 0.30 Multi-pinhole [10] 22.70 0.79 0.30 FlatCam [9] 20.68 0.74
0.40 Random Binary [11] 20.22 0.77 0.35 Random Speckle [7] 20.43 0.71 0.44
OpEnCam 22.90 0.80 0.32

TABLE I: Image Recovery via Keyed Decryption. The reconstruction quality of
OpEnCam is at par with state of the art lensless image reconstruction quality.
Best performance is made bold while the second best performance is underlined.

###  IV-D Image Reconstruction (Keyed Decryption)

The decryption module accepts three inputs, an encrypted image $Y$, a scaling
mask $S$, and a PSF $P$ due to the multiplexing mask $M$; and returns a
decrypted image $\hat{X}$. For an ideal system, decryption can be achieved by
simply applying the inverse of the two optical computations. However, in
practice, to account for noise, cropping, and imperfect calibration, we
instead employ a double-step estimation process to recover the scene $X$. In
the first step, we employ Tikhonov regularized least squares to recover an
initial estimate $\hat{X_{T}}$ of image $X$:

$\hat{X}_{T}=\operatorname*{arg\,min}_{X}||Y_{N}-P\ast
X||^{2}_{F}+\gamma||X||^{2}_{F}$ (8)

Here $Y_{N}$ is the measurement after scaling normalization i.e.
$Y_{N}=\frac{Y}{S+\epsilon}$. $\hat{X}_{T}$ has a closed-form solution and can
be implemented in the Fourier domain as Wiener filtering.

In the second step, we feed our initial estimate $\hat{X_{T}}$ plus
measurement $Y$, and multiplexing mask PSF $P$ to a Dense Prediction
Transformer (DPT) [28] $D:\mathbb{R}^{W_{1}\times H_{1}\times
7}\rightarrow\mathbb{R}^{W_{1}\times H_{1}\times 3}$, which we train via
stochastic gradient descent to produce a refined version $\hat{X}$, given
input $B=\\{Y,\hat{X_{T}},P\\}$. As in eq. 2, we use a combination of L1 loss
and VGG feature loss to learn the mapping with $w_{1}=0.5$ and $w_{2}=1.2$. We
found that this transformer-based refinement performed slightly better than
the U-Net-based refinement typically used in lensless imaging works like [34,
23, 35].

##  V Experimental Results

###  V-A Preliminaries

####  V-A1 Baselines

We compare the performance of OpEnCam design with other mask-based passive
optical encryption designs as well as existing lensless camera designs. Our
baselines for optical encryption include multi-pinhole mask generator (Ishii
et al.[10]), random binary mask generator (Wang et al.[11]) and random speckle
mask generator (Zang et al.[7]) We also compare against PhlatCam[2] lensless
camera mask generator, which uses a Perlin contour PSF multiplexing mask, and
FlatCam[1, 9], which uses separable Maximum Length Sequences (MLS) for mask
patterns.

####  V-A2 Dataset and Implementation

We validate the efficacy of our proposed framework via simulation and real-
world experiments. For training keyed-decryption and attack models described
in this section, we use the Places365 dataset [36]. We use a subset of
ImageNet[37] for quantitative evaluation of the keyed-decryption and attack
models. Given that each encryption approach (including OpEnCam) has its own
mask generator that randomly generates the corresponding mask patterns, we
simulate the captured measurement/ciphertext using these randomly generated
mask patterns using the above datasets. Once simulated, we feed these
measurements (along with the mask patterns for keyed-decryption) to the
corresponding neural networks for final predictions. During testing, we
generate predictions for a large number of mask patterns generated from each
mask generator and report the average performance.

Encryption Design PSNR $\downarrow$ SSIM $\downarrow$ LPIPS $\uparrow$
PhlatCam [2] 18.10 0.63 0.57 Multi-pinhole [10] 19.29 0.65 0.48 FlatCam [9]
19.00 0.65 0.54 Random Binary [11] 18.73 0.64 0.56 Random Speckle [7] 19.42
0.65 0.54 OpEnCam\- Single Mask 18.22 0.64 0.57 OpEnCam\- Double Mask 16.28
0.59 0.60

TABLE II: Transformer-based Blind COA. Double mask OpEnCam provides the best
security as seen from the reconstruction quality. Large gap between single and
double mask OpEnCam highlights the improvement offered by the scaling mask in
terms of privacy. Best performance is made bold while the second best
performance is underlined.

###  V-B Keyed Decryption

Given that high-quality scene recovery using the key or mask patterns defines
the usefulness of our system, we first experimentally validate the performance
of keyed decryption transformer. We compare the keyed-decryption performance
of OpEnCam against PhlatCam [2], Multi-pinhole [10], FlatCam [9], Random
Binary [11] and Random Speckle [7]. We use the same keyed-decryption strategy
as described in Section IV-D for a fair comparison. The keyed-decryption
network uses Tikhonov reconstructions described in Equation 8 along with the
measurement and the multiplexing mask as input and produces the final
reconstructions. We use DPT [28] as the transformer. We compute the average
Peak Signal to Noise ratio (PSNR), Structural Similarity Index Measure (SSIM)
and Learned Perceptual Image Patch Similarity (LPIPS) for a subset of ImageNet
dataset as described in [23]. Table I shows the average performance. Higher
PSNR, SSIM and lower LPIPS indicate better reconstruction quality. The
performance of OpEnCam is at par with state of the art lensless image
reconstruction performance [2, 34] while marginally outperforming existing
encryption methods like [10, 11, 7, 9]. A similar trend is observed in visual
results shown in fig. 5.

###  V-C COA: Ciphertext Only Attack

Figure 6: Blind COA.We show some visual reconstruction results for two scenes.
Double-mask OpEnCam reconstructions are barely recognizable, while existing
approaches reveal a significant amount of texture.

We train a Dense Prediction Transformer(DPT) [28] to learn a mapping from a
single ciphertext/measurement to the corresponding plaintext/scene for each of
the mask generators using the loss function described in 2. A secure optical
encryption system should be robust against these powerful data-driven
transformer-based attacks. We compare the scene reconstruction performance for
different optical encryption strategies. Among the existing works, we compare
against PhlatCam[2], random speckle[7], multi-pinhole[10], random binary[11]
and FlatCam[9]. We also compare against single-mask OpEnCam i.e. the OpEnCam
design without the scaling mask to highlight the importance of the double-mask
system. Lower PSNR, SSIM and higher LPIPS indicate lower reconstruction
quality and better security. In fig. 8, we show the privacy-utility tradeoff
plot. From the plot, it can be observed that OpEnCam clearly outperforms
existing optical encryption and lensless imaging designs. The reconstructions
from OpEnCam, shown in fig. 6, are barely recognizable, indicating that it is
extremely challenging to recover the original scene from a single encoded
measurement captured using OpEnCam when the key or mask patterns are not
available. Table II shows the quantitative comparison of the different
encryption approaches. Double Mask OpEnCam clearly outperforms all other
existing optical encryption and lensless imaging designs. The large gap
between OpEnCam single and double mask models suggests that the scaling mask
not only makes the system shift-variant but also adds another layer of privacy
to existing mask-based lensless camera designs.

Encryption Design PSNR $\downarrow$ SSIM $\downarrow$ LPIPS $\uparrow$
PhlatCam [2] 19.04 0.76 0.38 Multi-pinhole [10] 22.61 0.85 0.22 FlatCam [9]
18.79 0.66 0.52 Random Binary [11] 19.28 0.72 0.43 Random Speckle [7] 19.04
0.65 0.53 PhlatCam - Double Mask 18.28 0.67 0.50 OpEnCam\- Single Mask 18.02
0.66 0.57 OpEnCam\- Double Mask 16.28 0.59 0.60

TABLE III: Impulse - KPA. Attacker has access to an approximate PSF $P$
corresponding to a single bright source. Double-mask OpEnCam outperforms all
existing designs. PhlatCam even with a scaling mask (PhlatCam-Double) is
susceptible to this attack due to its sparse binary PSF. Best performance is
made bold while the second best performance is underlined.

###  V-D I-KPA: Impulse - Known Plaintext Attack

An optical encryption camera with a single multiplexing mask is susceptible to
an Impulse based Known Plaintext Attack (I-KPA). More specifically, if a
bright source appears in a scene, an approximate point spread function (PSF)
due to the multiplexing mask is revealed. Since for single mask-based
encryption, the PSF is the key, the attacker can use a bright-source
measurement to decrypt ciphertexts. Moreover, as pointed out in Section IV-B,
binary masks are more susceptible to these bright source attacks even with a
scaling mask. In this experiment, we show that our double-mask encryption of
OpEnCam is robust against these bright source attacks. To perform this
experiment, we first simulate a measurement of a scene with a bright source.
The relative intensity of the bright source with respect to the rest of the
scene is $10^{3}$. Assuming this measurement as the approximate PSF, we apply
our keyed decryption pipeline of Section IV-D. We compare against the
baselines described above: PhlatCam [2], multi-pinhole [10], random binary
[11], FlatCam MLS [9]. To highlight the importance of our double-mask system
we also compare it against single-mask version of OpEnCam without the scaling
mask. Furthermore, to highlight the disadvantage of binary masks like that of
PhlatCam under point source attacks, we compare against a modified version of
PhlatCam where we add an additional scaling mask. We show the privacy-utility
tradeoff plot for the I-KPA in fig. 8. Visual reconstruction results are shown
in fig. 7. It can be seen that OpEnCam with scaling mask significantly
outperforms all existing mask approaches. Moreover, despite the use of scaling
mask, PhlatCam-Double Mask is unable to provide any privacy due to the binary
nature of its multiplexing mask. Corresponding visual results for a few
important approaches are also compared in the same figure. We show the
quantitative results for the I-KPA in table III.

Figure 7: Impulse - KPA. We show some visual reconstruction results for two
scenes. Double-mask OpEnCam reconstructions are barely recognizable, while
other designs are completely susceptible to this attack, including Double-
masked PhlatCam. Figure 8: Privacy-Utility Plots. We show the privacy-utility
plot for blind COA and Impulse-KPA. For both attacks, double-masked OpEnCam,
outperforms other designs and occupies the desired top-left corner.

###  V-E U-KPA: Uniform - Known Plaintext Attack

####  V-E1 Attack due to uniform scene

To verify that a uniform scene response (USR) doesn’t reveal the scaling mask,
we approximated the scaling mask $S$ as $S_{USR}=Y_{USR}$, where $Y_{USR}$ is
a measurement for an all-ones scene. We used $S_{USR}$ to remove the scaling
mask component from the test measurement and decrypted it using the blind COA
transformer trained for single-mask OpEnCam. We show the decrypted result and
the error $|S_{USR}-S|$ in fig. 9(a). It can be seen that $S_{USR}$ is not an
accurate approximation of $S$.

Figure 9: Uniform - KPA. (a) From a USR, (b) from ciphertexts. It can be seen
that both the USR and the average of measurements fail to reveal the
underlying hidden scaling mask leading to poor plain text estimations.

####  V-E2 Attack due to multiple measurements

Another question that may arise is whether an average of diverse measurements
reveal the scaling mask. To verify this, we assume access to $N=12000$ diverse
OpEnCam measurements of scenes from the Places365 dataset and estimate the
scaling mask $S$ as $S_{avg}=(1/N)\sum_{j=1}^{N}Y^{(j)}$. We then used
$S_{avg}$ to remove the scaling mask component from the measurement and
decrypted it using the transformer trained for single-mask OpEnCam. We show
the results in fig. 9(b). It can be seen that $S_{avg}$ is not an accurate
approximation of $S$. This attack assumes that under the condition
$N\rightarrow\infty$, $S_{avg}$ will be approximately a USR.

###  V-F UI-KPA: Uniform-Impulse Known Plaintext Attack

In this attack, the attacker tries to jointly estimate the scaling mask $S$
and PSF $P$ and then uses them to decrypt the ciphertext. We assume the
attacker has access to the approximate impulse response for a scene with a
dominating bright source, and additional information in the form of a white-
wall-measurement. We then solve the optimization,

$\\{\hat{S},\hat{P}\\}=\operatorname*{arg\,min}_{S,P}\sum_{i=1}^{2}||Y_{i}-S\odot(P*X_{i})||_{F}^{2},$
(9)

to estimate the keys. $Y_{1}$ is a USR, and $Y_{2}$ is the captured impulse
response, while $X_{1}$ and $X_{2}$ are all-ones images and impulse images
respectively. $Y_{1}$ can also be an average measurement corresponding to a
large diverse set of natural scenes. Once $S$ and $P$ are estimated, the
attacker uses the keyed-decryption framework of Section IV-D to decrypt the
ciphertexts/measurements. We show the results for different relative
intensities of the point source with respect to the rest of the scene in fig.
10. For cases where the relative intensity of the bright source is less than
$10^{4}$, it is not possible to accurately estimate the scaling mask and the
PSF, and as a result, the quality of the decrypted scene is poor. For relative
intensities beyond $10^{4}$, it is possible to reconstruct the scene to a
reasonable extent. However, such large relative intensities of point sources
imply dark-room-like environments which are less likely to occur in the attack
model we have assumed i.e. the attacker doesn’t have access to the physical
location of the camera. Nevertheless, we acknowledge that this could be a
limitation of OpEnCam and leave improving security for such scenarios as
future work.

Figure 10:  Uniform-Impulse - KPA. (Top Rows) We show the decrypted scenes for
varying relative intensities of the bright source. (Bottom Row) We show the
scene with a point source (encircled) of varying relative intensity. OpEnCam
offers privacy for a wide range of practical relative point source
intensities.

###  V-G Real Experimental Results

For the real-world experiments, we construct a working prototype OpEnCam using
Basler Ace4024-29uc machine vision camera. To implement the multiplexing mask,
we use an amplitude mask printed using a conventional office printer and place
it over a Perlin contour phase mask. The scaling mask $S$ essentially
implements a spatially varying exposure pattern for each sensor pixel in
OpEnCam. Given that placing the scaling mask flush against the bare sensor
requires breaking the sensor’s protective covering, we implement its effect
through exposure bracketing. That is, we quantize the scaling mask into 16
levels and then capture the measurement due to the multiplexing mask at 16
different relative exposure values given by the 16 levels of quantization.
Finally, using the quantized scaling mask, we blend the measurement. We
capture the experimental scaling mask ($S_{exp}$) by displaying a white image
on a monitor and blending the measurements captured at 16 different exposure
values without the multiplexing mask. We then capture the experimental PSF
($P_{exp}$) in the same way, but with the multiplexing mask, i.e., we capture
the measurements for a point source at 16 different exposure values and blend
them. The true PSF is then the normalized version of $P_{exp}$ i.e.,
$P_{true}=P_{exp}/(S_{exp}+\epsilon)$. Finally, we display images of natural
scenes on the monitor, capture the corresponding exposure stack, and blend it
to obtain the experimental measurements. For keyed decryption on these
experimental data, we use the Tikhonov regularized reconstruction (Eq. (8)) to
obtain the intermediate scene estimates. Since these estimates are generally
noisy, we then use a denoising U-Net to refine them. For blind transformer-
based attack, we use the DPT trained in Section V-C on the experimental
measurements. We show the visual results in fig. 11. These results demonstrate
the efficacy of our approach – recovery of the original scene from an OpEnCam
measurement is possible only if one has access to the camera’s encryption key.

Figure 11: Real Results from OpEnCam Prototype. The top row shows OpEnCam
measurements. The second row shows the reconstructions from a DPT-based blind
COA. The third row shows keyed reconstructions. The bottom row shows the
ground-truth scene.

##  VI Conclusion and Discussion

We propose a novel design for lensless cameras called OpEnCam to perform
optical encryption. At the core of our method is the double-mask imaging model
which implements shift-variant forward encryption process. Moreover, using
ideas derived from signal processing, we design the mask patterns of our
system to not only allow high-quality lensless imaging but also be robust
against various forms of attacks. Through extensive experiments, we showed
that our proposed OpEnCam design is robust against powerful learning-based
ciphertext-only attacks and three special cases of known plaintext attacks.
Finally, we build an OpEnCam prototype that allows us to validate our claims
through promising preliminary results. Using a photolithography-based phase
mask printing process[2], our prototype’s performance can significantly
improve, and we plan to implement this in the future.

Although preliminary results from OpEnCam show promising results, it should be
noted that, like existing optical encryption systems, OpEnCam also has its
limitations. Accessing the hardware or allowing extensive control of the scene
has the potential to reveal the key. However, as pointed out through our
experiments, the OpEnCam design is still an improvement over the existing
optical encryption system and is a step towards more secure imaging.

In future, it would be interesting to look into the data-driven design of
OpEnCam systems and using multiple-masked lensless cameras for other
applications.

## References

  * [1] M. S. Asif, A. Ayremlou, A. Sankaranarayanan, A. Veeraraghavan, and R. G. Baraniuk, “Flatcam: Thin, lensless cameras using coded aperture and computation,” _IEEE Transactions on Computational Imaging_ , vol. 3, no. 3, pp. 384–397, 2016.
  * [2] V. Boominathan, J. K. Adams, J. T. Robinson, and A. Veeraraghavan, “Phlatcam: Designed phase-mask based thin lensless camera,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 42, no. 7, pp. 1618–1629, 2020.
  * [3] N. Antipa, G. Kuo, R. Heckel, B. Mildenhall, E. Bostan, R. Ng, and L. Waller, “Diffusercam: lensless single-exposure 3d imaging,” _Optica_ , vol. 5, no. 1, pp. 1–9, 2018.
  * [4] N. Antipa, P. Oare, E. Bostan, R. Ng, and L. Waller, “Video from stills: Lensless imaging with rolling shutter,” in _2019 IEEE International Conference on Computational Photography (ICCP)_. IEEE, 2019, pp. 1–8.
  * [5] K. Monakhova, K. Yanny, N. Aggarwal, and L. Waller, “Spectral diffusercam: lensless snapshot hyperspectral imaging with a spectral filter array,” _Optica_ , vol. 7, no. 10, pp. 1298–1307, 2020.
  * [6] P. Refregier and B. Javidi, “Optical image encryption based on input plane and fourier plane random encoding,” _Optics letters_ , vol. 20, no. 7, pp. 767–769, 1995.
  * [7] J. Zang, Z. Xie, and Y. Zhang, “Optical image encryption with spatially incoherent illumination,” _Opt. Lett._ , vol. 38, no. 8, pp. 1289–1291, Apr 2013. [Online]. Available: https://opg.optica.org/ol/abstract.cfm?URI=ol-38-8-1289
  * [8] J. Byrne, B. Decann, and S. Bloom, “Key-nets: Optical transformation convolutional networks for privacy preserving vision sensors,” in _British Machine Vision Conference (BMVC)_ , 2020.
  * [9] T. Nguyen Canh and H. Nagahara, “Deep compressive sensing for visual privacy protection in flatcam imaging,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision Workshops_ , 2019, pp. 0–0.
  * [10] Y. Ishii, S. Sato, and T. Yamashita, “Privacy-aware face recognition with lensless multi-pinhole camera,” in _European Conference on Computer Vision_. Springer, 2020, pp. 476–493.
  * [11] Z. W. Wang, V. Vineet, F. Pittaluga, S. N. Sinha, O. Cossairt, and S. Bing Kang, “Privacy-preserving action recognition using coded aperture videos,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops_ , 2019.
  * [12] X. Peng, H. Wei, and P. Zhang, “Chosen-plaintext attack on lensless double-random phase encoding in the fresnel domain,” _Optics letters_ , vol. 31, no. 22, pp. 3261–3263, 2006.
  * [13] M. Liao, D. Lu, W. He, and X. Peng, “Optical cryptanalysis method using wavefront shaping,” _IEEE Photonics Journal_ , vol. 9, no. 1, pp. 1–13, 2017.
  * [14] X. Peng, P. Zhang, H. Wei, and B. Yu, “Known-plaintext attack on optical encryption based on double random phase keys,” _optics letters_ , vol. 31, no. 8, pp. 1044–1046, 2006.
  * [15] U. Gopinathan, D. S. Monaghan, T. J. Naughton, and J. T. Sheridan, “A known-plaintext heuristic attack on the fourier plane encryption algorithm,” _Optics Express_ , vol. 14, no. 8, pp. 3181–3186, 2006.
  * [16] X. Peng, H.-Q. Tang, and J.-D. Tian, “Ciphertext-only attack on double random phase encoding optical encryption system,” 2007.
  * [17] C. Zhang, M. Liao, W. He, and X. Peng, “Ciphertext-only attack on a joint transform correlator encryption system,” _Optics express_ , vol. 21, no. 23, pp. 28 523–28 530, 2013.
  * [18] X. Liu, J. Wu, W. He, M. Liao, C. Zhang, and X. Peng, “Vulnerability to ciphertext-only attack of optical encryption scheme based on double random phase encoding,” _Optics express_ , vol. 23, no. 15, pp. 18 955–18 968, 2015.
  * [19] M. Liao, S. Zheng, S. Pan, D. Lu, W. He, G. Situ, X. Peng _et al._ , “Deep-learning-based ciphertext-only attack on optical double random phase encryption,” _Opto-Electronic Advances_ , vol. 4, no. 5, p. 05200016, 2021.
  * [20] M. Liao, W. He, D. Lu, and X. Peng, “Ciphertext-only attack on optical cryptosystem with spatially incoherent illumination: from the view of imaging through scattering medium,” _Scientific Reports_ , vol. 7, no. 1, pp. 1–9, 2017.
  * [21] Y. Wang, J. Gawedzinski, M. E. Pawlowski, and T. S. Tkaczyk, “3d printed fiber optic faceplates by custom controlled fused deposition modeling,” _Optics Express_ , vol. 26, no. 12, pp. 15 362–15 376, 2018.
  * [22] M. F. Duarte, M. A. Davenport, D. Takhar, J. N. Laska, T. Sun, K. F. Kelly, and R. G. Baraniuk, “Single-pixel imaging via compressive sampling,” _IEEE signal processing magazine_ , vol. 25, no. 2, pp. 83–91, 2008.
  * [23] S. S. Khan, V. Sundar, V. Boominathan, A. Veeraraghavan, and K. Mitra, “Flatnet: Towards photorealistic scene reconstruction from lensless measurements,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2020.
  * [24] D. Bagadthey, S. Prabhu, S. S. Khan, D. T. Fredrick, V. Boominathan, A. Veeraraghavan, and K. Mitra, “Flatnet3d: intensity and absolute depth from single-shot lensless capture,” _JOSA A_ , vol. 39, no. 10, pp. 1903–1912, 2022.
  * [25] X. Pan, X. Chen, S. Takeyama, and M. Yamaguchi, “Image reconstruction with transformer for mask-based lensless imaging,” _Optics Letters_ , vol. 47, no. 7, pp. 1843–1846, 2022.
  * [26] J. D. Rego, K. Kulkarni, and S. Jayasuriya, “Robust lensless image reconstruction via psf estimation,” in _Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision_ , 2021, pp. 403–412.
  * [27] G. Wetzstein, A. Ozcan, S. Gigan, S. Fan, D. Englund, M. Soljačić, C. Denz, D. A. Miller, and D. Psaltis, “Inference in artificial intelligence with deep optics and photonics,” _Nature_ , vol. 588, no. 7836, pp. 39–47, 2020.
  * [28] R. Ranftl, A. Bochkovskiy, and V. Koltun, “Vision transformers for dense prediction,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2021, pp. 12 179–12 188.
  * [29] C. Ledig, L. Theis, F. Huszár, J. Caballero, A. Cunningham, A. Acosta, A. Aitken, A. Tejani, J. Totz, Z. Wang _et al._ , “Photo-realistic single image super-resolution using a generative adversarial network,” in _CVPR_ , 2017, pp. 4681–4690.
  * [30] A. Dosovitskiy and T. Brox, “Generating images with perceptual similarity metrics based on deep networks,” in _Advances in Neural Information Processing Systems_ , 2016, pp. 658–666.
  * [31] K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” in _ICLR_ , 2015.
  * [32] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei, “ImageNet: A large-scale hierarchical image database,” in _CVPR_ , 2009, pp. 248–255.
  * [33] K. Perlin, “Improving noise,” in _Proceedings of the 29th annual conference on Computer graphics and interactive techniques_ , 2002, pp. 681–682.
  * [34] S. S. Khan, V. Adarsh, V. Boominathan, J. Tan, A. Veeraraghavan, and K. Mitra, “Towards photorealistic reconstruction of highly multiplexed lensless images,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2019, pp. 7860–7869.
  * [35] K. Monakhova, J. Yurtsever, G. Kuo, N. Antipa, K. Yanny, and L. Waller, “Learned reconstructions for practical mask-based lensless imaging,” _Optics express_ , vol. 27, no. 20, pp. 28 075–28 090, 2019.
  * [36] B. Zhou, A. Lapedriza, A. Khosla, A. Oliva, and A. Torralba, “Places: A 10 million image database for scene recognition,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2017.
  * [37] O. Russakovsky, J. Deng, H. Su, J. Krause, S. Satheesh, S. Ma, Z. Huang, A. Karpathy, A. Khosla, M. Bernstein, A. C. Berg, and L. Fei-Fei, “ImageNet Large Scale Visual Recognition Challenge,” _International Journal of Computer Vision (IJCV)_ , vol. 115, no. 3, pp. 211–252, 2015.


Robust approximate linear regression without correspondence

Amin Nejatbakhsh Department of Neurobiology and Behavior Columbia University
Erdem Varol Department of Statistics Columbia University

###### Abstract

We propose methods for estimating correspondence between two point sets under
the presence of outliers in both the source and target sets. The proposed
algorithms expand upon the theory of the regression without correspondence
problem to estimate transformation coefficients using unordered multisets of
covariates and responses. Previous theoretical analysis of the problem has
been done in a setting where the responses are a complete permutation of the
regressed covariates. This paper expands the problem setting by analyzing the
cases where only a subset of the responses is a permutation of the regressed
covariates in addition to some covariates being outliers. We term this problem
robust regression without correspondence and provide several algorithms based
on random sample consensus for exact and approximate recovery in a noiseless
and noisy one-dimensional setting as well as an approximation algorithm for
multiple dimensions. The theoretical guarantees of the algorithms are verified
in simulated data. We demonstrate an important computational neuroscience
application of the proposed framework by demonstrating its effectiveness in a
Caenorhabditis elegans neuron matching problem where the presence of outliers
in both the source and target nematodes is a natural tendency.

##  1 Introduction

Point set registration is one of the central problems in computer vision that
involves the optimization of a transformation that aligns two sets of point
clouds [1, 2]. Point set registration have been applied in numerous fields
including but not limited to robotics [3], medical imaging [4], object
recognition [5], panorama stitching [6] and computational neuroscience [7].
The types of allowable transformations and energy functions utilized in the
cost function have differentiated varying methods [8, 9, 10, 11, 12, 13, 14,
2, 15, 16]. In general, point set registration methods employ an iterative
strategy of solving the transformation and updating the matching which works
well in practice but there are no guarantees for reaching the global optima
[17]. Only a few methods have provided approximate globally optimal solutions
[18, 10]. These methods rely on severe constraints of the transformation
domains, such as the 3D rotation group SO(3), in order to employ branch and
bound techniques on discretizations.

Theoretical analysis of the recovery guarantees of point set registration has
not been performed for a general number of dimensions until recently when it
was termed as unlabelled sensing by [19] as a problem with duality connections
with the well-known problem of compressed sensing [20]. In this problem,
similar to linear regression, the response signal is modeled as a linear
combination of a set of covariates. However, the correspondence of the
responses to the covariates is modeled as having been shuffled by an unknown
permutation matrix. For this reason, the problem has also been termed as
linear regression with shuffled labels [21], linear regression with an unknown
permutation [22], homomorphic sensing [23] or linear regression without
correspondence (RWOC) [24], the latter of which will be used to refer to the
problem herein. Although RWOC is, in general, an NP-hard problem [22], there
have been several advances in recent years to propose signal to noise ratio
(SNR) bounds for recovery of the permutation matrix and the regression
coefficients [22, 25]. Conversely, the same works have also analyzed the SNR
and sampling regime by which no recovery is possible.

Nevertheless, the computer vision community has attempted to solve the point
set registration problem through consideration of outliers and missing
correspondences, which are typically encountered in real-world applications. A
common technique used in point set registration to robustify the optimization
against outliers is to employ random sampling consensus (RANSAC) subroutines
[26, 27, 28]. The main advantages of RANSAC are that the randomization
procedure employed can severely reduce the computational cost of an otherwise
combinatorial search.

Motivated by applications in computational neuroscience such as matching the
neuronal populations of Caenorhabditis elegans (C. elegans) across different
nematodes, we aim to unify the ideas presented in RWOC literature and robust
point set registration methods to provide provably approximate solutions to
the RWOC problem in the presence of outliers and missing measurements commonly
encountered in fluorescence microscopy data. Robustly and automatically
matching and identifying neurons in C. elegans could expedite the post-
experimental data analysis and hypothesis testing cycle [7, 29, 30, 31].

###  1.1 Main contributions

The main contributions presented in this paper are the introduction of
randomized algorithms for the recovery of the regression coefficients in the
RWOC problem that takes into account noise, missing data, and outliers. Hsu et
al. [24] provide algorithms for the noisy case without generative assumptions;
their algorithm takes into account square permutation matrices, which assumes
that the entire signal is captured in the responses and does not take into
account any missing correspondences or outliers. Unnikrishnan et al. [19, 25]
provide combinatorial existence arguments. Tsakiris et al. [23] provide an
algorithm that takes into account missing correspondences or outliers but not
both. Our method is designed for the practical purpose of matching point
clouds that may have noisy measurements, missing correspondences, and
outliers. This is undoubtedly the case in the application domain of neuron
tracking and matching in biological applications. Specifically, we demonstrate
the efficacy of the proposed method in the identification and tracking of in-
vivo (C. elegans) neurons. In summary, our contributions are four-fold:

  1. 1.

We introduce the notion of "robust" regression without correspondence (rRWOC)
that models missing correspondences between responses and covariates as well
as completely missed associations in the form of outliers and missing data.

  2. 2.

We introduce a polynomial-time algorithm to find the exact solution for the
one-dimensional noiseless rRWOC and the approximate solution in the noisy
regime.

  3. 3.

We introduce a randomized approximately correct algorithm that is more
efficient than pure-brute force approaches in multiple dimensional rRWOC.

  4. 4.

We demonstrate the computational neuroscience application of our approach to
point-set registration problems in the context of automatically matching and
identification of the cellular layout of the nervous system of the nematode C.
elegans.

Figure 1: Geometric intuition of the proposed algorithms. Top: Solving the
linear regression problem on the correct d-tuple correspondence yields
regression coefficients that result in a transformed set with maximal
alignment. Bottom: Incorrect correspondences of the sampled $d$-tuples results
in sub-maximal alignment.

###  1.2 Paper organization

In section 2, we introduce our statistical regression model (rRWOC) that
accounts for permuted correspondences, outliers, and noise. We then
demonstrate the added computational complexity of recovery of rRWOC in
contrast with simple linear regression and RWOC in a one-dimensional case in
section 3.1. In section 3.2, we provide a randomized algorithm for the rRWOC
problem in multiple dimensions with convergence bounds. Lastly, in section
4.1, we verify the theoretical recovery guarantees in simulated experiments
and in section 4.2 show the neuroscience application of the proposed
algorithms in the C. elegans neuron matching problem.

##  2 Regression model

First, we introduce notation. Let
$\mathbf{X}=[\mathbf{x}_{1}|\mathbf{x}_{2}|\ldots|\mathbf{x}_{m}]^{T}\in\mathbb{R}^{m\times
d}$ and
$\mathbf{Y}=[\mathbf{y}_{1}|\mathbf{y}_{2}|\ldots|\mathbf{y}_{n}]^{T}\in\mathbb{R}^{n\times
d}$ denote two d-dimensional point sets consisting of $m$ and $n$ points,
respectively. Let us call $\mathbf{X}$ the reference or source set. Let
$\mathbf{Y}$ denote the target set which may contain outliers and missing
correspondences. Note that the points in $\mathbf{X}$ that are missing
correspondences in $\mathbf{Y}$ can be seen as outliers in the source set,
hence justifying our claim that we model outliers in both the source and
target sets.

Let the set of indices
$\mathcal{I}=\\{i_{1},\ldots,i_{|\mathcal{I}|}\\}\subseteq[n]$ denote the
indices of $\mathbf{y}_{j}$ which are inliers. Conversely, let
$\mathcal{O}=\\{o_{1},\ldots,o_{|\mathcal{O}|}\\}\subseteq[n]$ denote set of
indices of $\mathbf{y}_{j}$ which are outliers. By construction, these sets
are a disjoint partition of the entire index set of target points:
$\mathcal{I}\bigcup\mathcal{O}=[n]$ and
$\mathcal{I}\bigcap\mathcal{O}=\emptyset$. Let
$\boldsymbol{\Pi}\in\mathcal{P}^{n\times m}$ denote a possibly unbalanced
permutation matrix where there are at most $\min\\{n,m\\}$ ones placed such
that no row or column has more than a single one. All other entries are
zeroes. Let $\pi(i)$ denote the location of the one in the $i$th row of the
permutation matrix $\boldsymbol{\Pi}$. Next, let
$\boldsymbol{\beta}\in\mathbb{R}^{d\times d}$ denote the regression
coefficients and $\epsilon\sim\mathcal{N}(0,\nu\boldsymbol{I})$ denote zero-
mean Gaussian noise. Lastly, let $\text{U}[\mathcal{C}]$ denote the uniform
distribution within some closed convex set $\mathcal{C}$. Given these
definitions, we can define the robust regression without correspondence
(rRWOC) model as

$\displaystyle\mathbf{y}_{i_{j}}$
$\displaystyle=\mathbf{x}_{\pi(i_{j})}\boldsymbol{\beta}+\epsilon~{}\quad~{}$
$\displaystyle\text{for}~{}i_{j}\in\mathcal{I}$
$\displaystyle\mathbf{y}_{o_{l}}$
$\displaystyle\sim\text{U}[\mathcal{C}]~{}\quad~{}$
$\displaystyle\text{for}~{}o_{l}\in\mathcal{O}$ (1)

Note that the bias terms in the regression can be modeled by padding
$\mathbf{x}$ with a constant column of ones.

In contrast with linear regression, where the sole objective is to recover the
coefficients $\boldsymbol{\beta}$, the two-fold objective of RWOC is to
recover the correct permutation matrix $\boldsymbol{\Pi}$, and the regression
coefficients $\boldsymbol{\beta}$. To add to the complexity of the problem,
the three-fold objective of rRWOC is to recover the inlier set $\mathcal{I}$,
the permutation $\boldsymbol{\Pi}$, and the coefficients $\boldsymbol{\beta}$.

##  3 Algorithms

To aid in the recovery of the solution in rRWOC, we introduce the following
assumption.

######  Assumption 1 (Maximal inlier set).

For point sets $\mathbf{X}$, $\mathbf{Y}$, there exists a triple
$\\{\mathcal{I}^{*},\boldsymbol{\beta}^{*},\boldsymbol{\Pi}^{*}\\}$ that is
maximal in the sense that $n\geq|\mathcal{I}^{*}|\geq|\mathcal{I}^{\prime}|$
such that any other triple
$\\{\mathcal{I}^{\prime},\boldsymbol{\beta}^{\prime},\boldsymbol{\Pi}^{\prime}\\}$
is not considered to be the underlying regression model.

Assumption 1 allows the identifiability of whether a given hypothetical index
set can be considered to be the true underlying inlier set or not. In
practical terms, suppose we generate simulated data with $n$ points in
$\mathbf{Y}$ of which $k>n/2$ are outliers generated uniformly and the
remainder generated with respect to a coefficient
$\boldsymbol{\beta}^{\mathcal{I}}$ such that
$\mathbf{Y}_{[\mathcal{I}]}=\mathbf{X}_{\pi(\mathcal{I})}\boldsymbol{\beta}^{\mathcal{I}}+\epsilon^{\mathcal{I}}$.
There may be cases such that uniformly generated "outliers",
$\mathbf{Y}_{[\mathcal{O}]}$, are structured such that there exists a
coefficient $\boldsymbol{\beta}^{\mathcal{O}}$ and permutation
$\boldsymbol{\Pi}^{\mathcal{O}}$ such that
$\mathbf{Y}_{[\mathcal{O}]}=\mathbf{X}_{\pi(\mathcal{O})}\boldsymbol{\beta}^{\mathcal{O}}+\epsilon^{\mathcal{O}}$
where
$\text{Var}(\epsilon^{\mathcal{I}})\geq\text{Var}(\epsilon^{\mathcal{O}})$. In
this case, $\boldsymbol{\beta}^{\mathcal{O}}$ is identifiable but not
verifiable as "correct."

Equipped with the rRWOC model and assumption 1, we now demonstrate the
progressive increase in the complexity of recovery of ordinary linear
regression, RWOC, and rRWOC in one-dimension.

###  3.1 Optimal regression in $d=1$

Linear regression in one-dimension with known correspondences, no offset term
and no outliers can be obtained in $O(n)$ time using the univariate normal
equation:
$\beta_{OLS}=\frac{\sum_{i}^{n}y_{i}x_{\pi(i)}}{\sum_{i}^{n}x_{\pi{i}}^{2}}.$
On the other hand, RWOC in the one-dimensional case with no noise can be
solved in $O(n\log(n))$ steps via the method of moments and a simple sorting
operation. Namely, first, the regressor $\beta_{RWOC}$ can be estimated using
the ratio of the first moments of the covariates to the responses:

$\displaystyle\beta_{RWOC}=\frac{\sum_{i=1}^{n}y_{i}}{\sum_{i=1}^{n}x_{i}}$
(2)

and then the permutation can be recovered using the re-arrangement inequality
[32],

$\displaystyle\min_{\boldsymbol{\Pi}}\sum_{i=1}^{n}(y_{i}-\hat{y}_{\pi(i)})^{2}=\sum_{i=1}^{n}(y_{(i)}-\hat{y}_{(i)})^{2}=$
(3)
$\displaystyle\|\boldsymbol{\Pi}_{y}\mathbf{y}-\boldsymbol{\Pi}_{\hat{y}}\hat{\mathbf{y}}\|_{2}^{2}\longrightarrow\boldsymbol{\Pi}_{RWOC}=\boldsymbol{\Pi}_{y}^{T}\boldsymbol{\Pi}_{\hat{y}}$

where $y_{(i)}$ denotes sorted $y_{i}$ and $\hat{y}_{(i)}$ denotes sorted
$x_{i}\beta_{RWOC}$ and $\boldsymbol{\Pi}_{y}$ and
$\boldsymbol{\Pi}_{\hat{y}}$ denote the permutation matrices that capture the
sorting operations.

In the case with outlier elements in $\mathbf{y}$, the problem is non-trivial,
even in one dimension, since sorting does not allow the identification of
outliers. To solve the one dimensional rRWOC, we introduce algorithm 1 which
recovers the triplet
$\\{\mathcal{I}^{*},\boldsymbol{\beta}^{*},\boldsymbol{\Pi}^{*}\\}$ in an
exhaustive fashion.

Algorithm 1 One dimensional robust regression without correspondence -
Exhaustive approach

Input:  Reference set: $\\{x_{1},\ldots,x_{m}\\}$, target set:
$\\{y_{1},\ldots,y_{n}\\}$, outlier margin: $\nu$

Require:  $k<\frac{n}{2}$ (number of outliers)

for $i=1,\ldots,n$ do

for $j=1,\ldots,m$ do

Compute $\beta^{i,j}=y_{i}/x_{j}$

Compute linear assignment [33]:

$\boldsymbol{\Pi}^{i,j}\leftarrow\underset{\boldsymbol{\Pi}\in\mathcal{P}^{n\times
m}}{\arg\min}\|\mathbf{x}\beta^{i,j}-\boldsymbol{\Pi}^{T}\mathbf{y}\|_{2}^{2}$

Compute hypothetical inliers:

$\mathcal{I}^{i,j}=\\{l:|x_{\pi^{i,j}(l)}\beta^{i,j}-y_{l}|\leq\nu\\}$

end for

end for

return $(i^{*},j^{*})=\underset{(i,j)}{\arg\max}|\mathcal{I}^{i,j}|$ ,
$\mathcal{I}^{*}=\mathcal{I}^{i^{*},j^{*}}$,
$\boldsymbol{\Pi}^{*}=\boldsymbol{\Pi}^{i^{*},j^{*}}$,
$\beta^{*}\leftarrow\frac{\sum_{l\in\mathcal{I}^{*}}y_{l}x_{\pi^{*}(l)}}{\sum_{l\in\mathcal{I}^{*}}x_{\pi^{*}(l)^{2}}}$

######  Proposition 1 (Correctness of Algorithm 1).

Suppose there exist $n-k$ inliers in $\mathbf{y}$ and that $k<n/2$. Then
algorithm 1 yields the correct regression coefficient $\beta^{*}=\beta$ with
probability 1 for noiseless data and with high probability for noisy data with
an appropriately selected margin parameter $\nu$.

###### Proof.

(The full proof is included in supplementary material) The overview of the
proof is as follows. In the noiseless case, if $j=\pi(i)$ then
$\beta^{i,j}=\frac{y_{i}}{x_{j}}=\beta^{*}$. The projection
$\mathbf{x}\beta^{i,j}$ maps all reference points to their exact corresponding
reference points. Thus the Hungarian algorithm will yield these as the
assignments since they incur minimal cost. Therefore, we will have
$|\mathcal{I}^{i,j}|\geq n-k$. The cardinality of inliers is lower bounded and
not equal to $n-k$ since outlier points may by chance be transformed to points
in $\mathbf{y}$ as well. Contrarily, suppose the transformation $\beta^{i,l}$
for $l\neq\pi(i)$ yields a larger hypothesized inlier set $\mathcal{I}^{i,l}$,
such that $|\mathcal{I}^{i,l}|>|\mathcal{I}^{i,j}|$ then this means that there
are more points in $\mathbf{x}\beta^{i,l}$ that are closer to $\mathbf{y}$
than $\mathbf{x}\beta^{i,j}$, contradicting the assumption that $n-k$ is the
maximal inlier set. ∎

The time complexity of algorithm 1 can be analyzed as follows. The main
computational cost is due to linear assignment which incurs a cost of
$O(\max\\{m,n\\}^{3})$ if [34] variant is used. Linear assignment is repeated
$mn$ times. If $m$ and $n$ are of the same order, then algorithm 1 has
complexity $O(n^{5})$.

However, if the ratio of inliers to outliers is relatively high, then it is
possible to use randomization procedures like RANSAC [26, 27] to speed up the
algorithm to yield the correct regression coefficient with high probability.
This is demonstrated in algorithm 2.

######  Proposition 2 (Correctness of Algorithm 2).

Suppose there are $n-k$ inliers in $\mathbf{x}$ and that $k<n/2$. In
$q\geq\frac{\log(1-\delta)}{\log(1-\frac{n-k}{mn})}$ iterations, algorithm 2
yields the correct regression coefficient $\beta^{*}=\beta$ with probability
$\delta\in(0,1)$ for an appropriately selected margin parameter $\nu$.

###### Proof.

The success of algorithm 1 relies on the fact that the exhaustive search
eventually hits a tuple $(i,j)$ such that $j=\pi(i)$ which yields the correct
regression coefficient. Therefore, when randomly sampling
$(i,j)\sim[n]\times[m]$, the probability of choosing a corresponding pair is
$\frac{n-k}{n}\frac{1}{m}$. The probability of iterating $q$ times such hat no
correct correspondence is selected is $(1-(n-k)/(nm))^{q}=(1-\delta)$ where
$\delta$ is the desired success rate. Taking logs yields,
$q=\frac{\log(1-\delta)}{\log(1-(n-k)/(nm))}$ ∎

Algorithm 2 One dimensional robust regression without correspondence -
Randomized approach

Input:  Reference set: $\\{x_{1},\ldots,x_{m}\\}$, target set:
$\\{y_{1},\ldots,y_{n}\\}$, $\delta$ (probability of success), outlier margin:
$\nu$

Require:  $k<\frac{n}{2}$ (number of outliers)

for $t=1,\ldots,q$ do

Sample $i\sim[n]$ and sample $j\sim[m]$

Compute $\beta^{t}=y_{i}/x_{j}$

Compute linear assignment [33]:

$\boldsymbol{\Pi}^{t}\leftarrow\underset{\boldsymbol{\Pi}\in\mathcal{P}^{n\times
m}}{\arg\min}\|\mathbf{x}\beta^{t}-\boldsymbol{\Pi}^{T}\mathbf{y}\|_{2}^{2}$

Compute hypothetical inliers:

$\mathcal{I}^{t}=\\{l:|x_{\pi^{t}(l)}\beta^{t}-y_{l}|\leq\nu\\}$

end for

return $t^{*}=\underset{t}{\arg\max}|\mathcal{I}^{t}|$ ,
$\mathcal{I}^{*}=\mathcal{I}^{t^{*}}$,

$\boldsymbol{\Pi}^{*}=\boldsymbol{\Pi}^{i^{*},j^{*}}$,
$\beta^{*}\leftarrow\frac{\sum_{l\in\mathcal{I}^{*}}y_{l}x_{\pi^{*}(l)}}{\sum_{l\in\mathcal{I}^{*}}x_{\pi^{*}(l)^{2}}}$

The time complexity of randomized algorithm 2 is
$O\bigg{(}\frac{\log(1-\delta)}{\log(1-(n-k)/n^{2})}n^{3}\bigg{)}$.

Figure 2: Left: Recovery rate (colorbar) vs. missing data ratio (y-axis) vs.
SNR (x-axis), Middle: Recovery rate (colorbar) vs. outlier ratio (y-axis) vs.
SNR (x-axis), Right: Recovery rate (y-axis) vs. outlier ratio (x-axis), blue:
rRWOC, red: ICP, yellow: randomized rRWOC ($\delta=0.9)$), purple: randomized
rRWOC ($\delta=0.6$)

###  3.2 Randomized approximation algorithm for $d\geq 2$

The exhaustive approach for the $d\geq 2$ dimensional case requires
$\binom{n}{d}\binom{m}{d}$ $d$-subset comparisons of $\mathbf{X}$,$\mathbf{Y}$
in order to guarantee hitting correct (in the noiseless case) or approximately
correct (in the noisy case) regression coefficients, with complexity
$O(m^{d}n^{d})$. However, especially in higher dimensions, the randomized
procedure enables a substantial reduction of iterations to yield a high
probability correct triplet of inlier set, permutation, and regression
coefficients. The randomized algorithm for rRWOC in $d\geq 2$ is demonstrated
in algorithm 3. Conceptually, the idea of the algorithm is illustrated in
figure 1. Random ordered $d$-tuples of reference and target point sets are
sampled and are used to align the remainder of the point set. The number of
hypothetical inliers for each hypothetical correspondence is assessed by
checking whether the transformed reference points are arbitrarily close to a
target point. With high probability, if correct a $d$-tuple correspondence is
captured, the number of transformed reference points matching a target point
will be high (Figure 1 top), otherwise it will result in a partial coverage
(Figure 1 bottom).

Algorithm 3 Robust regression without correspondence - Randomized approach

Input:
$\mathbf{X}=[\mathbf{x}_{1}|\ldots|\mathbf{x}_{m}]^{T}\in\mathbb{R}^{m\times
d}$ (reference points),
$\mathbf{Y}=[\mathbf{y}_{1}|\ldots|\mathbf{y}_{n}]^{T}\in\mathbb{R}^{n\times
d}$ (target points), $\delta$ (probability of success), $\nu$ (outlier margin)

Require:  $k<\frac{n}{2}$ (number of outliers)

Output:  $\mathcal{I}\subseteq[n]$ (index of inliers),
$\hat{\boldsymbol{\Pi}}\in\mathcal{P}^{m\times|\mathcal{I}|}$ (permutation
matrix), $\hat{\boldsymbol{\beta}}\in\mathbb{R}^{d\times d}$ (regression
coefficients)

for $t=1,\ldots,q$ do

Sample $\boldsymbol{i}=(i_{1},\ldots,i_{d})\sim[n]^{d}$ w/o replacement

Sample $\boldsymbol{j}=(j_{1},\ldots,j_{d})\sim[m]^{d}$ w/o replacement

Compute
$\boldsymbol{\beta}^{t}=\underset{\boldsymbol{\beta}}{\arg\min}\|\mathbf{X}_{[\boldsymbol{j}]}\boldsymbol{\beta}-\mathbf{Y}_{[\boldsymbol{i}]}\|_{F}^{2}$

Compute linear assignment via [33]:

$\boldsymbol{\Pi}^{t}\leftarrow\underset{\boldsymbol{\Pi}\in\mathcal{P}^{m\times
n}}{\arg\min}\|\mathbf{X}\boldsymbol{\beta}^{t}-\boldsymbol{\Pi}\mathbf{Y}\|_{F}^{2}$

Compute hypothetical inliers:

$\mathcal{I}^{t}=\\{l:\|\mathbf{x}_{\pi^{t}(l)}\boldsymbol{\beta}^{t}-\mathbf{y}_{l}\|_{2}\leq\nu\\}$

end for

return $t^{*}=\arg\max_{t}|\mathcal{I}^{t}|$,
$\mathcal{I}^{*}=\mathcal{I}^{t^{*}}$,

$\boldsymbol{\Pi}^{*}=\boldsymbol{\Pi}^{t^{*}}_{\mathcal{I}^{*}}$,$\boldsymbol{\beta}^{*}\leftarrow\arg\min_{\boldsymbol{\beta}}\|\mathbf{X}_{\pi^{*}(\mathcal{I^{*}})}\boldsymbol{\beta}-\mathbf{Y}_{\mathcal{I^{*}}}\|_{F}^{2}$

######  Proposition 3.

For
$q\geq\frac{\log(1-\delta)}{\log\bigg{(}1-\frac{\binom{m-k}{d}}{\binom{m}{d}\binom{n}{d}}\bigg{)}}$,
algorithm 3 recovers $\boldsymbol{\beta}^{*}$ and $\boldsymbol{\Pi}^{*}$ and
the set of inliers for the noiseless case with probability $(1-\delta)$ using
arbibrarily small $\nu$. For sufficiently small noise variance and
appropriately chosen $\nu$, algorithm 3 recovers approximate
$\boldsymbol{\beta}^{*}$ with high probability.

###### Proof.

Analogous to the analysis of algorithm 2, the probability of drawing $d$
inliers out of $n$ points with k outliers in $\mathbf{Y}$ is
$\frac{\binom{n-k}{d}}{\binom{n}{d}}$. The probability of matching the drawn
inliers with the $d$ corresponding sampled reference points in $\mathbf{X}$ is
$\frac{1}{\binom{m}{d}}$. Probability that any draw is not going to match is
$1-\frac{\binom{n-k}{d}}{\binom{m}{d}\binom{n}{d}}$. The probability that $q$
draws will be incorrect is
$\bigg{(}1-\frac{\binom{m-k}{d}}{\binom{m}{d}\binom{n}{d}}\bigg{)}^{q}$. If we
set this to be the probability of failure $(1-\delta)$, we then have the
estimate for the number of draws we need to make as
$q(\delta,n,m,k)\geq\log(1-\delta)/\log\bigg{(}1-\frac{\binom{m-k}{d}}{\binom{m}{d}\binom{n}{d}}\bigg{)}$
∎

The complexity of algorithm 3 can be analyzed as follows. In each inner loop,
the regression coefficient solution requires $O(d^{3})$ time, the Hungarian
algorithm requires $O(nmd)$ to compute the input distance matrix and then
$O(\max\\{n,m\\}^{3})$ to optimize the permutation matrix. The rest of the
operations are $O(d)$. Therefore, the overall time complexity is

$\displaystyle
O\bigg{(}\frac{\log(1-\delta)}{\log\bigg{(}1-\frac{\binom{m-k}{d}}{\binom{m}{d}\binom{n}{d}}\bigg{)}}(d^{3}+nmd+\max\\{n,m\\}^{3})\bigg{)}.$
(4)

Figure 3: 2D projection of 3D fluorescence microscopy image of C.elegans head
in Yemini et al. dataset. Superimposed annotation points denote neuron
locations. Outliers are detections that do not correspond to neurons and
missing data are undetected neurons.

####  3.2.1 Margin parameter ($\nu$) selection

Both of the proofs of the noiseless and the noisy cases of proposition 1 rely
on knowledge of the true regression coefficient and the noise variance in
order to estimate the margin coefficient $\nu$ and output the optimal
regression coefficient with high probability. However, in practice, as in many
RANSAC-like robust regression settings, these parameters cannot be known
apriori, and $\nu$ is typically determined via empirical heuristics and or
cross-validation [26].

In the noiseless case, an appropriate heuristic is choosing $\nu$ arbitrarily
small since the correct regression should yield zero residual. However, for
the noisy case, if available, supervised data should be used with known
correspondences to estimate the actual dispersion of point correspondences.

##  4 Numerical results

To verify the theoretical guarantees of the proposed algorithms, simulated
data in 3 dimensions was generated in both noisy and noiseless regimes.
Furthermore, iterative solutions of $\boldsymbol{\beta}$ and
$\boldsymbol{\Pi}$ were obtained to demonstrate the suboptimality of local
minima found using block coordinate descent for this non-convex problem.

The neuroscience application of rRWOC was demonstrated in the context of point
set matching of neurons of C. elegans worms recorded using fluorescence
microscopy imaging. The matching accuracy with respect to ground truth was
assessed for rRWOC as well as a robust variant of the iterative closest point
(ICP) algorithm [8] known as trimmed ICP [17].

Computational setup and code: All experiments were performed on an Intel
i5-7500 CPU at 3.40GHz with 32GB RAM. MATLAB code for 3D versions of algorithm
3 are included in supplementary material along with sample C. elegans neuron
point clouds.

Figure 4: A: Unaligned point sets of reference C. elegans neuron positions
(red) and target neuron positions (green) B: Alignment with coherent point
drift algorithm [9], C: Alignment with iterative closest point algorithm [17],
D: Alignment with proposed algorithm 3.

###  4.1 Simulated data

Three dimensional source point set $\mathbf{X}$ was generated by sampling
$\mathbf{x}_{j}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I}_{3})$ for
$j=1,\ldots,J$ where $J\in[20,\ldots,40]$. A random transformation
$\boldsymbol{\beta}$ was obtained by computing the QR factorization of a
$3\times 3$ random gaussian matrix $\mathbf{M}$ such that
$\mathbf{Q}\mathbf{R}=\mathbf{M}$, taking the orthonormal rotation component
$\mathbf{Q}$. This was randomly scaled by a factor between $s=[0.5,1.5]$ so
that $\boldsymbol{\beta}=s\mathbf{Q}$. For $k\in[1,\ldots,19]$, $20-k$ inlier
target points were generated by transforming a random $20-k$ subset of
$\mathbf{X}$ by $\boldsymbol{\beta}$ and adding gaussian noise with varying
$\sigma^{2}$:
$\mathbf{Y}_{\mathcal{I}}=\mathbf{X}_{\pi(\mathcal{I})}\boldsymbol{\beta}+\epsilon$.
Furthermore, $k$ points in $\mathbf{Y}$ were randomly uniformly sampled from
the convex hull of the $20-k$ inlier points:
$\mathbf{Y}_{\mathcal{O}}\sim\text{U}[\mathcal{C}(\mathbf{Y}_{\mathcal{I}})]$.
This procedure yielded two unordered multisets,
$\mathbf{X}\in\mathbb{R}^{J\times 3}$ and $\mathbf{Y}\in\mathbb{R}^{20\times
3}$. Using these unordered multisets as input to rRWOC, the regression
coefficients $\hat{\boldsymbol{\beta}}$ were estimated. If
$\|\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}\|_{F}\leq 1e-3$, the event was
considered a correct recovery, otherwise a failure. The margin parameter $\nu$
was set to be $\nu=\sigma$. Also, using the randomized algorithm 3, the
success probability parameter was set to $\delta=0.9$.

This procedure was repeated 100 times for varying $k=1,\ldots,19$, varying
$\sigma^{2}$ and varying $J=20,\dots,40$ to assess the empirical recovery rate
as a function of outlier amount, SNR and missing correspondences in the
target, respectively. The recovery rates vs. outlier ratio, and SNR can be
seen in figure 2-middle. The recovery rates vs. missing data ratio and SNR can
be seen in figure 2-left. Lastly, the comparison of the recovery rate of
exhaustive and randomized rRWOC versus iterative closest point can be seen in
figure 2-right.

These empirical results demonstrate that for a sufficiently high SNR and
outlier ratio of less than $50\%$, the proposed algorithm yields almost
perfect recovery rates. Furthermore, the comparisons with iterative closest
point algorithm (ICP) show that rRWOC is much more robust to outliers than ICP
since the inclusion of any outliers results in failure of ICP to recover the
true transformation.

###  4.2 Neuron matching of C. elegans

For this application, we have used the publicly available C. elegans
fluorescence imaging dataset of Nguyen et al. [30] found at
http://dx.doi.org/10.21227/H2901H as well as the neuronal position dataset
provided in [31]. The worm C. elegans is a widely known model organism for
studying the nervous system due to the known structural connectome of the 302
neurons it contains. The data provided 3D z-stack images of the head of 14
worms that each consists of approximately 185 to 200 neurons captured under
confocal microscopy using florescent tagged protein GFP. In figure 3, the
depth-colored 2D projection of an image frame can be seen superimposed with
annotation points delineating the locations of neurons. Figure 3 also
highlights the need for a method of matching and aligning worm point clouds
that is robust to outliers or missing associations. Here, we define outliers
as points where there is no neuron present and define missing data as neurons
with no detection present.

Of the 14 datasets of the head neurons of C.elegans worms, random pairs were
drawn to be the source and target point sets. From the remaining worms, the
positional covariance of each neuron was estimated using the supervised
alignment method of [35]. Since the positional variance of each neuron was
uniquely identified using training data, we used variable margin parameters
for rRWOC such that $\nu_{l}=\underset{i=1,2,3}{\max}\lambda_{i}(\Sigma_{l})$
where $\Sigma_{l}$ is the covariance matrix of the $l$th neuron and
$\lambda_{i}(\cdot)$ denotes the $i$th eigenvalue. Randomized RWOC (algorithm
3) was deployed with $\delta=0.9$. The results were compared with iterative
closest point(ICP) [8] as well as coherent point drift (CPD) [9] algorithms.

The recovery rates in terms of recovering the transformation
$\boldsymbol{\beta}^{*}$ as well as the permutation $\boldsymbol{\Pi}^{*}$,
are summarized in table 1. In general, rRWOC was able to recover both the
transformation and permutation better than ICP and CPD, which tend to be
initialization-dependent. In all of the experiments, ICP and CPD were
initialized with random rotation. rRWOC is invariant to initialization since
it is not a descent based method.

Method | $TP$ | $FP$ | $TN$ | $FN$ | $ACC$ | $F1$ | $PREC$ | $REC$ | $MD$
---|---|---|---|---|---|---|---|---|---
rRWOC | 155 | 41 | 0 | 40 | 0.65 | 0.79 | 0.79 | 0.79 | 2.83
ICP [17] | 10 | 186 | 0 | 185 | 0.02 | 0.055 | 0.05 | 0.05 | 3.4
CPD [9] | 20 | 176 | 0 | 175 | 0.05 | 0.10 | 0.10 | 0.10 | 10.89
Table 1: Transformation recovery and permutation recovery by rRWOC, ICP, and
CPD in the C.elegans dataset. TP = true positive, FP = false positive, TN =
true negative, FN = false negative, ACC = accuracy, F1 = F1 score , PREC =
precision, REC = recall, MD = mean distance

###  4.3 Discussion

In this paper, we expanded on the linear regression without correspondence
model [25, 21, 24, 22] to account for missing data and outliers. Furthermore,
we provided several exact and approximate algorithms for the recovery of
regression coefficients under noiseless and noisy regimes. The proposed
algorithms are combinatorial at worst with variable dimension. However,
randomization procedures make the average case complexity in constant
dimension tractable given enough tolerance for failure. We provided several
theoretical guarantees for exact recovery and running time complexity.
Furthermore, we empirically demonstrated the recovery rates of the proposed
algorithms in simulated and biological data. A future algorithmic direction is
to employ branch and bound techniques found in [23] to reduce computational
complexity of the brute force nature of the algorithms.

The proposed methods can be thought of as a general framework for dissociating
the outliers from a model-based data transformation process. The same
principles can apply for the cases where either the generative noise is non-
Gaussian, or some prior information exists about the structure of the
outliers. Case-specific noise analysis is required for a particular model
selection. Future work can focus on finding theoretical bounds on the
robustness of the inlier recovery as a function of the number of outliers and
the statistics of the generative and outlier distributions.

### Acknowledgements

The authors would like to acknowledge the sources of funding: NSF NeuroNex
Award DBI-1707398, and The Gatsby Charitable Foundation.

##  Appendix A Proof of proposition 1

###  A.1 Noiseless case

Among the $mn$ hypothetical regression coefficients obtained through all
possible pairs of $x_{i}$ and $y_{j}$, if a correct correspondence is
encountered (i.e. $j=\pi^{*}(i))$, we have $y_{\pi(i)}=x_{i}\beta^{*}$ where
$\beta^{*}$ is the true coefficient. Therefore if we let
$\beta^{i}=\frac{y_{\pi(i)}}{x_{i}}$ then $\beta^{i}=\beta^{*}$. Using this
estimate, the distances of the remaining covariates regressed to their
corresponding responses is

$\displaystyle x_{l}\beta^{i}-y_{\pi(l)}=x_{l}\beta^{*}-\mathbf{y}_{\pi(l)}=0$

Therefore, when computing
$\underset{\boldsymbol{\Pi}\in\mathcal{P}}{\min}\|\mathbf{x}\beta^{i}-\boldsymbol{\Pi}\mathbf{y}\|_{2}^{2}$
via the Hungarian algorithm [33], each column of the distance matrix
$[D]_{p,q}=|x_{p}\beta^{i}-y_{q}|$ corresponding to inlier points in
$\mathbf{y}$ (i.e. $q\in\mathcal{I}^{*}$) will have at least one zero element.
Thus, the optimal assignment $\boldsymbol{\Pi}^{i}$ will include all of the
permutations $\pi^{i}(l)=l$ since they incur zero cost. Since there are $m-k$
of them by assumption 1, then
$\sum_{l}\boldsymbol{1}(|x_{l}-y_{\pi^{i}(l)}|\leq\epsilon/2))\geq m-k$. This
is inequality because there might be additional outlier points that are by
chance close to the regressed points.

Conversely, for a pair $(x_{i},y_{\pi(k)})$ where $k\neq i$, we have the
estimated coefficient
$\beta^{i,k}=\frac{y_{\pi(k)}}{x_{i}}=\frac{x_{k}\beta^{*}}{x_{i}}$. The
distances of the remaining covariates regressed with this estimate to their
corresponding responses are

$\displaystyle
x_{l}\beta^{i,k}-y_{\pi(l)}=\frac{x_{l}x_{k}\beta^{*}}{x_{i}}-y_{\pi(l)}=x_{l}\beta^{*}\bigg{(}\frac{x_{k}}{x_{i}}-1\bigg{)}$

Therefore, without loss of generality, assuming $x_{l}\neq 0$ (if $x_{l}=0$
the correspondence $(x_{l},y_{\pi(l)})$ can be automatically inferred by
choosing any $y_{\pi(l)}=0$. If there aren’t any $y_{j}=0$, then this implies
$x_{l}$ is a point without correspondence in $\mathbf{y}$), we have

$\displaystyle|x_{l}\beta^{i,k}-y_{\pi(l)}|\geq\epsilon$

for some $\epsilon>0$. $\epsilon$ can be explicitly stated as

$\displaystyle\epsilon=\underset{i,l,k,~{}~{}i\neq
k}{\min}x_{l}\beta^{*}\bigg{(}\frac{x_{k}}{x_{i}}-1\bigg{)}$

On the other hand,

$\displaystyle x_{i}\beta^{i,k}-y_{\pi(k)}=0$

by construction.

Therefore, when computing
$\underset{\boldsymbol{\Pi}\in\mathcal{P}}{\min}\|\mathbf{x}\beta^{i,k}-\boldsymbol{\Pi}\mathbf{y}\|_{2}^{2}$
via Hungarian algorithm, there will less than $m-k$ assignments in the optimal
assignment $\boldsymbol{\Pi}^{i,k}$ such that
$|x_{l}-y_{\pi^{i,k}(l)}|\leq\epsilon/2$. Otherwise, this would imply the
coefficient $\beta^{i,k}$ is a coefficient that explains the inliers, which by
assumption 1 cannot be the case. Thus,
$\sum_{l}\boldsymbol{1}(|x_{l}-y_{\pi^{i,k}(l)}|\leq\epsilon/2))<m-k$.

This shows that the maximal cardinality of a hypothetical inlier set is at
least $m-k$, and it is only achieved for a coefficient that is obtained by a
correct correspondence pair. This is sufficient to show that algorithm 1
recovers the true coefficient $B^{*}$ under the noiseless regime.

###  A.2 Noisy case

Let the noise model of the inlier regression be
$\epsilon\sim\mathcal{N}(0,\sigma^{2})$. Therefore, if a correct
correspondence is encountered, we have $y_{\pi(i)}=x_{i}\beta^{*}+\epsilon$
where $\beta^{*}$ is the true coefficient. The coefficient estimated from this
pairing is
$\beta^{i}=\frac{y_{\pi}(i)}{x_{i}}=\beta^{*}+\frac{\epsilon}{x_{i}}$. When
this coefficient is applied to $\mathbf{x}$ we see that

$\displaystyle\text{E}(x_{l}\beta^{i}-y_{\pi(l)})=0,~{}\text{Var}(x_{l}\beta^{i}-y_{\pi(l)})=\bigg{(}\frac{x_{l}^{2}}{x_{i}^{2}}+1\bigg{)}\sigma^{2}$
$\displaystyle\text{E}(x_{l}\beta^{i}-y_{\pi(k)})=(x_{l}-x_{k})\beta^{*}$
$\displaystyle\text{Var}(x_{l}\beta^{i}-y_{\pi(k)})=\bigg{(}\frac{x_{l}^{2}}{x_{i}^{2}}+1\bigg{)}\sigma^{2}$

Therefore, if $\sigma^{2}$ is small (i.e. in the SNR regime of [22]), we have
$|x_{l}\beta^{i}-y_{\pi(l)}|<|x_{l}\beta^{i}-y_{\pi(k)}|$ for $l\neq k$ with
high probability. Thus the row-wise minimal cost assignment in the Hungarian
algorithm will be $\pi^{i}(l)=l$ with high probability. However, even if
$\pi^{i}(l)\neq l$, if we set margin $\nu$ such that
$\nu=\frac{1}{2}\underset{l,k~{}l\neq k}{\min}|(x_{l}-x_{k})\beta^{*}|$, with
high probability we will have that

$\displaystyle\sum_{l}\boldsymbol{1}(|x_{l}\beta^{i}-y_{\pi^{i}(l)}|\leq\nu)\geq\sum_{l}\boldsymbol{1}(|x_{l}\beta^{i,k}-y_{\pi^{i}(l)}|\leq\nu)$
(5)

where $\beta^{i,k}$ denotes the regression coefficient obtained via incorrect
correspondence $\beta^{i,k}=\frac{y_{pi(k)}}{x_{i}}$. Therefore, if
$\sigma^{2}$ is sufficiently small, with high probability, algorithm 1
recovers the coeffient $\beta^{i}=\beta^{*}+\frac{\epsilon}{x_{i}}$ for some
$i\in\mathcal{I}$ where $\mathcal{I}$ denotes the set of inliers.

## References

  * [1] Oliver Van Kaick, Hao Zhang, Ghassan Hamarneh, and Daniel Cohen-Or.  A survey on shape correspondence.  In Computer Graphics Forum, volume 30, pages 1681–1707. Wiley Online Library, 2011.
  * [2] Gary KL Tam, Zhi-Quan Cheng, Yu-Kun Lai, Frank C Langbein, Yonghuai Liu, David Marshall, Ralph R Martin, Xian-Fang Sun, and Paul L Rosin.  Registration of 3d point clouds and meshes: a survey from rigid to nonrigid.  IEEE transactions on visualization and computer graphics, 19(7):1199–1217, 2013.
  * [3] Ji Zhang and Sanjiv Singh.  Visual-lidar odometry and mapping: Low-drift, robust, and fast.  In 2015 IEEE International Conference on Robotics and Automation (ICRA), pages 2174–2181. IEEE, 2015.
  * [4] Michel A Audette, Frank P Ferrie, and Terry M Peters.  An algorithmic overview of surface registration techniques for medical imaging.  Medical image analysis, 4(3):201–217, 2000.
  * [5] Bertram Drost, Markus Ulrich, Nassir Navab, and Slobodan Ilic.  Model globally, match locally: Efficient and robust 3d object recognition.  In 2010 IEEE computer society conference on computer vision and pattern recognition, pages 998–1005. Ieee, 2010.
  * [6] Jean-Charles Bazin, Yongduek Seo, Richard Hartley, and Marc Pollefeys.  Globally optimal inlier set maximization with unknown rotation and focal length.  In European Conference on Computer Vision, pages 803–817. Springer, 2014.
  * [7] Greg Bubnis, Steven Ban, Matthew D DiFranco, and Saul Kato.  A probabilistic atlas for cell identification.  arXiv preprint arXiv:1903.09227, 2019.
  * [8] Paul J Besl and Neil D McKay.  Method for registration of 3-d shapes.  In Sensor Fusion IV: Control Paradigms and Data Structures, volume 1611, pages 586–607. International Society for Optics and Photonics, 1992\.
  * [9] Andriy Myronenko and Xubo Song.  Point set registration: Coherent point drift.  IEEE transactions on pattern analysis and machine intelligence, 32(12):2262–2275, 2010.
  * [10] Qian-Yi Zhou, Jaesik Park, and Vladlen Koltun.  Fast global registration.  In European Conference on Computer Vision, pages 766–782. Springer, 2016.
  * [11] Anders Hast, Johan Nysjö, and Andrea Marchetti.  Optimal ransac-towards a repeatable algorithm for finding the optimal set.  2013\.
  * [12] Sandy Irani and Prabhakar Raghavan.  Combinatorial and experimental results for randomized point matching algorithms.  Computational Geometry, 12(1-2):17–31, 1999.
  * [13] Dror Aiger, Niloy J Mitra, and Daniel Cohen-Or.  4-points congruent sets for robust pairwise surface registration.  In ACM transactions on graphics (TOG), volume 27, page 85. Acm, 2008\.
  * [14] David M Mount, Nathan S Netanyahu, and Jacqueline Le Moigne.  Efficient algorithms for robust feature matching.  Pattern recognition, 32(1):17–38, 1999.
  * [15] Piotr Indyk, Rajeev Motwani, and Suresh Venkatasubramanian.  Geometric matching under noise: Combinatorial bounds and algorithms.  In SODA, pages 457–465, 1999.
  * [16] Jonathan Pokrass, Alexander M Bronstein, Michael M Bronstein, Pablo Sprechmann, and Guillermo Sapiro.  Sparse modeling of intrinsic correspondences.  In Computer Graphics Forum, volume 32, pages 459–468. Wiley Online Library, 2013.
  * [17] Dmitry Chetverikov, Dmitry Svirko, Dmitry Stepanov, and Pavel Krsek.  The trimmed iterative closest point algorithm.  In Object recognition supported by user interaction for service robots, volume 3, pages 545–548. IEEE, 2002.
  * [18] Jiaolong Yang, Hongdong Li, Dylan Campbell, and Yunde Jia.  Go-icp: A globally optimal solution to 3d icp point-set registration.  IEEE transactions on pattern analysis and machine intelligence, 38(11):2241–2254, 2016.
  * [19] Jayakrishnan Unnikrishnan, Saeid Haghighatshoar, and Martin Vetterli.  Unlabeled sensing: Solving a linear system with unordered measurements.  In 2015 53rd Annual Allerton Conference on Communication, Control, and Computing (Allerton), pages 786–793. IEEE, 2015.
  * [20] David L Donoho et al.  Compressed sensing.  IEEE Transactions on information theory, 52(4):1289–1306, 2006\.
  * [21] Abubakar Abid, Ada Poon, and James Zou.  Linear regression with shuffled labels.  arXiv preprint arXiv:1705.01342, 2017.
  * [22] Ashwin Pananjady, Martin J Wainwright, and Thomas A Courtade.  Linear regression with an unknown permutation: Statistical and computational limits.  In 2016 54th Annual Allerton Conference on Communication, Control, and Computing (Allerton), pages 417–424. IEEE, 2016.
  * [23] Manolis Tsakiris and Liangzu Peng.  Homomorphic sensing.  In International Conference on Machine Learning, pages 6335–6344, 2019.
  * [24] Daniel J Hsu, Kevin Shi, and Xiaorui Sun.  Linear regression without correspondence.  In Advances in Neural Information Processing Systems, pages 1531–1540, 2017.
  * [25] Jayakrishnan Unnikrishnan, Saeid Haghighatshoar, and Martin Vetterli.  Unlabeled sensing with random linear measurements.  IEEE Transactions on Information Theory, 64(5):3237–3253, 2018\.
  * [26] Martin A Fischler and Robert C Bolles.  Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography.  Communications of the ACM, 24(6):381–395, 1981.
  * [27] Philip HS Torr and Andrew Zisserman.  Mlesac: A new robust estimator with application to estimating image geometry.  Computer vision and image understanding, 78(1):138–156, 2000.
  * [28] Heng Yang and Luca Carlone.  A polynomial-time solution for robust registration with extreme outlier rates.  arXiv preprint arXiv:1903.08588, 2019.
  * [29] Dagmar Kainmueller, Florian Jug, Carsten Rother, and Gene Myers.  Active graph matching for automatic joint segmentation and annotation of c. elegans.  In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 81–88. Springer, 2014.
  * [30] Jeffrey P Nguyen, Ashley N Linder, George S Plummer, Joshua W Shaevitz, and Andrew M Leifer.  Automatically tracking neurons in a moving and deforming brain.  PLoS computational biology, 13(5):e1005517, 2017.
  * [31] Eviatar Yemini, Albert Lin, Amin Nejatbakhsh, Erdem Varol, Ruoxi Sun, Gonzalo E Mena, Aravinthan DT Samuel, Liam Paninski, Vivek Venkatachalam, and Oliver Hobert.  Neuropal: A neuronal polychromatic atlas of landmarks for whole-brain imaging in c. elegans.  bioRxiv, page 676312, 2019.
  * [32] Edwin F Beckenbach and Richard Bellman.  Inequalities, volume 30.  Springer Science & Business Media, 2012.
  * [33] Harold W Kuhn.  The hungarian method for the assignment problem.  Naval research logistics quarterly, 2(1-2):83–97, 1955.
  * [34] Roy Jonker and Ton Volgenant.  Improving the hungarian assignment algorithm.  Operations Research Letters, 5(4):171–175, 1986.
  * [35] Georgios Dimitrios Evangelidis and Radu Horaud.  Joint alignment of multiple point sets with batch and incremental expectation-maximization.  IEEE transactions on pattern analysis and machine intelligence, 40(6):1397–1410, 2018.