Create a custom Raspbian image with pi-gen: part 2

In my previous post, I discussed how to setup user accounts and locales in a custom Raspbian image using pi-gen. In this follow-up post, I will discuss the main problems that I want to solve: automatically configuring the wireless network and ssh on a new Raspberry Pi without a terminal attached directly to the Pi.

Set up the wireless network

The WPA supplicant

The Pi's wireless credentials are configured in stage 2 in the file stage2/02-net-tweaks/files/wpa_supplicant.conf. Here's how it looks by default:

ctrl_interface=DIR=/var/run/wpa_supplicant GROUP=netdev
update_config=1

According to the blogs Learn Think Solve Create and the Raspberry Spy, the first thing we should do is add our country code to the top of this file with the line country=CH. (Use your own country code for this.) Next, we want to enter the details for our wireless network, which includes its name and the password. For security reasons that I hope are obvious, we should not store the password in this file. Instead, we create a hash of the password and put that inside the file. The command to create the password hash is

wpa_passphrase ESSID PASSWORD > psk

where ESSID is our wireless network's name. Note that I also used a space before the wpa_passphrase command to prevent the password being written to my .bash_history file. Now, we copy and paste the contents of the file psk into the wpa_supplicant.conf and remove the comment that contains the actual password:

country=CH
ctrl_interface=DIR=/var/run/wpa_supplicant GROUP=netdev
update_config=1
network={
        ssid=YOUR_ESSID_HERE
        psk=YOUR_PSK_HASH_HERE
}

Configure the wireless network interfaces

After having configured the supplicant, we next move on to configuring the network interfaces used by Raspbian. The appropriate file is found in stage1/02-net-tweaks/files/interfaces. In my post Connecting a Raspberry Pi to a Home Linux Network I described how to set up the network interfaces by editing /etc/network/interfaces. Much of the information presented in that post has now been superseded in Raspbian by the DHCP daemon. For now, we will use the interfaces file to instruct our Pi to use DHCP and will use /etc/dhcpcd.conf at a later time to set up a static IP address when provisioning the Pi.

We first need to make a few changes to make the interfaces file aware of the credentials in the wpa supplicant configuration file. According to the blog kerneldriver, we need modify the /etc/networe/interfaces file as such:

auto lo

iface lo inet loopback
iface eth0 inet dhcp

auto wlan0
iface wlan0 inet manual
     wpa-roam /etc/wpa_supplicant/wpa_supplicant.conf
iface default inet dhcp

In the first modification, I specify that I want the wirless interface wlan0 started automatically with auto wlan0. Next, I specify that the wlan0 interface should use the manual inet address family with the line iface wlan0 inet manual.

According to the man pages, "[the manual] method may be used to define interfaces for which no configuration is done by default." After this we use the wpa-roam command to specify the location of the wpa_supplicant.conf file that we previously modified. The wireless ESSID and password are therefore not defined in interfaces, but rather reference them inside wpa_supplicant.conf.

In case you noticed that wpa-roam doesn't appear as an option in documentation on the interfaces file and were wondering why, it's because other programs like wpasupplicant may provide additional options to the interfaces file. A similar command is wpa-conf, but I do not quite yet understand the difference between these two commands.

Following the wpa-roam command, we configure the default options for all networks in our wpa_supplicant.conf file with the line iface default inet dhcp. At this point, we save the setup of the static IP address for a later time.

For more information, see the interfaces man page for Debian Stretch.

Change the hostname

Our Pi's hostname may be changed from the default (raspberrypi) by modifying the line in stage1/02-net-tweaks/files/hostname. See RFC 1178 for tips on choosing a hostname.

In addition to modifying the hostname file, we need to update stage1/02-net-tweaks/00-patches/01-hosts.diff and change raspberrypi to the new hostname:

Index: jessie-stage1/rootfs/etc/hosts
===================================================================
--- jessie-stage1.orig/rootfs/etc/hosts
+++ jessie-stage1/rootfs/etc/hosts
@@ -3,3 +3,4 @@
 ff02::1                ip6-allnodes
 ff02::2                ip6-allrouters

+127.0.1.1      NEW_HOSTNAME_HERE

Set the DNS servers

DNS servers are configured in export-image/02-network/files/resolv.conf. By default, mine was already configured to use one of Google's DNS servers (8.8.8.8). I added a secondary Google DNS address as well:

nameserver 8.8.8.8
nameserver 8.8.4.4

Enable SSH

Enabling SSH is simple. Open stage2/01-sys-tweaks/01-run.sh and change the line systemctl disable ssh to systemctl enable ssh.

(I later learned that we can also enable ssh on a headless pi by adding an empty file named ssh to the boot partition of a standard Raspbian image. See here for more details: https://www.raspberrypi.org/documentation/remote-access/ssh/)

Configuring SSH keys (or not)

I decided after writing much of this tutorial that pi-gen was not necessarily the best tool for adding my public SSH keys. So long as I have network access and SSH enabled, I can easily add my keys using ssh-copy-id. Furthermore, after following this tutorial, there still remains a lot of setup and customization steps. These can more easily be performed manually or by server automation tools like Fabric or Ansible.

Therefore, I think that at this point we can stop with our customization of the image with pi-gen and move to a different tool. We have a basic Raspbian image that is already configured for our home network and that serves as a starting point for more complete customization.

Conclusion

This tutorial and my previous post demonstrated how to create a custom Raspbian image that is pre-configured for

  • our home wireless network
  • our locale information
  • ssh

Of course, we can do much, much more with pi-gen, but other tools exist for the purpose of configuring a server. These tutorials at least allow you to setup a new Raspberry Pi without having to manually configure its most basic functionality. Happy Pi'ing!

Create a custom Raspbian image with pi-gen: part 1

Docker has been an amazing tool for improving my development efficiency on the Raspberry Pi. For example, I recently used it to cross-compile a large C++ and Python library for the Pi's ARM architecture on my x86_64 laptop. However, in that post I took it for granted that I had already set up my Raspberry Pi with user accounts, packages, ssh keys, etc. Performing these steps manually on a fresh install of the Pi's Raspbian operating system can become tedious, especially because ssh needs to be manually enabled before doing any remote work.

Fortunately, the Raspberry Pi developers have provided us with pi-gen, a useful collection of Shell scripts and a Docker container for creating custom Raspbian images. In this post, I will summarize the steps that I take in using pi-gen to create my own, personalized Raspbian image.

After I wrote this post, I found a set of posts at Learn Think Solve Create that describe many of the tasks I explain here. Be sure to check them out for another take on modifying Raspbian images.

Clone the pi-gen repository

This is as easy as cloning the git repository.

git clone git@github.com:RPi-Distro/pi-gen.git

Alternatively, you can use the https address instead of ssh, which is https://github.com/RPi-Distro/pi-gen.git.

From now on, all directories in this post will be relative to the root pi-gen directory.

Build the official Raspbian images

By default, the pi-gen repository will build the official Raspbian images. Doing this once before making any modifications is probably a good idea; if you can't build the official images, how will you be able to build a custom image?

There are two main scripts that you can use to do this: build.sh and build-docker.sh. build.sh requires you to install the packages that are listed in the repository's README.md file, whereas build-docker.sh requires only that you have Docker already installed on your computer. I'm going to be using the Docker-based build script for the rest of this post. If you don't have Docker installed on your system, you can follow the instructions here for the community edition.

Name your image

First we need to give a name to our image, even if we use the default build. To do this, we assign a name to a variable called IMG_NAME inside a file called config that is located inside the root pi-gen folder.

echo "IMG_NAME=my_name" > config

Build the default image

Once we've named our image, we can go ahead and run the build script.

./build-docker.sh

Be prepared to wait a while when running this script; the full build took well over half an hour on my laptop and with the Docker volume located on a SSD. It also consumed several GB of space on the SSD.

Resuming a failed build

The first time I used pi-gen the build failed twice. Once, it hung without doing anything for several minutes, so I canceled it with a Ctrl-C command. The other time I encountered a hash error when installing a Debian package.

We can resume a failed build from the point of failure by assigning the value 1 to the CONTINUE variable when calling build-docker.sh again.

CONTINUE=1 ./build-docker.sh

If we don't want to run previously built stages, we can simply place a file inside the corresponding folder named SKIP. For example, if our build fails at stage2, we can place SKIP files inside the stage0 and stage1 folders, then rerun the build-docker.sh script with CONTINUE=1.

Unfortunately, I have sometimes noticed that I have to also rebuild the stage prior to the one where the build failed. In the worst case, I had to rebuild all the stages because the fixes I applied to a file in stage2 were not accounted for when I tried to skip building stages 0 and 1. YMMV with this; I have no idea how well the SKIP mechanism works for the normal build.sh script.

After a successful build, we can find our custom images located inside the deploy folder of the pi-gen directory. These may then be written onto a SD card and used as a standard Raspbian image.

We can ensure that the build container is preserved even after successful builds using

PRESERVE_CONTAINER=1 ./build-docker.sh

Custom Raspbian images

Now that we've got the default build working, let's start by customizing the build process. For this post, I have the following goals:

  • Build only the lite version of the Raspbian images
  • Add a custom user account and delete the default pi account
  • Set the Pi's locale information

In a follow-up post, I will discuss the following:

  • Setup the WiFi for a home network
  • Setup ssh so that we can log on to the Pi remotely on its first startup

Building just Raspbian Lite

Raspbian Lite is a minimal Raspbian image without the X windows server and speciality modules that would otherwise make Raspbian more user friendly. It's an ideal starting point for projects that are highly specialized, require only a few packages, and do not require a GUI.

pi-gen creates Raspbian images in sequential steps called stages. At the time of this writing, there were five stages, with stages 2, 4, and 5 producing images of the operating system. Building everything from stage 0 up to and including stage 2 produces a Raspbian Lite image. We can speed up the build process and save harddrive space by disabling all the later stages.

To disable the build for a particular a stage, we add an empty file called SKIP inside the corresponding stage folder of the pi-gen root directory, just as we did above when skipping previously built stages. We also disable the explicit creation of images by adding an empty file called SKIP_IMAGES to stages 4 and 5. (We don't need to add a SKIP_IMAGES file to the stage3 folder because no image is produced at this stage.)

touch ./stage3/SKIP ./stage4/SKIP ./stage5/SKIP
touch ./stage4/SKIP_IMAGES ./stage5/SKIP_IMAGES

Now, when we run build-docker.sh, pi-gen will only build and produce one image for Raspbian Lite in the deploy directory.

Add a custom user account

The default user in Raspbian is called pi. This account is created in stage1 in the the script stage1/01-sys-tweaks/00-run.sh. This account is not very secure because it and its password, raspberry, are the well-known defaults in Raspbian. Let's go ahead and change them.

The relevant lines in the script look like this:

on_chroot << EOF
if ! id -u pi >/dev/null 2>&1; then
     adduser --disabled-password --gecos "" pi
fi
echo "pi:raspberry" | chpasswd
echo "root:root" | chpasswd
EOF

The user pi is created with the line adduser --disabled-password --gecos "" pi if it doesn't already exist. According to the adduser man pages The --disabled-password flag prevents the program passwd from setting the account's password when adduser is run, but remote logins without password authentication to the pi account are still allowed. the --gecos "" flag simply adds an empty string to the /etc/passwd file for the pi account.

After the user is created, raspberry is set as pi's password and root is set as the root password in the lines echo "pi:raspberry" | chpasswd and echo "root:root" | chpasswd.

Let's start by modifying the pi account. For the sake of this example, let's change its name to alphapi. For the password, we will generate a temporary, random password and write it to a file in the deploy directory. We'll do the same for root. The modifications look like the following:

user_passwd=$(< /dev/urandom tr -dc _A-Z-a-z-0-9 | head -c${1:-8})
root_passwd=$(< /dev/urandom tr -dc _A-Z-a-z-0-9 | head -c${1:-8})

# Write passwords to a file.
cat <<EOF > /pi-gen/deploy/users
${user_passwd}
${root_passwd}
EOF

on_chroot << EOF
if ! id -u alphapi >/dev/null 2>&1; then
     adduser --disabled-password --gecos "" alphapi
fi
echo "alphapi:${user_passwd}" | chpasswd
echo "root:${root_passwd}" | chpasswd
EOF

The first two lines create random alphanumeric passwords for the users alphapi and root. They should be changed immediately when the image is first run.

user_passwd=$(< /dev/urandom tr -dc _A-Z-a-z-0-9 | head -c${1:-8})
root_passwd=$(< /dev/urandom tr -dc _A-Z-a-z-0-9 | head -c${1:-8})

This way of password generation works by reading random bytes from /dev/urandom and redirecting them to the standard input of the tr command, which filters the input so only alphanumeric characters remain. Next, the output is piped to the head command, which outputs only the first eight alphanumeric characters produced in this fashion.

The passwords are then written to a file named users inside the deploy directory where the outputs will eventually be placed.

# Write passwords to a file.
cat <<EOF > /pi-gen/deploy/users
${user_passwd}
${root_passwd}
EOF

The remaining parts of the script are more-or-less the same as before, except I changed pi to alphapi and used variable substitution for the passwords.

Running ./build-docker.sh at this point will raise an error in stage02 because it's at this stage where the user pi is added to the various groups on the system. We therefore need to open stage2/01-sys-tweaks/01-run.sh and modify the following lines, replacing pi with alphapi.

for GRP in adm dialout cdrom audio users sudo video games plugdev input gpio spi i2c netdev; do
    adduser alphapi $GRP
done

Set the locale information

The locale information used by your operating system may be modified as follows. Open stage0/01-locale/00-debconf. I personally changed every occurence of en_GB.UTF-8 to en_US.UTF-8, but you can set your locale accordingly.

# Locales to be generated:
# Choices: All locales, aa_DJ ISO-8859-1, aa_DJ.UTF-8 UTF-8, ...
locales locales/locales_to_be_generated multiselect en_US.UTF-8 UTF-8
# Default locale for the system environment:
# Choices: None, C.UTF-8, en_US.UTF-8
locales locales/default_environment_locale   select  en_US.UTF-8

Next, we open stage2/01-sys-tweaks/00-debconf. I currently live in Europe, so I made the following changes:

tzdata        tzdata/Areas    select  Europe

I also made the following changes to switch from the default British English to American English:

keyboard-configuration keyboard-configuration/xkb-keymap select us
keyboard-configuration keyboard-configuration/fvariant  select  English (US) - English (US\, international with dead keys)

Note that the comment in 00-debconf above the keyboard-configuration/xkb-keymap line erroneously states that American English is an option, but it's not. You need to change it from "gb" to "us" if you want the American layout.

Using the custom image

With all these changes, we can build our new image by running ./build-docker.sh and, if successful, find a .zip file inside the deploy directory with the image name and date.

To use this image, we unzip the file to extract the .img file inside it. Next, we need to copy it onto a SD card that will plug into the pi. I have a SD card reader/writer on my laptop for which I check for its Linux device name by running lsblk before and after plugging in the card. (The device that appears in the output of lsblk after plugging it in is its name, which is /dev/mmcblk0 on my laptop). Once I get its device name, I use the Linux dd command to copy the contents of the image onto the card. (Be sure to change /dev/mccblk0 to match the name that your system gives to your SD card device.)

sudo dd if=2018-07-21-my_name-lite.img of=/dev/mmcblk0 bs=4096; sync

Please be EXTREMELY careful that you get the device name right. It's not very difficult to write the contents of the image file over your root partition or other important data.

After writing the image, we can plug the SD card into our pi, boot it up, and try logging in as alphapi with the random password that was created in the users file. Be sure at this point to change your user's and root's password. We can also verify that the keyboard was set to US English by typing Shift-3 and observing whether we get a hashtag (#) symbol and not the symbol for the British pound currency.

In a follow-up post, I will describe how to setup the network and SSH so I can continue to setup my Raspberry Pi without ever needing a terminal.

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 cross-compilation workflow

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.

The final prerequisite is to register QEMU with the Docker build agent. First, install a few packages for QEMU. On Ubuntu, this looks like

$ sudo apt update
$ sudo install qemu qemu-user-static qemu-user binfmt-support

Finally, register the build agent with the command:

$ docker run --rm --privileged multiarch/qemu-user-static:register --reset

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.
RUN useradd -ms /bin/bash micro-manager
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.
RUN useradd -ms /bin/bash micro-manager
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;

# Radial coordinates, image space
r = res_lateral * np.arange(0, oversampling * max_radius) / oversampling

# Radial coordinates, pupil space
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).

The Raspberry Pi Camera Module V2

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.

Remove the main ribbon cable from the control board

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.

The small ribbon cable from the sensor is to the right.

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.

Pull the sensor off the control board.

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.

Cut carefully on one side of the housing.

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

Once the edge is cut, pop the housing off.

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.

The resistor near the top of the board was damaged when cutting the housing.

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

The IR blocking filter.

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.

Optical system with pupils

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)
cax = divider.append_axes("right", size="5%", pad=0.05)
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)
cax = divider.append_axes("right", size="5%", pad=0.05)
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
Sensitivity 5.88 ADU / e-
Temporal dark noise 2.29 e-
Pixel size 3.45 microns
Analog-to-digital converter (ADC) 12-bit

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)
cax = divider.append_axes("right", size="5%", pad=0.05)
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)
cax = divider.append_axes("right", size="5%", pad=0.05)
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)
cax = divider.append_axes("right", size="5%", pad=0.05)
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)
cax = divider.append_axes("right", size="5%", pad=0.05)
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
max_adu     = np.int(2**bitdepth - 1)
adu         = (electrons_out * sensitivity).astype(np.int)
adu[adu > max_adu] = bitdepth # models pixel saturation

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

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

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
adu      += baseline

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

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

plt.show()