Bilevel Inverse Problems in Neuromorphic Imaging

Harbir Antil H. Antil. The Center for Mathematics and Artificial Intelligence and Department of Mathematical Sciences, George Mason University, Fairfax, VA 22030, USA. [email protected] and David Sayre D. Sayre. The Center for Mathematics and Artificial Intelligence and Department of Mathematical Sciences, George Mason University, Fairfax, VA 22030, USA. [email protected]

Abstract.

Event or Neuromorphic cameras are novel biologically inspired sensors that record data based on the change in light intensity at each pixel asynchronously. They have a temporal resolution of microseconds. This is useful for scenes with fast moving objects that can cause motion blur in traditional cameras, which record the average light intensity over an exposure time for each pixel synchronously. This paper presents a bilevel inverse problem framework for neuromorphic imaging. Existence of solution to the inverse problem is established. Second order sufficient conditions are derived under special situations for this nonconvex problem. A second order Newton type solver is derived to solve the problem. The efficacy of the approach is shown on several examples.

Key words and phrases:

Bilevel optimization, neuromorphic imaging, existence of solution, second order sufficient conditions

2010 Mathematics Subject Classification:

65K05, 90C26, 90C46, 49J20

This work is partially supported by NSF grants DMS-2110263, DMS-1913004, the Air Force Office of Scientific Research (AFOSR) under Award NOs: FA9550-22-1-0248 and FA9550-19-1-0036.

1. Introduction

Event (Neuromorphic) cameras are novel biologically inspired sensors that record data based on the change in light intensity at each pixel asynchronously [9]. If the change in light intensity at a given pixel is greater than a preset threshold then an event is recorded at that pixel. For this reason if there is no change to the scene, be that movement or brightening/dimming of a light source, no events will be recorded. In contrast if a scene is dynamic from camera movement or from the movement of an object in the scene then each pixel of the event camera will record intensity changes with a temporal resolution on the order of microseconds [9]. This logging of events results in a non-redundant stream of events through the time dimension for each pixel [16]. The stream of data is exceptionally useful for scenes with fast moving objects that can cause motion blur in traditional cameras, which record the average light intensity over an exposure time for each pixel synchronously [13]. Figure 1 gives an illustration of the differences between traditional and event cameras.

Refer to caption — Figure 1. Comparison of the output of a standard frame-based camera and an event camera when facing a black dot on a rotating disk. The standard camera outputs frames at a fixed rate, thus sending redundant information when there is no motion in the scene [10].

The prior work, [13], proposed the Event-based Double Integral (EDI) and multiple Event-based Double integral (mEDI) algorithms to address motion blur for the underlying inverse problem. The EDI model utilizes a single blurry standard camera image, along with associated event data, to generate one or more clear images through the process of energy minimization. The mEDI model, an extension of the EDI model, utilizes multiple images in conjunction with event data to produce unblurred images. The EDI and mEDI models were an extension of the work from [16] where a continuous-time model was initially introduced. While these algorithms represent a leap forward over earlier work, questions remain about the validity of the models in various scenarios (see below), regarding what situations solutions exist and what optimization procedures are appropriate in order to find such solutions. In this paper we provide a mathematical foundation and further extensions to the previous methods by extending the model and providing a bilevel inverse problem framework. Indeed, as shown below, the EDI and mEDI models correspond to inverse problems in image deblurring [3, 4]. In particular,

(a)

We introduce a bilevel optimization problem that allows for the simultaneous optimization over a desired number of frames. We also illustrate that this optimization problem is indeed an inverse problem. Therefore, we use the terms ‘optimization’ and ‘inverse’, interchangeably.
(b)

We establish existence of solution to this problem. Under certain conditions, we derive the second order sufficient conditions and provide local uniqueness of solution to these nonconvex problems.
(c)

A fully implementable framework based on second order methods has been developed. The benefits of the proposed approach are illustrated using multiple numerical examples.

For completeness, we emphasize that the bilevel approaches to search for the hyperparameters in traditional imaging is not new, see for instance [2, 8]. However, the problem considered in this paper is naturally bilevel without having to search for the hyperparameters.

Outline: The remainder of the paper is organized as follows. Section 2 focuses on some notation and preliminary results. In Section 2.1, we discuss the basics of event based cameras. Section 2.3 focuses on existing models from [13] which serves as a foundation for the proposed bilevel optimization based variational model in Section 3 . Existence of solution to the proposed optimization problem is shown in Theorem 3.1. Moreover, local convexity of our reduced functional is shown in Theorem 3.5. Section 4 first focuses on how to prepare the data for the algorithm. It then discuss implementation details. The problem itself is solved using a second order Newton based method. Subsequently, several illustrative numerical examples are provided in Section 5 .

2. Notation and Preliminaries

Let $n_{t}$ represent the number of images and $n_{x}\times n_{y}$ denotes the number of pixels per image. We use $\bm{U}\in\mathbb{R}^{n_{t}\times n_{x}\times n_{y}}$ to denote a tensor. We use the vector

\bm{u}_{xy}\in\mathbb{R}^{n_{t}},

to represent the $(x,y)$ pixel value of $\bm{U}$ for all times. Moreover, we use the matrix

\bm{u}^{i}\in\mathbb{R}^{n_{x}\times n_{y}}.

to represent an image of size $\mathbb{R}^{n_{x}\times n_{y}}$ at a fixed time instance $i\in\{1,\dots,n_{t}\}$ . We use the scalar

\bm{u}^{i}_{xy}\in\mathbb{R}

to denote the value of $\bm{U}$ at a pixel location $(x,y)$ and time $i\in\{1,\dots,n_{t}\}$ . Graphical representations are given in Figure 2 . Finally, given quantities $\bm{u}_{xy}$ and $\bm{M}\in\mathbb{R}^{{n_{t}}\times{\eta}}$ for some $\eta\in\mathbb{N}$ , we define the operation $\langle\cdot,\cdot\rangle$ as:

\langle\bm{u}_{xy},\bm{M}\rangle\mathrel{\mathop{\mathchar 58\relax}}=(\bm{u}_{xy})^{\top}\bm{M}\in\mathbb{R}^{\eta}.

The remainder of the section is organized as follows. First in Section 2.1, we discuss the basic working principle behind the neuromorphic cameras. The main ideas behind the algorithms presented in [13] are described in Section 2.3. We briefly discuss some potential limitations of this approach which motivates our algorithm in Section 3.

2.1. Neuromorphic (Event Based) Cameras

Neuromorphic Camera Basics

Neuromorphic cameras are composed of independent pixels that detect light intensity changes in the environment as they occur. Light intensity is sampled $\mathcal{O}(\mu s)$ and events are logged if the intensity of the light is beyond a preset hardware defined threshold, $c$ . We denote the instantaneous light intensity at pixel location $(x,y)$ at time $s$ by $\bm{q}^{s}_{xy}\in\mathbb{R}$ . When the light intensity detected by the camera exceeds the threshold for a given pixel located at $(x,y)$ , at a given time $s$ , an event is logged. Then the reference intensity, $\bm{q}^{s_{ref}}_{xy}$ , for the pixel located at $(x,y)$ is updated. Due to the high rate of sampling, and independent nature of the camera’s pixels neuromorphic cameras are less susceptible to exposure issues, as well as image blurring. Figure 1 illustrates the difference in the data recorded between traditional and event based cameras.

Data representation

The output of an event based camera over some time interval, $\mathcal{I}$ , is of the form $(s,x,y,p)$ , with the following definitions:

$\bullet$

$s\in\mathcal{I}=\{s_{1},s_{2},\dots,s_{n}\,\mathrel{\mathop{\mathchar 58\relax}}\,s_{1}<s_{2}<\dots<s_{n}\}$ , where $\mathcal{I}$ represents the exposure interval, each $s_{j}$ represents a discrete time at which an event occurred and $n$ represents the total number of events recorded during the interval, $\mathcal{I}$ . When discussing methods that involve multiple images, we will denote the exposure interval for the $i$ - $th$ image as $\mathcal{I}_{i}$ . We will also denote the set of events associated to each image $i$ as $\{s^{i}_{j}\}_{j=1}^{n_{i}}$ . With this notation we have $i=1,\dots,n_{t}$ and $j=1,\dots,n_{i}$ with $n_{i}$ representing the total number of events associated to image $i$ .
$\bullet$

$x,y$ represent pixel coordinates.
$\bullet$

$p_{x,y}^{s_{j}}\mathrel{\mathop{\mathchar 58\relax}}=\begin{cases}+1,&\log\left(\frac{\bm{q}^{s_{j}}_{xy}}{\bm{q}^{s_{ref}}_{xy}}\right)\geq c\\ -1,&\log\left(\frac{\bm{q}^{s_{j}}_{xy}}{\bm{q}^{s_{ref}}_{xy}}\right)\leq-c\end{cases}$ is the polarity of the event. An increase in light intensity above a threshold, $c>0$ , we regard as a polarity 1 and a decrease of intensity as a polarity shift of $-1$ . No event is logged if $p_{x,y}^{s_{j}}\in(-c,c)$ .

We can reformulate these events into a datacube with the following definition:

{\bm{EC}}_{xy}^{s_{j}}=p_{x,y}^{s_{j}},\quad j=1,\dots,n,

(2.1)

which results in a sparse 3-dimensional matrix with entries of $0,1,-1$ . We refer to Figure 5 for 2D and 3D examples of the datacube. Representing the data in this way allows us to track intensity changes over time on a per pixel basis.

2.2. Inverse problems in Image Deblurring

Image deblurring is a classic example of an inverse problem, where the goal is to reconstruct the unknown original image from a degraded, blurred observation. In image deblurring, the degradation process can be modeled as a interaction between the original image and a blur kernel which represents the blur caused by various factors such as movement [19, 20].

In the context of image deblurring, the field can be divided into two categories, blind and non-blind, based on the amount of information known about the blur kernel [1, 15, 17, 18, 19]. In the non-blind case, the blur kernel is known which leads to a simpler ableit ill-posed inverse problem [17, 18]. In the blind case, no information about the blur kernel is available, making the deblurring problem more challenging as the blur kernel must be estimated from the degraded image in addition to the original image [15, 19]. The semi-blind case is a particular instance of a blind image deblurring problem in that some but not all information about the blur kernel is available [5, 12, 14]. The problems under consideration in this project are of semi-blind type.

2.3. Existing Model, Algorithm, and Limitations

The approach discussed in this paper is motivated by [13]. The article [13] seeks to find latent un-blurred images using the Event Based Double Integral (EDI) and multi-Event Based Double Integral (mEDI) models. We briefly discuss this next. This will be followed by our new proposed models.

EDI Model

In [13], the discrete events outlined in Section 2.1 are used to generate a continuous function, $\bm{e}_{xy}$ , for each $(x,y)$ pixel location. This is done by generating a series of continuous unit bump functions centered at each $s_{j}$ : $\phi_{s_{j}}(t)$ . Then we can define $\bm{e}_{xy}\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}\rightarrow\mathbb{R}$ as:

\bm{e}_{xy}(t)\mathrel{\mathop{\mathchar 58\relax}}=\sum_{s_{j}\in\;\mathcal{I}}{\bm{EC}}_{xy}^{s_{j}}\cdot\phi_{s_{j}}(t)\in\mathbb{R}.

(2.2)

We graphically illustrate building this function for a single pixel below with the following string of events:

[1,0,1,-1,1].

In Figure 3 we see the function $\bm{e}_{xy}$ overlayed onto a stem plot of events.

From 2.2, we define the ‘sum of events’ from time of interest $t_{i}$ to arbitrary time $t$ as:

\bm{E}^{i}_{xy}(t,t_{i})=\int_{t_{i}}^{t}\bm{e}_{xy}(r)\;dr,

(2.3)

where $\bm{E}^{i}_{xy}(t,t_{i})\in\mathbb{R}$ corresponds to one pixel and $\bm{E}^{i}(t,t_{i})=\{\bm{E}^{i}_{xy}(t,t_{i})\}_{xy}\in\mathbb{R}^{{n_{x}}\times{n_{y}}}$ consists of all the pixels. We note that $t_{i}$ , $i=1,\dots,n_{t}$ are manually chosen and represent the time at which we wish to generate an image. In practice we will set each $t_{i}$ to correspond to the timestamps of the standard camera images we are using in the model. For convenience we will omit $t_{i}$ from $\bm{E}^{i}_{xy}(t,t_{i})$ unless needed for clarity. Introducing a scalar variable $z$ gives us the EDI Model:

\bm{B}^{i}_{xy}=\frac{1}{|\mathcal{I}_{i}|}\int_{t_{i}-|\mathcal{I}_{i}|/2}^{t_{i}+|\mathcal{I}_{i}|/2}\bm{L}^{i}_{xy}\cdot\exp\big{(}z\cdot\bm{E}^{i}_{xy}(t)\big{)}dt\,.

(2.4)

Where $\bm{B}^{i}_{xy}\in\mathbb{R}$ and $\bm{B}^{i}=\{\bm{B}^{i}_{xy}\}_{xy}\in\mathbb{R}^{n_{x}\times n_{y}}$ is a standard camera image. Furthermore, $\bm{L}^{i}_{xy}\in\mathbb{R}$ is the pixel value of the unblurred image at location $(x,y)$ , and $|\mathcal{I}_{i}|$ is the exposure time for the $i^{th}$ image. We note that the superscript $i$ in $\bm{B}^{i}_{xy}$ and $\bm{L}^{i}_{xy}$ are simply placeholders for notation since the EDI model only considers a single image at a time. Taking the log and re-arranging terms of (2.4) leads to:

\bm{v}_{xy}^{i}=\bm{d}_{xy}^{i}-\bm{g}_{xy}^{i}(z),

(2.5)

where

\displaystyle\begin{aligned} \bm{v}_{xy}^{i}&=\log(\bm{L}_{xy}^{i}),\quad\bm{d}_{xy}^{i}=\log\bm{B}_{xy}^{i},\\ \bm{g}^{i}_{xy}(z)&=\log\bigg{(}\frac{1}{|\mathcal{I}_{i}|}\int_{t_{i}-|\mathcal{I}_{i}|/2}^{t_{i}+|\mathcal{I}_{i}|/2}\exp\big{(}z\cdot\bm{E}_{xy}^{i}(t)\big{)}dt\bigg{)}.\end{aligned}

(2.6)

The formulation in (2.5) gives us the correlation between the EDI model and other semi-blind inverse problems in image deblurring. Here $\bm{B^{i}}$ is our blurry image, $\bm{L^{i}}$ is our latent unblurred image, and $\bm{g^{i}}$ is our kernel function. As shown in [1, 3, 4, 6, 7] this formulation is common for many inverse problems. We note the EDI/mEDI models would fall into the category of semi-blind inverse problems mentioned in Section 2.2 since the kernel, $\bm{g^{i}}$ , is known up to the parameter $z$ .

mEDI Model

If we consider two latent images $\bm{v}^{i}$ and $\bm{v}^{i+1}$ we know that event data will give us the per pixel changes that have occurred between the two images. We can mathematically describe this change as $\int_{t_{i}}^{t_{i+1}}\bm{e}_{xy}(s)ds$ where $\bm{e}_{xy}$ is as given in (2.2) and $t_{i}$ are user chosen times of interest. Incorporating the threshold variable $z$ allows us to formulate the following equations

\bm{v}^{i+1}_{xy}-\bm{v}^{i}_{xy}=z\int_{t_{i}}^{t_{i+1}}\bm{e}_{xy}(s)ds=\mathrel{\mathop{\mathchar 58\relax}}\bm{b}_{xy}^{i}(z).

(2.7)

Combining (2.5) and (2.7) gives us the following system of equations:

\left(\begin{array}[]{rrrrrrrr}-1&1&&&&&&\\ &-1&1&&&&&\\ &&&\ddots&\ddots&&&\\ &&&&&-1&1&\\ &&&&&&-1&1\\ \hline\cr 1&&&&&&&\\ &1&&&&&&\\ &&1&&&&&\\ &&&\ddots&&&&\\ &&&&1&&&\\ &&&&&1&&\\ &&&&&&1&\end{array}\right){\bm{v}}_{xy}=\left(\begin{array}[]{rrr}\bm{b}_{xy}\\ \bm{d}_{xy}-\bm{g}_{xy}\end{array}\right)

(2.8)

Notice that, (2.8) gives us the mEDI model described by [13]. The size the upper half of the submatrix in (2.8) is $\mathbb{R}^{(n_{t}-1)\times n_{t}}$ and the lower half matrix is $\mathbb{R}^{{n_{t}}\times{n_{t}}}$ .

mEDI Model - Optimizing for z

The variable $z$ is found by solving the following minimization problem for each pixel $(x,y)$ of an image of interest $i$ :

\min_{z}\frac{1}{2}\|{\bm{v}}^{i}_{xy}({z})+{\bm{g}}^{i}_{xy}(z)-{\bm{d}}^{i}_{xy}\|^{2}_{2}.

(2.9)

Here $\bm{v}_{xy}$ , ${\bm{d}}^{i}_{xy}$ , and ${\bm{g}}^{i}_{xy}(z)$ are as given in (2.6). We refer to [13] for their solution approach. Notice that (2.7) is not used as constraint to the minimization problem (2.9), but (2.9) is used to generate the image reconstructions.

Limitations of the mEDI and EDI Models

While the mEDI model uses multiple frames in order to accurately represent de-blurred images some limitations for the model remain.

(a)

Performance of the model can degrade for images that contain large variations of light intensity within a single scene. Some examples of this are having very dark objects in a bright room with a highly dynamic scene. In these instances we can see “shadowing” effects where shadows appear across multiple frames. We refer to Section 5.3 for specific examples.
(b)

The optimization problem introduced (2.9) to identify $z$ only seeks to optimize the pixels of one image at a time and neglects the constraints (2.7). Then uses (2.7) to reconstruct multiple images with the same $z$ value.
(c)

Recall from 2.7 that $\bm{b}^{i}_{xy}\mathrel{\mathop{\mathchar 58\relax}}=z\int_{t_{i}}^{t_{i+1}}\bm{e}_{xy}(s)\;ds$ . This definition does not include any information from the standard images $\bm{B}^{i}$ . Due to this the mEDI model is not well suited for generating image reconstructions across time-spans that include multiple standard images.

3. Proposed Model and Algorithm

3.1. Model

For each pixel $(x,y)$ , we consider the following bilevel inverse problem

	$\min_{\bm{z}_{xy}}J(\bm{u}_{xy},\bm{z}_{xy})\mathrel{\mathop{\mathchar 58\relax}}=\frac{1}{2}\\|\bm{u}_{xy}(\bm{z}_{xy})+\bm{g}_{xy}(\bm{z}_{xy})-\bm{d}_{xy}\\|_{2}^{2}+\frac{\lambda_{1}}{2}\\|\bm{z}_{xy}\\|_{2}^{2}$	(3.1a)
subject to constraints
	$\min_{\bm{u}_{xy}}\frac{1}{2}\\|{\bm{A}}\bm{u}_{xy}-\bm{b}_{xy}(\bm{z}_{xy})\\|_{2}^{2}+\frac{\lambda_{2}}{2}\\|\bm{u}_{xy}\\|_{2}^{2}.$	(3.1b)

Notice, that this is an optimization problem for one pixel across all images. By $0<\lambda_{1},\lambda_{2}$ , we denote the regularization parameters. Moreover,

\bm{A}=\left(\begin{array}[]{rrrrrrrr}-1&1&&&&&&\\ &-1&1&&&&&\\ &&&\ddots&\ddots&&&\\ &&&&&-1&1&\\ &&&&&&-1&1\\ \end{array}\right)\in\mathbb{R}^{(n_{t}-1)\times n_{t}}

is a submatrix of the matrix given in (2.8). Recall in Equation (2.9) that $z$ is taken as a scalar and the optimization is done over a single image at a time. We consider $\bm{z}_{xy}$ as size $n_{t}\times 1$ so that we can optimize over the same pixel $(x,y)$ for all images simultaneously. Also, in contrast to [13] since we have defined $\bm{z}_{xy}$ as a vector over all images, we also modify the variable $\bm{b}$ with the following definition:

\bm{b}^{i}_{xy}(\bm{z}_{xy})\mathrel{\mathop{\mathchar 58\relax}}=(\bm{d}_{xy}^{i+1}-\bm{g}_{xy}^{i+1}(\bm{z}_{xy}))-(\bm{d}_{xy}^{i}-\bm{g}_{xy}^{i}(\bm{z}_{xy}))\in\mathbb{R}^{n_{t}-1}.

(3.2)

Where $\bm{z}^{i}_{xy}\in\mathbb{R}$ and is the $i$ - $th$ component of the vector $\bm{z}_{xy}\in\mathbb{R}^{n_{t}}$ . This is done so the events associated to the $i$ - $th$ image are being optimized by the $i$ - $th$ $\bm{z}_{xy}$ component.

Notice that, $\bm{A}$ above is not a square matrix. For every $\lambda_{2}>0$ , the lower level problem is uniquely solvable and the solution is given by

\bm{u}_{xy}(\bm{z}_{xy})=(\bm{A}^{\top}\bm{A}+\lambda_{2}\bm{I})^{-1}\bm{A}^{\top}\bm{b}_{xy}(\bm{z}_{xy})=\bm{K}^{-1}\bm{A^{\top}}\bm{b}_{xy}(\bm{z}_{xy})

(3.3)

with

\bm{K}\mathrel{\mathop{\mathchar 58\relax}}=(\bm{A}^{\top}\bm{A}+\lambda_{2}\bm{I})\in\mathbb{R}^{n_{t}\times n_{t}}.

Substituting, this in the upper level problem, we obtain the following reduced problem

\min_{\bm{z}_{xy}\in Z_{ad}}\mathcal{J}(\bm{z}_{xy})\mathrel{\mathop{\mathchar 58\relax}}=J(\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}(\bm{z}_{xy}),\bm{z}_{xy})

(3.4)

which is now just a problem in $\bm{z}_{xy}$ . Notice, that the reduced problem (3.4) is equivalent to the full space formulation

	$\min_{\bm{z}_{xy}}J(\bm{u}_{xy},\bm{z}_{xy})$	(3.5a)
subject to constraints
	$\bm{K}\bm{u}_{xy}=\bm{A}^{\top}\bm{b}_{xy}(\bm{z}_{xy}).$	(3.5b)

The next result establishes existence of solution to the reduced problem (3.4) and equivalently to (3.5).

Theorem 3.1.

There exists a solution to the reduced problem (3.4).

Proof.

The proof follows from the standard Direct method of calculus of variations. We will omit the subscript $xy$ for notation simplicity. We begin by noticing that $\mathcal{J}(\cdot)$ is bounded below by 0. Therefore, there exists a minimizing sequence $\{\bm{z}_{n}\}_{n\in\mathbb{N}}$ such that

\lim_{n\rightarrow\infty}\mathcal{J}(\bm{z}_{n})=\inf_{\bm{z}}\mathcal{J}(\bm{z}).

It then follows that

\frac{\lambda_{1}}{2}\lim_{n\rightarrow\infty}\|\bm{z}_{n}\|^{2}\leq\lim_{n\rightarrow\infty}\mathcal{J}(\bm{z}_{n})=\inf_{\bm{z}}\mathcal{J}(\bm{z})\leq\mathcal{J}(\bm{0})=C<\infty,

where the constant $C$ is independent of $n\in\mathbb{N}$ . Thus, $\{\bm{z}_{n}\}_{n\in\mathbb{N}}$ is a bounded sequence. Therefore, it has a convergent subsequence, still denoted by $\{\bm{z}_{n}\}_{n\in\mathbb{N}}$ , such that $\bm{z}_{n}\rightarrow\bar{\bm{z}}$ . It them remains to show that $\bar{\bm{z}}$ is the minimizer. Since, $\bm{u}_{n}\mathrel{\mathop{\mathchar 58\relax}}={\bm{u}}(\bm{z}_{n})$ solves (3.3), therefore, we have that $\bm{u}_{n}\rightarrow{\bm{u}}(\bar{\bm{z}})$ and $\bm{g}(\bm{z}_{n})\rightarrow\bm{g}(\bm{\bar{z}})$ as $n\rightarrow\infty$ . Subsequently, from the definition of infimum and these convergence estimates, we readily obtain that

\lim_{n\rightarrow\infty}\mathcal{J}(\bm{z})\geq\liminf_{n\rightarrow\infty}\frac{1}{2}\|\bm{u}(\bm{z}_{n})+\bm{g}(\bm{z}_{n})-\bm{d}\|_{2}^{2}+\liminf_{n\rightarrow\infty}\frac{\lambda_{1}}{2}\|\bm{z}_{n}\|_{2}^{2}=\mathcal{J}(\bar{\bm{z}}).

Thus $\bar{\bm{z}}$ is the minimizer and the proof is complete. ∎

In order to develop a gradient based solver for (3.4), we next write down the expression of gradient of $\mathcal{J}$ .

Lemma 3.2.

The reduced objective $\mathcal{J}$ is continuously differentiable and the derivative is given by

\displaystyle\nabla\mathcal{J}(\bm{z})=\left\langle\bm{K}^{-1}\bm{A}^{\top}\bm{b}(\bm{z}_{xy})+\bm{g}(\bm{z}_{xy})-\bm{d}_{xy}\,,\,\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}^{\prime}(\bm{z}_{xy})+\bm{g}_{xy}^{\prime}(\bm{z_{xy}})\right\rangle+\lambda\bm{z}_{xy},

where $\bm{K}^{-1}\bm{A}^{\top}\in\mathbb{R}^{n_{t}\times(n_{t}-1)}$ ,

\displaystyle\bm{b}_{xy}^{\prime}(\bm{z}_{xy})=\begin{pmatrix}(\bm{g}_{xy}^{1})^{\prime}&-(\bm{g}_{xy}^{2})^{\prime}&0&\cdots&0&0&0\\ 0&(\bm{g}_{xy}^{2})^{\prime}&-(\bm{g}_{xy}^{3})^{\prime}&0&\cdots&0&0\\ 0&0&(\bm{g}_{xy}^{3})^{\prime}&-(\bm{g}_{xy}^{4})^{\prime}&0&\cdots&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 0&0&0&0&0&(\bm{g}_{xy}^{n_{t}-1})^{\prime}&-(\bm{g}_{xy}^{n_{t}})^{\prime}\end{pmatrix}\in\mathbb{R}^{(n_{t}-1)\times n_{t}},

and $\bm{g}_{xy}^{\prime}(\bm{z_{xy}})$ is a diagonal matrix given by

	$\displaystyle\bm{g}_{xy}^{\prime}(\bm{z_{xy}})=$		(3.6)
	$\displaystyle\begin{pmatrix}\frac{\int_{\mathcal{I}_{1}}\bm{E}^{1}_{xy}(t)\exp(\bm{z}^{1}_{xy}\cdot\bm{E}^{1}_{xy}(t))dt}{\int_{\mathcal{I}_{1}}\exp(\bm{z}^{1}_{xy}\cdot\bm{E}^{1}_{xy}(t))dt}&0&0&0&0\\ 0&\frac{\int_{\mathcal{I}_{2}}\bm{E}_{xy}^{2}(t)\exp(\bm{z}^{2}_{xy}\cdot\bm{E}_{xy}^{2}(t))dt}{\int_{\mathcal{I}_{2}}\exp(\bm{z}^{2}_{xy}\cdot\bm{E}_{xy}^{2}(t))dt}&0&0&0\\ 0&0&0&\ddots&0\\ 0&0&0&0&\frac{\int_{\mathcal{I}_{n_{t}}}\bm{E}_{xy}^{n_{t}}(t)\exp(\bm{z}^{n_{t}}_{xy}\cdot\bm{E}_{xy}^{n_{t}}(t))dt}{\int_{\mathcal{I}_{n_{t}}}\exp(\bm{z}^{n_{t}}_{xy}\cdot\bm{E}_{xy}^{n_{t}}(t))dt}\end{pmatrix}.$		(3.7)

Proof.

The proof follows by simple calculations. ∎

The next result establishes that $(\bm{g}^{i}_{xy}(\bm{z}_{xy}))^{\prime\prime}$ is bounded and strictly positive for each $i$ , with $\mathcal{I}_{i}$ denoting the $i$ - $th$ exposure time interval.

Theorem 3.3.

Let $\mathcal{I}_{i}$ denote the $i$ - $th$ exposure time interval with $i=1,\dots,n_{t}$ . If there exists $t\in\mathcal{I}_{i}$ , for all $i$ , with ${\bm{E}}^{i}_{xy}(t)$ being nonzero then $(\bm{g}^{i}_{xy}(\bm{z}_{xy}))^{\prime\prime}\in\mathbb{R}$ is bounded and strictly positive for all $i$ .

Proof.

First we will calculate $(\bm{g}^{i}_{xy}(\bm{z}_{xy}))^{\prime\prime}$ . Recall:

\displaystyle\bm{g}^{i}_{xy}(\bm{z}_{xy})=\log\bigg{(}\frac{1}{|\mathcal{I}_{i}|}\int_{\mathcal{I}_{i}}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))dt\bigg{)}\implies(\bm{g}^{i}_{xy}(\bm{z}_{xy}))^{\prime}=\frac{\int_{\mathcal{I}_{i}}\bm{E}^{i}_{xy}(t)\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))dt}{\int_{\mathcal{I}_{i}}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))dt},

yielding

\displaystyle\begin{aligned} (\bm{g}^{i}_{xy}(\bm{z}_{xy}))^{\prime\prime}=\frac{\boxed{\textrm{I}}-\bigg{(}\int_{\mathcal{I}_{i}}\bm{E}^{i}_{xy}(t)\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))dt\bigg{)}^{2}}{\big{[}\int_{\mathcal{I}_{i}}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))dt\big{]}^{2}},\end{aligned}

(3.8)

where $\boxed{\textrm{I}}=\bigg{(}\int_{\mathcal{I}_{i}}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))dt\bigg{)}\bigg{(}\int_{\mathcal{I}_{i}}\big{[}\bm{E}^{i}_{xy}(t)\big{]}^{2}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))dt\bigg{)}$ . Let $\bigg{(}\int_{\mathcal{I}_{i}}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))dt\bigg{)}=\mathrel{\mathop{\mathchar 58\relax}}\alpha$ . Then we have the following:

	$\displaystyle\boxed{\textrm{I}}$	$\displaystyle=\bigg{(}\int_{\mathcal{I}_{i}}\big{[}\bm{E}^{i}_{xy}(t)\big{]}^{2}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))\cdot\alpha\;dt\bigg{)}$
		$\displaystyle=\bigg{(}\int_{\mathcal{I}_{i}}\big{[}\bm{E}^{i}_{xy}(t)\big{]}^{2}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))\cdot\bigg{(}\int_{\mathcal{I}_{i}}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(s))ds\bigg{)}\;dt\bigg{)}$
		$\displaystyle=\bigg{(}\int_{\mathcal{I}_{i}}\int_{\mathcal{I}_{i}}\big{[}\bm{E}^{i}_{xy}(t)\big{]}^{2}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))\cdot\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(s))ds\;dt\bigg{)}.$

Using Fubini’s Theorem and following the procedure above for the two integral products in the numerator of (3.8) we get the following equalities:

$\displaystyle\eqref{eq:int}$	$\displaystyle=\frac{\boxed{\textrm{I}}-\bigg{(}\iint\limits_{\mathcal{I}_{i}\times\mathcal{I}_{i}}\bm{E}^{i}_{xy}(t)\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))\cdot\;\bm{E}^{i}_{xy}(s)\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(s))\;ds\;dt\bigg{)}}{\big{[}\int_{\mathcal{I}_{i}}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))dt\big{]}^{2}}$
	$\displaystyle=\frac{\bigg{(}\iint\limits_{\mathcal{I}_{i}\times\mathcal{I}_{i}}\big{[}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(s))\big{]}\bigg{(}\big{[}\bm{E}^{i}_{xy}(t)\big{]}^{2}-\bm{E}^{i}_{xy}(s)\bm{E}^{i}_{xy}(t)\bigg{)}\;ds\;dt\bigg{)}}{\big{[}\int_{\mathcal{I}_{i}}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))dt\big{]}^{2}}$
	$\displaystyle=\frac{\bigg{(}\iint\limits_{\mathcal{I}_{i}\times\mathcal{I}_{i}}\big{[}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(s))\big{]}\bigg{(}\frac{1}{2}\big{[}\bm{E}^{i}_{xy}(t)-\bm{E}^{i}_{xy}(s)\big{]}^{2}\bigg{)}\;ds\;dt\bigg{)}}{\big{[}\int_{\mathcal{I}_{1}}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))dt\big{]}^{2}}.$	(3.9)

We notice that Equation (3.9) is positive for $\bm{E}^{i}_{xy}(t)\neq\bm{E}^{i}_{xy}(s)$ (since all terms are positive) and the expression is $0$ for $t=s$ . By assumption there exists some $t\in\mathcal{I}_{i}$ such that $\bm{E}^{i}_{xy}(t)\neq 0$ which implies that an event has occurred at the pixel ${xy}$ during the interval $\mathcal{I}_{i}.$ Without loss of generality suppose a single event occurred at time $\bar{t}$ . Recall from (2.3) that $\bm{E}^{i}_{xy}(t,t_{i})=\int_{t_{i}}^{t}\bm{e}_{xy}(r)\;dr$ . Then $0=\bm{E}^{i}_{xy}(t_{i},t_{i})\neq\bm{E}^{i}_{xy}(\bar{t},t_{i})$ which implies $\bm{E}^{i}_{xy}(t_{i})\neq\bm{E}^{i}_{xy}(\bar{t})$ . Which further impiles that $\eqref{eq:pos}>0$ . Thus we conclude $(\bm{g}^{i}_{xy}(\bm{z}_{xy}))^{\prime\prime}>0$ . Next, using (3.9), we obtain the following upper bound:

	$\displaystyle(\bm{g}^{i}_{xy}(\bm{z}_{xy}))^{\prime\prime}$	$\displaystyle\leq\frac{\bigg{(}\iint\limits_{\mathcal{I}_{i}\times\mathcal{I}_{i}}\big{[}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(s))\big{]}\max\limits_{t,s\in\mathcal{I}_{i}}\bigg{\{}\frac{1}{2}\big{[}\bm{E}^{i}_{xy}(t)-\bm{E}^{i}_{xy}(s)\big{]}^{2}\bigg{\}}\;ds\;dt\bigg{)}}{\big{[}\int_{\mathcal{I}_{i}}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}_{xy}(t))dt\big{]}^{2}}$
		$\displaystyle=\max\limits_{t,s\in\mathcal{I}_{i}}\bigg{\{}\frac{1}{2}\big{[}\bm{E}^{i}_{xy}(t)-\bm{E}^{i}_{xy}(s)\big{]}^{2}\bigg{\}}\frac{\bigg{(}\iint\limits_{\mathcal{I}_{i}\times\mathcal{I}_{i}}\big{[}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(s))\big{]}\;ds\;dt\bigg{)}}{\big{[}\int_{\mathcal{I}_{i}}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))dt\big{]}^{2}}$
		$\displaystyle=\max\limits_{t,s\in\mathcal{I}_{i}}\bigg{\{}\frac{1}{2}\big{[}\bm{E}^{i}_{xy}(t)-\bm{E}^{i}_{xy}(s)\big{]}^{2}\bigg{\}}.$

Hence

0<(\bm{g}^{i}_{xy}(\bm{z}_{xy}))^{\prime\prime}\leq\max\limits_{t,s\in\mathcal{I}_{i}}\bigg{\{}\frac{1}{2}\big{[}\bm{E}^{i}_{xy}(t)-\bm{E}^{i}_{xy}(s)\big{]}^{2}\bigg{\}}.

We now conclude that $(\bm{g}^{i}_{xy}(\bm{z}_{xy}))^{\prime\prime}$ is bounded and strictly positive. ∎

In order to show the local convexity of $\mathcal{J}$ , we study the structure of the Hessian of the first term in the definition of $\mathcal{J}$ .

Lemma 3.4.

The matrix $\big{\langle}\text{{Hess}}(\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}(\bm{z}_{xy}))+\text{{Hess}}(\bm{g}_{xy}(\bm{z}_{xy})),\bm{u}_{xy}(\bm{z}_{xy})+\bm{g}_{xy}(\bm{z}_{xy})-\bm{d}_{xy}\big{\rangle}$ is a diagonal matrix.

Proof.

Without loss of generality, suppose that $n_{t}=3$ . Recall by our definition:

\displaystyle\bm{g}_{xy}(\bm{z}_{xy})

\displaystyle=\begin{pmatrix}\log(\frac{1}{|\mathcal{I}_{1}|}\int_{\mathcal{I}_{1}}\exp(\bm{z}^{1}_{xy}\cdot\bm{E}^{1}_{xy}(t))dt)\\ \log(\frac{1}{|\mathcal{I}_{2}|}\int_{\mathcal{I}_{2}}\exp(\bm{z}^{2}_{xy}\cdot\bm{E}^{2}_{xy}(t))dt)\\ \log(\frac{1}{|\mathcal{I}_{3}|}\int_{\mathcal{I}_{3}}\exp(\bm{z}^{3}_{xy}\cdot\bm{E}^{3}_{xy}(t))dt)\\ \end{pmatrix}.

Then, recall the expression of $\bm{g}^{\prime}_{xy}(\bm{z}_{xy})$ from (3.7). Now to calculate the Hessian of $\bm{g}$ , we need to take the derivative with respect to each variable once again. This will result in $\text{{Hess}}(\bm{g}_{xy}(\bm{z}_{xy}))\in\mathbb{R}^{n_{t}\times n_{t}\times n_{t}}$ and will be of the form:

\displaystyle\text{{Hess}}(\bm{g}_{xy}(\bm{z}_{xy}))=\begin{Bmatrix}\begin{pmatrix}(\bm{g}^{1}_{xy}(\bm{z}_{xy}))^{\prime\prime}&0&0\\ 0&0&0\\ 0&0&0\end{pmatrix},\;\begin{pmatrix}0&0&0\\ 0&(\bm{g}^{2}_{xy}(\bm{z}_{xy}))^{\prime\prime}&0\\ 0&0&0\end{pmatrix},\;\begin{pmatrix}0&0&0\\ 0&0&0\\ 0&0&(\bm{g}_{xy}^{3}(\bm{z}_{xy}))^{\prime\prime}\end{pmatrix}\end{Bmatrix},

where $(\bm{g}^{i}(\bm{z}_{xy}))^{\prime\prime}$ is defined in Theorem 3.3 . Next we consider $\text{{Hess}}(\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}(\bm{z}_{xy}))$ . Observe:

\displaystyle\text{{Hess}}(\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}(\bm{z}_{xy}))

\displaystyle=\bm{K}^{-1}\bm{A}^{\top}\cdot\text{{Hess}}(\bm{b}_{xy}(\bm{z}_{xy}))

where

\displaystyle\text{{Hess}}(\bm{b}_{xy}(\bm{z}_{xy}))=\begin{Bmatrix}\begin{pmatrix}(\bm{g}^{1}_{xy}(\bm{z}_{xy}))^{\prime\prime}&0&0\\ 0&0&0\end{pmatrix},\begin{pmatrix}0&-(\bm{g}^{2}_{xy}(\bm{z}_{xy}))^{\prime\prime}&0\\ 0&(\bm{g}^{2}_{xy}(\bm{z}_{xy}))^{\prime\prime}&0\end{pmatrix},\begin{pmatrix}(0&0&0\\ 0&0&\bm{g}^{3}_{xy}(\bm{z}_{xy}))^{\prime\prime}\end{pmatrix}\end{Bmatrix},

\displaystyle\text{{Hess}}(\bm{b}_{xy}(\bm{z}_{xy}))\in\mathbb{R}^{{n_{t}-1}\times n_{t}\times n_{t}},

and

\displaystyle\bm{K}^{-1}\bm{A}^{\top}\mathrel{\mathop{\mathchar 58\relax}}=\begin{pmatrix}|&|\\ \bm{\kappa}_{1}&\bm{\kappa}_{2}\\ |&|\end{pmatrix}\in\mathbb{R}^{n_{t}\times(n_{t}-1)}.

Here $\bm{\kappa}_{i}$ , $i=1...n_{t}$ represent the columns of $\bm{K}^{-1}\bm{A}^{\top}$ . Then we have the following:

	$\displaystyle\text{{Hess}}(\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}(\bm{z}_{xy}))$
	$\displaystyle=\begin{Bmatrix}(\bm{g}^{1}_{xy}(\bm{z}_{xy}))^{\prime\prime}\cdot\begin{pmatrix}\|&0&0\\ \bm{\kappa}_{1}&0&0\\ \|&0&0\end{pmatrix},(\bm{g}^{2}_{xy}(\bm{z}_{xy}))^{\prime\prime}\cdot\begin{pmatrix}0&\|&0\\ 0&\bm{\kappa}_{2}-\bm{\kappa}_{1}&0\\ 0&\|&0\end{pmatrix},(\bm{g}^{3}_{xy}(\bm{z}_{xy}))^{\prime\prime}\cdot\begin{pmatrix}0&0&\|\\ 0&0&\bm{\kappa}_{3}\\ 0&0&\|\end{pmatrix}\end{Bmatrix},$

with $\text{{Hess}}(\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}(\bm{z}_{xy})\in\mathbb{R}^{n_{t}\times n_{t}\times n_{t}}$ . Next we can add $\text{{Hess}}(\bm{g}_{xy}(\bm{z}_{xy}))$ and $\text{{Hess}}(\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}(\bm{z}_{xy}))$ to get:

	$\displaystyle\text{{Hess}}(\bm{g}_{xy}(\bm{z}_{xy}))+\text{{Hess}}(\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}(\bm{z}_{xy}))$
	$\displaystyle=\begin{Bmatrix}(\bm{g}^{1}_{xy}(\bm{z}_{xy}))^{\prime\prime}\cdot\begin{pmatrix}\|&0&0\\ \bm{\gamma}_{1}&0&0\\ \|&0&0\end{pmatrix},(\bm{g}^{2}_{xy}(\bm{z}_{xy}))^{\prime\prime}\cdot\begin{pmatrix}0&\|&0\\ 0&\bm{\gamma}_{2}&0\\ 0&\|&0\end{pmatrix},(\bm{g}^{3}_{xy}(\bm{z}_{xy}))^{\prime\prime}\cdot\begin{pmatrix}0&0&\|\\ 0&0&\bm{\gamma}_{3}\\ 0&0&\|\end{pmatrix}\end{Bmatrix},$

where $\bm{\gamma}_{i}\in\mathbb{R}^{n_{t}},i=1,\dots,n_{t}$ . Lastly let $\bm{u}_{xy}(\bm{z}_{xy})+\bm{g}_{xy}(\bm{z}_{xy})-\bm{d}_{xy}=\bm{\alpha}_{xy}\in\mathbb{R}^{n_{t}}$ . Then

	$\displaystyle\big{\langle}\text{{Hess}}(\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}(\bm{z}_{xy}))+\text{{Hess}}(\bm{g}_{xy}(\bm{z}_{xy})),\bm{u}_{xy}(\bm{z}_{xy})+\bm{g}_{xy}(\bm{z}_{xy})-\bm{d}_{xy}\big{\rangle}$
	$\displaystyle=\begin{Bmatrix}(\bm{g}^{1}_{xy}(\bm{z}_{xy}))^{\prime\prime}\cdot\begin{pmatrix}\bm{\gamma}_{1}^{\top}\bm{\alpha}_{xy}\\ 0\\ 0\end{pmatrix},(\bm{g}^{2}_{xy}(\bm{z}_{xy}))^{\prime\prime}\cdot\begin{pmatrix}0\\ \bm{\gamma}_{2}^{\top}\bm{\alpha}_{xy}\\ 0\end{pmatrix},(\bm{g}^{3}_{xy}(\bm{z}_{xy}))^{\prime\prime}\cdot\begin{pmatrix}0\\ 0\\ \bm{\gamma}_{3}^{\top}\bm{\alpha}_{xy}\end{pmatrix}\end{Bmatrix}$
	$\displaystyle\cong\begin{pmatrix}(\bm{g}^{1}_{xy}(\bm{z}_{xy}))^{\prime\prime}\cdot\bm{\gamma}_{1}^{\top}\bm{\alpha}_{xy}&0&0\\ 0&(\bm{g}^{2}_{xy}(\bm{z}_{xy}))^{\prime\prime}\cdot\bm{\gamma}_{2}^{\top}\bm{\alpha}_{xy}&0\\ 0&0&(\bm{g}^{3}_{xy}(\bm{z}_{xy}))^{\prime\prime}\cdot\bm{\gamma}_{3}^{\top}\bm{\alpha}_{xy}\end{pmatrix}.$		(3.9)

Here $\cong$ denotes the mode 1 unfolding of the $\big{\langle}\text{{Hess}}(\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}(\bm{z}_{xy}))+\text{{Hess}}(\bm{g}_{xy}(\bm{z}_{xy})),\bm{u}_{xy}(\bm{z}_{xy})+\bm{g}_{xy}(\bm{z}_{xy})-\bm{d}_{xy}\big{\rangle}$ tensor. Further details of tensor matrix multiplication can be found in [11]. The proof is complete. ∎

Next, we establish that $\mathcal{J}$ is strictly locally convex on any finite interval when $\lambda_{1}$ is chosen appropriately.

Theorem 3.5.

Let $\mathcal{I}_{i}$ denotes the $i$ - $th$ exposure time interval with $i=1,\dots,n_{t}$ . If there exists $t\in\mathcal{I}_{i}$ , for all $i$ , with ${\bm{E}}^{i}_{xy}(t)$ being nonzero then for any closed ball $\Omega=\overline{\mathcal{B}_{\delta}(0)}$ centered at 0 and radius $\delta$ , there exists $\lambda_{1}>0$ such that for $\bm{z}_{xy}\in\Omega$ , $\mathcal{J}$ is strictly convex.

Proof.

First we calculate Hess( $\mathcal{J}$ ) as:

	$\displaystyle\text{Hess}(\mathcal{J})$	$\displaystyle=\big{\langle}\;\bm{u}_{xy}^{\prime}(\bm{z}_{xy})+\bm{g}_{xy}^{\prime}(\bm{z}_{xy})),\bm{u}_{xy}^{\prime}(\bm{z}_{xy})+\bm{g}_{xy}^{\prime}(\bm{z}_{xy}))\;\big{\rangle}$
		$\displaystyle\quad+\big{\langle}\text{{Hess}}(\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}(\bm{z}_{xy}))+\textbf{Hess}(\bm{g}_{xy}(\bm{z}_{xy})),\bm{u}_{xy}(\bm{z}_{xy})+\bm{g}_{xy}(\bm{z}_{xy})-\bm{d}_{xy}\big{\rangle}+\lambda_{1}\bm{I}.$

We notice that the term $\big{\langle}\;\bm{u}_{xy}^{\prime}(\bm{z}_{xy})+\bm{g}_{xy}^{\prime}(\bm{z}_{xy})),\bm{u}_{xy}^{\prime}(\bm{z}_{xy})+\bm{g}_{xy}^{\prime}(\bm{z}_{xy}))\;\big{\rangle}$ is of them form $\bm{M}^{\top}\bm{M}$ , which implies that it is symmetric positive semi-definite. Next, from Lemma 3.4 we recall that the term $\big{\langle}\text{{Hess}}(\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}(\bm{z}_{xy}))+\textbf{Hess}(\bm{g}_{xy}(\bm{z}_{xy})),\bm{u}_{xy}(\bm{z}_{xy})+\bm{g}_{xy}(\bm{z}_{xy})-\bm{d}_{xy}\big{\rangle}$ is a diagonal matrix. In the case where the diagonal entries (eigenvalues) are greater than 0, then we are done. Suppose some eigenvalue is less than 0 and let $\bm{z}_{xy}\in\Omega\mathrel{\mathop{\mathchar 58\relax}}=\overline{\mathcal{B}_{\delta}(0)}$ which implies $\|\bm{z}_{xy}\|_{\infty}\leq\delta$ . Given $\delta>0$ , choose $\lambda_{1}>0$ such that

0<\delta<\frac{\lambda_{1}-(\|\bm{\gamma}\|_{\infty}\cdot M_{1}\cdot\|\bm{d}_{xy}\|_{\infty})(2\|\bm{K}^{-1}\bm{A}^{\top}\|_{\infty}+1)}{(\|\bm{\gamma}\|_{\infty}\cdot M_{1}\cdot M_{2})(2\|\bm{K}^{-1}\bm{A}^{\top}\|_{\infty})+1)}.

Where $M_{2}=\max\limits_{i}\bigg{\{}\big{|}\min\limits_{t\in\mathcal{I}_{i}}\{\bm{E}^{i}_{xy}(t)\}\big{|},\max\limits_{t\in\mathcal{I}_{i}}\big{\{}\bm{E}^{i}_{xy}(t)\big{\}}\bigg{\}}$ and $M_{1}=\max\limits_{i}\biggl{\{}\max\limits_{t,s\in\mathcal{I}_{i}}\bigg{\{}\frac{1}{2}\big{[}\bm{E}^{i}_{xy}(t)-\bm{E}^{i}_{xy}(s)\big{]}^{2}\bigg{\}}\biggr{\}}.$ Recall from Theorem 3.3 that $0<(\bm{g}^{i}_{xy}(\bm{z}_{xy}))^{\prime\prime}\leq\max\limits_{t,s\in\mathcal{I}_{i}}\bigg{\{}\frac{1}{2}\big{[}\bm{E}^{i}_{xy}(t)-\bm{E}^{i}_{xy}(s)\big{]}^{2}\bigg{\}}$ , which implies $(\bm{g}^{i}_{xy}(\bm{z}_{xy}))^{\prime\prime}\leq M_{1}$ for any $i$ . Then using 3.9 we have the following:

	$\displaystyle\\|\big{\langle}\text{{Hess}}(\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}(\bm{z}_{xy}))+\textbf{Hess}(\bm{g}_{xy}(\bm{z}_{xy})),\bm{u}_{xy}(\bm{z}_{xy})+\bm{g}_{xy}(\bm{z}_{xy})-\bm{d}_{xy}\big{\rangle}\\|_{\infty}$
	$\displaystyle\leq\begin{Vmatrix}\begin{pmatrix}(\bm{g}^{1}_{xy}(\bm{z}_{xy}))^{\prime\prime}\cdot\bm{\gamma}_{1}^{\top}\bm{\alpha}_{xy}&0&0&0\\ 0&(\bm{g}^{2}_{xy}(\bm{z}_{xy}))^{\prime\prime}\cdot\bm{\gamma}_{2}^{\top}\bm{\alpha}_{xy}&0&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\dots&(\bm{g}^{n_{t}}_{xy}(\bm{z}_{xy}))^{\prime\prime}\cdot\bm{\gamma}_{n_{t}}^{\top}\bm{\alpha}_{xy}\end{pmatrix}\\ \end{Vmatrix}_{\infty}$
	$\displaystyle\leq M_{1}\cdot\\|\bm{\gamma}\\|_{\infty}\bigg{(}\\|\bm{u}_{xy}(\bm{z}_{xy})\\|_{\infty}+\\|\bm{g}_{xy}(\bm{z}_{xy})\\|_{\infty}+\\|\bm{d}_{xy}\\|_{\infty}\bigg{)}$
	$\displaystyle\leq M_{1}\cdot\\|\bm{\gamma}\\|_{\infty}\bigg{(}\\|\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}\\|_{\infty}+M_{2}\\|\bm{z}_{xy}\\|_{\infty}+\\|\bm{d}_{xy}\\|_{\infty}\bigg{)}$
	$\displaystyle\leq M_{1}\cdot\\|\bm{\gamma}\\|_{\infty}\bigg{(}\\|\bm{K}^{-1}\bm{A}^{\top}\\|_{\infty}\\|\bm{b}_{xy}\\|_{\infty}+M_{2}\\|\bm{z}_{xy}\\|_{\infty}+\\|\bm{d}_{xy}\\|_{\infty}\bigg{)}$
	$\displaystyle\leq M_{1}\cdot\\|\bm{\gamma}\\|_{\infty}\bigg{(}\\|\bm{K}^{-1}\bm{A}^{\top}\\|_{\infty}\;(2\\|\bm{d}_{xy}\\|_{\infty}+2M_{2}\\|\bm{z}_{xy}\\|_{\infty})+M_{2}\\|\bm{z}_{xy}\\|_{\infty}+\\|\bm{d}_{xy}\\|_{\infty}\bigg{)}$
	$\displaystyle<\lambda_{1}.$

Here

\bm{\gamma}\mathrel{\mathop{\mathchar 58\relax}}=\begin{pmatrix}|&|&\dots&|\\ \bm{\gamma}_{1}&\bm{\gamma}_{2}&\dots&\bm{\gamma}_{n_{t}}\\ |&|&\dots&|\end{pmatrix}^{\top},

$\bm{\gamma}_{i}$ and $\bm{\alpha}_{xy}$ are defined in Equation 3.9 . Notice that $M_{2}$ is derived from bounds on $\bm{g}_{xy}(\bm{z}_{xy})$ as follows:

\displaystyle\bm{g}^{i}_{xy}(\bm{z}_{xy})=\log\bigg{(}\frac{1}{|\mathcal{I}_{i}|}\int_{\mathcal{I}_{i}}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))dt\bigg{)}\leq\bm{z}^{i}_{xy}\cdot\max\big{\{}\bm{E}^{i}_{xy}(t)\big{\}},\text{ where }\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t)\geq 0

and

	$\displaystyle\bm{g}^{i}_{xy}(\bm{z}_{xy})$	$\displaystyle=\log\bigg{(}\frac{1}{\|\mathcal{I}_{i}\|}\int_{\mathcal{I}_{i}}\exp(\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t))dt\bigg{)}\geq\bm{z}^{i}_{xy}\cdot\min\big{\{}\bm{E}^{i}_{xy}(t)\big{\}},\text{ where }\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(t)\leq 0$
		$\displaystyle\implies\|\bm{g}^{i}_{xy}(\bm{z}_{xy})\|\leq\|\bm{z}^{i}_{xy}\cdot M_{2}\|.$

Recall $\big{\langle}\text{{Hess}}(\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}(\bm{z}_{xy}))+\textbf{Hess}(\bm{g}_{xy}(\bm{z}_{xy})),\bm{u}_{xy}(\bm{z}_{xy})+\bm{g}_{xy}(\bm{z}_{xy})-\bm{d}_{xy}\big{\rangle}$ is diagonal from 3.9 . Denote the diagonal entries as $\sigma_{i}$ , $i=1,\dots,n_{t}$ . Without loss of generality let $|\sigma_{1}|\geq|\sigma_{i}|$ , $i=2,\dots,n_{t}$ . Then we have the following:

|\sigma_{i}|\leq|{\sigma_{1}}|=\|\big{\langle}\text{{Hess}}(\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}(\bm{z}_{xy}))+\textbf{Hess}(\bm{g}_{xy}(\bm{z}_{xy})),\bm{u}_{xy}(\bm{z}_{xy})+\bm{g}_{xy}(\bm{z}_{xy})-\bm{d}_{xy}\big{\rangle}\|_{\infty}<\lambda_{1}.

Thus $\sigma_{i}<\lambda_{1}$ for all $i$ . Furthermore $\sigma_{i}+\lambda_{1}>0$ for all $i$ . We notice

\big{\langle}\text{{Hess}}(\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}(\bm{z}_{xy}))+\textbf{Hess}(\bm{g}_{xy}(\bm{z}_{xy})),\bm{u}_{xy}(\bm{z}_{xy})+\bm{g}_{xy}(\bm{z}_{xy})-\bm{d}_{xy}\big{\rangle}+\lambda_{1}\bm{I}

is symmetric positive definite since it is a diagonal matrix with entries: $\sigma_{i}+\lambda_{1}>0$ . Hence Hess( $\mathcal{J}$ ) is symmetric positive definite due to it being the summation of a symmetric positive definite matrix and a symmetric positive semi-definite matrix. We conclude that for any convex set $\Omega$ there exists $\lambda_{1}$ such that $\mathcal{J}$ is strictly convex over $\Omega$ . ∎

We further note that from Theorem 3.5 once an interval is chosen a lower bound for $\lambda_{1}$ can be calculated by:

\lambda_{1}>\|\bm{\gamma}\|_{\infty}\cdot M_{1}\bigg{(}\|\bm{K}^{-1}\bm{A}^{\top}\|_{\infty}\;(2\|\bm{d}_{xy}\|_{\infty}+2M_{2}\|\bm{z}_{xy}\|_{\infty})+M_{2}\|\bm{z}_{xy}\|_{\infty}+\|\bm{d}_{xy}\|_{\infty}\bigg{)}.

4. Experimental Setup

This section focuses on how to prepare the dataset to be used in our computations.

Standardizing $\bm{B^{i}}$

Depending on the event camera used, the entries of the matrix $\bm{B^{i}}$ from (2.4), where $i$ identifies a unique standard camera image in the sequence being considered, may take on various ranges to include $[0,255]$ or $[-255,255]$ . We standardize $\bm{B^{i}}$ using MATLAB function mat2gray. This function will map the values of $\bm{B^{i}}$ to the range $[0,1]$ on a log scale. That being said since we define the matrix $\bm{d}^{i}=\log(\bm{B^{i}})$ we have to be careful not to allow $\bm{B^{i}}$ to have any zero values. We can easily protect against this by inserting a manual range into the mat2gray function of $[\min{\bm{B^{i}}}-\epsilon,\max{\bm{B^{i}}}+\epsilon]$ , for a sufficiently small $\epsilon$ . In our computations we set $\epsilon=1e-3$ . We note $\min\bm{B}^{i}$ and $\max\bm{B}^{i}$ are the minimum and maximum values over all values in the matrix $\bm{B}^{i}$ , respectively.

Building the data cube

In order to view the problem in three dimensions, we create a data cube as shown in Figure $\eqref{f:ex1}$ . As stated in (2.1) we define this cube by:

{\bm{EC}}_{xy}^{s_{j}}=p_{x,y}^{s_{j}}.

Since the event data is generated on a nearly continuous basis, the data cube will be of size $\mathbb{R}^{\nu\times{n_{x}}\times{n_{y}}}$ where $\nu\mathrel{\mathop{\mathchar 58\relax}}=\sum\limits_{i}n_{i}$ is the total number of events used in the reconstruction over all images. Note that $n_{i}$ is explained in further detail in 2.1 . On average our image reconstructions use approximately $\nu\approx 25,000$ events. It is impractical numerically to use this data cube due to its size. This issue is resolved by reducing the size of the data cube via time compression. Choose some $r,k\in\mathbb{N}$ such that $\nu\approx r\cdot k$ and define a new data cube as follows:

\displaystyle\bm{dc}^{\pi}=\sum_{\omega=1+(\pi-1)\cdot k}^{\pi\cdot k}\bm{EC}^{\omega}

Then we have $\bm{dc}\in\mathbb{R}^{r\times{n_{x}}\times{n_{y}}}$ where $r\approx\frac{\nu}{k}$ , and $\pi=1,\dots,r$ . For our examples we have chosen $k\approx 200$ .

Calculate $\bm{g}(\bm{z})$

Once we have built the data cube then we can calculate the function $\bm{g}(\bm{z})$ which is defined componentwise as:

\bm{g}_{xy}^{i}(\bm{z}_{xy})=\log\bigg{(}\frac{1}{|\mathcal{I}_{i}|}\int_{\mathcal{I}_{i}}\exp\big{(}\bm{z}^{i}_{xy}\cdot\bm{E}^{i}_{xy}(\tau,t_{i})\big{)}dt\bigg{)},\quad i=1,\dots,n_{t}.

Numerically we approximate this using the trapezoid rule as:

\bm{g}_{xy}^{i}(\bm{z}_{xy})\approx\log\bigg{(}\frac{1}{|\mathcal{I}_{i}|}\sum_{l=2}^{r}\big{(}\frac{\triangle\tau_{l-1}}{2}\big{)}\big{(}\exp\big{(}\bm{z}^{i}_{xy}\cdot\bm{E}^{i}(\tau_{l-1},t_{i})\big{)}+\exp\big{(}\bm{z}^{i}_{xy}\cdot\bm{E}^{i}t_{xy}(\tau_{l},t_{i})\big{)}\big{)}\bigg{)},

where $t_{i}$ is defined in (2.3) and $\Delta\tau_{l-1}=\tau_{l}-\tau_{l-1}$ is the step length. We evaluate the function $\bm{b}(\bm{z})$ as well as the derivatives of $\bm{b}(\bm{z})$ and $\bm{g}(\bm{z})$ in a similar way.

Omitting unnecessary calculations for $\bm{z}$

We recall that $\bm{z}$ must be calculated on a per pixel basis. To decrease computation time we only perform an optimization for $\bm{z}_{xy}$ over the pixels for which event data has been captured.

Image Reconstruction

As shown in Figure 5 our method is capable of generating two image representations. Given some value ${\bm{z}}$ the $\bm{u}$ representation is given as: $\bm{u}({\bm{z}})$ defined in $\eqref{eq:u}$ . The $\bm{u}$ representation reconstructs only the dynamical part of the image. This is due to our definition of $\bm{b}$ in (3.2). In order to reconstruct the full image including the portion of the image with no dynamics, we use

\bm{v}_{xy}^{i}=\bm{d}_{xy}^{i}-\bm{g}_{xy}^{i}({\bm{z}}),

a variation of the definition given in (2.5). In general we will refer to the latter formulation as the reconstruction, and will specifically denote the $\bm{u}$ representation when presented.

5. Numerical Examples

Next, we present a series of examples which establishes that the proposed approach could be beneficial in practice.

5.1. Benchmark Problem: Unit Bump

The first example is a synthetic case where we consider a blurred unit bump, see Figure 4 . In order to perform this analysis, we require a baseline image for comparison, a series of consecutive images as input for the ESIM generator [10] and our proposed method, a blurry image, and the corresponding event data. To begin, we will construct a series of consecutive images by first creating a white circle on a black background (not shown). Using this initial image, we move the white circle two pixels to the right and two pixels down eight times. This set of nine images completes our series of images. We now generate a tenth image, the blurry image (middle in Figure 4), by averaging the sequence of nine images. The fifth image of the set is designated as our baseline image (left in Figure 4). Lastly, to generate the associated event data, we use the ESIM generator from [10]. The ESIM generator takes a series of images as input, and will output event data. Recall that our model requires a series of images as well as event data in order to generate a reconstructed image. We set $\bm{B}^{1}$ to be the fourth image of the nine image set, $\bm{B}^{2}$ to be the blurry image, and $\bm{B}^{3}$ to be the sixth image. The result of our method (right in Figure 4) is compared to the baseline image. As shown, we obtain a high quality reconstruction with SSIM of 0.96 and PSNR of 27.3. This example serves a benchmark for us to apply our approach on the event based dataset.

Figure 4. From left to right we have our baseline image, a simulated blurry image, and finally our reconstruction. Events for this reconstruction were produced using the ESIM generator from [10]. Using parameters

\lambda_{1}=1

and

\lambda_{2}=1e-3

, we achieve an SSIM value of 0.96 and a PSNR value of 27.3.

5.2. Multiple Event Based Dataset Examples

Figures 5-7 show the application of our method to various examples. In each of these examples three standard camera images are used to define the $\bm{B}$ variable, while the $\bm{g}$ and $\bm{b}$ variables are calculated using the associated event data. In each of these examples the $\bm{B}^{2}$ image was chosen to be a blurry image. We will refer to the $\bm{B}^{2}$ image as the blurry image for the remainder of this section.

In Figure 5 for the chosen three image sequence $\approx 39,000$ events were recorded. The blurry image is shown (left) while our reconstruction is shown (right). Our reconstruction deblurs the fence-line as well as captures the grass and tree texture that is not apparent in the blurred image. The blurred image (left) in Figure 6 shows a person swinging a pillow. Approximately 26,000 events were captured for this example. Our model has successfully reconstructed the boundary of the pillow as well as the persons hand and thumb. Figure 7 is an example of a person jumping and the three image sequence used has $\approx$ 16000 associated events. Notice in the blurred image (left) the persons arm is barely visible. In contrast our model is able to reconstruct the persons arm (right).

5.3. Event Based Problem: Shadowing Effects

As noted in Section 2.3 one of the limitations to the mEDI model is the appearance of shadows in a high contrasting environment. Examples of this shadowing effect can be seen in [13] and Figures 8 and 9. As shown in Figures 8 and 9 our model is able to substantially reduce the shadowing effects of the EDI and mEDI models.

Figure 10 shows the behavior of five randomly chosen pixels. We observe the typical convergence behavior of a Newton’s method.

Conclusion, Limitation, and Future Work

In this study, we introduce a novel approach to image deblurring by combining event data and multiple standard camera images within a variational framework. Our method demonstrates the existence of solutions for a locally convex problem and we apply a second order Newton solver to several examples, showcasing the effectiveness of our approach.

Though our method does offer a systematic way to construct de-blurred images across time there are limitations that still exist in the model.

(a)

Our method requires a pixel by pixel optimization and depending on the camera resolution can run into scalability issues. This issue would be most evident with images where dynamics occur at a majority of pixels.
(b)

Manual tuning may be required of certain parameters as described in Section 4 .

To mitigate the scalability risk, we propose to filter out pixels that have significantly lower event counts compared to other pixels. In future work, we also plan to test the model’s robustness in a domain shift scenario, such as a camera moving in and out of water. Moreover, we aim to develop a parameter learning framework to optimize the regularization parameters for better performance.

Acknowledgement

The authors are grateful to Dr. Patrick O’Neil and Dr. Diego Torrejon (BlackSky) and Dr. Noor Qadri (US Naval Research Lab, Washington DC), for bringing the neuromorphic imaging topic to their attention and for several fruitful discussions.

References

[1] Mariana SC Almeida and Luis B Almeida. Blind and semi-blind deblurring of natural images. IEEE Transactions on Image Processing, 19(1):36–52, 2009.
[2] Harbir Antil, Zichao Wendy Di, and Ratna Khatri. Bilevel optimization, deep learning and fractional laplacian regularization with applications in tomography. Inverse Problems, 36(6):064001, May 2020.
[3] Simon Arridge, Peter Maass, Ozan Öktem, and Carola-Bibiane Schönlieb. Solving inverse problems using data-driven models. Acta Numerica, 28:1–174, 2019.
[4] Mario Bertero, Patrizia Boccacci, and Christine De Mol. Introduction to inverse problems in imaging. CRC press, 2021.
[5] Alessandro Buccini, Marco Donatelli, and Ronny Ramlau. A semiblind regularization algorithm for inverse problems with application to image deblurring. SIAM Journal on Scientific Computing, 40(1):A452–A483, 2018.
[6] Leon Bungert and Matthias J Ehrhardt. Robust image reconstruction with misaligned structural information. IEEE Access, 8:222944–222955, 2020.
[7] Julianne Chung, Matthias Chung, and Dianne P O’Leary. Designing optimal spectral filters for inverse problems. SIAM Journal on Scientific Computing, 33(6):3132–3152, 2011.
[8] Juan Carlos De los Reyes, Carola-Bibiane Schönlieb, and Tuomo Valkonen. Bilevel parameter learning for higher-order total variation regularisation models. J. Math. Imaging Vision, 57(1):1–25, 2017.
[9] Guillermo Gallego, Tobi Delbrück, Garrick Orchard, Chiara Bartolozzi, Brian Taba, Andrea Censi, Stefan Leutenegger, Andrew J. Davison, Jörg Conradt, Kostas Daniilidis, and Davide Scaramuzza. Event-based vision: A survey. IEEE Transactions on Pattern Analysis and Machine Intelligence, 44(1):154–180, 2022.
[10] Daniel Gehrig, Henri Rebecq, Guillermo Gallego, and Davide Scaramuzza. Asynchronous, photometric feature tracking using events and frames. In Proceedings of the European Conference on Computer Vision (ECCV), pages 750–765, 2018.
[11] Tamara G. Kolda and Brett W. Bader. Tensor decompositions and applications. SIAM Review, 51(3):455–500, 2009.
[12] Renaud Morin, Stéphanie Bidon, Adrian Basarab, and Denis Kouamé. Semi-blind deconvolution for resolution enhancement in ultrasound imaging. In 2013 IEEE International Conference on Image Processing, pages 1413–1417. IEEE, 2013.
[13] Liyuan Pan, Cedric Scheerlinck, Xin Yu, Richard Hartley, Miaomiao Liu, and Yuchao Dai. Bringing a blurry frame alive at high frame-rate with an event camera. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 6820–6829, 2019.
[14] Se Un Park, Nicolas Dobigeon, and Alfred O Hero. Semi-blind sparse image reconstruction with application to mrfm. IEEE Transactions on Image Processing, 21(9):3838–3849, 2012.
[15] Pooja Satish, Mallikarjunaswamy Srikantaswamy, and Nataraj Kanathur Ramaswamy. A comprehensive review of blind deconvolution techniques for image deblurring. Traitement du Signal, 37(3), 2020.
[16] Cedric Scheerlinck, Nick Barnes, and Robert E. Mahony. Continuous-time intensity estimation using event cameras. CoRR, abs/1811.00386, 2018.
[17] Shu Tang, Weiguo Gong, Weihong Li, and Wenzhi Wang. Non-blind image deblurring method by local and nonlocal total variation models. Signal processing, 94:339–349, 2014.
[18] Shipeng Xie, Xinyu Zheng, Wen-Ze Shao, Yu-Dong Zhang, Tianxiang Lv, and Haibo Li. Non-blind image deblurring method by the total variation deep network. IEEE Access, 7:37536–37544, 2019.
[19] Zhengrong Xue. Blind image deblurring: a review. CoRR, abs/2201.10522, 2022.
[20] Min Zhang, Geoffrey S Young, Yanmei Tie, Xianfeng Gu, and Xiaoyin Xu. A new framework of designing iterative techniques for image deblurring. Pattern Recognition, 124:108463, 2022.

	$\displaystyle\\|\big{\langle}\text{{Hess}}(\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}(\bm{z}_{xy}))+\textbf{Hess}(\bm{g}_{xy}(\bm{z}_{xy})),\bm{u}_{xy}(\bm{z}_{xy})+\bm{g}_{xy}(\bm{z}_{xy})-\bm{d}_{xy}\big{\rangle}\\|_{\infty}$
	$\displaystyle\leq\begin{Vmatrix}\begin{pmatrix}(\bm{g}^{1}_{xy}(\bm{z}_{xy}))^{\prime\prime}\cdot\bm{\gamma}_{1}^{\top}\bm{\alpha}_{xy}&0&0&0\\ 0&(\bm{g}^{2}_{xy}(\bm{z}_{xy}))^{\prime\prime}\cdot\bm{\gamma}_{2}^{\top}\bm{\alpha}_{xy}&0&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\dots&(\bm{g}^{n_{t}}_{xy}(\bm{z}_{xy}))^{\prime\prime}\cdot\bm{\gamma}_{n_{t}}^{\top}\bm{\alpha}_{xy}\end{pmatrix}\\ \end{Vmatrix}_{\infty}$
	$\displaystyle\leq M_{1}\cdot\\|\bm{\gamma}\\|_{\infty}\bigg{(}\\|\bm{u}_{xy}(\bm{z}_{xy})\\|_{\infty}+\\|\bm{g}_{xy}(\bm{z}_{xy})\\|_{\infty}+\\|\bm{d}_{xy}\\|_{\infty}\bigg{)}$
	$\displaystyle\leq M_{1}\cdot\\|\bm{\gamma}\\|_{\infty}\bigg{(}\\|\bm{K}^{-1}\bm{A}^{\top}\bm{b}_{xy}\\|_{\infty}+M_{2}\\|\bm{z}_{xy}\\|_{\infty}+\\|\bm{d}_{xy}\\|_{\infty}\bigg{)}$
	$\displaystyle\leq M_{1}\cdot\\|\bm{\gamma}\\|_{\infty}\bigg{(}\\|\bm{K}^{-1}\bm{A}^{\top}\\|_{\infty}\\|\bm{b}_{xy}\\|_{\infty}+M_{2}\\|\bm{z}_{xy}\\|_{\infty}+\\|\bm{d}_{xy}\\|_{\infty}\bigg{)}$
	$\displaystyle\leq M_{1}\cdot\\|\bm{\gamma}\\|_{\infty}\bigg{(}\\|\bm{K}^{-1}\bm{A}^{\top}\\|_{\infty}\;(2\\|\bm{d}_{xy}\\|_{\infty}+2M_{2}\\|\bm{z}_{xy}\\|_{\infty})+M_{2}\\|\bm{z}_{xy}\\|_{\infty}+\\|\bm{d}_{xy}\\|_{\infty}\bigg{)}$
	$\displaystyle<\lambda_{1}.$

Bilevel Inverse Problems in Neuromorphic Imaging

Abstract.

Key words and phrases:

2010 Mathematics Subject Classification:

1. Introduction

2. Notation and Preliminaries

2.1. Neuromorphic (Event Based) Cameras

Neuromorphic Camera Basics

Data representation

2.2. Inverse problems in Image Deblurring

2.3. Existing Model, Algorithm, and Limitations

EDI Model

mEDI Model

mEDI Model - Optimizing for z

Limitations of the mEDI and EDI Models

3. Proposed Model and Algorithm

3.1. Model

Theorem 3.1.

Proof.

Lemma 3.2.

Proof.

Theorem 3.3.

Proof.

Lemma 3.4.

Proof.

Theorem 3.5.

Proof.

4. Experimental Setup

Standardizing 𝑩𝒊\bm{B^{i}}

Building the data cube

Calculate 𝒈​(𝒛)\bm{g}(\bm{z})

Omitting unnecessary calculations for 𝒛\bm{z}

Image Reconstruction

5. Numerical Examples

5.1. Benchmark Problem: Unit Bump

5.2. Multiple Event Based Dataset Examples

5.3. Event Based Problem: Shadowing Effects

Conclusion, Limitation, and Future Work

Acknowledgement

References

Standardizing $\bm{B^{i}}$

Calculate $\bm{g}(\bm{z})$

Omitting unnecessary calculations for $\bm{z}$