Skip to content

Commit

Permalink
code reorganization, and docs
Browse files Browse the repository at this point in the history
  • Loading branch information
tlambert03 committed Mar 5, 2019
1 parent 4d79a90 commit 4d7b86f
Show file tree
Hide file tree
Showing 11 changed files with 250 additions and 138 deletions.
8 changes: 1 addition & 7 deletions .appveyor.yml
Original file line number Diff line number Diff line change
Expand Up @@ -30,13 +30,7 @@ for:
environment:
LABEL: 'main'
BRANCH: 'master'
-
branches:
only:
- develop
environment:
LABEL: 'dev'
BRANCH: 'develop'


install:
# If there is a newer build queued for the same PR, cancel this one.
Expand Down
5 changes: 5 additions & 0 deletions docs/affine.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
Affine Transformations
======================

.. automodule:: pycudadecon
:members: affineGPU, deskewGPU, rotateGPU
23 changes: 20 additions & 3 deletions docs/deconvolution.rst
Original file line number Diff line number Diff line change
@@ -1,5 +1,22 @@
deconvolution functions
=======================
Deconvolution
=============

The primary function for performing deconvolution is :func:`pycudadecon.decon`.
This convenience function is capable of receiving a variety of input types (filenames, directory names, numpy arrays, or a list of any of those) and will handle setting up and breaking down the FFT plan on the GPU for all files being deconvolved. Keywords arguments will be passed internally to the :class:`pycudadecon.RLContext` context manager or the :func:`pycudadecon.makeotf` :func:`pycudadecon.rl_decon` functions.

The setup and breakdown for the GPU-deconvolution can also be performed manually:
1. call :func:`pycudadecon.rl_init` with the shape of the raw data and path to OTF file
2. perform deconvolution(s) with :func:`pycudadecon.rl_decon`
3. cleanup with :func:`pycudadecon.rl_cleanup`

As a convenience, the :class:`pycudadecon.RLContext` context manager will perform the setup and breakdown automatically:

.. code-block:: python
data = tiffile.imread('some_file.tif')
otf = 'path_to_otf.tif'
with RLContext(data.shape, otf) as ctx:
result = rl_decon(data, output_shape=ctx.out_shape)
.. automodule:: pycudadecon
:members: decon, TemporaryOTF
:members: decon, rl_init, rl_decon, rl_cleanup, RLContext,
11 changes: 6 additions & 5 deletions docs/index.rst
Original file line number Diff line number Diff line change
@@ -1,15 +1,16 @@
pyCUDAdecon: C++/CUDA-accelerated deconvolution with a python wrapper
=====================================================================
pyCUDAdecon
===========

This package provides a python wrapper and convenience functions for `cudaDeconv <https://github.com/dmilkie/cudaDecon>`_ (which is a CUDA/C++ implementation of an accelerated Richardson Lucy Deconvolution). cudaDeconv was originally written by `Lin Shao <https://github.com/linshaova>`_ and modified by `Dan Milkie <https://github.com/dmilkie>`_, at Janelia Research campus. This package makes use of a shared library interface that I wrote for cudaDecon while developing `LLSpy <https://github.com/tlambert03/LLSpy>`_, that adds a couple additional kernels for affine transformations and camera corrections. The code here is mostly extracted from that package to allow it to be used independently of LLSpy.

The primary function is
pyCUDAdecon is a python wrapper and set of convenience functions for `cudaDeconv <https://github.com/dmilkie/cudaDecon>`_ (which is a CUDA/C++ implementation of an accelerated Richardson Lucy Deconvolution). cudaDeconv was originally written by `Lin Shao <https://github.com/linshaova>`_ and modified by `Dan Milkie <https://github.com/dmilkie>`_, at Janelia Research campus. This package makes use of a shared library interface that I wrote for cudaDecon while developing `LLSpy <https://github.com/tlambert03/LLSpy>`_, that adds a couple additional kernels for affine transformations and camera corrections.

.. toctree::
:maxdepth: 2
:caption: Contents:

deconvolution
otf
affine




Expand Down
7 changes: 7 additions & 0 deletions docs/otf.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
Creating Optical Transfer Functions
===================================

These are functions for converting a 3D point spread function (PSF) volume into a radially averaged 2D complex Optical Transfer Function (OTF) that can be used for deconvolution. You can either write the OTF to file, for later use, using the :func:`pycudadecon.makeotf` function, or use the :class:`pycudadecon.TemporaryOTF` context manager to create and delete a temporary OTF from a 3D PSF input.

.. automodule:: pycudadecon
:members: makeotf, TemporaryOTF
6 changes: 4 additions & 2 deletions pycudadecon/__init__.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,9 @@
from .version import __version__
from .deconvolution import decon, TemporaryOTF

from .deconvolution import decon
from .otf import makeotf, TemporaryOTF
from .libcudawrapper import RLContext, rl_decon, rl_init, rl_cleanup
from .affine import affineGPU, deskewGPU, rotateGPU

__all__ = (decon, TemporaryOTF, RLContext, rl_decon, rl_init, rl_cleanup,
affineGPU, deskewGPU, rotateGPU)
affineGPU, deskewGPU, rotateGPU, makeotf, TemporaryOTF)
21 changes: 20 additions & 1 deletion pycudadecon/affine.py
Original file line number Diff line number Diff line change
Expand Up @@ -119,7 +119,26 @@ def affineGPU(im, tmat, dzyx=None):


def rotateGPU(im, dzdata, dxdata=0.1, angle=31.5, reverse=False):
# TODO: crop smarter
"""Rotate image around Y axis by some angle
This is a convenience function that will apply the appropriate affine
transformation for rotating a volume around the Y axis by some angle.
This is typically done with images acquired on inverted light sheet
microscopes where the image plane is not parallel to the coverslip
(such as lattice light sheet, or diSPIM microscopes), in order to change
the coordinate system of the image volume such that the Z axis is
orthogonal to the coverslip
Args:
im (np.ndarray): 3D volume to be rotated
dzdata (float): Z-step size in microns of the image volume
dxdata (float): XY pixel size of the volume (default: {0.1})
angle (float): Angle to rotate around Y axis (default: {31.5})
reverse (bool): Rotate in the opposite direction. (default: {False})
Returns:
np.ndarray: Rotated 3D volume
"""
angle = float(angle)
xzRatio = dxdata / (np.deg2rad(angle) * dzdata)

Expand Down
62 changes: 1 addition & 61 deletions pycudadecon/deconvolution.py
Original file line number Diff line number Diff line change
@@ -1,10 +1,8 @@
from .util import is_otf
from .libcudawrapper import RLContext, rl_decon
from .otf import makeotf
from .otf import TemporaryOTF
import os
import tifffile as tf
from fnmatch import fnmatch
import tempfile
import numpy as np


Expand Down Expand Up @@ -34,64 +32,6 @@ def _yield_arrays(images, fpattern='*.tif'):
yield from _yield_arrays(item) # noqa


class TemporaryOTF(object):
"""Context manager to read OTF file or generate a temporary OTF from a PSF.
Normalizes the input PSF to always provide the path to an OTF file,
converting the PSF to a temporary file if necessary.
``self.path`` can be used within the context to get the filepath to
the temporary OTF filepath.
Args:
psf (str, np.ndarray): 3D PSF numpy array, or a filepath to a 3D PSF
or 2D complex OTF file.
**kwargs: optional keyword arguments will be passed to the :func:`pycudadecon.otf.makeotf` function
Note:
OTF files cannot currently be provided directly as 2D complex np.ndarrays
Raises:
ValueError: If the PSF/OTF is an unexpected type
NotImplementedError: if the PSF/OTF is a complex 2D numpy array
Example:
>>> with TemporaryOTF(psf, **kwargs) as otf:
print(otf.path)
/tmp/...
"""
def __init__(self, psf, **kwargs):
self.psf = psf
self.kwargs = kwargs

def __enter__(self):
if not is_otf(self.psf):
self.temp = tempfile.NamedTemporaryFile()
if isinstance(self.psf, np.ndarray):
with tempfile.NamedTemporaryFile() as tpsf:
tf.imsave(tpsf.name, self.psf)
makeotf(tpsf.name, self.temp.name, **self.kwargs)
elif isinstance(self.psf, str) and os.path.isfile(self.psf):
makeotf(self.psf, self.temp.name, **self.kwargs)
else:
raise ValueError('Did not expect PSF file as {}'
.format(type(self.psf)))
self.path = self.temp.name
elif is_otf(self.psf) and os.path.isfile(self.psf):
self.path = self.psf
elif is_otf(self.psf) and isinstance(self.psf, np.ndarray):
raise NotImplementedError('cannot yet handle OTFs as numpy arrays')
else:
raise ValueError('Unrecognized input for otf')
return self

def __exit__(self, typ, val, traceback):
try:
self.temp.close()
except Exception:
pass


def decon(images, psf, fpattern='*.tif', **kwargs):
"""Deconvolve an image or images with a PSF or OTF file
Expand Down
144 changes: 99 additions & 45 deletions pycudadecon/libcudawrapper.py
Original file line number Diff line number Diff line change
Expand Up @@ -113,9 +113,103 @@ def camcor(imstack):
return result


def rl_init(rawdata_shape, otfpath, dzdata=0.5, dxdata=0.1, dzpsf=0.1,
dxpsf=0.1, deskew=0, rotate=0, width=0, **kwargs):
"""Initialize GPU for deconvolution
prepares cuFFT plan for deconvolution with a given data shape and OTF.
Must be used prior to :func:`pycudadecon.rl_decon`
Args:
rawdata_shape (tuple, list): 3-tuple of data shape [nz, ny, nx]
otfpath (str): Path to OTF TIF
dzdata (float): Z-step size of data (default: {0.5})
dxdata (float): XY pixel size of data (default: {0.1})
dzpsf (float): Z-step size of the OTF (default: {0.1})
dxpsf (float): XY pixel size of the OTF (default: {0.1})
deskew (float): Deskew angle. If not 0.0 then deskewing will be
performed before deconv (default: {0})
rotate (float): Rotation angle; if not 0.0 then rotation will be
performed around Y axis after deconvolution (default: {0})
width (int): If deskewed, the output image's width (default: do not crop)
Example:
>>> rl_init(im.shape, otfpath)
>>> decon_result = rl_decon(im)
>>> rl_cleanup()
"""
nz, ny, nx = rawdata_shape
RL_interface_init(nx, ny, nz, dxdata, dzdata, dxpsf, dzpsf, deskew, rotate,
width, otfpath.encode())


def rl_decon(im, background=80, n_iters=10, shift=0, save_deskewed=False,
output_shape=None, napodize=15, nz_blend=0,
pad_val=0.0, dup_rev_z=False, **kwargs):
"""Perform Richardson Lucy Deconvolution
Performs actual deconvolution, after GPU has been initialized with
:func:`pycudadecon.rl_init`
Args:
im (np.ndarray): 3D image volume to deconvolve
background: User-supplied background to subtract. If 'auto', the
median value of the last Z plane will be used as background.
(default: {80})
n_iters: Number of decon iterations (default: {10})
shift: If deskewed, the output image's extra shift in X
(positive->left). (default: {0})
save_deskewed: Save deskewed raw data as well as deconvolution
result (default: {False})
output_shape: Specify the output shape after deskewing. Usually this
is unnecessary and will be autodetected. Mostly intended for
use within a :class:`pycudadecon.RLContext` context. (default: autodetect)
napodize: [description] (default: {15})
nz_blend: [description] (default: {0})
pad_val: [description] (default: {0.0})
dup_rev_z: [description] (default: {False})
"""
nz, ny, nx = im.shape
if output_shape is None:
output_shape = (get_output_nz(), get_output_ny(), get_output_nx())
else:
assert len(output_shape) == 3, 'Decon output shape must have length==3'
decon_result = np.empty(tuple(output_shape), dtype=np.float32)

if save_deskewed:
deskew_result = np.empty_like(decon_result)
else:
deskew_result = np.empty(1, dtype=np.float32)

# must be 16 bit going in
if not np.issubdtype(im.dtype, np.uint16):
im = im.astype(np.uint16)

if isinstance(background, str) and background == 'auto':
background = np.median(im[-1])

rescale = False # not sure if this works yet...

if not im.flags['C_CONTIGUOUS']:
im = np.ascontiguousarray(im)
RL_interface(im, nx, ny, nz, decon_result, deskew_result,
background, rescale, save_deskewed, n_iters, shift,
napodize, nz_blend, pad_val, dup_rev_z)

if save_deskewed:
return decon_result, deskew_result
else:
return decon_result


def quickDecon(im, otfpath, save_deskewed=False, **kwargs):
"""Perform deconvolution of im with otf at otfpath
Not currently used...
kwargs can be:
dxdata float
dzdata float
Expand All @@ -139,21 +233,15 @@ def quickDecon(im, otfpath, save_deskewed=False, **kwargs):
return decon_result


def rl_init(rawdata_shape, otfpath, dzdata=0.5, dxdata=0.1, dxpsf=0.1,
dzpsf=0.1, deskew=0, rotate=0, width=0, **kwargs):
nz, ny, nx = rawdata_shape
RL_interface_init(nx, ny, nz, dxdata, dzdata, dxpsf, dzpsf, deskew, rotate,
width, otfpath.encode())


class RLContext(object):
""" Context manager to setup the GPU for RL decon
Takes care of handing the OTF to the GPU and cleaning up after decon
Takes care of handing the OTF to the GPU, preparing a cuFFT plane,
and cleaning up after decon
EXAMPLE:
with RLContext(data.shape, otfpath, dz) as ctx:
return rl_decon(data, ctx.out_shape)
Example:
>>> with RLContext(data.shape, otfpath, dz) as ctx:
result = rl_decon(data, ctx.out_shape)
"""
def __init__(self, shape, otfpath, **kwargs):
Expand All @@ -173,37 +261,3 @@ def __exit__(self, typ, val, traceback):
# exit receives a tuple with any exceptions raised during processing
# if __exit__ returns True, exceptions will be supressed
rl_cleanup()


def rl_decon(im, background=80, n_iters=10, shift=0, save_deskewed=False,
rescale=False, output_shape=None, napodize=15, nz_blend=0,
pad_val=0.0, dup_rev_z=False, **kwargs):
nz, ny, nx = im.shape
if output_shape is None:
output_shape = (get_output_nz(), get_output_ny(), get_output_nx())
else:
assert len(output_shape) == 3, 'Decon output shape must have length==3'
decon_result = np.empty(tuple(output_shape), dtype=np.float32)

if save_deskewed:
deskew_result = np.empty_like(decon_result)
else:
deskew_result = np.empty(1, dtype=np.float32)

# must be 16 bit going in
if not np.issubdtype(im.dtype, np.uint16):
im = im.astype(np.uint16)

if isinstance(background, str) and background == 'auto':
background = np.median(im[-1])

if not im.flags['C_CONTIGUOUS']:
im = np.ascontiguousarray(im)
RL_interface(im, nx, ny, nz, decon_result, deskew_result,
background, rescale, save_deskewed, n_iters, shift,
napodize, nz_blend, pad_val, dup_rev_z)

if save_deskewed:
return decon_result, deskew_result
else:
return decon_result

0 comments on commit 4d7b86f

Please sign in to comment.