# An extremely simple deconvolution algorithm

A number of sophisticated deconvolution algorithms have been proposed in the literature, most of which involve solving an $L_1$ optimization problem.

While trying to develop a fast, on-line spike detection algorithm for Scanbox I tried a simple idea that seems to work reasonably well and the resulting algorithm is trivial.

In all likelihood this must have been published somewhere.  If you happen to know a relevant reference, just let me know in the comments below.

Suppose the calcium signal, $y_n$, is the response of a linear system with an exponential impulse response to a train of impulses $u_n$.

Then, the response satisfies $y_n = \alpha y_{n-1} + u_n + \epsilon_n$, where we assume $\epsilon_n$ is Gaussian noise.

One simple way to estimate $u_n$ is to first find the value of $\alpha$ that minimizes the error $E \{ ( y_n - \alpha y_{n-1} - u_n )^2 \}$.

A few calculations yield $\alpha = ( \mu^2 - m_{12} ) / ( \mu^2 - m_{02} ),$ where $\mu = E\{ y_n \}$, $m_{12} = E\{y_n y_{n-1}\}$ and $m_{02} = E\{y_n^2\}$.

As all these quantities are easily estimated recursively as new data arrive, one can also maintain a recursive estimate of the optimal $\alpha$.

Then, the estimated input to the system can be estimated by computing $\hat{u}_n = y_n - \alpha y_{n-1}$.  So, all we need to do is to keep a recursive estimate of $\alpha$ and the last two samples of the signal.

Below is an example of how this simple algorithm performs.A simulated spike train, $u_n$ (second panel form top) is convolved with an exponential impulse response and corrupted by additive Gaussian white noise to yield a simulated signal $y_n$ (top panel).  The third panel shows the estimated $\hat{u}_n$. Finally, we threshold this signal to obtain the estimated spike train shown in the bottom panel (thresholding is done by means of Otsu’s algorithm, implemented by Matlab’s gray thresh).

Here is the code if you want to play with it:

% generate ca-trace

u = double(rand(1,10000)&amp;lt;0.1); % original spikes
h = exp(-(0:40)/10);       % exponential impulse response
y = conv(h,u);             % convolve
y = y+0.15*randn(size(y)); % additive gaussian noise

% estimate alpha, uhat and input spike train
m = mean(y);
m02 = mean(y.^2);
m12 = mean(y(2:end-1) .* y(1:end-2));
a = (m^2-m12) / (m^2-m02);
uhat = [0 y(2:end) - a * y(1:end-1)];    % uhat
shat = uhat&amp;gt;graythresh(uhat);         % estimated spike train


Here is another example with a higher noise level and some bursting:

Here is how it runs on a some of the cases of the Chen et al dataset that have varying degrees of difficulty.  The red dots represent the spikes detected as described in the original paper, the black dots represents the spikes as detected by the deconvolution algorithm above.

Given the simplicity of the algorithm, it performs remarkably well and can be run across many ROIs in real time.

But wait… there is more. One may expect the impulse response to be not purely exponential. In such a case, it may be helpful to generalize the above to a higher-order model. You can then estimate $\hat{u}_n$ in just two lines of Matlab as follows:

a = lpc(y,10); % order=10
uhat = y-filter([0 -a(2:end)],1,y);


The above should bring back memories to anyone who studied speech compression and recognition at some point in time.

The graph below shows what happens when the impulse response is a gamma function (left).  Now, the LPC coefficients are non-zero for lags larger than one.

Still, and simple thresholding of $\hat{u}_n$ still recovers the original spike train even under the presence of substantial noise:

Go ahead, load some of your calcium imaging data and plot both your signal $y_n$ along with $\hat{u}_n$ estimated by the two lines above. Here is are a couple of data segments, again, from the Chen et al dataset.

Note: There are other pre-processing steps I took (like projecting out the neuropil signal) that I have not described above.

The calculation:

We want to minimize $J = E\{( y_n - \alpha y_{n-1} - u_n )^2\}$.  Taking the derivative with respect to $\alpha$ and equating to zero leads to

$E\{y_n y_{n-1}\} - \alpha E\{y_{n-1}^2\} - E\{u_n y_{n-1}\} = 0$

or

$m_{12}- \alpha m_{02} - \lambda \mu = 0$

where $E\{ u_n\} = \lambda$ and $E\{ y_n \} = \mu$.

But taking expectations on both sides of $y_n = \alpha y_{n-1} + u_n$ leads to

$\mu = \alpha \mu + \lambda$ or $\lambda = \mu (1-\alpha)$.

Substituting $\lambda$ in the equation above and solving for $\alpha$ gives $\alpha = (\mu^2 - m_{12})/(\mu^2 - m_{02})$.

# A montage display for real-time display of volumetric data

By default, Scanbox displays the incoming image stream on its main window. Thus, during volumetric scanning, one sees the incoming images as depth is changing over time.  If one is imaging only a handful of optical planes, it is difficult to see what is really going on.

A different way to visualize the data in such recordings is to have a montage showing the different optical sections separately, each being updated as new data becomes available. Here, we offer a plug-in for Scanbox that implements this mode of visualization.

As discussed in previous examples, Scanbox shares date with other processing by means of memory-mapped files. A header at the beginning of the file provides a mechanism for exclusive access to the data via a semaphore and exposes basic information about the data, such as the size of the images, the frame # being shared at any one time, and the period (in frames) of the volumetric scanning waveform, among other information.

Such data, along with a readily available montage function from Matlab, allows one to easily display the data as separate optical planes during acquisition.

Here is the code:


% Plug-in Demo: Display optical sections separately

close all; % Close all open figs

% Open memory mapped file -- define just the header first

mmfile = memmapfile('scanbox.mmap','Writable',true, ...
'Format', { 'int16' [1 16] 'header' } , 'Repeat', 1);
flag = 1;

% Process all incoming frames until Scanbox stops

while(true)

while(mmfile.Data.header(1)<0) % wait for a new frame...
if(mmfile.Data.header(1) == -2) % exit if Scanbox stopped
return;
end
end

if(flag) % first time? Format chA according to lines/columns in data
mmfile.Format = {'int16' [1 16] 'header' ; ...
mchA = double(intmax('uint16')-mmfile.Data.chA);
flag = 0;
I = zeros([size(mchA) 1 nplanes]);
I(:,:,1,1) = mchA;
imaqmontage(I);
axis off;           % remove axis
colormap gray;      % use gray colormap
else
imaqmontage(I);
end

mmfile.Data.header(1) = -1; % signal Scanbox that frame has been consumed!
drawnow limitrate;

end

clear(mmfile); % close the memory mapped file
close all;     % close all figures



And here are, side-by-side, the result of viewing the incoming stream in the main Scanbox window and in the volumetric plug-in.

# A rolling average plug-in for scanbox

Scanbox offers a memory-mapped interface mechanism to expose incoming data to other processes. This facilitates a wide rage of customization by users that may want to display their data in different ways or do some on-line processing not currently supported by Scanbox.

We previously provided one example of how a real-time histogram of the data can be generated. Today we offer a simple modification of this code to show a rolling average display, where frames are weighted with an exponential window, can be implemented.

These additional processes can run concurrently with the Scanbox live display if your computer is fast enough. To lighten the load on the computer you can now choose to disable the Scanbox display in the “Image Display” panel.

Remember that plug-ins run on a separate Matlab process and from the yeti/mmap directory.

The code is self-explanatory:


% Simply rolling average plug-in for Scanbox

% Open memory mapped file -- define just the header first

mmfile = memmapfile('scanbox.mmap','Writable',true, ...
'Format', { 'int16' [1 16] 'header' } , 'Repeat', 1);
flag = 1;

% Define the forgetting factor  0 < delta <= 1

delta = 0.9;  % this will generate an exponential decaying memory window: lambda^n

% Process all incoming frames until Scanbox stops

while(true)

while(mmfile.Data.header(1)<0) % wait for a new frame...
if(mmfile.Data.header(1) == -2) % exit if Scanbox stopped
return;
end
end

display(sprintf('Frame %06d',mmfile.Data.header(1))); % print frame# being processed

if(flag) % first time? Format chA according to lines/columns in data
mmfile.Format = {'int16' [1 16] 'header' ; ...
mchA = double(intmax('uint16')-mmfile.Data.chA);
flag = 0;
ih = imagesc(mchA); % setup display
axis off;           % remove axis
colormap gray;      % use gray colormap
truesize            % true image size
else
mchA = delta*mchA + (1-delta)*double(intmax('uint16')-mmfile.Data.chA);
ih.CData = mchA;
end

mmfile.Data.header(1) = -1; % signal Scanbox that frame has been consumed!

drawnow limitrate;

end

clear(mmfile); % close the memory mapped file
close all;     % close all figures



If you have any questions just post them below. Homework: Create a plug-in that shows areas of the image that are saturated in red (due date: Aug 8, 2016)

# Reading scanbox files in Python

One quick way to get your data into Python is by using the Matlab API for Python.

First, install the MATLAB engine API for Python by following he instructions here.  You will need to install as administrator.

Then, from your Python code you can simply do, as an example:


import matlab.engine
import numpy as np
from matplotlib import pyplot as plt

eng = matlab.engine.start_matlab()

q = np.array(z._data).reshape(z.size[::-1]).T

m = np.squeeze(np.ndarray.mean(q,axis=3))
plt.imshow(m,interpolation='nearest',cmap='gray')
plt.show()


This will read 50 frames starting at index 0, make it into a numpy array, take the mean over time, and display the result.

Creating the communication to the Matlab engine may take a few seconds, but after that the calls to sbxread() should be relatively fast.  This is, of course, not the most efficient way to read the data into Python.

An alternative, offered by Francisco Luongo at Caltech, is to use the following native Python function:

import os
import scipy.io as spio

'''
as it cures the problem of not properly recovering python dictionaries
from mat files. It calls the function check keys to cure all entries
which are still mat-objects
'''
return _check_keys(data)

def _check_keys(dict):
'''
checks if entries in dictionary are mat-objects. If yes
todict is called to change them to nested dictionaries
'''

for key in dict:
if isinstance(dict[key], spio.matlab.mio5_params.mat_struct):
dict[key] = _todict(dict[key])
return dict

def _todict(matobj):
'''
A recursive function which constructs from matobjects nested dictionaries
'''

dict = {}
for strg in matobj._fieldnames:
elem = matobj.__dict__[strg]
if isinstance(elem, spio.matlab.mio5_params.mat_struct):
dict[strg] = _todict(elem)
else:
dict[strg] = elem
return dict

'''
Input: filename should be full path excluding .sbx
'''
# Check if contains .sbx and if so just truncate
if '.sbx' in filename:
filename = filename[:-4]

#print info.keys()

# Defining number of channels/size factor
if info['channels'] == 1:
info['nChan'] = 2; factor = 1
elif info['channels'] == 2:
info['nChan'] = 1; factor = 2
elif info['channels'] == 3:
info['nChan'] = 1; factor = 2

# Determine number of frames in whole file
max_idx = os.path.getsize(filename + '.sbx')/info['recordsPerBuffer']/info['sz'][1]*factor/4-1

# Paramters
k = 0; #First frame
N = max_idx; #Last frame

nSamples = info['sz'][1] * info['recordsPerBuffer'] * 2 * info['nChan']

# Open File
fo = open(filename + '.sbx')

# Note: There is a weird inversion that happns thus I am dding the negative sign....
fo.seek(k*nSamples, 0)
x = np.fromfile(fo, dtype = 'uint16',count = nSamples/2*N)

x = -x.reshape((info['nChan'], info['sz'][1], info['recordsPerBuffer'], N), order = 'F')

return x