Vectorfit

This code is distributed as accompanying software for Calibration-free estimation of field dependent aberrations for single molecule localization microscopy across large fields of view by Isabel Droste, Erik Schuitema, Sajjad A. Khan, Myron Hensgens, Stijn Heldens, Carlas S. Smith, Ben van Werkhoven, Hylkje Geertsema, Keith A. Lidke, Sjoerd Stallinga, and Bernd Rieger.

The paper is available here.

This code implements full vectorial point spread function (PSF) fitting for Single Molecule Localization Microscopy (SMLM) including field dependent aberrations that are estimated directly from the single molecule data.

The code is developed in MATLAB and includes a fast CPU/GPU implementation in C++ with CUDA support. This implementation can be accessed directly from MATLAB using MEX files.

Operating system

The CPU/GPU implementation has been tested on Linux and Windows.

Alternatively, the code can run entirely within MATLAB on both Linux and Windows without requiring compilation. However, this approach has significantly slower localization performance.

Requirements and installation

Linux CPU/GPU

The code was tested to work on Ubuntu 22.04 and Ubuntu 24.04 with Matlab R2024b, GCC 10 and CUDA 12.6.2, but may also work with other configurations. Make sure you use a version of GCC that is compatible with your version of Matlab. You can check the supported GCC compiler for your version of Matlab here.

Download or clone the repository

git clone https://gitlab.tudelft.nl/imphys/ci/vectorfit.git

The code requires Eigen3, FFTW, HDF5, and PkgConfig libraries. These are preinstalled in most high-performance computing environments. If you want to install them locally on Ubuntu you can use:

sudo apt install libeigen3-dev libfftw3-dev libhdf5-dev pkg-config

When running the code on an HPC environment, you need to load Matlab before compilation, together with a compatible version of GCC.

module load matlab/R2024b
module load compiler/gcc/10
module load cuda/12.6.2

omit the last line when running without GPU support.

To compile, please run the following commands inside the path of the source code. First, move to the cpp folder

cd cpp

Then, run the following commands to compile the code

mkdir build
cd build
cmake -DCMAKE_BUILD_TYPE=Release -DUSECUDA=1 ..
make

Omit the variable -DUSECUDA=1 to disable GPU support.

In case you get an error that looks like this

version 'GLIBCXX_3.4.32' not found (required by .../xxx.mexa64)

make sure you are using the correct version of GCC that is compatible with your version of matlab.

When running the code on an HPC environment with multiple versions of GCC installed, it can be the case that the systems picks up on the default GCC version, even when you loaded the correct version. In that case, you can force the correct GCC version using

export CC=/usr/bin/gcc-10
export CXX=/usr/bin/g++-10
export PATH=/usr/bin:$PATH

If that does not work, try:

export LD_PRELOAD=/usr/lib/x86_64-linux-gnu/libstdc++.so.6

To enable the CPU/GPU mode, set params.cpp = true in the parameter file. For the example code, the parameter file can be found in set_parameters_scripts/set_parameters_microtubules_3D. In addition, to use the GPU mode, set params.cpp_fitmode = "gpu-lowaccuracy" to perform the Fourier transforms in single precision and params.cpp_fitmode = "gpu" to perform the Fourier transforms in double precision. To use the CPU mode, set params.cpp_fitmode = "cpu".

The DIPImage toolbox for MATLAB is required, please see http://www.diplib.org for installation instructions.

Windows CPU

The code was tested to work on Windows 11 with Matlab R2024b, but may also work with other configurations. Unfortunately, however, it does not work for GPU on Windows.

Download or clone the repository

git clone https://gitlab.tudelft.nl/imphys/ci/vectorfit.git

The code requires Eigen3, FFTW, HDF5, and PkgConfig libraries, which can be installed using vcpkg.'

First, install vcpkg using

git clone https://github.com/microsoft/vcpkg.git
cd vcpkg
.\bootstrap-vcpkg.bat

Then install required libraries

vcpkg install eigen3 fftw3 hdf5 pkgconf

Integrate vcpkg with Visual Studio

vcpkg integrate install

To compile, please run the following commands inside the path of the source code. First, move to the cpp folder in the command line

cd cpp

Run the following commands (be sure to have modified the path to vcpkg)

mkdir build
cd build
cmake -G "Visual Studio 17 2022" -A x64 -DCMAKE_BUILD_TYPE=Release -DCMAKE_TOOLCHAIN_FILE=path/to/vcpkg/scripts/buildsystems/vcpkg.cmake ..

Then, in a Powershell, move to the cpp folder and run

cd build
cmake --build . --config Release

Matlab

The code was tested to work with Matlab R2024b but may work with other versions as well. Download the code or clone the repository using

git clone https://gitlab.tudelft.nl/imphys/ci/vectorfit.git

The DIPImage toolbox for MATLAB is required, please see http://www.diplib.org for installation instructions. It is also necessary to install the Parallel Computing Toolbox, Image Processing Toolbox and Bioinformatics Toolbox that can be found in the MATLAB Add-ons.

No further compilation of the code is needed. To enable the purely Matlab mode, set params.cpp = false; in the parameter file. For the example code, the parameter file can be found in set_parameters_scripts/set_parameters_microtubules_3D.

Software Demo

Follow the steps below to apply the SMLM pipeline to a test dataset.

• Download the test data and place the file microtubules_3D_data0001.tif in the folder vectorfit/data/0_raw_data/. This data contains the first 800 frames from the 3D microtuble dataset from Figure 2 from the article (This is 1% of the full dataset).

  When you are running the code on linux machine via ssh, you can copy the test data from your local machine to the linux machine via

  scp /path/to/local/file/microtubules_3D_data0001.tif username@remote_host:/path/to/vectorfit/data/0_raw_data/
  

• The parameter file is located at set_parameter_scripts/set_parameters_microtubules_3D.m. It contains all the imaging and fitting parameters. When using your own data, you need adjust the parameters to your own microscope parameters. When you change the name of the parameter file, make sure to declare it at the top of each of the three scripts segmentation.m, estimate_aberrations.m and localize.m that we will run below.

• When you are running the code on a high performance computing environment, you can load the required modules via

  module load matlab/R2024b cuda/12.6.2 diplib
  

• When you are running the code in an environment that supports graphical output (such as the Matlab user interface, or when running Matlab from the terminal with matlab -nodesktop), make sure that params.show_plots = true in the parameter file to enable plot visualisation. Otherwise, set params.show_plots = false.

• Start Matlab and add all scripts in the code directory to Matlab's path. You can do this by running

  addpath(genpath('rootDir'))
  

  where you replace rootDir with to path to the code directory, so for example: addpath(genpath('/home/idroste/vectorfit'))

• Within Matlab, navigate into the root directory of the code, and run the scripts below:

1. Run the segmentation script

   run('segmentation.m')
   

   This will segment the raw SMLM data into regions of iterest (ROIs). The expected output is

   Starting parallel pool (parpool) using the 'Threads' profile ...
   Connected to parallel pool with 12 workers.
   Found 1 files
   Input file = microtubules_3D_data0001.tif
   Found 800 frames
   Load 800 frames
   Finished loading 800 frames
   Elapsed time is 6.020052 seconds.
   Skip first 0 frames
   Found 169332 ROIs in file 1
   Saved segmentation in data/1_segmentation/segmentation_microtubules_3D_data0001.mat
   Elapsed time is 9.060499 seconds.
   
   Finished segmentation
   Elapsed time is 17.018615 seconds.
   Parallel pool using the 'Threads' profile is shutting down.
   

   When your output gives an error and shows Found 0 files on the third line, make sure that within Matlab you are inside the root directory of the code.

   Depending on your system, you may also get the error A minimum pool size of 10 was requested. The maximum thread-based pool size is currently 2. This is because MATLAB can use 10 compute threads internally, but thread-based parallel pool are limited to 2 workers on your system (in this case). To correct this, change params.nr_of_parallel_workers = feature('numcores'); to params.nr_of_parallel_workers = 2; in set_parameter_scripts/set_parameters_microtubules_3D.m.

   When plotting is enabled, the following 4 plots will appear:

   

   The upper left plot shows that raw SMLM data. You can step trough all of the frames by selecting Actions>Step through slices and then clicking the plot. The upper right plot can be used to determine the segmentation threshold. An optimal value lies slightly above the noise level. For the test data, the segmentation threshold is set to 15.5. You can change the segmentation threshold by changing the parameter params.segmentation_ROI_thr in the parameter file. The bottom left plot shows the centers of all ROIs in the field of view. The bottom right plot shows all segmented ROIs: you can step through them by selecting Actions>Step through slices and then clicking the plot.

   The segmented ROIs will be saved in the folder data/1_segmentation.

2. Next, run the aberration estimation script

   run('estimate_aberrations.m')
   

   This will load the segmented ROIs and estimate the field dependent aberrations from a subset of 5000 ROIs. When plotting is enabled, the following 4 plots will appear:

   

   The first plot shows the estimated aberration surfaces. The second plot shows the selected measured ROIs and the corresponding modeled PSFs, which you can step through by selecting Actions>Step through slices and then clicking the plot. The third plot shows the measured and modeled PSFs, averaged and upsampled on a 3x3 grid of the FOV. The fourth plot shows the evolution of gamma coefficients depending on the number of iterations.

   The estimated NAT coefficients can be found in theta_full.global. The output will be saved in data/2_estimate_aberrations/.

   [To do: add instruction on how to use the code for 2D data]

3. Finally, run the localization script

   run('localize.m')
   

   This will localize all segmented ROIs using the estimated field dependent aberrations. When plotting is enabled, the final image will be rendered:

   

   You can change the super-resolution pixel size and the rendered area of the FOV by adusting pixelsize, xlim and ylim at the bottom of the script localize.m.

   Use params.FlagOTF = true to use the precomputation of OTFs. This will only work in the CPU/GPU mode. The output localizations can be found in the array localizations and represent the parameters: [ID, x (nm), y (nm), z(nm), framenumber, CRLBx, CRLBy, CRLBz, Nph, bg]. The output will be saved in data/3_localize/.

Bead fitting

When you want to estimate aberrations from a z-stack of beads, use the script bead_fitting/fit_bead_aberrations.m. This code segments the raw data into 3D ROIs, and estimated the Zernike coefficient aberrations for each ROI separately. The parameters file set_parameters/set_parameters_nuclearlamin_240819_beads_crop.m can serve as an example parameter file for bead fitting.

Helpful notes

When using you own segmentation, make sure that your input to the rest of the code contains the following, where Ncfg is the total number of emitters:

allspots: Array of size [Mx,My,Mz,Ncfg] that contains the segmented ROIs

framelist: Array of size [1,Ncfg] that contains the framenumber of each ROI. NB: the first frame is 1 and not 0

ID: Array of size [1,Ncfg] that contains a unique number for each ROI. You can use 1,2,...,Ncfg

roixy: Array of size [2,Ncfg] that contains the pixel coordinates of the segmented ROIs within the FOV.

Troubleshooting and contributing

If you run into issues installing or using the software, please create an issue on the issue tracker, or contact the authors: Isabel Droste i.e.a.c.droste@tudelft.nl, Sjoerd Stallinga s.stallinga@tudelft.nl, Bernd Rieger b.rieger@tudelft.nl.