# How I built a cross-compilation workflow for the Raspberry Pi

Some of you may know I tinker with the Raspberry Pi in my free time and that one of my current projects is to build a lensless microscope with the Pi as the brains. To control the microscope, I decided a while ago that I would use Micro-Manager, an open-source software package for microscope control. I made this decision for a few reasons:

1. I already knew the Micro-Manager codebase since I use it frequently at work.
2. The Micro-Manager core provides a device-independent interface to hardware.
3. I've contributed to the project in the past and feel a sense of loyalty to the project and the people involved. Expanding Micro-Manager into embedded microscopy would be a great way for me to give back to the community.

Building Micro-Manager from source code presents its own set of challenges. After ensuring that you have the correct build environment, you need to actually compile it, and here's where things get tricky in Raspyberry Pi development. The Pi has an ARM processor, whereas most laptops and workstations use a x86_64 processor. This means that code compiled on a typical desktop PC will not work on the Pi. As I showed in my earlier post, you can compile the code directly on the Pi to circumvent this, but this unfortunately is quite cumbersome because the code base and dependencies are quite large. (They are nearly 8 GB in total). Furthermore, compiling the project on the Pi is slow and requires connecting to it via ssh or working directly on a TV screen or monitor.

These problems extend beyond Micro-Manager to other large-scale projects that require code compilation for a specific processor architecture. In this post, I'll describe the workflow that I developed for cross-compiling projects for the Raspberry Pi.

## Previous attempts

Prior to the workflow that is the main topic of this post, I managed to cross-compile Micro-Manager using a chroot environment and the QEMU emulator. chroot is a Linux command that changes the apparent root (or '/') directory for a running process. With this approach, I mount an image of the Raspbian operating system that contains the gcc and g++ compilers and libraries for the ARM architecture. Then, I chroot into the image and run a setup script that builds the software. During execution of this script, the QEMU static libraries run the ARM compilers from within the chroot environment to build the project. The compiled code remains inside the image, which I then burn onto a micro SD card to insert into the Pi. I uploaded a gist of the bash script which orchestrates all this, and my inspiration for this approach came from a great series of blog posts from Disconnected Systems.

Ultimately this approach is a huge amount of work. As you can see in the gist, it's fairly complicated bash scripting that's not easy to debug. Furthermore, the setup script that is run inside the image needs to do a lot of work beyond cross-compiling, like setting up the user, permissions, network, etc. Debugging the final product is also a challenge because you need to verify that it's working on the Pi, which requires burning the image to a micro SD card.

## Cross-compiling with Docker

After a bit of research I decided I would try instead to use Docker for cross-compilation and deployment to the Pi. I had just started using Docker at work to build reproducible environments for scientific computing research. In particular, and unlike my chroot script, I had learned that a Docker container that built the project could work on nearly any system that had Docker installed. Furthermore, deploying updates can be done on any Raspberry Pi that's running Docker.

I liked the idea of a portable cross-compilation workflow, so I dove into the Docker documentation and managed to get everything working in a few weeks of tinkering at home.

### An overview of Docker

You can find many resources online about Docker, so I won't go into the details here. The main thing you need to know is that Docker is a system for creating, running, and sharing containers, which are something like light weight virtual machines. Containers solve the problem in software development of how to build and deploy programs that have a complex set of dependencies. It does this by isolating the environment in which a program runs from the rest of the operating system. For example, if you have a computer that has a certain version of gcc (the GNU C compiler) installed, but your application requires a different version, then you can install the required gcc along with your application inside a container and they will not interfere with the version of gcc that belongs to your operating system. This also means that you can send your container to any machine that has Docker installed and it should just run without having to do any setup.

Other important things to know about Docker are:

• There are two main types of objects: images and containers. Images are sort of like blueprints that define what is inside a container, whereas containers are like the actual buildings specified by the blueprints. There can be many containers that come from a single image.
• Containers are meant to be immutable. When you stop them and restart them, they always restart in the same state as when they were first created.
• Since containers are immutable, some of your application data may need to be placed in a volume, which is either another container or a folder on the host system. A volume gets connected to your application container and exists even when your application container is not running.

## The cross-compilation workflow

Now that we have established the essential background to this project, let's look at the cross-compilation workflow. Below is a picture that provides a sense of the entire process, moving in general from left-to-right.

The process involves two Docker containers: one for building Micro-Manager and the other for running the application. The build dependencies and the QEMU emulator are both located inside the build container, having been specified when its image was created. These allow us to compile Micro-Manager for the ARM architecture. The source code is connected to the build container as a bind mount, which is a folder from the host workstation that is mounted inside the build container when it is run.

Once the libraries are compiled, they are installed into a folder inside the bind mount so that the host system will have access to them after the build container closes. Next, the compiled libraries are copied directly into an image that defines the application container. This image defines only the essential run-time requirements for running Micro-Manager and nothing else. The application image is stored on the registry server which I set up on my local network. This makes it easy for the Raspberry Pi to download the latest image and run the Micro-Manager application container whenever I make changes.

An important aspect of this workflow is how the data is passed between the systems and containers. Unlike what you will find in many introductory tutorials on Docker, I do not add the Micro-Manager source code directly to the build image/containers but instead use a bind mount. The reason for this is that the source code and 3rd party libraries are quite large, about 8 GB in total. By using a bind mount, I avoid needless copying of this data. Another reason for using a bind mount is that the source code will change frequently during development. If I add the source code to the image, then I will have to recreate the image every time the source code changes.

Once the libraries are built, I directly copy them into the application image because they are much, much smaller than the source code. I also want the code stored directly in the image so that the application image is all the Raspberry Pi needs to run the program. The image is stored in my local Docker registry server so that once I push an updated image to the server, the Raspberry Pi can download it and use it immediately.

### Step 0: Prerequisites

I am going to assume that you already have installed Docker. (If not, follow these directions.) I am also going to assume that you are somewhat familiar with how to work on a Linux system. The Raspberry Pi runs Linux, so you probably wouldn't be here if you didn't already know at least a little.

For this article, I am working with these versions of Docker and Ubuntu on my host workstation.:

kmdouglass@xxxxx:~$uname -a Linux xxxxx 4.13.0-39-generic #44~16.04.1-Ubuntu SMP Thu Apr 5 16:43:10 UTC 2018 x86_64 x86_64 x86_64 GNU/Linux kmdouglass@xxxxx:~$ docker version
Client:
Version:      18.03.1-ce
API version:  1.37
Go version:   go1.9.5
Git commit:   9ee9f40
Built:        Thu Apr 26 07:17:20 2018
OS/Arch:      linux/amd64
Experimental: false
Orchestrator: swarm

Server:
Engine:
Version:      18.03.1-ce
API version:  1.37 (minimum version 1.12)
Go version:   go1.9.5
Git commit:   9ee9f40
Built:        Thu Apr 26 07:15:30 2018
OS/Arch:      linux/amd64
Experimental: false


Finally, below is how my project directory structure is laid out.:

kmdouglass@xxxxx:~/src/alphapi/docker$tree -L 2 . └── rpi-micromanager ├── 2.0-python │ ├── build │ └── Dockerfile └── build ├── build ├── Dockerfile ├── run └── setup  I have two folders; build, which contains the files for the build container, and 2.0-python, which contains the files for creating the Micro-Manager application container. (In my case, I am going to build the Python wrapper for Micro-Manager 2.0.) Inside each folder are the scripts and Dockerfiles that execute the various steps of the workflow. ### Step 1: Create the build image Inside the build folder, I have a file called Dockerfile. Here are its contents: # Copyright (C) 2018 Kyle M. Douglass # # Defines a build environment for Micro-Manager on the Raspberry Pi. # # Usage: docker build \ # -t NAME:TAG \ # . # FROM resin/raspberrypi3-debian:stretch MAINTAINER Kyle M. Douglass <kyle.m.douglass@gmail.com> RUN [ "cross-build-start" ] # Get the build dependencies. RUN apt-get update && apt-get -y install --no-install-recommends \ autoconf \ automake \ build-essential \ git \ libatlas-base-dev \ libboost-dev \ libboost-all-dev \ libtool \ patch \ pkg-config \ python3-dev \ python3-pip \ python3-setuptools \ python3-wheel \ swig \ && apt-get clean && rm -rf /var/lib/apt/lists/* \ && pip3 install numpy RUN [ "cross-build-end" ] # Set up the mount point for the source files and setup script. ADD setup /micro-manager/ VOLUME /micro-manager/src WORKDIR /micro-manager/src ENTRYPOINT [ "/sbin/tini", "-s", "--" ] CMD [ "/micro-manager/setup" ]  A Dockerfile defines the steps in building an image -- in this case, the build image. Let's break this file down into pieces. In the first two lines that follow the comments, I specify that my image is based on the resin/raspberrypi3-debian:stretch image and that I am the maintainer.: FROM resin/raspberrypi3-debian:stretch MAINTAINER Kyle M. Douglass <kyle.m.douglass@gmail.com>  Images from Resin are freely available and already have the QEMU emulator installed. Next, I specify what commands should be run for the ARM architecture. Any commands located between RUN [ "cross-build-start" ] and RUN [ "cross-build-end" ] will be run using the emulator. Inside these two commands, I install the build dependencies for Micro-Manager using apt-get and pip. (These are just standard commands for installing software on Debian/Ubuntu Linux machines and from PyPI, respectively.) After the installation of the requirements completes, I add the setup script to the folder /micro-manager inside the image with the ADD setup /micro-manager/ command. The setup script contains the commands that will actually compile Micro-Manager. I then define a mount point for the source code with VOLUME /micro-manager/src. It's important to realize here that you do not mount volumes inside images, you mount volumes inside containers. This command is just telling the image to expect a folder to be mounted at this location when the container is run. The last three lines set the working directory, the entrypoint and the default container command, respectively.: WORKDIR /micro-manager/src ENTRYPOINT [ "/sbin/tini", "-s", "--" ] CMD [ "/micro-manager/setup" ]  This specific entrypoint tells Docker that any containers built from this image should first run Tini, which is a lightweight init system for Docker containers. If you do not specify Tini as the entry point, then it will not be able to reap zombies. (I don't know what this means exactly, but it sounds cool and you can read about it here: https://github.com/krallin/tini) By default, the container will run the setup script, but, since I used the CMD directive, this can be overriden in case we need to perform some manual steps. Roughly speaking, you can think of the entrypoint as the command that can not be overridden and the CMD command as the one that can be. In other words, Tini will always be executed when containers created from this image are launched, whereas you can choose not to run the setup script but instead to enter the container through a Bash shell, for example. To build the image, I use the following build script located in the same directory as the Dockerfile for convenience: #!/bin/bash # Copyright (C) 2018 Kyle M. Douglass # # Usage: ./build # docker build \ -t localhost:5000/rpi-micromanager:build \ .  By using -t localhost:5000/rpi-micromanager:build argument I am giving the image a name of rpi-micromanager, a tag of build, and specifying that I will eventually host this image on my local registry server (localhost) on port 5000. In case you are wondering about the contents of the setup script, don't worry. I'll explain it in the next section. ### Step 2: Compile Micro-Manager After the image is built, I create a container and use it to compile Micro-Manager. For this, I use the run script in the build directory: #!/bin/bash # Copyright (C) 2018 Kyle M. Douglass # # Usage: ./run DIR CONFIGURE # # DIR is the parent folder containing the micro-manager Git # repository, the 3rdpartypublic Subversion repository, and any # additional build resources. # # If CONFIGURE=true, the build system is remade and the configure # script is rerun before running 'make' and 'make install'. If # CONFIGURE=false, only 'make' and 'make install' are run. # # The compiled program files are stored in a bind mount volume so that # they may be copied into the deployment container. # src_dir=$1
cmd="/micro-manager/setup $2" # Remove the build artifacts from previous builds. if [ "$2" == true ] || [ "$2" == false ]; then rm -rf${src_dir}/build || true
fi

docker run --rm \
-v ${src_dir}:/micro-manager/src \ --name mm-build \ localhost:5000/rpi-micromanager:build \${cmd}


The script takes two arguments. The first is the path to the folder containing all the source code (see below for details). The second argument is either true or false. (It can actually be anything, but it will only compile Micro-Manager if either true or false are provided.) If true, the full build process is run, including setting up the configure script; if false, only make and make install are run, which should recompile and install only recently updated files.

The run script uses the -v argument to docker run to mount the source directory into the container at the point specified by the VOLUME command in the Dockerfile. The directory layout on my host file system for the source directory looks like this:

kmdouglass@xxxxx:/media/kmdouglass/Data/micro-manager$tree -L 1 . ├── 3rdpartypublic ├── micro-manager └── patches  The patches folder is not necessary and only there to fix a bug in the WieneckeSinscke device adapter. (This bug may be fixed by now.) 3rdpartypublic is the large Subversion repository of all the required software to build Micro-Manager, and micro-manager is the cloned GitHub repository. Prior to building, I checkout the mm2 branch because I am interested in developing my application for Micro-Manager 2.0. The setup script that is run inside the container and mentioned in the previous section looks like this: #!/bin/bash # # # Copyright (C) 2018 Kyle M. Douglass # # Builds Micro-Manager. # # Usage: ./setup CONFIGURE # # If CONFIGURE=true, the build system is remade and the configure # script is rerun before running 'make' and 'make install'. If # CONFIGURE=false, only 'make' and 'make install' or run. # # Kyle M. Douglass, 2018 # # Move into the source directory. cd micro-manager # Undo any previous patches. git checkout -- DeviceAdapters/WieneckeSinske/CAN29.cpp git checkout -- DeviceAdapters/WieneckeSinske/WieneckeSinske.cpp # Patch the broken WieneckeSinske device adapter. patch DeviceAdapters/WieneckeSinske/CAN29.cpp < ../patches/CAN29.cpp.diff \ && patch DeviceAdapters/WieneckeSinske/WieneckeSinske.cpp < ../patches/WieneckeSinske.cpp.diff # Compile MM2. if [ "$1" = true ]; then
# Remake the entire build system, then compile from scratch.
./autogen.sh
PYTHON="/usr/bin/python3" ./configure \
--prefix="/micro-manager/src/build" \
--with-python="/usr/include/python3.5" \
--with-boost-libdir="/usr/lib/arm-linux-gnueabihf" \
--with-boost="/usr/include/boost" \
--disable-java-app \
--disable-install-dependency-jars \
--with-java="no"
make
make install
chmod -R a+w /micro-manager/src/build
elif [ "$1" = false ]; then # Only recompile changed source files. make make install chmod -R a+w /micro-manager/src/build else echo "$1 : Unrecognized argument."
echo "Pass \"true\" to run the full build process."
echo "Pass \"false\" to run only \"make\" and \"make install\"."
fi


Most important in this script is the call to configure. You can see that the compiled libraries and Python wrapper will be written to the build folder inside the mounted directory. This gives the host file system access to the compiled artifacts after the container has stopped.

### Step 3: Build the application image

Once the libraries are compiled, we can add them to an application image that contains only the essentials for running Micro-Manager.

For this, I use a separate Dockerfile inside the 2.0-python directory:

# Copyright (C) 2018 Kyle M. Douglass
#
# Builds the Micro-Manager 2.0 Python wrapper for the Raspberry Pi.
#
# Usage: docker build \
#          -t NAME:TAG \
#          .
#

FROM resin/raspberrypi3-debian:stretch
MAINTAINER Kyle M. Douglass <kyle.m.douglass@gmail.com>

RUN [ "cross-build-start" ]

# Install the run-time dependencies.
RUN apt-get update && apt-get -y install --no-install-recommends \
libatlas-base-dev \
libboost-all-dev \
python3-pip \
python3-setuptools \
python3-wheel \
&& pip3 install numpy \
&& apt-get clean && rm -rf /var/lib/apt/lists/*

# Copy in the Micro-Manager source files.
WORKDIR /home/micro-manager/app
COPY --chown=micro-manager:micro-manager . .

RUN [ "cross-build-end" ]

# Final environment configuration.
USER micro-manager:micro-manager
ENV PYTHONPATH /home/micro-manager/app/lib/micro-manager
ENTRYPOINT ["/sbin/tini", "-s", "--"]
CMD ["/usr/bin/python3"]


As before, I use a clean resin base image. However, this time I only install the essential software to run Micro-Manager.

After apt-getting and pip-installing everything, I create a new user called micro-manager and a new folder called app inside this user's home directory:

# Copy in the Micro-Manager source files.
WORKDIR /home/micro-manager/app


Next, I directly copy the compiled libraries into the image with the COPY command:

COPY --chown=micro-manager:micro-manager . .


The two periods (.) mean that I copy the current host directory's contents into the container's current working directory (/home/micro-manager/app). What is the current host directory? Well, as I explain below, I actually run this Dockerfile from inside the build folder that was created to hold the compiled libraries in the previous step. But first, I'll end my explanation of the Dockerfile by saying that I switch the USER so that I do not run the container as root, add the library to the PYTHONPATH environment variable, and setup the default command as the python3 interpreter.

To build this image, I use the following build script:

#!/bin/bash
# Copyright (C) 2018 Kyle M. Douglass
#
# Usage: ./build DIR
#
# DIR is the root directory containing the Micro-Manager build
# artifacts. These artifacts will be added to the Docker image.
#

src_dir=$1 cp Dockerfile${src_dir}
cd ${src_dir} docker build \ -t localhost:5000/rpi-micromanager:2.0-python \ .  This script takes one argument, which is the build directory containing the compiled source code. The script first copies the Dockerfile into this directory and then changes into it with the cd command. (This explains the two periods (.) in the COPY command in the Dockerfile.) Finally, I build the image and give it a name of localhost:5000/rpi-micromanager:2.0-python. ### Step 4: Add the image to the local registry server Now we need a way to get the image from the workstation onto the Raspberry Pi. Of course, I could manually transfer the file with a USB stick or possibly use ssh, but what if I have multiple Pi's? This process could become cumbersome. Docker provides a few ways to push and pull images across a network. The most obvious is Dockerhub, a site for freely sharing images. For the moment I don't want to use Dockerhub, though, because I have not yet checked all the software licenses and am unsure as to what my rights are for putting an image with Micro-Manager software on a public repository. A better option, especially for testing, is to use a local registry server. This server operates only on my home network and already allows my workstation and Pi's to communicate with one another. Following the official registry documentation and this blog post by Zachary Keeton, I managed to setup the registry as follows. #### Host setup First, we need to setup a transport layer security (TLS) certificate. It's possible to run the server without one if you don't expect your network to be attacked, but it's good practice so let's create one. To do this, I edit the /etc/ssl/openssl.cnf file and add the following to the top of the [ v3_ca ] section.: subjectAltName = IP:192.168.XXX.XXX  where the IP address is the address of the workstation on the network. Next, I actually create the certificate. I make a directory called certs inside my workstation home directory and then use openssl to make the cerficate. During the prompts, I press ENTER at every step except the FQDN (fully qualified domain name). For the FQDN, I enter the same IP address as above.: mkdir certs openssl req -newkey rsa:4096 -nodes -sha256 \ -keyout certs/domain.key -x509 -days 365 \ -config /etc/ssl/openssl.cnf -out certs/domain.crt  I had to add the -config /etc/ssl/openssl.cnf argument for the subject alternative name to be added to the certificate. This part was tricky, because if this argument is not included, then the key generation step will use some other .cnf file (I am not sure which). This results in the following SAN error when attemptingt to connect to the registry.: cannot validate certificate for 192.168.XXX.XXX because it doesn't contain any IP SANs  After the domain.key and domain.crt files have been created, I run the official registry server container. (See how handy Docker containers are? There's no messy installation beyond grabbing the container.): docker run -d -p 5000:5000 \ --restart=always \ --name registry \ -v$(pwd)/certs:/certs \
-e REGISTRY_HTTP_TLS_CERTIFICATE=/certs/domain.crt \
-e REGISTRY_HTTP_TLS_KEY=/certs/domain.key \
registry:2


If the registry:2 image is not already downloaded, then it will be downloaded for automatically when running the container. Note that the -p 5000:5000 argument indicates that the server is using port 5000 on both the host system and inside the container. Note also that the certs directory is relative to the current directory because I use the ($pwd) command. You can change this to an absolute path if you wish on your setup. Let's go ahead and push the application image to the server now that it's running.: docker push localhost:5000/rpi-micromanager:2.0-python  #### Setup the Pi Now, startup the Pi. I will assume that you have already installed Docker on it and know how to communicate with it via ssh and copy files to it using scp. I copy the certificate from the host with scp: sudo mkdir -p /etc/docker/certs.d/192.168.XXX.XXX:5000/ sudo scp kmdouglass@192.168.XXX.XXX:/home/kmdouglass/certs/domain.crt /etc/docker/certs.d/192.168.XXX.XXX:5000/ca.crt  The IP address that I am using is the one to the machine where the registry server is running. After this step, I make the operating system trust the certificate: sudo scp kmdouglass@192.168.XXX.XXX:/home/kmdouglass/certs/domain.crt /usr/local/share/ca-certificates/192.168.XXX.XXX.crt sudo update-ca-certifications  Finally, I restart the Docker daemon.: sudo service docker restart  If everything is working, then I should be able to pull the image from your network's registry server: docker pull 192.168.XXX.XXX:5000/rpi-micromanager:python2.0  ### Step 5: Run Micro-Manager! And now the moment of truth: running the application container. Since it's setup to run Python automatically, I use a pretty simple docker run command.: docker run -it --rm \ --name micro-manager \ 192.168.XXX.XXX:5000/rpi-micromanager:2.0-python  I verify that the Micro-Manager Python wrapper is working by trying to import it and run a few basic commands: >>> import MMCorePy >>> mmc = MMCorePy.CMMCore() >>> mmc.getVersionInfo()  If these work without error, then congratulations! You're now ready to start building your embedded microscopy system ;) ### Step 6: Running the whole process The beauty of having scripted all these steps is that the full workflow may be executed quite simply. From the host system's build folder, run: kmdouglass@xxxxx:~/src/alphapi/docker/rpi-micromanager/build$ ./build
kmdouglass@xxxxx:~/src/alphapi/docker/rpi-micromanager/build$./run /path/to/source true  From the 2.0-python folder: kmdouglass@xxxxx:~/src/alphapi/docker/rpi-micromanager/2.0-python ./build /path/to/source/artifacts kmdouglass@xxxxx:~$ docker push localhost:5000/rpi-micromanager:2.0-python


And from the Raspberry Pi:

pi@yyyyy:~$docker pull 192.168.XXX.XXX:5000/rpi-micromanager:2.0-python pi@yyyyy:~$ docker run -it --rm \
--name micro-manager \
192.168.XXX.XXX:5000/rpi-micromanager:2.0-python


Hopefully this is enough to get you started building Micro-Manager for the Raspberry Pi with Docker. Though I focused on Micro-Manager, the workflow should be generally applicable to any large scale project in which you want to isolate the build environment from the host machine.

If you have any questions, just leave them in the comments. Happy programming!

# Implementing a fast Gibson-Lanni PSF solver in Python

Many techniques in modern fluorescence microscopy--such as deconvolution and single molecule localization microscopy--rely on the computation of accurate microscope point spread functions (PSFs). Several models exist for these purposes with varying degrees of accuracy and generality. One such model, the so-called Gibson-Lanni model, approximates the PSF of an immersion microscope objective by accounting for the different refractive indexes of the immersion medium, coverslip, and sample, as well as variations in the thicknesses of the different layers. This post describes a Python implementation of a fast computational approach to the G&L PSF that was presented in Fast and accurate three-dimensional point spread function computation for fluorescence microscopy by Li, Xue, and Blu. The approach is fast because it avoids the need to perform any explicit integration of the Kirchhoff diffraction integral, such as is done in PSF Generator. The prototype is based on the manuscript and the associated MATLAB code at http://www.ee.cuhk.edu.hk/~tblu/monsite/phps/fastPSF.php

## References¶

In [1]:
import sys
%pylab inline
import scipy.special
from scipy.interpolate import interp1d
from matplotlib import animation
plt.style.use('dark_background')
plt.rcParams['figure.facecolor'] = '#272b30'
plt.rcParams['image.cmap'] = 'viridis'

print('Python {}\n'.format(sys.version))
print('NumPy\t\t{}'.format(np.__version__))
print('matplotlib\t{}'.format(matplotlib.__version__))
print('SciPy\t\t{}'.format(scipy.__version__))

Populating the interactive namespace from numpy and matplotlib
Python 3.5.2 |Continuum Analytics, Inc.| (default, Jul  2 2016, 17:53:06)
[GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]

NumPy		1.13.3
matplotlib	2.1.0
SciPy		0.19.1


## Simulation setup¶

### Define the simulation parameters¶

The Gibson-Lanni PSF model accounts for differences between design parameters and real parameters that describe the fluorescence microscope, such as the refractive indexes and thicknesses of the immersion medium, coverslip, and sample. For more information, please see the above references or, if you do not have access to the publications, the description of the model here: http://bigwww.epfl.ch/algorithms/psfgenerator/

In [2]:
# Image properties
# Size of the PSF array, pixels
size_x = 256
size_y = 256
size_z = 128

# Precision control
num_basis    = 100  # Number of rescaled Bessels that approximate the phase function
num_samples  = 1000 # Number of pupil samples along radial direction
oversampling = 2    # Defines the upsampling ratio on the image space grid for computations

# Microscope parameters
NA          = 1.4
wavelength  = 0.610 # microns
M           = 100   # magnification
ns          = 1.33  # specimen refractive index (RI)
ng0         = 1.5   # coverslip RI design value
ng          = 1.5   # coverslip RI experimental value
ni0         = 1.5   # immersion medium RI design value
ni          = 1.5   # immersion medium RI experimental value
ti0         = 150   # microns, working distance (immersion medium thickness) design value
tg0         = 170   # microns, coverslip thickness design value
tg          = 170   # microns, coverslip thickness experimental value
res_lateral = 0.1   # microns
res_axial   = 0.25  # microns
pZ          = 2     # microns, particle distance from coverslip

# Scaling factors for the Fourier-Bessel series expansion
min_wavelength = 0.436 # microns
scaling_factor = NA * (3 * np.arange(1, num_basis + 1) - 2) * min_wavelength / wavelength


### Create the coordinate systems¶

In [3]:
# Place the origin at the center of the final PSF array
x0 = (size_x - 1) / 2
y0 = (size_y - 1) / 2

# Find the maximum possible radius coordinate of the PSF array by finding the distance
# from the center of the array to a corner
max_radius = round(sqrt((size_x - x0) * (size_x - x0) + (size_y - y0) * (size_y - y0))) + 1;

r = res_lateral * np.arange(0, oversampling * max_radius) / oversampling

a = min([NA, ns, ni, ni0, ng, ng0]) / NA
rho = np.linspace(0, a, num_samples)

# Stage displacements away from best focus
z = res_axial * np.arange(-size_z / 2, size_z /2) + res_axial / 2


## Step 1: Approximate the pupil phase with a Fourier-Bessel series¶

z.reshape(-1,1) flips z from a row array to a column array so that it may be broadcast across rho.

The coefficients C are found by a least-squares solution to the equation

\begin{equation*} \mathbf{\phi} \left( \rho , z \right)= \mathbf{J} \left( \rho \right) \mathbf{c} \left( z \right) \end{equation*}

$\mathbf{c}$ has dimensions num_basis $\times$ len(z). The J array has dimensions num_basis $\times$ len(rho) and the phase array has dimensions len(z) $\times$ len(rho). The J and phase arrays are therefore transposed to get the dimensions right in the call to np.linalg.lstsq.

In [4]:
# Define the wavefront aberration
OPDs = pZ * np.sqrt(ns * ns - NA * NA * rho * rho) # OPD in the sample
OPDi = (z.reshape(-1,1) + ti0) * np.sqrt(ni * ni - NA * NA * rho * rho) - ti0 * np.sqrt(ni0 * ni0 - NA * NA * rho * rho) # OPD in the immersion medium
OPDg = tg * np.sqrt(ng * ng - NA * NA * rho * rho) - tg0 * np.sqrt(ng0 * ng0 - NA * NA * rho * rho) # OPD in the coverslip
W    = 2 * np.pi / wavelength * (OPDs + OPDi + OPDg)

# Sample the phase
# Shape is (number of z samples by number of rho samples)
phase = np.cos(W) + 1j * np.sin(W)

# Define the basis of Bessel functions
# Shape is (number of basis functions by number of rho samples)
J = scipy.special.jv(0, scaling_factor.reshape(-1, 1) * rho)

# Compute the approximation to the sampled pupil phase by finding the least squares
# solution to the complex coefficients of the Fourier-Bessel expansion.
# Shape of C is (number of basis functions by number of z samples).
# Note the matrix transposes to get the dimensions correct.
C, residuals, _, _ = np.linalg.lstsq(J.T, phase.T)


### Compare the function to its Fourier-Bessel expansion¶

In [5]:
# Which z-plane to compute
z0 = 24

# The Fourier-Bessel approximation
est = J.T.dot(C[:,z0])

fig, ax = plt.subplots(ncols=2, sharey=True, figsize=(10,5))
ax[0].plot(rho, np.real(phase[z0, :]), label=r'$\exp{ \left[ jW \left( \rho \right) \right] }$')
ax[0].plot(rho, np.real(est), '--', label='Fourier-Bessel')
ax[0].set_xlabel(r'$\rho$')
ax[0].set_title('Real')
ax[0].legend(loc='upper left')

ax[1].plot(rho, np.imag(phase[z0, :]))
ax[1].plot(rho, np.imag(est), '--')
ax[1].set_xlabel(r'$\rho$')
ax[1].set_title('Imaginary')
ax[1].set_ylim((-1.5, 1.5))

plt.show()


## Step 2: Compute the PSF¶

Here, we use the Fourier-Bessel series expansion of the phase function and a Bessel integral identity to compute the approximate PSF. Each coefficient $c_{m} \left( z \right)$ needs to be multiplied by

\begin{equation*} R \left(r; \mathbf{p} \right) = \frac{\sigma_m J_1 \left( \sigma_m a \right) J_0 \left( \beta a \right)a - \beta J_0 \left( \sigma_m a \right) J_1 \left( \beta a \right)a }{\sigma_m^2 - \beta^2} \end{equation*}

and the resulting products summed over the number of basis functions. $\mathbf{p}$ is the parameter vector for the Gibson-Lanni model, $\sigma_m$ is the scaling factor for the argument to the $m'th$ Bessel basis function, and $\beta = kr\text{NA}$.

b is defined such that R has dimensions of len(r) $\times$ len(rho).

In [6]:
b = 2 * np. pi * r.reshape(-1, 1) * NA / wavelength

# Convenience functions for J0 and J1 Bessel functions
J0 = lambda x: scipy.special.jv(0, x)
J1 = lambda x: scipy.special.jv(1, x)

# See equation 5 in Li, Xue, and Blu
denom = scaling_factor * scaling_factor - b * b
R = (scaling_factor * J1(scaling_factor * a) * J0(b * a) * a - b * J0(scaling_factor * a) * J1(b * a) * a)
R /= denom


Now compute the point-spread function via

\begin{equation*} PSF \left( r, z; z_p, \mathbf{p} \right) = \left| \mathbf{R} \left( r; \mathbf{p} \right) \mathbf{c} \left( z \right) \right|^2 \end{equation*}

In [7]:
# The transpose places the axial direction along the first dimension of the array, i.e. rows
# This is only for convenience.
PSF_rz = (np.abs(R.dot(C))**2).T

# Normalize to the maximum value
PSF_rz /= np.max(PSF_rz)


### Profile in r and z¶

In [8]:
fig, ax = plt.subplots()

ax.imshow(PSF_rz, extent=(r.min(), r.max(), z.max(), z.min()))
ax.set_xlim((0,2.5))
ax.set_ylim((-6, 0))
ax.set_xlabel(r'r, $\mu m$')
ax.set_ylabel(r'z, $\mu m$')

plt.show()


z is the stage displacement in an inverted microscope. So, negative z's correspond to moving "up" in the sample.

## Step 3: Resample the PSF onto a rotationally-symmetric Cartesian grid¶

Here we generate a two dimensional grid where the value at each grid point is the distance of the point from the center of the grid. These values are supplied to an interpolation function computed from PSF_rz to produce a rotationally-symmetric 2D PSF at each z-position.

In [9]:
# Create the fleshed-out xy grid of radial distances from the center
xy      = np.mgrid[0:size_y, 0:size_x]
r_pixel = np.sqrt((xy[1] - x0) * (xy[1] - x0) + (xy[0] - y0) * (xy[0] - y0)) * res_lateral

PSF = np.zeros((size_y, size_x, size_z))

for z_index in range(PSF.shape[2]):
# Interpolate the radial PSF function
PSF_interp = interp1d(r, PSF_rz[z_index, :])

# Evaluate the PSF at each value of r_pixel
PSF[:,:, z_index] = PSF_interp(r_pixel.ravel()).reshape(size_y, size_x)


## Save z-stack to a movie¶

In [10]:
fig, ax    = plt.subplots(nrows=1, ncols=1)

# Rescale the PSF values logarithmically for better viewing
PSF_log = np.log(PSF)
vmin, vmax = PSF_log.min(), PSF_log.max()

img = ax.imshow(PSF[:,:,64], vmin=vmin, vmax=vmax)
img.set_extent([-r.max(), r.max(), -r.max(), r.max()])

txt = ax.text(3, 7, 'z = {:.1f} $\mu m$'.format(z[64]))
ax.set_xlim((-8, 8))
ax.set_xticks((-5, 0, 5))
ax.set_xlabel(r'x, $\mu m$')
ax.set_ylim((-8, 8))
ax.set_yticks((-5, 0, 5))
ax.set_ylabel(r'y, $\mu m$')
ax.spines['top'].set_visible(False)
ax.spines['bottom'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['right'].set_visible(False)

# Initialize the figure with an empty frame
def init():
img.set_data(np.zeros((size_y, size_x)))
return img,

# This function is called for each frame of the animation
def animate(frame):
img.set_data(np.log(PSF[:, :, frame]))

txt.set_text('z = {:.1f} $\mu m$'.format(z[frame]))
return img,

# Create the animation
anim = animation.FuncAnimation(fig, animate, frames=len(z),
init_func=init, interval=20, blit=True)

# Save the animation
myWriter = animation.FFMpegWriter(fps=10, extra_args=['-vcodec', 'libx264'])
anim.save('gibson-lanni.mp4', writer = myWriter)

plt.close()


## Summary¶

In summary, this method for calculating the Gibson-Lanni PSF is fast because it reduces the problem to a numer of matrix operations, rather than a numerical integration of the Kirchhoff diffraction integral. I especially find the Fourier-Bessel approximation of the phase function very clever, since this enables us to apply the integral identity described above. Moreover, it seems like the approach can extend to computing PSFs from other wavefront aberrations, though at the moment I do not know how it might handle non-rotationally symmetric pupils. (I am referring to the W variable above.)

I would like to thank the first author of the work, Jizhou Li, for encouraging me to write this as a blog post and for providing feedback on the implementation presented here.

# Accessing the Raspberry Pi camera image sensor

WARNING You can easily ruin your Raspberry Pi camera module by following the steps in this post. Proceed with caution.

Over the past half year I have been making slow but steady progress on my lensless imager project. The purpose, aside from having a bit of fun, is to create an imaging system for basic cell biology that doesn't use an expensive microscope objective.

A lensless imager works just as its name implies: an image sensor records the scattered light from a microscopic, transparent object and computationally reconstructs an image of that object, all without a lens. The best resolutions are achieved when the object is relatively close to the image sensor. Ideally, the separation between the object and the sensor's pixels would be at most about one millimeter. This limit is partly determined by the light source's spatial coherence, but also by the fact that high resolution is achieved by recording the scattered light at very large angles, which is possible only when the sample is close to the sensor.

Today I had a bit of free time so I decided to see whether I could remove the housing that surrounds the Raspberry Pi camera's sensor. The Raspberry Pi Camera Module version 2 sensor is a Sony IMX219 color chip. Directly above the sensor is a filter, a lens, and the housing for both of these that prevent me from placing anything closer than about half a centimeter from the sensor plane. If I would want to use this camera for the lensless imager, then the housing would have to go.

Now, even without the housing the Raspberry Pi camera is not necessarily the best option for the project because the IMX219 is a color sensor. This means that there is a Bayer filter over its pixels, which would cut the resolution of the imager since I would be using a very narrow band light source. Effectively, only a quarter of the pixels would be receiving any light. Regardless, I had a spare second camera and it interfaces well with the Raspberry Pi, so I figured it would make for a good prototype.

As you will see, I did slightly damage my sensor, though it seems to still work. You can easily ruin your camera module or cut your finger by following these steps, so proceed at your own risk.

## Step 0: The Raspberry Pi Camera Module V2

In the picture below you see my Raspberry Pi Camera Module. From above, you can see the small circular aperture with the lens immediately behind it. The sensor is inside the small gray rectangular housing that is attached to the control board by a small ribbon cable (to its right) and a bit of two-sided sticky foam tape (underneath the sensor; not visible).

## Step 1: Remove the main ribbon cable

To make working on the board a bit easier, I removed the white ribbon cable that attaches the module to the Pi. I did this by pulling on the black tabs on the two ends of the connecter until the cable is easily removed. I labeled the sides of the ribbon cable just in case.

## Step 2: Detach the sensor's ribbon cable

Next, I used my finger and thumb nail to remove the small ribbon cable that attaches the sensor to the control board. I essentially applied a small torque to the bottom edge of the connector until it just "popped" up, as seen in the second image below.

## Step 3: Remove the sensor from the control board

In the third step, I used my thumbnail to gently pry the sensor from the control board. The sensor is attached with some two-sided sticky tape and may need a few minutes of work to come free.

## Step 4: Remove the rectangular housing

In this step you risk cutting your finger, so please be careful.

The housing around the sensor is glued. To remove it, you will need to gently work a knife (or, better yet, a thin screw driver) between the housing and the sensor board, taking care not to let the blade go too far into the housing and possibly ruining one of the resistors or wire bonds.

Once you get a knife between the two, try popping the housing off of the sensor.

When I did this I cut on three sides of the housing, but in retrospect I should have only cut on the side opposite the ribbon cable and pried the other sides loose. This is because I damaged a small resistor when the knife blade went too far into the housing. You can see this below and, at the same time, get an idea of the layout of the sensor board so you know where you can and can't cut.

If you have the normal version of the camera, then you can also find the IR blocking filter inside the housing.

Fortunately for me the camera still works, despite the damaged resistor. I can now place samples directly on the sensor if I wanted to, though the wire bonds from the sensor to its control board appear quite fragile. For this reason, it may make more sense to build a slide holder that holds a sample just above the surface without touching it. For now, I can use this exposed sensor to prototype different methods for mounting the sample.

# Modeling noise for image simulations

In my previous post, I described a method for simulating inline holography images. The ultimate purpose of doing this was to generate holography data that could be used to test the performance of different reconstruction algorithms. Since I would know the input object, I could apply the reconstruction algorithm to the holograms and see how closely the reconstructed object matched the known input.

My modeling however ignored an important aspect about the holographic process, namely the acquisition of digital images and the assoicated noise that would accompany the raw data. This is important because holographic reconstruction algorithms need to be robust against the uncertainty that will inevitably be introduced by recording the hologram with a digital camera.

In [1]:
%pylab inline
import scipy
from numpy.fft import fft2, fftshift, ifftshift
from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib import animation
plt.style.use('dark_background')
plt.rcParams['figure.facecolor'] = '#272b30'
plt.rcParams['image.cmap'] = 'viridis'

print('Python version:\n{}\n'.format(sys.version))
print('Numpy version:\t\t{}'.format(np.__version__))
print('matplotlib version:\t{}'.format(matplotlib.__version__))
print('Scipy version:\t\t{}'.format(scipy.__version__))

Populating the interactive namespace from numpy and matplotlib
Python version:
3.5.2 |Continuum Analytics, Inc.| (default, Jul  2 2016, 17:53:06)
[GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]

Numpy version:		1.11.3
matplotlib version:	2.0.0
Scipy version:		0.19.0


## A Mathematical Model for Camera Noise¶

One of the best descriptions of camera noise that I have ever found may be the EMVA 1288 Standard. This document details a standard set of models and procedures to characterize CMOS and CCD cameras and light sensors. From its website:

EMVA launched the initiative to define a unified method to measure, compute and present specification parameters for cameras and image sensors used for machine vision applications... [The EMVA Standard] creates transparency by defining reliable and exact measurement procedures as well as data presentation guidelines and makes the comparison of cameras and image sensors much easier.

In particular, chapters 2 and 3 of the Standard provide a very clear and concise description the mathematical model for CCD/CMOS camera noise. Below, I have adapted the image of the mathematical model presented at the beginning of chapter 2.

In this model, the camera/sensor is a single pixel that takes an input (a number of photons) and transforms this number to an output (a digital value in analog-to-digital units, or ADU). In performing this transformation, various sources of noise are added to the signal so that, given a camera's output in ADU, we can only make probabilisitic inferences about the actual number of photons impinging upon it. There are also a few assumptions in this model that, when taken together, constitute what's called an ideal image sensor. They are:

1. the numer of photons collected by the pixel depends on the product of irradiance $E \, \left[ W / m^2 \right]$ and the exposure time $t_{exp} \, \left[ s \right]$;
2. the sensor is linear, i.e. the digital signal $y$ increases linearly with the number of photons received;
3. all noise sources are constant in time and space;
4. only the total quantum efficiency of the camera is dependent on the wavelength of the incident irradiation;
5. and only the dark current depends on the temperature.

From these five assumptions, the laws of quantum mechanics, and the laws of probability, expressions for the mean output, $\mu_y$, the variance of the output $\sigma_y^2$, and the signal-to-noise ratio ($SNR$) may be derived.

\begin{equation*} \mu_y = K \mu_d + K \eta \mu_p \\ \sigma_y^2 = K^2 \sigma_d^2 + K \left( \mu_y - K \mu_d \right) \\ SNR = \frac{\eta \mu_p}{\sqrt{\sigma_d^2 + \sigma_q^2/K^2 + \eta \mu_p}} \end{equation*}

In these expressions, $\mu$ and $\sigma^2$ represent the mean and variance of the number of photons hitting the camera ($p$), the dark noise ($d$), and the output signal in ADU's ($y$). The number of electrons $\mu_e$ generated by $\mu_p$ photons hitting the active area of the sensor is obtained from the expression for the quantum efficiency $\eta$, an engineered property of the camera that in general depends on the wavelength.

\begin{equation*} \eta \left( \lambda \right) = \frac{\mu_e}{\mu_p} \end{equation*}

The sensitivity of the camera $K \, \left[ ADU / e^- \right]$ represents the amplification of the voltage in the pixel from the photoelectrons and is also a property of the camera. Finally, $\sigma_q^2$ represents the noise introduced by digitzing the continous voltage signal into a digital one.

The EMVA 1288 Standard explains how to derive these expressions from the assumptions and known physical laws, so I will not go into their derivations in this post. Instead, I will focus on explaining how we can use this model to simulate the effects of camera noise in our holograms.

## From Model to Simulation¶

The goal of the simulation will be to take a simulated hologram as input and convert it into an image that might be captured by a digital camera. This means that the final image should contain noise that is realistic and follows the rules of the model.

The model will depend on five free parameters that vary from camera-to-camera. They are:

1. the total quantum efficiency $\eta$;
2. the read noise variance $\sigma_{r}^2$ (or its standard deviation);
3. the dark current $\mu_I$;
4. the sensitivity $K$;
5. and the bit-depth of the camera $k$.

The dark noise $\sigma_d^2$ is actually comprised of two separate noise sources, i.e. the read noise $\sigma_r^2$ and the dark current $\mu_I$. Because the noise sources are independent, we can write

\begin{equation*} \sigma_d^2 = \sigma_r^2 + \sigma_I^2 \\ \sigma_d^2 = \sigma_r^2 + \mu_I t_{exp} \end{equation*}

where the last step relies on the fact that the generation of thermal electrons is a Poisson process whose mean is equal to its variance. The units of $\mu_I$ are $e^- / s$ and values for $\mu_I$ are often provided on camera spec. sheets. Likewise, the read noise $\sigma_r^2$ is also provided on spec. sheets, although in my experience the units are likely to be in root-mean-square electrons or median electrons and not electrons-squared.

Other than the inputs, we do not need to provide values for the other numbers in the model because they are related to the free parameters by various physical and mathematical laws.

Finally, let's ignore spatial noise sources in these simulations because it's difficult to characterize them with a single number. Spatial noise is a variation in the free parameters of the camera from pixel to pixel. An example of a common type of spatial noise is photon response non-uniformity and is manifest in a sensitivity $K$ that varies between pixels. Ignoring spatial noise in the simulations simply means that values for the read noise, dark current, quantum efficiency, and sensitivty are the same for all pixels.

### Step 0: Convert the electric field to photons¶

Before doing anything, we need to take the electric field from our Fourier optics simulations and convert it to units of photons. If it isn't clear how you would first get an electric field, you can refer to my previous post on calculating the hologram of a thin, circular phase-only object. Typically, the calcuation of the electric field will follow from the laws of diffraction and serve as an input for the camera model.

To get the number of photons, we first compute the irradiance of the field using the magnitude of the time-averaged Poynting vector of a linearly-polarzed electromagnetic wave in free space, which is

\begin{equation*} E \left(x, y \right) = \frac{1}{2Z_o} \left| U\left(x, y \right) \right|^2 \end{equation*}

Here, $E \left( x, y \right)$ is the irradiance at point $\left( x, y \right)$, $U \left( x, y \right)$ is the scalar electric field at the same point, and $Z_0 \\approx 377 \, \Omega$ is the impedance of free space. We would change the value for the impedance if we were imaging in a medium other than vacuum. The units of the field are $V / m$ which, when squared and divided by the units of the impedance (Ohms), gives units of power per area, or $W / m^2$. It's also worth mentioning that we are implicitly assuming the field varies smoothly on the length scale of a pixel. This is why we can use the expression for the Poynting vector of a plane wave.

The mean number of photons is obtained from the irradiance by the expression

\begin{align*} \mu_p &= \frac{\lambda A}{hc}Et_{exp} \\ &= \frac{\lambda A}{2hcZ_0} \left| U \right|^2 t_{exp} \end{align*}

where $A$ is the pixel's area and $\frac{hc}{\lambda}$ is the energy per photon in units of Joules, with $hc \approx 2 \times 10^{-25} \, J \cdot m$.

This step is a bit cumbersome, and in the rest of what follows I will skip this step and will use a set number of photons directly as the input to the simulations, rather than a scalar field.

### Step 1: Compute the photon shot noise¶

For the sake of simplicity, let's start by assuming that we have a camera with $256 \times 256$ pixels and a mean incident photon flux of 500 photons on each pixel. This would produce an image that looks like this:

In [2]:
num_photons = 500
num_pixels  = 256
mu_p        = num_photons * np.ones((num_pixels, num_pixels))

fig, ax = plt.subplots()
img = ax.imshow(mu_p, vmin=400, vmax=600)
ax.set_xticks([])
ax.set_yticks([])
ax.set_title('No noise')

cb = plt.colorbar(img)
cb.set_label('Photons')

plt.show()


Photon shot noise arises from the fact that the incident number of photons is not deterministic but rather a random process whose statistics are Poissonian. In a Poisson process, the variance of the noise $\sigma_p^2$ is equal to the mean number of photons $\mu_p$.

Poisson random numbers are easy to generate with NumPy. First, we set a seed so that we can reproducibly generate the same random numbers each time and initialize a RandomState instance with it. Using this RandomState instance, we call the poisson method with a mean of num_photons and a size (num_pixels, num_pixels). Finally, we plot the image and compare it to the one where no shot noise is present.

In [3]:
seed       = 42
rs         = np.random.RandomState(seed)
shot_noise = rs.poisson(num_photons, (num_pixels, num_pixels))

fig, (ax0, ax1) = plt.subplots(ncols=2)
img0 = ax0.imshow(mu_p, vmin=400, vmax=600)
ax0.set_xticks([])
ax0.set_yticks([])
ax0.set_title('No shot noise')

divider = make_axes_locatable(ax0)
cb0 = plt.colorbar(img0, cax=cax)
cb0.set_ticks([])

img1 = ax1.imshow(shot_noise, vmin=400, vmax=600)
ax1.set_xticks([])
ax1.set_yticks([])
ax1.set_title('Shot noise')

divider = make_axes_locatable(ax1)
cb1 = plt.colorbar(img1, cax=cax)
cb1.set_label('Photons')

plt.show()

In [4]:
# It's also interesting to look at the distribution of photons hitting each pixel.
plt.hist(shot_noise.ravel(), bins=np.arange(350, 650))
plt.xlabel('Number of photons per pixel')
plt.ylabel('Frequency')
plt.show()


### Step 2: Compute the number of photoelectrons¶

At this point we're going to need some numbers for our camera. For this example, I chose a camera using the popular Sony IMX264 CMOS chip, the FLIR Chamleon 3 (part number CM3-U3-50S5M-CS). Here are some of the key specs for our modeling:

Spec Value
Quantum efficiency (at 525 nm) 0.69
Temporal dark noise 2.29 e-
Pixel size 3.45 microns

Assuming our light has a wavelength of 525 nm, then the number of photoelectrons is simply the product of the quantum efficiency with the realization of the field.

In [5]:
quantum_efficiency = 0.69
# Round the result to ensure that we have a discrete number of electrons
electrons = np.round(quantum_efficiency * shot_noise)

fig, (ax0, ax1) = plt.subplots(ncols=2)
img0 = ax0.imshow(shot_noise, vmin=200, vmax=600)
ax0.set_xticks([])
ax0.set_yticks([])
ax0.set_title('Photons')

divider = make_axes_locatable(ax0)
cb0 = plt.colorbar(img0, cax=cax)
cb0.set_ticks([])

img1 = ax1.imshow(electrons, vmin=200, vmax=600)
ax1.set_xticks([])
ax1.set_yticks([])
ax1.set_title('Electrons')

divider = make_axes_locatable(ax1)
cb = plt.colorbar(img1, cax=cax)

plt.show()


### Step 3: Simulate read noise and dark current¶

Following the schematic of the camera model above, we see that we next need to add the dark noise. Dark noise (not to be confused with dark current) refers to the noise in the pixels when there is no light on the camera. Typically we deal with two sources of dark noise: readout (or just "read") noise and dark current. The camera specs in our example only provide us with the total dark noise, which is not surprising since dark current is not very significant in microscopy or machine vision. (It is however significant in astronomy because of the very long exposure times of several seconds or more, which allows appreciable dark current electrons to accumulate during the exposure.)

To the best of my knowledge, when dark noise is primarily read noise we may well model it as a Gaussian distribution whose standard deviation is equivalent to the dark noise spec of the camera. This modeling approach is taken, for example, in the 2016 SMLMS Challenge, a community-driven competition to identify the best reconstruction algorithms for single molecule localization microscopy. Please note that the noise model employed in the challenge was for an electron multiplying CCD (EMCCD), so you may find that it differs slightly from what I am discussing here.

In [6]:
dark_noise = 2.29 # electrons
electrons_out = np.round(rs.normal(scale=dark_noise, size=electrons.shape) + electrons)

fig, (ax0, ax1) = plt.subplots(ncols=2)
img0 = ax0.imshow(electrons, vmin=250, vmax=450)
ax0.set_xticks([])
ax0.set_yticks([])
ax0.set_title('Electrons In')

divider = make_axes_locatable(ax0)
cb0 = plt.colorbar(img0, cax=cax)
cb0.set_ticks([])

img1 = ax1.imshow(electrons_out, vmin=250, vmax=450)
ax1.set_xticks([])
ax1.set_yticks([])
ax1.set_title('Electrons Out')

divider = make_axes_locatable(ax1)
cb = plt.colorbar(img1, cax=cax)
cb.set_label('Electrons')

plt.show()

In [7]:
# Plot the difference between the two
fig, ax = plt.subplots()
img = ax.imshow(electrons - electrons_out, vmin=-10, vmax=10)
ax.set_xticks([])
ax.set_yticks([])
ax.set_title('Difference')

cb = plt.colorbar(img)
cb.set_label('Electrons')

plt.show()


### Step 4: Convert from $e^-$ to ADU¶

Now we need to convert the value of each pixel from electrons to ADU, or analog-to-digital units. ADU's are the units of grey scale that are output by monochrome cameras. To do this, we multiply the number of electrons after the addition of read noise by the sensitivity. We also need to ensure that the final ADU count is discrete and to set the maximum upper value of ADU's to the $2^k - 1$ where $k$ is the camera's bit-depth. The latter operation ensures that saturation of the pixels is correctly modeled.

In [8]:
sensitivity = 5.88 # ADU/e-
bitdepth    = 12

fig, ax = plt.subplots()
ax.set_xticks([])
ax.set_yticks([])

cb = plt.colorbar(img)

plt.show()


### Step 5: Add a basline¶

Many cameras have a built-in baseline ADU value. This prevents the number of ADU's from becoming negative at low input signal. The Andor Zyla 4.2, for example, has a baseline value of 100 ADU. If you take a large number of dark frames and compute the average of each pixel in a Zyla, the average values will very nearly be 100.

We don't have baseline information for this camera, but 100 seems like a reasonable assumption. Here's what our final image will look like if we had a constant average of 500 photons incident on each pixel.

In [9]:
baseline  = 100 # ADU