elgw/deconwolf

Deconvolution of widefield microscopy images and generation of point spread functions

deconwolf v0.4.0

deconwolf is a software for 3-D deconvolution of fluorescent wide-field images:

• The deconvolved images shows very mild boundary effects which means that relatively small regions of interest can be used.
• RAM usage can be reduced at the cost of slightly longer computation times by tiling. That makes it possible to deconvolve large images on small machines.
• Critical parts run on separate threads (as many as you would like). However, for maximal throughput (and if you have enough RAM), run several instances of dw in parallel.
• It is shipped with program to generate Point Spread Functions (PSFs) according to the Born and Wolf model. Our program is the only one we know of that actually integrate the PSF over each pixel.
• The programs are developed under Linux but can be compiled on Mac and Windows as well.
• Fully open source. And we embrace contributions and suggestions.

Deconwolf does not:

• Have a full-featured Graphical User Interface (GUI), however, if you like to click buttons there is a limited GUI.
• Show your images, for that, try ImageJ
• Diagnose your imaging system and correct for typical CMOS or CCD artifacts.
• Estimate your PSF based on real images.
• Tell you how many iterations that you should use although there are three ways to make it stop i) after a fixed number of iterations ii) when the progress is slow, or iii) at a fixed difference between the input image and the current guess (convolved with the PSF).
• If you miss one or more of these features, or just want something else, there is an (incomplete) list of alternative software.

Except for this README.me there is also a short USAGE.md, a CHANGELOG.md and a TODO.md.

After building and installing you will find two binaries:

• dw -- for deconvolution, man page
• dw_bw -- to create PSFs using the Born-Wolf model man page

At the moment we don't provide pre-built packages. You will have to build deconwolf from the source code, instructions follows below.

Method

The deconvolution algorithm is based on the Richardson-Lucy (RL) method ¹². The scaled heavy ball (SHB) ³ is used for acceleration, i.e. to do more work per iteration. Compared to the "Biggs" acceleration, it use less memory, it also has a proven convergence without any modifications to the algorithm.

The default boundary handling method is based on ⁴ , although extended to 3D. In our experiments it cause the least artifacts close to the image borders, especially important in the axial direction. To disable boundary handling, i.e. use periodic deconvolution, the option --periodic can be used.

The PSF calculations in dw_bw implements the "Born and Wolf" model ⁵ which is also decribed well in ⁶. It is possible to speed up the calculations using the method from ⁷ using the --li argument. At some point (after some more testing) it might become the default method. An important difference to the PSFGenerator is that each pixel in the PSF is integrated full, rather than just using one sample at the pixel centre. This modification gives measureable improvements, especially when the pixel size to resolution limit is low. dw psf also contains a few other PSF models which eventually will be documented as well.

For GPU calculations dw use VkFFT ⁸ via OpenCL. Until version 0.3.6 libclfft2 was used for this purpose. Since OpenCL is vendor neutral, this option should work for all major GPU brands. On the CPU side FFTW is used, although it can of course be swapped with API compatible alternatives like MKL.

Build and install

Deconwolf runs on 64-bit machines with, both aarch64 and x86_64, and require no special hardware. To compile and install deconwolf should take less than a minute on a Linux machine but might be more cumbersome on MacOS and Windows. For platform specific build instructions, see INSTALL.md. We hope to provide pre-compiled binaries in the future.

Dependencies

Deconwolf uses:

• FFTW for Fourier Transforms. Or alternatively Intel MKL.
• libtiff to read and write TIFF files.
• GNU Scientific Library for integration and special functions.
• OpenMP for automatic parallelization of code, used in dw.

If these libraries are available for your platform, chances are that deconwolf can be built as well.

To enable GPU acceleration, deconwolf also requires a GPU and OpenCL drivers installed.

Performance hints

The performance depends on the system and the image size. Except for some initialization (read images, padding, cropping ... ) and writing the output, the time is linear with the number of iterations.

Here is a table with performance figures for some real images. Run on a system with: a 4-core Intel i7-6700K CPU, 64 GB RAM, NVIDIA GeForce RTX 3090, using the --gpu flag and 50 iterations (--iter 50).

software	image size	job size	time (s)	sys-mem (Mb)
dw 1.3.7	2048x2048x35	2228x2228x103	37	28,748
dw 1.3.7	1024x1024x35	1204x1204x103	11	9,686
dw 1.3.7	512x512x35	692x692x103	4	7,564

Running only on the CPU the throughput drops drastically on this machine, having only 4 cores:

software	image size	job size	time (s)	sys-mem (Mb)
dw 1.3.7	2048x2048x35	2228x2228x103	1,079	10,802
dw 1.3.7	1024x1024x35	1204x1204x103	222	3,300
dw 1.3.7	512x512x35	692x692x103	78	1,238

Yes, deconvolution is time consuming without a GPU :(

Below there will be some hints on the performance relative to other software. For that purpose synthetic data will be used where periodic boundary conditions (--periodic) can be used. dw should not try to crop the PSF (which is usually done automatically to save some memory/gain some speed) in this case, hence the flag --xyfactor 0 is added as well.

Benchmarking is performed on the microtubules image using the accompanying PSF. Please note that it does not simulate "real" wide field data very well since it is created by periodic convolution.

System: Ubuntu 22.04.4 LTS, AMD Ryzen 7 3700X 8-Core Processor, 64 GB RAM, 12 GB RX 6700 XT GPU. Iterations: 115.

software	time (s)	self-mem (Mb)	sys-mem (Mb)
DeconvolutionLab2	1,025	1,582	48,511
DeconvolutionLab2 + FFTW2	862	1,353	47,387
MATLAB/deconvlucy	104		5,270
dw 1.3.7 --threads 1	52		344
dw 1.3.7 --threads 2	32		419
dw 1.3.7 --threads 4	21		566
dw 1.3.7 --threads 8	18		1,124
dw 1.3.7 --gpu	3		5,085

Notes:

• sys-mem is measured by parsing the VmPeak value from /proc/pid/status. In the case of DeconvolutionLab2 the values does not necessarily reflect the required memory since it is written in Java which is garbage collected. For MATLAB/deconvlucy the memory includes the full MATLAB environment.

• self-mem is the memory usage reported by the software if available.

• DeconlutionLab2 use vanilla RL, i.e. without any acceleration which means that more iterations will be needed before convergence.

• MATLAB/deconvlucy use "Biggs" acceleration. Matlab version R2020b was used in this case.

• For "real" data, when --periodic is not used, the input image is padded automagically during processing and the relevant part is cropped out at the end.

• Also works on Raspberry PI5 :) In that case it took 146 s.

Installation

These instructions should work under Linux, BSD, Windows (via WSL), and MacOS.

1. Get the dependencies. The required libraries should be found on most platforms, however, the installation process differs slightly. See INSTALL.md for per-platform advice.

2. Compile and install

mkdir builddir
cd builddir
cmake ..
cmake --build .
sudo make install

Once again, see INSTALL.md for more options.

Minimal usage example

To generate a parametric PSF and deconvolve an image, all you need is something like this:

# generate PSF.tif
dw_bw --resxy 130 \
--resz 250 \
--NA 1.46 \
--ni 1.518 \
--lambda 460 PSF.tiff
# Deconvolve image.tiff -> dw_image.tiff
dw --iter 50 image.tiff PSF.tiff

For available options, please see

dw --help
dw_bw --help

To validate that dw does not create random garbage, run it on /demo/dapi_001.tif

cd demo
make
imagej dapi_001.tif dw_dapi_001.tif

The run time on an AMD 3700x was 8s. To use GPU accelerated deconvolution, test

dw --iter 20 dapi_001.tif PSF_dapi.tif --gpu --prefix dwgpu

The results should be visually identical.

For more documentation see the short usage guide, and the manual pages for both binaries, man dw man dw_bw.

Bugs

Most likely there are bugs and they can only be fixed when they are known. Please open a new ticket if you have any issues with the program.

Alternatives

This is a non-complete list of alternative deconvolution software:

Deconvolution

Free and open source:

• Deconvolution Lab2

Commercial:

• Huygens
• Microvolution

Point spread functions

• PSF Generator

Others

• nd2tool Convert nd2 files (Nikon) to tif files. Can also generate scripts for deconvolution with dw.

References

Footnotes

1. Richardson, William Hadley (1972). "Bayesian-Based Iterative Method of Image Restoration". JOSA. 62 (1): 55–59. doi ↩

2. Lucy, L. B. (1974). "An iterative technique for the rectification of observed distributions". Astronomical Journal. 79 (6): 745–754. doi ↩

3. Wang H, et al. Scaled Heavy-Ball Acceleration of the Richardson-Lucy Algorithm for 3D Microscopy Image Restoration. IEEE Trans Image Process. 2014 doi. ↩

4. M. Bertero and P. Boccacci, A simple method for the reduction of boundary effects in the Richardson-Lucy approach to image deconvolution, A&A 437, 369-374 (2005). doi ↩

5. Max Born. Principles of optics : electromagnetic theory of propagation, interference, and diffraction of light. Cambridge: Cambridge University Press, 2019. ISBN: 978-1-108-47743-7. ↩

6. F. Aguet. “Super-Resolution Fluorescence Microscopy Based on Physical Models”. EPFL Thesis no. 4418 (2009), 209 p. Swiss Federal Institute of Technology Lausanne (EPFL), May 2009 url ↩

7. Jizhou Li, Feng Xue, and Thierry Blu. “Fast and accurate three-dimensional point spread function computation for fluorescence microscopy”. In: Journal of the Optical Society of America A 34.6 (May 2017), p. 1029. doi ↩

8. D. Tolmachev, "VkFFT-A Performant, Cross-Platform and Open-Source GPU FFT Library," in IEEE Access, vol. 11, pp. 12039-12058, 2023, doi ↩

deconwolf's Issues

include checksum/hash of the input/output files

For example XXH64

dw maxproj -- only works for files in the current folder

dw maxproj /some/dir/file.tif

does not work, while

cd /some/dir/
dw file.tif

does.

'cl_complex_mul.h' file not found when using GPU

Failed when using --gpu

(base) akikawa@DESKTOP-90048F8:/image$ dw --gpu --iter 40 --out ./dw_GPU/dw_dapi.tif ./raw/dapi.tif ./dw/dw_dapi_psf.tif
Reading /image/raw/dapi.tif
Reading /image/dw/dw_dapi_psf.tif
Output: ./dw_GPU/dw_dapi.tif(.log.txt)
Deconvolving using shbcl2 (using inplace)
image: [2304x2304x40], psf: [181x181x79], job: [2484x2484x118]
1 error generated.

Sorry, a fatal error occurred and dw can't continue
File: /home/akikawa/sources/repos/deconwolf/src/cl_util.c, Function: clu_kernel_newa Line: 1777
This happened when trying to compile the kernel cl_complex_mul
Build log: 
error: /home/akikawa/.cache/pocl/kcache/tempfile_au06iD.cl:1:10: 'cl_complex_mul.h' file not found
Device pthread-11th Gen Intel(R) Core(TM) i7-11700F @ 2.50GHz failed to build the program, log: error: /home/akikawa/.cache/pocl/kcache/tempfile_au06iD.cl:1:10: 'cl_complex_mul.h' file not found

and when using cpu it works well:

(base) akikawa@DESKTOP-90048F8:/image$ dw --iter 40 --out ./dw_GPU/dw_dapi.tif ./raw/dapi.tif ./dw/dw_dapi_psf.tif
Reading /image/raw/dapi.tif
Reading /image/dw/dw_dapi_psf.tif
Output: ./dw_GPU/dw_dapi.tif(.log.txt)
image: [2304x2304x40], psf: [181x181x79], job: [2484x2484x118]
creating fftw3 plans ... 
Iterating .^C

I wonder what the problem should be and if someone can give me some debugging directions. Thanks

btw, I am running it on wsl2 and have installed pocl as described in microsoft/WSL#6951 so that my opencl can detect GPU:

(base) akikawa@DESKTOP-90048F8:/image$ clinfo
Number of platforms                               1
  Platform Name                                   Portable Computing Language
  Platform Vendor                                 The pocl project
  Platform Version                                OpenCL 2.0 pocl 1.8  Linux, None+Asserts, RELOC, LLVM 11.1.0, SLEEF, DISTRO, POCL_DEBUG
  Platform Profile                                FULL_PROFILE
  Platform Extensions                             cl_khr_icd cl_pocl_content_size
  Platform Extensions function suffix             POCL

  Platform Name                                   Portable Computing Language
Number of devices                                 1
  Device Name                                     pthread-11th Gen Intel(R) Core(TM) i7-11700F @ 2.50GHz
  Device Vendor                                   GenuineIntel
  Device Vendor ID                                0x6c636f70
  Device Version                                  OpenCL 1.2 pocl HSTR: pthread-x86_64-pc-linux-gnu-icelake-client
  Driver Version                                  1.8
  Device OpenCL C Version                         OpenCL C 1.2 pocl
  Device Type                                     CPU
  Device Profile                                  FULL_PROFILE
  Device Available                                Yes
  Compiler Available                              Yes
  Linker Available                                Yes
  Max compute units                               16
  Max clock frequency                             2496MHz
  Device Partition                                (core)
    Max number of sub-devices                     16
    Supported partition types                     equally, by counts
    Supported affinity domains                    (n/a)
  Max work item dimensions                        3
  Max work item sizes                             4096x4096x4096
  Max work group size                             4096
  Preferred work group size multiple (kernel)     8
  Preferred / native vector sizes                 
    char                                                16 / 16      
    short                                               16 / 16      
    int                                                 16 / 16      
    long                                                 8 / 8       
    half                                                 0 / 0        (n/a)
    float                                               16 / 16      
    double                                               8 / 8        (cl_khr_fp64)
  Half-precision Floating-point support           (n/a)
  Single-precision Floating-point support         (core)
    Denormals                                     Yes
    Infinity and NANs                             Yes
    Round to nearest                              Yes
    Round to zero                                 Yes
    Round to infinity                             Yes
    IEEE754-2008 fused multiply-add               Yes
    Support is emulated in software               No
    Correctly-rounded divide and sqrt operations  Yes
  Double-precision Floating-point support         (cl_khr_fp64)
    Denormals                                     Yes
    Infinity and NANs                             Yes
    Round to nearest                              Yes
    Round to zero                                 Yes
    Round to infinity                             Yes
    IEEE754-2008 fused multiply-add               Yes
    Support is emulated in software               No
  Address bits                                    64, Little-Endian
  Global memory size                              65205055488 (60.73GiB)
  Error Correction support                        No
  Max memory allocation                           17179869184 (16GiB)
  Unified memory for Host and Device              Yes
  Minimum alignment for any data type             128 bytes
  Alignment of base address                       1024 bits (128 bytes)
  Global Memory cache type                        Read/Write
  Global Memory cache size                        16777216 (16MiB)
  Global Memory cache line size                   64 bytes
  Image support                                   Yes
    Max number of samplers per kernel             16
    Max size for 1D images from buffer            1073741824 pixels
    Max 1D or 2D image array size                 2048 images
    Max 2D image size                             32768x32768 pixels
    Max 3D image size                             2048x2048x2048 pixels
    Max number of read image args                 128
    Max number of write image args                128
  Local memory type                               Global
  Local memory size                               524288 (512KiB)
  Max number of constant args                     8
  Max constant buffer size                        524288 (512KiB)
  Max size of kernel argument                     1024
  Queue properties                                
    Out-of-order execution                        Yes
    Profiling                                     Yes
  Prefer user sync for interop                    Yes
  Profiling timer resolution                      1ns
  Execution capabilities                          
    Run OpenCL kernels                            Yes
    Run native kernels                            Yes
  printf() buffer size                            16777216 (16MiB)
  Built-in kernels                                (n/a)
  Device Extensions                               cl_khr_byte_addressable_store cl_khr_global_int32_base_atomics cl_khr_global_int32_extended_atomics cl_khr_local_int32_base_atomics cl_khr_local_int32_extended_atomics cl_khr_3d_image_writes cl_khr_fp64 cl_khr_int64_base_atomics cl_khr_int64_extended_atomics cl_khr_fp64

NULL platform behavior
  clGetPlatformInfo(NULL, CL_PLATFORM_NAME, ...)  No platform
  clGetDeviceIDs(NULL, CL_DEVICE_TYPE_ALL, ...)   No platform
  clCreateContext(NULL, ...) [default]            No platform
  clCreateContext(NULL, ...) [other]              Success [POCL]
  clCreateContextFromType(NULL, CL_DEVICE_TYPE_DEFAULT)  Success (1)
    Platform Name                                 Portable Computing Language
    Device Name                                   pthread-11th Gen Intel(R) Core(TM) i7-11700F @ 2.50GHz
  clCreateContextFromType(NULL, CL_DEVICE_TYPE_CPU)  Success (1)
    Platform Name                                 Portable Computing Language
    Device Name                                   pthread-11th Gen Intel(R) Core(TM) i7-11700F @ 2.50GHz
  clCreateContextFromType(NULL, CL_DEVICE_TYPE_GPU)  No devices found in platform
  clCreateContextFromType(NULL, CL_DEVICE_TYPE_ACCELERATOR)  No devices found in platform
  clCreateContextFromType(NULL, CL_DEVICE_TYPE_CUSTOM)  No devices found in platform
  clCreateContextFromType(NULL, CL_DEVICE_TYPE_ALL)  Success (1)
    Platform Name                                 Portable Computing Language
    Device Name                                   pthread-11th Gen Intel(R) Core(TM) i7-11700F @ 2.50GHz

Environment build error on windows10

I have installed the dependencies mentioned in INSTALL.md like:

git clone https://github.com/microsoft/vcpkg
.\vcpkg\bootstrap-vcpkg.bat
.\vcpkg integrate install
.\vcpkg\vcpkg.exe install fftw3
.\vcpkg\vcpkg.exe install tiff
.\vcpkg\vcpkg.exe install gsl
.\vcpkg\vcpkg.exe install opencl

and the result here:

(base) PS E:\TMC\sources> .\vcpkg\vcpkg.exe install fftw3
Computing installation plan...
The following packages are already installed:
    fftw3:[email protected]#8
fftw3:x64-windows is already installed
Total install time: 120 us
fftw3 provides CMake targets:

  # this is heuristically generated, and may not be correct
  find_package(FFTW3 CONFIG REQUIRED)
  target_link_libraries(main PRIVATE FFTW3::fftw3)

  find_package(FFTW3f CONFIG REQUIRED)
  target_link_libraries(main PRIVATE FFTW3::fftw3f)        

  find_package(FFTW3l CONFIG REQUIRED)
  target_link_libraries(main PRIVATE FFTW3::fftw3l)        

fftw3 provides pkg-config modules:

  # fast Fourier transform library
  fftw3

  # fast Fourier transform library
  fftw3f

  # fast Fourier transform library
  fftw3l

(base) PS E:\TMC\sources> .\vcpkg\vcpkg.exe install tiff
Computing installation plan...
The following packages are already installed:     
    tiff[core,zip,lzma,jpeg]:[email protected]#4  
tiff:x64-windows is already installed
Total install time: 102 us
tiff is compatible with built-in CMake targets:   

    find_package(TIFF REQUIRED)
    target_link_libraries(main PRIVATE TIFF::TIFF)

tiff provides pkg-config modules:

    #  Tag Image File Format (TIFF) library.      
    libtiff-4

(base) PS E:\TMC\sources> .\vcpkg\vcpkg.exe install gsl
Computing installation plan...
The following packages are already installed:
    gsl:[email protected]
gsl:x64-windows is already installed
Total install time: 184 us
The package gsl is compatible with built-in CMake targets:    

    find_package(GSL REQUIRED)
    target_link_libraries(main PRIVATE GSL::gsl GSL::gslcblas)

(base) PS E:\TMC\sources> .\vcpkg\vcpkg.exe install opencl
Computing installation plan...
The following packages are already installed:
    opencl:[email protected]#2
opencl:x64-windows is already installed
Total install time: 82.1 us
The package opencl is compatible with built-in CMake targets via CMake v3.6 and prior syntax

    find_package(OpenCL REQUIRED)
    target_link_libraries(main PRIVATE ${OpenCL_LIBRARIES})
    target_include_directories(main PRIVATE ${OpenCL_INCLUDE_DIRS})

and the CMake v3.7 and beyond imported target syntax

    find_package(OpenCL REQUIRED)
    target_link_libraries(main PRIVATE OpenCL::OpenCL)

This package is only an OpenCL SDK. To actually run OpenCL code you also need to install an implementation.

    WINDOWS: implementations typically ship with the drivers of you CPU/GPU vendors.
    LINUX: implementations may be installed from your distro's repo or manually. There are too many to count.       
    APPLE: consult your distribution vendor on the state of OpenCL support: https://support.apple.com/en-us/HT202823

but error occurs when i tried to build deconwolf:

(base) PS E:\TMC\sources\deconwolf\winbuild> cmake .. -G "Visual Studio 17 2022" -T ClangCL -A x64
-- Selecting Windows SDK version 10.0.22621.0 to target Windows 10.0.19045.
-- The C compiler identification is Clang 17.0.3 with MSVC-like command-line
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Check for working C compiler: D:/Microsoft Visual Studio/2022/Community/VC/Tools/Llvm/x64/bin/clang-cl.exe - skipped
-- Detecting C compile features
-- Detecting C compile features - done
-- Build type: Release
-- The CXX compiler identification is Clang 17.0.3 with MSVC-like command-line
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Check for working CXX compiler: D:/Microsoft Visual Studio/2022/Community/VC/Tools/Llvm/x64/bin/clang-cl.exe - skipped
-- Detecting CXX compile features
-- Detecting CXX compile features - done
Enabling IPO
-- Could NOT find PkgConfig (missing: PKG_CONFIG_EXECUTABLE) 
CMake Error at D:/CMake/share/cmake-3.29/Modules/FindPackageHandleStandardArgs.cmake:230 (message):
  Could NOT find FFTW (missing: FFTW_INCLUDE_DIRS FLOAT_LIB FLOAT_OPENMP_LIB)
Call Stack (most recent call first):
  D:/CMake/share/cmake-3.29/Modules/FindPackageHandleStandardArgs.cmake:600 (_FPHSA_FAILURE_MESSAGE)
  util/findFFTW/FindFFTW.cmake:399 (find_package_handle_standard_args)
  CMakeLists.txt:175 (find_package)


-- Configuring incomplete, errors occurred!

https://no-color.org/

support this.

write to log when GPU runs out of memory

As it is now the message appear to screen but not to the log file.

Cannot properly update parameters in dw_guide

When running dw_guide I cannot change the parameters properly. If I select threads I can only input single-character options. Any additional keystroke overwrites the first one... Let me know if you need additional details.

use quoting when saving the command line to the log files

So that the command can be copied/pasted directly without editing if you want to run it again.

iterdump usage

1. Change output names to:

"itd%05d_%s", iter, file

that would be more clear than using x_...

2. Don't write any

dw_<file> 

at the end. This file is already likely to exist (or to be generated later on).

Scaling of 16 bit images

• Currently images are scaled to have a max of 2^16-1. It would be better to scale them up to 2^16-2 in order to not get warnings about saturation downstream in the analysis.
• It would make sense to have --float as the default and use a flag (--16bit?) to save as 16-bit.

Missing header files in src/kernels ?

Hello,
When trying to install deconwolf :

[ 61%] Building C object CMakeFiles/dw.dir/src/cl_util.c.o
/home/s_montap01/src/deconwolf-0.3.8/src/cl_util.c:4:10: fatal error: kernels/cl_complex_mul.h: No such file or directory
#include "kernels/cl_complex_mul.h"

In the github repo, there are no .h files under "kernels/", and I cannot find cl_complex_mul.h anywhere. Did I miss something

Gaussian low pass filter for maxprojections

Add an option --sigma or --lsigma (lateral sigma) for Gaussian low pass during max projections.