Repository: PyAbel/PyAbel (language: python)

Issue:
Adding Tikhonov regularization 
Hi,

I started using the PyAbel package and it works very well. I was just wondering if you considered adding Tikhonov regularization to the other Abel transforms available in the PyAbel package. My understanding is that currently only BASEX uses Tikhonov regularization. I previously applied Tikhonov regularization to Onion Peeling following the description by Daun et al.:
https://www.osapublishing.org/ao/abstract.cfm?uri=ao-45-19-4638

For the flame images I work with, this significantly improved the reconstruction (i.e., less noisy result). For me, BASEX works fine but for others Tikhonov regularization of a faster transform might be useful.

Thanks for the great package and regards,
Jochen
Need a help to save the output file 
In Example: Hansen–Law xenon i want to save the  output_text 
i ran the code and i seems to be working properly and i was able to save the output plot
is there way to save this dat file 


Discussion / hints:
Now BASEX is among the fastest transforms (see Fig. 5 in https://arxiv.org/abs/1902.09007). Its only problem is the relatively slow generation of the basis sets, but this needs to be done only once.

I need to read the Daun's article carefully, but from a quick glance it seems that they are using a Dasch-like approach, but with the regularized inversion of the forward-transform matrix instead of directly computing the inverse-transform matrix. In principle, adding this to the Dasch methods would not be difficult, but the resulting speed will be the same as in BASEX (except the faster basis generation).

Hi there! Are you basically trying to save a text file of the transformed image?

You can try adding something like:

```python
np.savetxt('abel_transform.txt', recon)
```

towards the end of the example script.

Here is a modified example script that saves the Abel transformed image as a text file using the line above:

https://gist.github.com/DanHickstein/1a63ac94621027784603ec59ce636854
I think that this issue has been addressed, so closing this issue. 

Fix patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 47a47a69..1ba6e85b 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -9,11 +9,12 @@ Unreleased
   nonintuitive behavior (PR #318, PR #319).
 * Avoid unnecessary calculations in transform.Transform() for the
   symmetry_axis=(0, 1) case (PR #324).
+* New method by Daun et al. and its extensions (PR #326).
 
 v0.8.4 (2021-04-15)
 -------------------
 * Added odd angular orders to tools.vmi.Distributions (PR #266).
-* Important! Some "center" functions/parameters are renamed to "origin" or
+* **Important!** Some "center" functions/parameters are renamed to "origin" or
   "method"; using old names still works but will print deprecation warnings,
   please update your code accordingly. Image origin is now always in the
   (row, column) format for consistency within PyAbel and with NumPy/SciPy; this
diff --git a/README.rst b/README.rst
index 12922bc6..533e6b96 100644
--- a/README.rst
+++ b/README.rst
@@ -29,23 +29,25 @@ Transform Methods
 
 The outcome of the numerical Abel transform depends on the exact method used. So far, PyAbel includes the following `transform methods <https://pyabel.readthedocs.io/en/latest/transform_methods.html>`__:
 
-1. ``basex`` - Gaussian basis set expansion of Dribinski and co-workers.
+1. ``basex`` – Gaussian basis set expansion of Dribinski and co-workers.
 
-2. ``hansenlaw`` - recursive method of Hansen and Law.
+2. ``hansenlaw`` – recursive method of Hansen and Law.
 
-3. ``direct`` - numerical integration of the analytical Abel transform equations.
+3. ``direct`` – numerical integration of the analytical Abel transform equations.
 
-4. ``two_point`` - the "two point" method of Dasch and co-workers.
+4. ``two_point`` – the "two point" method of Dasch and co-workers.
 
-5. ``three_point`` - the "three point" method of Dasch and co-workers.
+5. ``three_point`` – the "three point" method of Dasch and co-workers.
 
-6. ``onion_peeling`` - the "onion peeling" deconvolution method of Dasch and co-workers.
+6. ``onion_peeling`` – the "onion peeling" deconvolution method of Dasch and co-workers.
 
-7. ``onion_bordas`` - "onion peeling" or "back projection" method of Bordas *et al.* based on the MatLab code by Rallis and Wells *et al.*
+7. ``onion_bordas`` – "onion peeling" or "back projection" method of Bordas *et al.* based on the MatLab code by Rallis and Wells *et al.*
 
-8. ``linbasex`` - the 1D-spherical basis set expansion of Gerber *et al.*
+8. ``linbasex`` – the 1D-spherical basis set expansion of Gerber *et al.*
 
-9. ``rbasex`` - a pBasex-like method formulated in terms of radial distributions.
+9. ``rbasex`` – a pBasex-like method formulated in terms of radial distributions.
+
+10. ``daun`` – the regularized deconvolution method by Daun and co-workers, with additional capabilities.
 
 
 Installation
diff --git a/abel/benchmark.py b/abel/benchmark.py
index e5b61b51..9d414572 100644
--- a/abel/benchmark.py
+++ b/abel/benchmark.py
@@ -8,6 +8,7 @@
 import abel
 from . import basex
 from . import dasch
+from . import daun
 from . import direct
 from . import hansenlaw
 from . import linbasex
@@ -194,6 +195,7 @@ def __init__(self, n=[301, 501], select='all',
         # which methods need half and whole images
         need_half = frozenset([
             'basex',
+            'daun',
             'direct_C',
             'direct_Python',
             'hansenlaw',
@@ -398,6 +400,40 @@ def basex_var():
         # discard all caches
         basex.cache_cleanup()
 
+    @_skip(('bs', 'daun'),
+           (['inverse', 'forward'], 'daun'),
+           ('inverse', 'daun(var.reg.)'))
+    def _time_daun(self):
+        # benchmark the basis generation (default parameters)
+        def gen_basis():
+            daun.cache_cleanup()
+            daun.get_bs_cached(self.w)
+        self._benchmark('bs', 'daun',
+                        gen_basis)
+
+        # benchmark all transforms
+        for direction in ['inverse', 'forward']:  # (default first)
+            # cache the transform matrix
+            daun.get_bs_cached(self.w, direction=direction)
+            # benchmark the transform itself
+            self._benchmark(direction, 'daun',
+                            daun.daun_transform,
+                            self.half_image, direction=direction,
+                            verbose=False)
+            if direction == 'inverse':
+                # benchmark the transform with variable regularization
+                def daun_var():
+                    daun.daun_transform(self.half_image,
+                                        reg=1.0 + np.random.random(),
+                                        verbose=False)
+                self._benchmark('inverse', 'daun(var.reg.)',
+                                daun_var)
+            # discard the unneeded transform matrix
+            daun.cache_cleanup(direction)
+
+        # discard all caches
+        daun.cache_cleanup()
+
     @_skip((['inverse', 'forward'], 'direct_C'))
     def _time_direct_C(self):
         for direction in ['inverse', 'forward']:
diff --git a/abel/daun.py b/abel/daun.py
new file mode 100644
index 00000000..a958d466
--- /dev/null
+++ b/abel/daun.py
@@ -0,0 +1,461 @@
+# -*- coding: utf-8 -*-
+
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+import sys
+import os
+import glob
+
+import numpy as np
+from scipy.linalg import inv, toeplitz, solve_banded, solve_triangular
+from scipy.optimize import nnls
+
+
+def daun_transform(data, reg=0.0, degree=0, dr=1.0, direction='inverse',
+                   basis_dir=None, verbose=True):
+    """
+    Forward and inverse Abel transforms based on onion-peeling deconvolution
+    using Tikhonov regularization described in
+
+    K. J. Daun, K. A. Thomson, F. Liu, G. J. Smallwood,
+    "Deconvolution of axisymmetric flame properties using Tikhonov
+    regularization",
+    `Appl. Opt. 45, 4638–4646 (2006)
+    <https://doi.org/10.1364/AO.45.004638>`__.
+
+    with additional basis-function types and regularization methods.
+
+    This function operates on the “right side” of an image, that it, just one
+    half of a cylindrically symmetric image, with the axial pixels located in
+    the 0-th column.
+
+    Parameters
+    ----------
+    data : m × n numpy array
+        the image to be transformed.
+        ``data[:, 0]`` should correspond to the central column of the image.
+    reg : float or tuple or str
+        regularization for the inverse transform:
+
+        ``strength``:
+            same as ``('diff', strength)``
+        ``('diff', strength)``:
+            Tikhonov regularization using the first-order difference operator
+            (first-derivative approximation), as described in the original
+            article
+        ``('L2', strength)``:
+            Tikhonov regularization using the :math:`L_2` norm, like in
+            :ref:`BASEX`
+        ``('L2c', strength)``:
+            same as ``('L2', strength)``, but with an intensity correction
+            applied to compensate the drop near the symmetry axis
+        ``'nonneg'``:
+            non-negative least-squares solution.
+
+            `Warning: this regularization method is very slow, typically taking
+            up to a minute for a megapixel image.`
+    degree : int
+        degree of basis-function polynomials:
+
+        0:
+            rectangular functions (step-function approximation), corresponding
+            to “onion peeling” from the original article
+        1:
+            triangular functions (piecewise linear approximation)
+        2:
+            piecewise quadratic functions (smooth approximation)
+        3:
+            piecewise cubic functions (cubic-spline approximation)
+    dr : float
+        pixel size in the radial direction. This only affects the absolute
+        scaling of the transformed image.
+    direction : str: ``'forward'`` or ``'inverse'``
+        type of Abel transform to be performed
+    basis_dir : str, optional
+        path to the directory for saving / loading the basis set (potentially
+        useful only for **degree** = 3 and transform without regularization;
+        time savings in other cases are small and might be negated by the
+        disk-access overhead). If ``None`` (default), the basis set will not be
+        loaded from or saved to disk.
+    verbose : bool
+        determines whether progress report should be printed
+
+    Returns
+    -------
+    recon : m × n numpy array
+        the transformed (half) image
+    """
+    # make sure that the data has the right shape (1D must be converted to 2D)
+    # and type:
+    dim = len(data.shape)
+    data = np.atleast_2d(data).astype(float)
+    h, w = data.shape
+
+    if reg in [None, 0]:
+        reg_type, strength = None, 0
+    elif reg == 'nonneg':
+        reg_type, strength = reg, None
+    elif isinstance(reg, (float, int)):
+        reg_type, strength = 'diff', reg
+    elif np.ndim(reg) == 0 and reg != 'nonneg':
+        raise ValueError('Wrong regularization format "{}"'.format(reg))
+    else:
+        reg_type, strength = reg
+
+    # load the basis sets and compute the transform matrix
+    M = get_bs_cached(w, degree, reg_type, strength, direction, basis_dir,
+                      verbose)
+
+    if reg == 'nonneg':
+        if verbose:
+            print('Solving NNLS equations...')
+            sys.stdout.flush()
+        recon = np.empty_like(data)
+        # transform row by row
+        for i in range(h):
+            if verbose:
+                print('\r{}/{}'.format(1 + i, h), end='')
+                sys.stdout.flush()
+            recon[i] = nnls(M.T, data[i])[0]
+        if verbose:
+            print('\nDone!')
+    else:
+        # do the linear transform
+        if direction == 'inverse' and strength == 0 and degree != 3:
+            # (this is faster than general-purpose multiplication by inverse)
+            recon = solve_triangular(M.T, data.T).T  # (here M is forward)
+        else:
+            recon = data.dot(M)
+
+    # apply pixel scaling
+    if dr != 1.0:
+        if direction == 'forward':
+            recon *= dr
+        else:  # 'inverse'
+            recon /= dr
+
+    if dim == 1:
+        return recon[0]
+    else:
+        return recon
+
+
+# Caches and their parameters
+_bs = None  # basis set
+_bs_prm = None  # [size, degree]
+_tr = None  # inverse-transform matrix
+_tr_prm = None  # [size, type, strength]
+
+
+def get_bs_cached(n, degree=0, reg_type='diff', strength=0,
+                  direction='inverse', basis_dir=None, verbose=False):
+    """
+    Internal function.
+
+    Gets the basis set and calculates the necessary transform matrix (notice
+    that inverse direction with ``'nonneg'`` regularization, as well as with
+    **strength** = 0 for **degree** ≠ 3, gives the forward (triangular)
+    matrix, to be used in solvers).
+
+    Parameters
+    ----------
+    n : int
+        half-width of the image in pixels, must include the axial pixel
+    degree : int
+        polynomial degree for basis functions (0–3)
+    reg_type: None or str
+        regularization type (``None``, ``'diff'``, ``'L2'``, ``'L2c'``,
+        ``'nonneg'``)
+    strength : float
+        Tikhonov regularization parameter (for **reg_type** = ``'diff'`` and
+        ``'L2'``/``'L2c'``, ignored otherwise)
+    direction : str: ``'forward'`` or ``'inverse'``
+        type of Abel transform to be performed
+    basis_dir : str
+        path to the directory for saving / loading the basis set.
+        If ``None``, the basis sets will not be saved to disk.
+    verbose : bool
+        print some debug information
+
+    Returns
+    -------
+    M : n × n numpy array
+        matrix of the Abel transform (forward or inverse)
+    """
+    global _bs, _bs_prm, _tr, _tr_prm
+
+    bs_OK = (_bs is not None and  # cached
+             _bs_prm[1] == degree and  # correct degree
+             (_bs_prm[0] == n if degree == 3 else _bs_prm[0] >= n))  # good size
+    if not bs_OK:
+        # try to load
+        _bs = _load_bs(basis_dir, n, degree, verbose)
+        if _bs is None:
+            # generate and cache
+            _bs = _bs_daun(n, degree, verbose)
+            _save_bs(basis_dir, n, degree, _bs, verbose)
+            # (does nothing for basis_dir == None)
+        _bs_prm = [n, degree]
+        # reset cached inverse-transform matrix
+        _tr = None
+        _tr_prm = None
+
+    if direction == 'forward' or reg_type == 'nonneg':
+        return _bs[:n, :n]  # (cropping if too large)
+
+    if reg_type is None:
+        strength = 0
+
+    if strength == 0:
+        if _tr is None or _tr_prm[2] != 0:  # not cached or for strength != 0
+            # compute transform matrix (full-sized, just in case)
+            if degree == 3:
+                # general-purpose inverse of non-triangular
+                _tr = inv(_bs)
+            else:
+                # triangular matrix as is — will be used in solve_triangular,
+                # which is even faster than multiplication by cached inverse
+                _tr = _bs
+            _tr_prm = [_bs_prm[0], reg_type, strength]
+        return _tr[:n, :n]  # (cropping if too large — safe for triangular)
+
+    if _tr_prm != [n, reg_type, strength]:
+        # square of Tikhonov matrix
+        if reg_type == 'diff':
+            # of difference operator (approx. derivative operator)
+            LTL = toeplitz([2, -1] + [0] * (n - 2))
+            LTL[0, 0] = LTL[-1, -1] = 1
+        else:
+            # for L2 norm
+            LTL = np.eye(n)
+
+        # regularized inverse
+        # (transposed compared to the Daun article, since our data are in rows)
+        A = _bs[:n, :n]
+        _tr = A.T.dot(inv(A.dot(A.T) + strength * LTL))
+        if reg_type == 'L2c':
+            # apply correction: divide by regularized inverse of
+            # forward-transformed uniform distribution
+            _tr /= A.sum(axis=0).dot(_tr)
+        _tr_prm = [n, reg_type, strength]
+    return _tr
+
+
+def _load_bs(basis_dir, n, degree, verbose=False):
+    """
+    Internal function.
+
+    Try to load the basis set from the best suitable file.
+
+    Returns None on failure.
+    """
+    if basis_dir is None:
+        return None
+
+    file_mask = 'daun_basis_*_{}.npy'.format(degree)
+    best_file = None
+    best_size = np.inf
+    for f in glob.glob(os.path.join(basis_dir, file_mask)):
+        size = int(f.split('_')[-2])  # (from '...basis_<size>_<degree>.npy')
+        if n <= size < best_size:
+            best_size = size
+            best_file = f
+
+    if best_file is None:
+        return None
+
+    if verbose:
+        print('Loading basis set from', best_file)
+    try:
+        bs = np.load(best_file)
+    except ValueError:
+        print('Incompatible cached basis-set file!')
+        return None
+
+    if size > n:
+        bs = bs[:n, :n]
+        if verbose:
+            print('(cropped to {})'.format(n))
+
+    return bs
+
+
+def _save_bs(basis_dir, n, degree, bs, verbose=False):
+    """
+    Internal function.
+
+    Try to save the basis set.
+    """
+    if basis_dir is None:
+        return
+
+    file_name = 'daun_basis_{}_{}.npy'.format(n, degree)
+    if verbose:
+        print('Saving basis set to disk as', file_name)
+    np.save(os.path.join(basis_dir, file_name), bs)
+
+
+def _bs_daun(n, degree=0, verbose=False):
+    """
+    Internal function.
+
+    Generate the projected basis set.
+
+    Parameters
+    ----------
+    n : int
+        half-width of the image in pixels, must include the axial pixel
+    degree : int
+        polynomial degree for basis functions (0–3)
+
+    Returns
+    -------
+    A : n × n numpy array
+        coefficient matrix (transposed projected basis set)
+    """
+    # pixel coordinates
+    x = np.arange(float(n))
+    # (and common subexpressions for projections)
+    x2 = x**2
+    if degree > 0:
+        x2logx = x2 * np.log(x, np.zeros_like(x), where=x > 0)
+
+    # projections (Abel transforms) of given-degree basis functions
+    if degree == 0:
+        def p(j):
+            o = np.zeros(n)
+            o[:j + 1] += np.sqrt((j + 1/2)**2 - x2[:j + 1])
+            if j > 0:
+                o[:j] -= np.sqrt((j - 1/2)**2 - x2[:j])
+            return 2 * o
+    elif degree == 1:
+        def p(j):
+            def P(R):
+                y = np.sqrt(R**2 - x2[:R])
+                return y * R - x2[:R] * np.log(y + R)
+            o = np.zeros(n)
+            o[:j + 1] += P(j + 1)
+            o[:j] -= 2 * P(j)
+            o[j] += x2logx[j]
+            if j > 0:
+                o[:j - 1] += P(j - 1)
+                o[j - 1] -= x2logx[j - 1]
+            return o
+    elif degree == 2:
+        def p(j):
+            def P(R, a, b, c):
+                x2_R = x2[:int(R + 1/2)]
+                y = np.sqrt(R**2 - x2_R)
+                return y * (a * 2 + b * R + c * 4/3 * (R**2 / 2 + x2_R)) + \
+                       b * x2_R * np.log(y + R)
+            o = np.zeros(n)
+            o[:j + 1] += P(j + 1, 2 * (j + 1)**2, -4 * (j + 1), 2) - \
+                         P(j + 1/2, (2 * j + 1)**2, -4 * (2 * j + 1), 4)
+            if j > 0:
+                o[:j] += P(j - 1/2, (2 * j - 1)**2, -4 * (2 * j - 1), 4)
+                o[j] -= 4 * j * x2logx[j]
+                o[:j - 1] -= P(j - 1, 2 * (j - 1)**2, -4 * (j - 1), 2)
+                o[j - 1] += 4 * (j - 1) * x2logx[j - 1]
+            return o
+    elif degree == 3:
+        def p(j):
+            def P(R, a, b, c, d):
+                y = np.sqrt(R**2 - x2[:R])
+                return y * (a * 2 + (b + (c * 2/3 + d * R / 2) * R) * R +
+                            (d * 3/4 * R + c * 4/3) * x2[:R]) + \
+                       (b + d * 3/4 * x2[:R]) * x2[:R] * np.log(y + R)
+            o = np.zeros(n)
+            o[:j + 1] += P(j + 1,
+                           -j**2 * (2 * j + 3) + 1, 6 * j * (j + 1),
+                           -3 * (2 * j + 1),  2)
+            o[:j] += P(j, 4 * j**3, -12 * j**2, 12 * j, -4)
+            o[j] -= (6 * j * (j + 1) + 3/2 * x2[j]) * x2logx[j]
+            if j > 0:
+                o[:j - 1] -= P(j - 1,
+                               j**2 * (2 * j - 3) + 1, -6 * j * (j - 1),
+                               3 * (2 * j - 1), -2)
+                o[j - 1] += 6 * (j * (j - 1) + x2[j - 1] / 4) * x2logx[j - 1]
+            return o
+
+        # derivative basis functions
+        def q(j):
+            def P(R, a, b, c, d):
+                y = np.sqrt(R**2 - x2[:R])
+                return y * (a * 2 + (b + (c * 2/3 + d * R / 2) * R) * R +
+                            (d * 3/4 * R + c * 4/3) * x2[:R]) + \
+                       (b + d * 3/4 * x2[:R]) * x2[:R] * np.log(y + R)
+            o = np.zeros(n)
+            o[:j + 1] += P(j + 1,
+                           -j * (j * (j + 2) + 1), j * (3 * j + 4) + 1,
+                           -3 * j - 2, 1)
+            o[:j] += P(j, 4 * j**2, -8 * j, 4, 0)
+            o[j] -= (j * (3 * j + 4) + 1 + 3/4 * x2[j]) * x2logx[j]
+            if j > 0:
+                o[:j - 1] -= P(j - 1,
+                               -j * (j * (j - 2) + 1), j * (3 * j - 4) + 1,
+                               -3 * j + 2, 1)
+                o[j - 1] -= (j * (3 * j - 4) + 1 + 3/4 * x2[j - 1]) * \
+                            x2logx[j - 1]
+            return o
+    else:
+        raise ValueError('Wrong degree={} (must be 0, 1, 2 or 3).'.
+                         format(repr(degree)))
+
+    if verbose:
+        print('Generating basis projections for '
+              'n = {}, degree = {}...'.format(n, degree))
+
+    # fill the coefficient matrix
+    # (transposed compared to the Daun article, since our data are in rows)
+    A = np.empty((n, n))
+    for j in range(n):
+        A[j] = p(j)
+
+    if degree == 3:
+        # coefficient matrix for derivative functions
+        B = np.empty((n, n))
+        for j in range(n):
+            B[j] = q(j)
+        # solve for smooth derivative and modify A accordingly
+        C = solve_banded((1, 1), ([0] + [1] * (n - 2) + [0],
+                                  [4] * n,
+                                  [0] + [1] * (n - 2) + [0]),
+                         3 * B)[1:-1, 1:-1]
+        A[2:,  1:-1] += C
+        A[:-2, 1:-1] -= C
+
+    return A
+
+
+def cache_cleanup(select='all'):
+    """
+    Utility function.
+
+    Frees the memory caches created by :func:`get_bs_cached`.
+    This is usually pointless, but might be required after working
+    with very large images, if more RAM is needed for further tasks.
+
+    Parameters
+    ----------
+    select : str
+        selects which caches to clean:
+
+        ``'all'`` (default)
+            everything, including basis set
+        ``'inverse'``
+            only inverse transform
+
+    Returns
+    -------
+    None
+    """
+    global _bs, _bs_prm, _tr, _tr_prm
+
+    if select == 'all':
+        _bs = None
+        _bs_prm = None
+    _tr = None
+    _tr_prm = None
diff --git a/abel/rbasex.py b/abel/rbasex.py
index dd8892be..f444ccca 100644
--- a/abel/rbasex.py
+++ b/abel/rbasex.py
@@ -35,7 +35,7 @@
 #    hydroxymethyl radical (CH₂OH and CD₂OH)”,
 #   Ph.D. dissertation, University of Southern California, 2012.
 #   https://search.proquest.com/docview/1289069738
-#   http://digitallibrary.usc.edu/cdm/ref/collection/p15799coll3/id/112619
+#   https://digitallibrary.usc.edu/asset-management/2A3BF169XWB4
 #
 ###############################################################################
 
diff --git a/abel/tools/polynomial.py b/abel/tools/polynomial.py
index 03aeb915..31f364a0 100644
--- a/abel/tools/polynomial.py
+++ b/abel/tools/polynomial.py
@@ -31,7 +31,9 @@ class Polynomial(object):
         the function is defined as the polynomial on [**r_min**, **r_max**]
         and zero outside it;
         0 ≤ **r_min** < **r_max** ≲ **max r**
-        (**r_max** might exceed maximal **r**, but usually by < 1 pixel)
+        (**r_max** might exceed maximal **r**, but usually by < 1 pixel;
+        negative **r_min** or **r_max** are allowed for convenience but are
+        interpreted as 0)
     c: numpy array
         polynomial coefficients in order of increasing degree:
         [c₀, c₁, c₂] means c₀ + c₁ *r* + c₂ *r*\ ²
@@ -47,13 +49,24 @@ class Polynomial(object):
         at large r (and might improve numeric accuracy)
     """
     def __init__(self, r, r_min, r_max, c, r_0=0.0, s=1.0, reduced=False):
+        n = r.shape[0]
+
+        # trim negative r limits
+        if r_max <= 0:
+            # both func and abel must be zero everywhere
+            self.func = np.zeros(n)
+            self.abel = np.zeros(n)
+            return
+        if r_min < 0:
+            r_min = 0
+
         # remove zero high-order terms
         c = np.array(np.trim_zeros(c, 'b'), float)
         # if all coefficients are zero
         if len(c) == 0:
             # then both func and abel are also zero everywhere
-            self.func = np.zeros_like(r)
-            self.abel = self.func
+            self.func = np.zeros(n)
+            self.abel = np.zeros(n)
             return
         # polynomial degree
         K = len(c) - 1
@@ -83,7 +96,6 @@ def __init__(self, r, r_min, r_max, c, r_0=0.0, s=1.0, reduced=False):
         odd = np.any(c[1::2])
 
         # index limits
-        n = r.shape[0]
         i_min = np.searchsorted(r, r_min)
         i_max = np.searchsorted(r, r_max)
 
diff --git a/abel/transform.py b/abel/transform.py
index 910675bf..244833bf 100644
--- a/abel/transform.py
+++ b/abel/transform.py
@@ -11,6 +11,7 @@
 
 from . import basex
 from . import dasch
+from . import daun
 from . import direct
 from . import hansenlaw
 from . import linbasex
@@ -53,6 +54,9 @@ class Transform(object):
         ``basex``
             the Gaussian "basis set expansion" method of Dribinski et al.
             (2002).
+        ``daun``
+            the deconvolution method with Tikhonov regularization of Daun et
+            al. and its extensions.
         ``direct``
             a naive implementation of the analytical formula by Roman Yurchuk.
         ``hansenlaw``
@@ -231,6 +235,22 @@ class Transform(object):
         `Rev. Sci. Instrum. 73, 2634–2642 (2002)
         <https://doi.org/10.1063/1.1482156>`__.
 
+    ``daun`` *
+        Methods based on onion-peeling deconvolution using Tikhonov
+        regularization described in
+
+        K. J. Daun, K. A. Thomson, F. Liu, G. J. Smallwood,
+        "Deconvolution of axisymmetric flame properties using Tikhonov
+        regularization",
+        `Appl. Opt. 45, 4638–4646 (2006)
+        <https://doi.org/10.1364/AO.45.004638>`__.
+
+        In addition to the original implicit step-functions basis (“onion
+        peeling”) and the derivative regularization, linear and quadratic basis
+        functions are implemented, as well as the :math:`L_2`-norm Tikhonov
+        regularization (like in ``basex``) and non-negative least-squares
+        solution.
+
     ``direct``
         This method attempts a direct integration of the Abel-transform
         integral. It makes no assumptions about the data (apart from
@@ -314,8 +334,7 @@ class Transform(object):
         CD\ :sub:`2`\ OH)”,
         Ph.D. dissertation, University of Southern California, 2012
         (`ProQuest <https://search.proquest.com/docview/1289069738>`__,
-        `USC <http://digitallibrary.usc.edu/cdm/ref/collection/p15799coll3/id/
-        112619>`__).
+        `USC <https://digitallibrary.usc.edu/asset-management/2A3BF169XWB4>`__).
 
     ``three_point`` *
         The "Three Point" Abel transform method exploits the observation that
@@ -491,6 +510,7 @@ def _abel_transform_image_by_quadrant(self, **transform_options):
 
         abel_transform = {
             "basex": basex.basex_transform,
+            "daun": daun.daun_transform,
             "direct": direct.direct_transform,
             "hansenlaw": hansenlaw.hansenlaw_transform,
             "linbasex": linbasex.linbasex_transform,
diff --git a/doc/abel.rst b/doc/abel.rst
index ea808f46..1feca1c2 100644
--- a/doc/abel.rst
+++ b/doc/abel.rst
@@ -18,18 +18,26 @@ abel.basex module
     :undoc-members:
     :show-inheritance:
 
-abel.linbasex module
---------------------
+abel.dasch module
+-----------------
 
-.. automodule:: abel.linbasex
+.. automodule:: abel.dasch
     :members:
     :undoc-members:
     :show-inheritance:
 
-abel.rbasex module
+abel.daun module
+----------------
+
+.. automodule:: abel.daun
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+abel.direct module
 ------------------
 
-.. automodule:: abel.rbasex
+.. automodule:: abel.direct
     :members:
     :undoc-members:
     :show-inheritance:
@@ -42,10 +50,10 @@ abel.hansenlaw module
     :undoc-members:
     :show-inheritance:
 
-abel.dasch module
------------------
+abel.linbasex module
+--------------------
 
-.. automodule:: abel.dasch
+.. automodule:: abel.linbasex
     :members:
     :undoc-members:
     :show-inheritance:
@@ -58,10 +66,10 @@ abel.onion_bordas module
     :undoc-members:
     :show-inheritance:
 
-abel.direct module
+abel.rbasex module
 ------------------
 
-.. automodule:: abel.direct
+.. automodule:: abel.rbasex
     :members:
     :undoc-members:
     :show-inheritance:
diff --git a/doc/conf.py b/doc/conf.py
index a41f1e94..34391637 100755
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -174,7 +174,9 @@ def __getattr__(cls, name):
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
 # so a file named "default.css" will overwrite the builtin "default.css".
-# html_static_path = ['_static']
+html_static_path = ['css']
+
+html_css_files = ['custom.css']
 
 # Add any extra paths that contain custom files (such as robots.txt or
 # .htaccess) here, relative to this directory. These files are copied
diff --git a/doc/css/custom.css b/doc/css/custom.css
new file mode 100644
index 00000000..faa72f5f
--- /dev/null
+++ b/doc/css/custom.css
@@ -0,0 +1,22 @@
+/*
+ * Necessary parts from
+ * Sphinx stylesheet -- basic theme
+ * absent from sphinx_rtd_theme
+ * (see https://github.com/readthedocs/sphinx_rtd_theme/issues/301)
+ */
+
+/* -- math display ---------------------------------------------------------- */
+
+span.eqno {
+    float: right;
+}
+
+span.eqno a.headerlink {
+    position: absolute;
+    z-index: 1;
+    visibility: hidden;
+}
+
+div.math:hover a.headerlink {
+    visibility: visible;
+}
diff --git a/doc/example_daun_degree.rst b/doc/example_daun_degree.rst
new file mode 100644
index 00000000..733b050e
--- /dev/null
+++ b/doc/example_daun_degree.rst
@@ -0,0 +1,6 @@
+Example: Daun degree
+====================
+
+
+.. plot:: ../examples/example_daun_degree.py
+   :include-source:
diff --git a/doc/example_daun_reg.rst b/doc/example_daun_reg.rst
new file mode 100644
index 00000000..c99d4904
--- /dev/null
+++ b/doc/example_daun_reg.rst
@@ -0,0 +1,6 @@
+Example: Daun regularizations
+=============================
+
+
+.. plot:: ../examples/example_daun_reg.py
+   :include-source:
diff --git a/doc/examples.rst b/doc/examples.rst
index f307e73a..fc11c4bb 100644
--- a/doc/examples.rst
+++ b/doc/examples.rst
@@ -18,5 +18,7 @@ Contents:
    example_linbasex
    example_rbasex_block
    example_rbasex_multimass
+   example_daun_reg
+   example_daun_degree
    example_anisotropy_parameter
    example_circularize_image
diff --git a/doc/tools/polynomial.rst b/doc/tools/polynomial.rst
index 4b4935b2..da2658ba 100644
--- a/doc/tools/polynomial.rst
+++ b/doc/tools/polynomial.rst
@@ -61,11 +61,28 @@ so
     \int r^k \,dy = y \sum_{m=1}^k C_m r^m x^{k-m} + C_1 x^{k+1} \ln (y + r),
     \qquad (\text{summing over odd}\ m)
 
-with the same expressions for :math:`C_m`.
-
-These sums are computed using Horner's method in :math:`x`, which requires only
-:math:`x^2`, :math:`y` (see above), :math:`\ln (y + r)` (for polynomials with
-odd degrees), and powers of :math:`r` up to :math:`K`.
+with the same expressions for :math:`C_m`. For example, here are explicit
+formulas for several low degrees:
+
+    ========= =====================
+    :math:`k` :math:`\int r^k \,dy`
+    ========= =====================
+    0         :math:`y`
+    1         :math:`\frac12 r y + \frac12 x^2 \ln(y + r)`
+    2         :math:`\left(\frac13 r^2 + \frac23 x^2\right) y`
+    3         :math:`\left(\frac14 r^3 + \frac38 r x^2\right) y +
+              \frac38 x^4 \ln(y + r)`
+    4         :math:`\left(\frac15 r^4 + \frac4{15} r^2 x^2 +
+              \frac8{15} x^4\right) y`
+    5         :math:`\left(\frac16 r^5 + \frac5{24} r^3 x^2 +
+              \frac5{16} r x^4\right) y + \frac5{16} x^6 \ln(y + r)`
+    ...       :math:`\dots`
+    ========= =====================
+
+The sums over :math:`m` are computed using `Horner's method
+<https://en.wikipedia.org/wiki/Horner's_method>`__ in :math:`x`, which requires
+only :math:`x^2`, :math:`y` (see above), :math:`\ln (y + r)` (for polynomials
+with odd degrees), and powers of :math:`r` up to :math:`K`.
 
 The sum of the integrals, however, is computed by direct addition. In
 particular, this means that an attempt to use this method for high-degree
diff --git a/doc/transform_methods.rst b/doc/transform_methods.rst
index f2e06d63..9cc95a8a 100644
--- a/doc/transform_methods.rst
+++ b/doc/transform_methods.rst
@@ -27,6 +27,8 @@ Transform Methods
 ..
 ..     9. ``*`` The :doc:`rbasex <transform_methods/rbasex>` method is based on the pBasex method of Garcia et al. and basis functions developed by Ryazanov. Evaluates radial distributions of velocity-map images and transforms them to radial distributions of the reconstructed 3D distributions. This method makes additional assumptions about the symmetry of the data is not applicable to all situations! :meth:`abel.rbasrx.rbasex_transform`
 ..
+..     10. ``*`` The :doc:`daun <transform_methods/daun>` method of Daun et al. applies Tikhonov regularization to the onion-peeling deconvolution approach. The implementation is extended by several additional basis sets and regularization methods. :meth:`abel.daun.daun_transform`
+..
 ..     ``*`` Methods marked with an asterisk generate basis sets and allow these basis sets to be saved to the disk to speed up future transforms.
 
 
@@ -34,10 +36,11 @@ Contents:
 
 .. toctree::
    :maxdepth: 2
-   
+
    transform_methods/comparison
 
    transform_methods/basex
+   transform_methods/daun
    transform_methods/direct
    transform_methods/hansenlaw
    transform_methods/linbasex
@@ -46,4 +49,3 @@ Contents:
    transform_methods/rbasex
    transform_methods/three_point
    transform_methods/two_point
-   
\ No newline at end of file
diff --git a/doc/transform_methods/comparison.rst b/doc/transform_methods/comparison.rst
index e29d381d..2fcff3e8 100644
--- a/doc/transform_methods/comparison.rst
+++ b/doc/transform_methods/comparison.rst
@@ -67,7 +67,9 @@ Below is a list that describes the basic approach and characteristics of all the
 
 - ``onion_bordas`` -- The onion-peeling method of Bordas et al. [bordas1996]_ is a Python adaptation of the MATLAB implementation of Rallis et al. [rallis2014]_. While it is conceptually similar to the ``onion_peeling`` method, the numerical implementation is significantly different. This method is reasonably slow, and is therefore not recommended for general use. See :doc:`onion_bordas` and :meth:`abel.onion_bordas.onion_bordas_transform`
 
-- ``rbasex`` (F, S) --  the rBasex method is based on the pBasex method of Garcia et al. [garcia2004]_, using basis functions developed by Ryazanov [ryazanov2012]_. This method evaluates radial distributions of velocity-map images and transforms them to radial distributions of the reconstructed 3D distributions. Similar to ``linbasex``, the ``rbasex`` method makes additional assumptions about the symmetry of the data is not applicable to all situations. See :doc:`rbasex` and :meth:`abel.rbasex.rbasex_transform`.
+- ``rbasex`` (F, S) -- The rBasex method is based on the pBasex method of Garcia et al. [garcia2004]_, using basis functions developed by Ryazanov [ryazanov2012]_. This method evaluates radial distributions of velocity-map images and transforms them to radial distributions of the reconstructed 3D distributions. Similar to ``linbasex``, the ``rbasex`` method makes additional assumptions about the symmetry of the data is not applicable to all situations. See :doc:`rbasex` and :meth:`abel.rbasex.rbasex_transform`.
+
+- ``daun`` (F, S) -- The method by Daun et al. [daun2006]_ applies Tikhonov regularization to onion-peeling deconvolution. It is conceptually similar to "BASEX" (``basex``), but instead of :math:`L_2` regularization uses the first-order difference operator (approximating the derivative operator) as the Tikhonov matrix to suppress high-frequency oscillations, making the transform less sensitive to perturbations in the projected data. The PyAbel implementation also includes several extensions to the original method. First, in addition to the rectangular basis functions implied in onion peeling, explicit basis sets of piecewise polynomials up to 3rd degree (cubic splines) can be chosen. Second, the :math:`L_2` regularization (as in BASEX) is implemented for comparison. And most importantly, the non-negative least-squares solution to the deconvolution problem can be obtained, which produces meaningful results in situations where the transformed intensities must not be negative, and at the same time greatly reduces the baseline noise. See :doc:`daun` and :meth:`abel.daun.daun_transform`.
 
 
 Implementation
@@ -105,7 +107,7 @@ A similar ambiguity exists with regards to the left--right and top--bottom symme
 
 In the following performance benchmarks, left--right symmetry is assumed, because this is the most common benchmark presented in other studies ([rallis2014]_, [harrison2018]_). However, the image size is listed as the width of a square image. For example, :math:`n=513` corresponds to the time for the transformation of a 513×513-pixel image with the axis of symmetry located in the center. Since the Abel transform makes the assumption of cylindrical symmetry, both sides of the image are identical, and it is sufficient to perform the Abel transform on only one side of the image, or on an average of the two sides. So, to complete an Abel transform of a typical 513×513-pixel image, it is only necessary to perform the Abel transform on a 513×257-pixel array.
 
-Another fundamental question about real-world Abel transforms is whether negative values are allowed in the transform result. In most situations, negative values are not physical, and some implementations set all negative values to zero. In contrast, PyAbel allows negative values, which enables its use in situations where negative values are physically reasonable. Moreover, maintaining negative values keeps the transform methods linear and gives users the option to average, smooth, or fit images either before or after the Abel transform without causing a systematic error in the baseline. Suppression of negative values in a transformed image ``im`` can easily be achieved by executing ``im[im<0] = 0``.
+Another fundamental question about real-world Abel transforms is whether negative values are allowed in the transform result. In most situations, negative values are not physical, and some implementations set all negative values to zero. In contrast, PyAbel allows negative values, which enables its use in situations where negative values are physically reasonable. Moreover, maintaining negative values keeps the transform methods linear and gives users the option to average, smooth, or fit images either before or after the Abel transform without causing a systematic error in the baseline. Suppression of negative values in a transformed image ``im`` can easily be achieved by executing ``im[im<0] = 0``. On the other hand, the ``daun`` and ``rbasex`` methods offer optional non-linear regularization methods specifically designed to produce non-negative values without systematically shifting the baseline. It is recommended to use them instead of artificially zeroing negative values in situations where negative values are undesirable but the transform speed is not essential.
 
 
 Comparison of Transform Results
@@ -113,7 +115,7 @@ Comparison of Transform Results
 
 Since PyAbel incorporates numerous Abel-transform methods into the same interface, it is straightforward to directly compare the results. Consequently, a good approach is to simply try several (or all!) of the transform methods and see which produces the best results or performance for a specific application. Nevertheless, the following provides a brief comparison of the various transform methods in several cases. First, the methods are applied to a simple Gaussian function (for which an analytical Abel transform exists) in order to assess the accuracy of each transform method. Second, each method is applied to a "comb" function constructed of narrow peaks with noise added in order to closely examine the fundamental resolution of each method and how noise accumulates. Third, each method is used to provide the inverse Abel transform a high-resolution photoelectron-spectroscopy image in order to examine the ability of each method to handle real-world data.
 
-The Abel transform of a Gaussian is simply a Gaussian, which allows a comparison of each numerical transform method with the analytical result in the case of a one-dimensional (1D) Gaussian (:numref:`fig_gaussian`). As expected, each transform method exhibits a small discrepancy compared with the analytical result. However, as the number of pixels is increased, the agreement between the transform and the analytical result improves. Even with only 70 points (the case shown in :numref:`fig_gaussian`), all of the method produce reasonable agreement. While all methods show a systematic error as :math:`r` approaches zero, the ``basex``, ``three_point``, and ``onion_peeling`` methods seem to provide the best agreement with the analytical result. The direct methods show fairly good agreement with the analytical curve, which is a result of the "correction" discussed above. We note that the results from the ``direct_Python`` and the ``direct_C`` methods produce identical results to within a factor of :math:`10^{-9}`.
+The Abel transform of a Gaussian is simply a Gaussian, which allows a comparison of each numerical transform method with the analytical result in the case of a one-dimensional (1D) Gaussian (:numref:`fig_gaussian`). As expected, each transform method exhibits a small discrepancy compared with the analytical result. However, as the number of pixels is increased, the agreement between the transform and the analytical result improves. Even with only 70 points (the case shown in :numref:`fig_gaussian`), all of the method produce reasonable agreement. While all methods show a systematic error as :math:`r` approaches zero, the ``basex``, ``daun`` (especially with 3rd-degree basis functions), ``three_point``, and ``onion_peeling`` methods seem to provide the best agreement with the analytical result. The direct methods show fairly good agreement with the analytical curve, which is a result of the "correction" discussed above. We note that the results from the ``direct_Python`` and the ``direct_C`` methods produce identical results to within a factor of :math:`10^{-9}`.
 
 
 .. plot:: transform_methods/comparison/fig_gaussian/gaussian.py
@@ -124,10 +126,10 @@ The Abel transform of a Gaussian is simply a Gaussian, which allows a comparison
     :alt: gaussian
     :figclass: align-center
 
-    Comparison of inverse Abel-transform methods for a 1D Gaussian function with 70 points. All of the inverse Abel transform methods show reasonable agreement for the inverse Abel transform of a Gaussian function. The root-mean-square error (RMSE) for each method is listed in the figure legend. In the limit of many pixels, the error trends to zero. However, when a small number of pixels is used, systematic errors are seen, especially near the origin (:math:`r=0`). The error near the origin is more pronounced in some methods than others. The lowest error seen from the ``basex``, ``three_point``, and ``onion_peeling`` methods. The ``linbasex`` and ``rBasex`` methods are not included in this figure because they are not applicable to 1D functions.
+    Comparison of inverse Abel-transform methods for a 1D Gaussian function with 70 points. All of the inverse Abel transform methods show reasonable agreement for the inverse Abel transform of a Gaussian function. The root-mean-square error (RMSE) for each method is listed in the figure legend. In the limit of many pixels, the error trends to zero. However, when a small number of pixels is used, systematic errors are seen, especially near the origin (:math:`r=0`). The error near the origin is more pronounced in some methods than others. The lowest error seen from the ``basex``, ``daun``, ``three_point``, and ``onion_peeling`` methods. The ``daun`` method with degree=0 is identical to ``onion_peeling`` and with degree=2 is slightly better (RMSE=0.05%). The ``linbasex`` and ``rbasex`` methods are not included in this figure because they are not applicable to 1D functions.
 
 
-Applying the various transform methods to a synthetic "comb" function that consists of triangular peaks with one-pixel halfwidth -- the sharpest features representable on the pixel grid -- allows the fundamental resolution of each method to be visualized (:numref:`fig_comb`). In order to provide an understanding of how each method responds to noise, the function transformed in :numref:`fig_comb` also has uniformly distributed random noise added to each pixel. The figure reveals that some methods (``basex``, ``hansenlaw``), ``onion_peeling``, and ``two_point``) are capable of faithfully reproducing the sharpest features, while other methods (``direct``, ``onion_bordas``, and ``three_point``) provide some degree of smoothing. In general, the methods that provide the highest resolution also produce the highest noise, which is most obvious at low *r* values. The exception is the ``basex`` method using a moderate regularization factor (:numref:`fig_comb` b), which exhibits low noise near the center, while still displaying good resolution. Thus, it seems that experiments that benefit from an optimal balance of noise suppression and resolution would benefit from inverse Abel-transform methods that incorporate regularization.
+Applying the various transform methods to a synthetic "comb" function that consists of triangular peaks with one-pixel halfwidth -- the sharpest features representable on the pixel grid -- allows the fundamental resolution of each method to be visualized (:numref:`fig_comb`). In order to provide an understanding of how each method responds to noise, the function transformed in :numref:`fig_comb` also has uniformly distributed random noise added to each pixel. The figure reveals that some methods (``basex``, ``daun``, ``hansenlaw``, ``onion_peeling``, and ``two_point``) are capable of faithfully reproducing the sharpest features, while other methods (``direct``, ``onion_bordas``, and ``three_point``) provide some degree of smoothing. In general, the methods that provide the highest resolution also produce the highest noise, which is most obvious at low *r* values. The exception is the ``basex`` and ``daun`` methods using a moderate regularization factor (:numref:`fig_comb` b, c), which exhibit low noise near the center, while still displaying good resolution. The ``daun`` method with non-negativity regularization (:numref:`fig_comb` d), besides producing no negative values, significantly suppresses the baseline noise without affecting the sharp features. Thus, it seems that experiments that benefit from an optimal balance of noise suppression and resolution would benefit from inverse Abel-transform methods that incorporate regularization.
 
 
 .. plot:: transform_methods/comparison/fig_comb/comb.py
@@ -137,7 +139,7 @@ Applying the various transform methods to a synthetic "comb" function that consi
 .. figure:: comparison/fig_comb/comb.svg
     :figclass: align-center
 
-    Inverse Abel-transform methods applied to a synthetic "comb" function of one-pixel-width peaks with noise added. The gray line represents the analytical inverse Abel transform in the absence of noise. Some methods reproduce the height of the peaks, while other methods reduce noise while somewhat smoothing the peaks. The regularization in the ``basex`` method provides strong noise suppression near the origin, while maintaining peak height at higher values of :math:`r`.
+    Inverse Abel-transform methods applied to a synthetic "comb" function of one-pixel-width peaks with noise added. The gray line represents the analytical inverse Abel transform in the absence of noise. Some methods reproduce the height of the peaks, while other methods reduce noise while somewhat smoothing the peaks. The regularization in the ``basex`` and ``daun`` methods provides strong noise suppression near the origin, while maintaining peak height at higher values of :math:`r`. The ``daun`` method without regularization is identical to ``onion_peeling``, and its :math:`L_2` regularization is very similar to ``basex`` regularization.
 
 
 Applying the various inverse Abel-transform methods to an experimental photoelectron-spectroscopy image (:numref:`fig_experiment`) provides a comparison of how the noise in the reconstructed image depends on the transform method.
@@ -150,17 +152,17 @@ Applying the various inverse Abel-transform methods to an experimental photoelec
 .. figure:: comparison/fig_experiment/experiment.svg
     :figclass: align-center
 
-    Comparison of inverse Abel-transform methods applied to an experimental photoelectron velocity-map image. While all methods provide a faithful reconstruction of the experimental image, some of them cause a greater amplification of the noise present in the original image. The ``linbasex`` and ``rbasex`` methods models the image using a basis set of functions that vary slowly as a function of angle, which strongly reduces the high-frequency noise seen in the other transform methods. Besides the ``basex`` method with adjustable regularization, the ``direct`` and ``three_point`` methods seem particularly suited for providing a low-noise transform. This dataset is the photoelectron spectrum of O\ :sub:`2`\ :sup:`−` photodetachment using a 455 nm laser, as described by Van Duzor et al. [vanduzor2010]_.
+    Comparison of inverse Abel-transform methods applied to an experimental photoelectron velocity-map image. While all methods provide a faithful reconstruction of the experimental image, some of them cause a greater amplification of the noise present in the original image. The ``linbasex`` and ``rbasex`` methods models the image using a basis set of functions that vary slowly as a function of angle, which strongly reduces the high-frequency noise seen in the other transform methods. Besides the ``basex`` and ``daun`` method with regularization, the ``direct`` and ``three_point`` methods seem particularly suited for providing a low-noise transform. The ``daun`` method without regularization is identical to ``onion_peeling``, and its :math:`L_2` regularization is very similar to ``basex`` regularization. This dataset is the photoelectron spectrum of O\ :sub:`2`\ :sup:`−` photodetachment using a 455 nm laser, as described by Van Duzor et al. [vanduzor2010]_.
 
 
-To a first approximation, the results of all the transform methods look similar. The ``rbasex`` and ``linbasex`` methods produces the "smoothest" image, which is a result of the fact that it models the projection using functions fitted to the image, that vary only slowly as a function of angle. The ``basex`` method incorporates a user-adjustable Tikhonov regularization factor, which tends to suppress noise, especially near the symmetry axis. Here, we set the regularization factor to 200, which provides significant noise suppression while providing no noticeable broadening of the narrow features. When the regularization factor is set to zero, the ``basex`` method provides a transform that appears very similar to the ``onion_peeling`` method. For the other transform methods, the ``direct`` and ``three_point`` methods appear to have the strongest noise-filtering properties.
+To a first approximation, the results of all the transform methods look similar. The ``rbasex`` and ``linbasex`` methods produces the "smoothest" image, which is a result of the fact that it models the projection using functions fitted to the image, that vary only slowly as a function of angle. The ``basex`` and ``daun`` methods incorporate a user-adjustable Tikhonov regularization factor, which tends to suppress noise, especially near the symmetry axis. Here, we set the regularization factor to 200 for ``basex`` and 100 for ``daun``, which provides significant noise suppression without noticeable broadening of the narrow features. When the regularization factor is set to zero, the ``basex`` and ``daun`` methods provide a transform that appears very similar to the ``onion_peeling`` method. For the other transform methods, the ``direct`` and ``three_point`` methods appear to have the strongest noise-filtering properties.
 
 
 .. _fig_integration:
 .. figure:: comparison/fig_experiment/integration.svg
     :figclass: align-center
 
-    Comparison of photoelectron spectra obtained by angular integration of the transformed images shown in :numref:`fig_experiment`, corresponding to various inverse Abel-transform methods applied to the same experimental velocity-map image. a) Looking at the entire photoelectron speed distribution, all of the transform methods appear to produce similar results. b) Closely examining two of the peaks shows that all of the methods produce similar results, but that some methods produce broader peaks than others. c) Examining the small peaks in the low-energy region reveals that some methods accumulate somewhat more noise than others.
+    Comparison of photoelectron spectra obtained by angular integration of the transformed images shown in :numref:`fig_experiment`, corresponding to various inverse Abel-transform methods applied to the same experimental velocity-map image. a) Looking at the entire photoelectron speed distribution, all of the transform methods appear to produce similar results. b) Closely examining two of the peaks shows that all of the methods produce similar results, but that some methods produce broader peaks than others. c) Examining the small peaks in the low-energy region reveals that some methods accumulate somewhat more noise than others. Notice the absence on negative intensities in the ``daun`` method with non-negativity regularization and the corresponding suppression of baseline oscillations.
 
 
 :numref:`fig_integration` uses the same dataset as :numref:`fig_experiment`, but with an angular integration performed to show the 1D photoelectron spectrum. Good agreement is seen between most of the methods, even on a one-pixel level. Small but noticeable differences can be seen in the broadness of the peaks (:numref:`fig_integration`\ b). The ``hansenlaw``, ``onion_peeling`` and ``two_point`` methods show the sharpest peaks, suggesting that they provide enhanced ability to resolve sharp features. Of course, the differences between the methods are emphasized by the very high resolution of this dataset. In most cases, more pixels per peak yield a much better agreement between the transform methods. Interestingly, the ``linbasex`` method shows more baseline noise than the other methods. :numref:`fig_integration`\ c shows a close examination of the two lowest-energy peaks in the image. The methods that produce that sharpest peaks (``hansenlaw``, ``onion_peeling``, and ``two_point``) also exhibit somewhat more noise than the rest (except ``linbasex``).
@@ -176,7 +178,7 @@ For many applications of the inverse Abel transform, the speed at which transfor
 
 While PyAbel offers improvements to the raw computational efficiency of each transform method, it also provides improvements to the efficiency of the overall workflow, which are likely to provide a significant improvements for most applications. For example, since PyAbel provides a straightforward interface to switch between different transform methods (using :meth:`abel.transform.Transform`), a comparison of the results from each method can easily be made and the fastest method that produces acceptable results can be selected. Additionally, PyAbel provides fast algorithms for angular and radial integration, which can be the rate-limiting step for some data-processing workflows.
 
-In addition, when the computational efficiency of the various Abel transform methods is evaluated, a distinction must be made between those methods that can pre-compute, save, and re-use information for a specific image size (``basex``, ``linbasex``, ``onion_peeling``, ``rbasex``, ``three_point``, ``two_point``) and those that do not (``direct``, ``hansenlaw``, ``onion_bordas``). Often, the time required for the pre-computation is orders of magnitude longer than the time required to complete the transform. One solution to this problem is to pre-compute information for a specific image size and provide this data as part of the software. Indeed, the popular BASEX application [dribinski2002]_ includes a "basis set" for transforming 1000×1000-pixel images. While this approach relieves the end user of the computational cost of generating basis sets, it often means that the ideal basis set for efficiently transforming an image of a specific size is not available. Thus, "padding" is necessary for smaller images, resulting in increased computational time, while larger higher-resolution images must be downsampled or cropped.
+In addition, when the computational efficiency of the various Abel transform methods is evaluated, a distinction must be made between those methods that can pre-compute, save, and re-use information for a specific image size (``basex``, ``daun``, ``linbasex``, ``onion_peeling``, ``rbasex``, ``three_point``, ``two_point``) and those that do not (``direct``, ``hansenlaw``, ``onion_bordas``). Often, the time required for the pre-computation is orders of magnitude longer than the time required to complete the transform. One solution to this problem is to pre-compute information for a specific image size and provide this data as part of the software. Indeed, the popular BASEX application [dribinski2002]_ includes a "basis set" for transforming 1000×1000-pixel images. While this approach relieves the end user of the computational cost of generating basis sets, it often means that the ideal basis set for efficiently transforming an image of a specific size is not available. Thus, "padding" is necessary for smaller images, resulting in increased computational time, while larger higher-resolution images must be downsampled or cropped.
 
 PyAbel provides the ability to pre-compute information for any image size and cache it to disk for future use. Moreover, a cached basis set intended for transforming a larger image can be automatically cropped for use on a smaller image, avoiding unnecessary computations. The ``basex`` algorithm in PyAbel also includes the ability to extend a basis set intended for transforming a smaller image for use on a larger image. This allows the ideal basis set to be efficiently generated for an arbitrary image size.
 
@@ -195,7 +197,6 @@ Among the widely available options, the `Intel Math Kernel Library <https://en.w
 Another important issue for modern Intel CPUs is that they suffer a performance degradation when `denormal numbers <https://en.wikipedia.org/wiki/Denormal_number>`_ are encountered, which sometimes happens in the intermediate calculations even if the input and output are “normal”. In this case, configuring the CPU to treat denormals as zeros does help. There is no official way to achieve this in NumPy/SciPy, but a third-party module `daz <https://github.com/chainer/daz>`_ can be used for this purpose. At least some modern AMD CPUs are less or not affected by this issue, although it's always better to run the tests to make sure.
 
 
-
 Speed benchmarks
 """"""""""""""""
 
@@ -208,7 +209,8 @@ The :class:`abel.benchmark.AbelTiming` class provides the ability to benchmark t
 
 A comparison of the time required to complete an inverse Abel transform versus the width of a square image is presented in :numref:`fig_transform_time`. All methods are benchmarked using their default parameters, with the following exceptions:
 
-* **basex(var)** means “variable regularization”, that is changing the regularization parameter for each transformed image.
+* **basex(var)** and **daun(var)** mean “variable regularization”, that is changing the regularization parameter for each transformed image.
+* **daun(nonneg)** shows only the result of transforming the O\ :sub:`2`\ :sup:`−` image (see :numref:`fig_experiment`) with non-negativity regularization. Since the time needed for this non-linear transform strongly depends on the data, it is impossible to provide "universal" benchmarks; however, the general scaling is also expected to be roughly cubic.
 * **direct_C** and **direct_Python** correspond to the “direct” method using its C (Cython) and Python backends respectively.
 * **linbasex** and **rbasex** show whole-image (*n* × *n*) transforms, while all other methods show half-image (*n* rows, (*n* + 1)/2 columns) transforms.
 * **rbasex(None)** means no output-image creation (only the transformed radial distributions).
@@ -224,7 +226,7 @@ A comparison of the time required to complete an inverse Abel transform versus t
     Computational efficiency of inverse Abel-transform methods. The time to complete an inverse Abel transform increases with the size of the image. Most of the methods display a roughly :math:`n^3` scaling (dashed gray line). The ``basex``, ``onion_peeling``, ``three_point``, and ``two_point`` methods all rely on similar matrix-algebra operations as their rate-limiting step, and consequently exhibit identical performance for typical experimental image sizes.
 
 
-:numref:`fig_transform_time` reveals the computational scaling of each method as the image size is increased. At image sizes below :math:`n=100`, most of the transform methods exhibit a fairly flat relationship between image size and transform time, suggesting that the calculation is limited by the computational overhead. For image sizes of 1000 pixels and above, all the methods show a steep increase in transform time with increasing image size. A direct interpretation of the integral for the inverse Abel transform involves three nested loops, one over :math:`z`, one over :math:`r`, and one over :math:`y`, and we should expect :math:`n^3` scaling. Indeed, the ``direct_C`` and ``direct_Python`` methods scale as nearly :math:`n^3`. Several of the fastest methods (``basex``, ``onion_peeling``, ``three_point``, and ``two_point``) rely on matrix multiplication. These methods also scale roughly as :math:`n^3`, which is approximately the expected scaling for matrix-multiplication operations [coppersmith1990]_. For typical image sizes (~500--1000 pixels width), ``basex`` and the methods of Dasch [dasch1992]_ consistently out-perform other methods, often by several orders of magnitude. Interestingly, the ``hansenlaw`` and ``rbasex`` algorithms exhibits a nearly :math:`n^2` scaling and should outperform other algorithms for large image sizes. While the ``linbasex`` method does not provide the fastest transform, we note that it analytically provides the angular-integrated intensity and anisotropy parameters. Thus, if those parameters are desired outcomes -- as they often are during the analysis of photoelectron spectroscopy datasets -- then ``linbasex`` may provide an efficient analysis. The ``rbasex`` method also provides the intensity and anisotropy distributions directly. Moreover, if only these qualities are needed, without the transformed image, the transform can be completed faster and starts to outperform the fastest general-purpose methods for image sizes of ≳1000 pixels (extracting the desired distributions from the results of these methods requires additional time, not included in their plotted transform times).
+:numref:`fig_transform_time` reveals the computational scaling of each method as the image size is increased. At image sizes below :math:`n=100`, most of the transform methods exhibit a fairly flat relationship between image size and transform time, suggesting that the calculation is limited by the computational overhead. For image sizes of 1000 pixels and above, all the methods show a steep increase in transform time with increasing image size. A direct interpretation of the integral for the inverse Abel transform involves three nested loops, one over :math:`z`, one over :math:`r`, and one over :math:`y`, and we should expect :math:`n^3` scaling. Indeed, the ``direct_C`` and ``direct_Python`` methods scale as nearly :math:`n^3`. Several of the fastest methods (``basex``, ``onion_peeling``, ``three_point``, and ``two_point``) rely on matrix multiplication (or back substitution in case of ``daun``). These methods also scale roughly as :math:`n^3`, which is approximately the expected scaling for matrix-multiplication operations [coppersmith1990]_. For typical image sizes (~500--1000 pixels width), ``basex``, ``daun`` and the methods of Dasch [dasch1992]_ consistently out-perform other methods, often by several orders of magnitude. Interestingly, the ``hansenlaw`` and ``rbasex`` algorithms exhibits a nearly :math:`n^2` scaling and should outperform other algorithms for large image sizes. While the ``linbasex`` method does not provide the fastest transform, we note that it analytically provides the angular-integrated intensity and anisotropy parameters. Thus, if those parameters are desired outcomes -- as they often are during the analysis of photoelectron spectroscopy datasets -- then ``linbasex`` may provide an efficient analysis. The ``rbasex`` method also provides the intensity and anisotropy distributions directly. Moreover, if only these qualities are needed, without the transformed image, the transform can be completed faster and starts to outperform the fastest general-purpose methods for image sizes of ≳1000 pixels (extracting the desired distributions from the results of these methods requires additional time, not included in their plotted transform times).
 
 
 .. plot:: transform_methods/comparison/fig_benchmarks/throughput.py
@@ -247,8 +249,7 @@ A comparison of the time required to complete an inverse Abel transform versus t
     Computational efficiency of the basis-set generation calculation.
 
 
-
-The ``basex``, ``onion_peeling``, ``three_point``, and ``two_point`` methods run much faster if appropriately sized basis sets have been pre-calculated. For the ``basex`` method, the time for this pre-calculation is orders of magnitude longer than the transform time (:numref:`fig_btime`). For the Dasch methods (``onion_peeling``, ``three_point``, and ``two_point``), the pre-calculation is significantly longer than the transform time for image sizes smaller than 2000 pixels. For larger image sizes, the pre-calculation of the basis sets approaches the same speed as the transform itself. In particular, for the ``two_point`` method, the pre-calculation of the basis sets actually becomes faster than the image transform for *n* ≳ 4000. For the ``linbasex`` method, the pre-calculation of the basis sets is consistently faster than the transform itself, suggesting that the pre-calculation of basis sets isn't necessary for this method.
+The ``basex``, ``onion_peeling``, ``three_point``, and ``two_point`` methods run much faster if appropriately sized basis sets have been pre-calculated. For the ``basex`` method, the time for this pre-calculation is orders of magnitude longer than the transform time (:numref:`fig_btime`). For the Dasch methods (``onion_peeling``, ``three_point``, and ``two_point``), the pre-calculation is significantly longer than the transform time for image sizes smaller than 2000 pixels. For larger image sizes, the pre-calculation of the basis sets approaches the same speed as the transform itself. In particular, for the ``two_point`` method, the pre-calculation of the basis sets actually becomes faster than the image transform for *n* ≳ 4000. For the ``daun`` and ``linbasex`` methods, the pre-calculation of the basis sets is consistently faster than the transform itself, suggesting that the pre-calculation of basis sets isn't necessary for these methods.
 
 
 Conclusion
@@ -259,6 +260,7 @@ The various Abel-transform methods in PyAbel provide advantages for different si
 Methods recommended for general-purpose Abel transforms:
 
 * ``basex``
+* ``daun``
 * ``hansenlaw``
 * ``three-point``
 * ``two-point``
@@ -288,6 +290,8 @@ References
 
 .. [dasch1992] Cameron J. Dasch. One-dimensional tomography: a comparison of Abel, onion-peeling, and filtered backprojection methods. Appl. Opt., **31**, 1146, 1992. DOI: `10.1364/AO.31.001146 <https://doi.org/10.1364/AO.31.001146>`_.
 
+.. [daun2006] Kyle J. Daun, Kevin A. Thomson, Fengshan Liu, Fengshan J. Smallwood, Deconvolution of axisymmetric flame properties using Tikhonov regularization. Appl. Opt., **45**, 4638, 2006. DOI: `10.1364/AO.45.004638 <https://doi.org/10.1364/AO.45.004638>`_.
+
 .. [demicheli2017] Enrico De Micheli. A fast algorithm for the inversion of Abel’s transform. Appl. Math. Comput., **301**, 12, 2017. DOI: `10.1016/j.amc.2016.12.009 <https://doi.org/10.1016/j.amc.2016.12.009>`_.
 
 .. [dick2014] Bernhard Dick. Inverting ion images without Abel inversion: maximum entropy reconstruction of velocity maps. Phys. Chem. Chem. Phys., **16**, 570, 2014. DOI: `10.1039/C3CP53673D <http://doi.org/10.1039/C3CP53673D>`_.
diff --git a/doc/transform_methods/comparison/fig_benchmarks/basis_time.py b/doc/transform_methods/comparison/fig_benchmarks/basis_time.py
index 4b221fa8..c920ef33 100644
--- a/doc/transform_methods/comparison/fig_benchmarks/basis_time.py
+++ b/doc/transform_methods/comparison/fig_benchmarks/basis_time.py
@@ -4,8 +4,10 @@
 from matplotlib.ticker import LogLocator
 
 transforms = [
-    ("basex",         '#880000', {}),
-    ("basex(var)",    '#880000', {'mfc': 'w'}),
+    ("basex",         '#006600', {}),
+    ("basex(var)",    '#006600', {'mfc': 'w'}),
+    ("daun",          '#880000', {}),
+    ("daun(var)",     '#880000', {'mfc': 'w'}),
     ("direct_C",      '#EE0000', {}),
     ("direct_Python", '#EE0000', {'mfc': 'w'}),
     ("hansenlaw",     '#CCAA00', {}),
diff --git a/doc/transform_methods/comparison/fig_benchmarks/benchmarks.rst b/doc/transform_methods/comparison/fig_benchmarks/benchmarks.rst
index 32e20c14..aae97a6d 100644
--- a/doc/transform_methods/comparison/fig_benchmarks/benchmarks.rst
+++ b/doc/transform_methods/comparison/fig_benchmarks/benchmarks.rst
@@ -8,7 +8,7 @@ Speed benchmarks
 
 All methods are benchmarked using their default parameters, with the following exceptions:
 
-* **basex(var)** means “variable regularization”, that is changing the regularization parameter for each transformed image.
+* **basex(var)** and **daun(var)** mean “variable regularization”, that is changing the regularization parameter for each transformed image.
 * **direct_C** and **direct_Python** correspond to the “direct” method using its C (Cython) and Python backends respectively.
 * **linbasex** and **rbasex** show whole-image (*n* × *n*) transforms, while all other methods show half-image (*n* rows, (*n* + 1)/2 columns) transforms.
 * **rbasex(None)** means no output-image creation (only the transformed radial distributions).
diff --git a/doc/transform_methods/comparison/fig_benchmarks/i7-6700_Linux/daun(var).dat b/doc/transform_methods/comparison/fig_benchmarks/i7-6700_Linux/daun(var).dat
new file mode 100644
index 00000000..38290ccf
--- /dev/null
+++ b/doc/transform_methods/comparison/fig_benchmarks/i7-6700_Linux/daun(var).dat
@@ -0,0 +1,15 @@
+# daun(var)       iabel(ms)           
+    5             0.04809
+    9             0.04917
+   17             0.05251
+   33             0.06064
+   65             0.09011
+  129             0.46200
+  257             0.96465
+  513             3.34389
+ 1025            16.89925
+ 2049            92.87464
+ 4097           654.97368
+ 8193          5630.99326
+16385         44244.44146
+32769        358215.71806
diff --git a/doc/transform_methods/comparison/fig_benchmarks/i7-6700_Linux/daun.dat b/doc/transform_methods/comparison/fig_benchmarks/i7-6700_Linux/daun.dat
new file mode 100644
index 00000000..71360c16
--- /dev/null
+++ b/doc/transform_methods/comparison/fig_benchmarks/i7-6700_Linux/daun.dat
@@ -0,0 +1,15 @@
+# daun            iabel(ms)           basis(ms)
+    5             0.01801             0.02375
+    9             0.03927             0.03721
+   17             0.04102             0.06304
+   33             0.04390             0.11677
+   65             0.05894             0.22311
+  129             0.08766             0.44174
+  257             0.20657             0.92827
+  513             0.84640             1.98533
+ 1025             4.85983             4.27869
+ 2049            35.90545             9.95490
+ 4097           201.18627            27.41542
+ 8193          1272.86124            76.51444
+16385          9042.48185           242.99418
+32769         69927.85735           857.24778
diff --git a/doc/transform_methods/comparison/fig_benchmarks/i7-9700_Linux/daun(var).dat b/doc/transform_methods/comparison/fig_benchmarks/i7-9700_Linux/daun(var).dat
new file mode 100644
index 00000000..2daf071a
--- /dev/null
+++ b/doc/transform_methods/comparison/fig_benchmarks/i7-9700_Linux/daun(var).dat
@@ -0,0 +1,15 @@
+# daun(var)       iabel(ms)           
+    5             0.04662
+    9             0.04763
+   17             0.04966
+   33             0.05632
+   65             0.08008
+  129             0.27905
+  257             0.60093
+  513             2.05494
+ 1025             9.34317
+ 2049            60.91675
+ 4097           404.64816
+ 8193          3576.25377
+16385         24801.12200
+32769        186857.49893
diff --git a/doc/transform_methods/comparison/fig_benchmarks/i7-9700_Linux/daun.dat b/doc/transform_methods/comparison/fig_benchmarks/i7-9700_Linux/daun.dat
new file mode 100644
index 00000000..b6791c9d
--- /dev/null
+++ b/doc/transform_methods/comparison/fig_benchmarks/i7-9700_Linux/daun.dat
@@ -0,0 +1,15 @@
+# daun            iabel(ms)           basis(ms)
+    5             0.01711             0.02373
+    9             0.01955             0.03707
+   17             0.02057             0.06232
+   33             0.02293             0.11695
+   65             0.02886             0.22029
+  129             0.04329             0.43562
+  257             0.10706             0.90285
+  513             0.42341             1.89876
+ 1025             4.23736             4.05734
+ 2049            19.78891             9.56412
+ 4097           109.93354            27.52181
+ 8193           727.86471            74.85047
+16385          4903.41616           227.30092
+32769         37612.45776           818.95221
diff --git a/doc/transform_methods/comparison/fig_benchmarks/pyabel-plot.py b/doc/transform_methods/comparison/fig_benchmarks/pyabel-plot.py
index 82ddfeb6..4bdbae4a 100644
--- a/doc/transform_methods/comparison/fig_benchmarks/pyabel-plot.py
+++ b/doc/transform_methods/comparison/fig_benchmarks/pyabel-plot.py
@@ -5,8 +5,10 @@
 benchmark_dir = 'working'
 
 transforms = [
-    ("basex",         '#880000', {}),
-    ("basex(var)",    '#880000', {'mfc': 'w'}),
+    ("basex",         '#006600', {}),
+    ("basex(var)",    '#006600', {'mfc': 'w'}),
+    ("daun",          '#880000', {}),
+    ("daun(var)",     '#880000', {'mfc': 'w'}),
     ("direct_C",      '#EE0000', {}),
     ("direct_Python", '#EE0000', {'mfc': 'w'}),
     ("hansenlaw",     '#CCAA00', {}),
diff --git a/doc/transform_methods/comparison/fig_benchmarks/pyabel-timing.py b/doc/transform_methods/comparison/fig_benchmarks/pyabel-timing.py
index 30ecfeba..3d6b6a19 100644
--- a/doc/transform_methods/comparison/fig_benchmarks/pyabel-timing.py
+++ b/doc/transform_methods/comparison/fig_benchmarks/pyabel-timing.py
@@ -56,8 +56,8 @@ def bench(method, repeats, min_time, max_time, benchmark_dir, append=True):
     header = f'# {method:15s} {"iabel(ms)":20s}'
     if has_basis:
         header += 'basis(ms)'
-    if method in ['basex', 'rbasex']:  # special cases
-        if method == 'basex':
+    if method in ['basex', 'daun', 'rbasex']:  # special cases
+        if method in ['basex', 'daun']:
             varnone = 'var'
             varlab = 'var.reg.'
         else:
@@ -89,8 +89,8 @@ def bench(method, repeats, min_time, max_time, benchmark_dir, append=True):
             print(t_iabel, end='')
             fp.write(t_iabel)
 
-            if method in ['basex', 'rbasex']:  # special cases
-                if method == 'basex':
+            if method in ['basex', 'daun', 'rbasex']:  # special cases
+                if method in ['basex', 'daun']:
                     varlab = method+'(var.reg.)'
                 else:
                     varlab = method+'(None)'
@@ -106,7 +106,7 @@ def bench(method, repeats, min_time, max_time, benchmark_dir, append=True):
             print(t_bs)
             fp.write(t_bs+'\n')
 
-    if method in ['basex', 'rbasex']:
+    if method in ['basex', 'daun', 'rbasex']:
         fpvar.close()
 
     # spacing between methods
@@ -172,6 +172,7 @@ def bench(method, repeats, min_time, max_time, benchmark_dir, append=True):
   'three_point',
   'onion_peeling',
   'hansenlaw',
+  'daun',
   'basex',
   'rbasex',
   'direct_C',
diff --git a/doc/transform_methods/comparison/fig_benchmarks/throughput.py b/doc/transform_methods/comparison/fig_benchmarks/throughput.py
index 5bec5673..fedf9cda 100644
--- a/doc/transform_methods/comparison/fig_benchmarks/throughput.py
+++ b/doc/transform_methods/comparison/fig_benchmarks/throughput.py
@@ -2,11 +2,11 @@
 import numpy as np
 import matplotlib.pyplot as plt
 
-directory = 'benchmarks_i7-9700_Linux_5.4.0-26-generic'
-
 transforms = [
-  ("basex",         '#880000', {}),
-  ("basex(var)",    '#880000', {'mfc': 'w'}),
+  ("basex",         '#006600', {}),
+  ("basex(var)",    '#006600', {'mfc': 'w'}),
+  ("daun",          '#880000', {}),
+  ("daun(var)",     '#880000', {'mfc': 'w'}),
   ("direct_C",      '#EE0000', {}),
   ("direct_Python", '#EE0000', {'mfc': 'w'}),
   ("hansenlaw",     '#CCAA00', {}),
@@ -43,6 +43,10 @@ def plot(directory, xlim, ylim, va):
         t = times[1] * 1e-3  # in ms
         plt.plot(n, n**2 / t, 'o-', label=meth, ms=ms, **pargs)
 
+        # add an empty entry to end column 1 for more logical grouping
+        if meth == 'daun(var)':
+            plt.plot(np.NaN, np.NaN, 'o-', color='none', label=' ')
+
     plt.xlabel('Image size ($n$, pixels)')
     plt.xscale('log')
     plt.xlim(xlim)
@@ -70,5 +74,5 @@ def plot(directory, xlim, ylim, va):
 
 
 if __name__ == "__main__":
-    plot('i7-9700_Linux', xlim=(5, 1e5), ylim=(4e3, 2e9), va='top')
+    plot('i7-9700_Linux', xlim=(5, 1e5), ylim=(3e3, 2e9), va='top')
     plt.savefig('throughput.svg')
diff --git a/doc/transform_methods/comparison/fig_benchmarks/transform_time.py b/doc/transform_methods/comparison/fig_benchmarks/transform_time.py
index cb3a3bb0..2d12a6e4 100644
--- a/doc/transform_methods/comparison/fig_benchmarks/transform_time.py
+++ b/doc/transform_methods/comparison/fig_benchmarks/transform_time.py
@@ -4,8 +4,10 @@
 from matplotlib.ticker import LogLocator
 
 transforms = [
-  ("basex",         '#880000', {}),
-  ("basex(var)",    '#880000', {'mfc': 'w'}),
+  ("basex",         '#006600', {}),
+  ("basex(var)",    '#006600', {'mfc': 'w'}),
+  ("daun",          '#880000', {}),
+  ("daun(var)",     '#880000', {'mfc': 'w'}),
   ("direct_C",      '#EE0000', {}),
   ("direct_Python", '#EE0000', {'mfc': 'w'}),
   ("hansenlaw",     '#CCAA00', {}),
@@ -65,6 +67,15 @@ def plot(directory, xlim, ylim, linex):
         t = times[1] * 1e-3  # in ms
         plt.plot(n, t, 'o-', label=meth, ms=ms, **pargs)
 
+    # Daun with reg='nonneg' for the O2-ANU1024 example
+    if __name__ == "__main__":
+        x = 1023
+        y = 46  # done on i5-4278U, please replace!
+        plt.plot(x, y, 's', c='#880000', mfc='w')
+        plt.annotate('daun(nonneg)\n\n', (x, y), ha='right', va='center',
+                     xytext=(5, 0), textcoords='offset points')
+        #                    5 = fontsize / 2
+
     plt.legend()
 
     # plt.show()
diff --git a/doc/transform_methods/comparison/fig_comb/comb.py b/doc/transform_methods/comparison/fig_comb/comb.py
index 2ba52772..1fdf955d 100644
--- a/doc/transform_methods/comparison/fig_comb/comb.py
+++ b/doc/transform_methods/comparison/fig_comb/comb.py
@@ -28,8 +28,10 @@ def peak(i):
 func = comb.abel + np.random.random(comb.abel.size)*1.2
 
 transforms = [
-  ("basex",          abel.basex.basex_transform,               '#880000'),
-  ("basex (reg=10)", abel.basex.basex_transform,               '#880000'),
+  ("basex",          abel.basex.basex_transform,               '#006600'),
+  ("basex (reg=10)", abel.basex.basex_transform,               '#006600'),
+  ("daun (reg=5)",   abel.daun.daun_transform,                 '#880000'),
+  ("daun (nonneg)",  abel.daun.daun_transform,                 '#880000'),
   ("direct",         abel.direct.direct_transform,             '#EE0000'),
   ("hansenlaw",      abel.hansenlaw.hansenlaw_transform,       '#CCAA00'),
   ("onion_bordas",   abel.onion_bordas.onion_bordas_transform, '#00AA00'),
@@ -40,7 +42,7 @@ def peak(i):
 
 ntrans = len(transforms)  # number of transforms
 
-fig, axs = plt.subplots(ntrans, 1, figsize=(5, 9),
+fig, axs = plt.subplots(ntrans, 1, figsize=(5, 10),
                         sharex=True, sharey=True)
 
 
@@ -53,7 +55,10 @@ def mysum(x,  dx=1, axis=1):
                                                           transforms)):
     print(label)
     if 'reg' in label:
-        targs = dict(reg=10)
+        reg = label[label.rfind('=') + 1:-1]
+        targs = dict(reg=float(reg))
+    elif 'nonneg' in label:
+        targs = dict(reg='nonneg')
     else:
         targs = dict()
 
@@ -69,10 +74,10 @@ def place_letter(letter, ax, color='k', offset=(0, 0)):
                 color=color, ha='left', va='top', weight='bold')
 
 
-for ax, letter in zip(axs.ravel(), 'abcdefgh'):
+for ax, letter in zip(axs.ravel(), 'abcdefghij'):
     ax.legend(loc='upper right', frameon=False, borderaxespad=0)
     ax.set_xlim(0, 60)
-    ax.set_ylim(-0.2, 1.3)
+    ax.set_ylim(-0.2, 1.4)
     ax.set_yticks([0, 0.5, 1])
 
     ax.grid(alpha=0.2)
@@ -80,7 +85,7 @@ def place_letter(letter, ax, color='k', offset=(0, 0)):
     place_letter(letter+')', ax)
 
 axs[-1].set_xlabel('$r$ (pixels)')
-axs[3].set_ylabel('$z$')
+axs[4].set_ylabel('$z$')
 
 
 fig.tight_layout(pad=0)
diff --git a/doc/transform_methods/comparison/fig_experiment/experiment.py b/doc/transform_methods/comparison/fig_experiment/experiment.py
index 0317e5eb..71c16698 100644
--- a/doc/transform_methods/comparison/fig_experiment/experiment.py
+++ b/doc/transform_methods/comparison/fig_experiment/experiment.py
@@ -4,15 +4,17 @@
 import bz2
 
 transforms = [
-  ("basex",         abel.basex.basex_transform,               '#880000'),
-  ("direct",        abel.direct.direct_transform,             '#EE0000'),
-  ("hansenlaw",     abel.hansenlaw.hansenlaw_transform,       '#CCAA00'),
-  ("onion_bordas",  abel.onion_bordas.onion_bordas_transform, '#00AA00'),
-  ("onion_peeling", abel.dasch.onion_peeling_transform,       '#00CCFF'),
-  ("three_point",   abel.dasch.three_point_transform,         '#0000FF'),
-  ("two_point",     abel.dasch.two_point_transform,           '#CC00FF'),
-  ("linbasex",      abel.linbasex.linbasex_transform,         '#AAAAAA'),
-  ("rbasex",        abel.rbasex.rbasex_transform,             '#AACC00'),
+  ("basex (reg=200)", abel.basex.basex_transform,               '#006600'),
+  ("daun (reg=100)",  abel.daun.daun_transform,                 '#880000'),
+  ("daun (nonneg)",   abel.daun.daun_transform,                 '#880000'),
+  ("direct",          abel.direct.direct_transform,             '#EE0000'),
+  ("hansenlaw",       abel.hansenlaw.hansenlaw_transform,       '#CCAA00'),
+  ("onion_bordas",    abel.onion_bordas.onion_bordas_transform, '#00AA00'),
+  ("onion_peeling",   abel.dasch.onion_peeling_transform,       '#00CCFF'),
+  ("three_point",     abel.dasch.three_point_transform,         '#0000FF'),
+  ("two_point",       abel.dasch.two_point_transform,           '#CC00FF'),
+  ("linbasex",        abel.linbasex.linbasex_transform,         '#AAAAAA'),
+  ("rbasex",          abel.rbasex.rbasex_transform,             '#AACC00'),
 ]
 
 infile = bz2.BZ2File('../../../../examples/data/O2-ANU1024.txt.bz2')
@@ -26,25 +28,29 @@
 
 h, w = np.shape(IM)
 
-fig, axs = plt.subplots(3, 3, figsize=(7, 6.5), sharex=True, sharey=True)
-fig1, axs1 = plt.subplots(3, 1, figsize=(6, 7))
+fig, axs = plt.subplots(4, 3, figsize=(7, 8.6), sharex=True, sharey=True)
+fig1, axs1 = plt.subplots(3, 1, figsize=(7, 8))
 
 for num, (ax, (label, transFunc, color), letter) in enumerate(zip(axs.ravel(),
                                                               transforms,
-                                                              'abcdefghij')):
+                                                              'abcdefghijk')):
     print(label)
 
-    if label == 'linbasex':
-        targs = dict(proj_angles=np.arange(0, np.pi, np.pi/10))
-    else:
-        targs = dict()
+    targs = dict()
+    pargs = dict(lw=1, color=color, zorder=-num)
 
-    if label == 'basex':
+    if label == 'basex (reg=200)':
         targs = dict(reg=200)
+    elif label == 'daun (reg=100)':
+        targs = dict(reg=100)
+    elif label == 'daun (nonneg)':
+        targs = dict(reg='nonneg')
+        pargs['ls'] = '--'
+    elif label == 'linbasex':
+        targs = dict(proj_angles=np.arange(0, np.pi, np.pi/10))
 
     if label == 'rbasex':
         trans = transFunc(Q0, direction="inverse", origin='ll')[0]
-
     else:
         trans = transFunc(Q0, direction="inverse", **targs)
 
@@ -64,14 +70,12 @@
                  x=0.05, y=0.93, ha='left', va='top',
                  weight='bold', color='k')
 
-    pargs = dict(lw=1, color=color, zorder=-num)
-
     axs1[0].plot(r, inten,              **pargs)
     axs1[1].plot(r, inten, label=label, **pargs)
     axs1[2].plot(r, inten,              **pargs)
 
 
-axc = fig.add_axes([0.93, 0.06, 0.01, 0.93])
+axc = fig.add_axes([0.93, 0.05, 0.01, 0.94])
 cbar = plt.colorbar(im, orientation="vertical", cax=axc, label='Intensity')
 cbar.ax.xaxis.set_ticks_position('top')
 cbar.ax.xaxis.set_label_position('top')
@@ -90,7 +94,7 @@
                        rotation='vertical', verticalalignment='center')
     ax.set_yticks(minor, minor=True)
 
-fig.subplots_adjust(left=0.05, bottom=0.06, right=0.92, top=0.99,
+fig.subplots_adjust(left=0.06, bottom=0.05, right=0.92, top=0.99,
                     wspace=0.04, hspace=0.04)
 
 for ax in axs[-1]:
diff --git a/doc/transform_methods/comparison/fig_gaussian/gaussian.py b/doc/transform_methods/comparison/fig_gaussian/gaussian.py
index 14f2e9c4..05b83f89 100644
--- a/doc/transform_methods/comparison/fig_gaussian/gaussian.py
+++ b/doc/transform_methods/comparison/fig_gaussian/gaussian.py
@@ -4,6 +4,10 @@
 
 transforms = [
   ("basex",         abel.basex.basex_transform),
+  (r"daun\ (degree=1)", lambda *args, **kwargs:
+                    abel.daun.daun_transform(*args, degree=1, **kwargs)),
+  (r"daun\ (degree=3)", lambda *args, **kwargs:
+                    abel.daun.daun_transform(*args, degree=3, **kwargs)),
   ("direct",        abel.direct.direct_transform),
   ("hansenlaw",     abel.hansenlaw.hansenlaw_transform),
   ("onion_bordas",  abel.onion_bordas.onion_bordas_transform),
@@ -42,7 +46,7 @@ def a(n, x):
     proj = 2*np.sqrt(1-r**2)
     dr = r[1] - r[0]
 
-fig, axs = plt.subplots(ntrans, 1, figsize=(5, 9), sharex=True, sharey=True)
+fig, axs = plt.subplots(ntrans, 1, figsize=(5, 10), sharex=True, sharey=True)
 
 for row, (label, transFunc) in enumerate(transforms):
     axs[row].plot(r, func, label='Analytical' if row == 0 else None, lw=1)
@@ -60,23 +64,25 @@ def a(n, x):
 
     axs[row].axhline(0, color='k', alpha=0.3, lw=1)
 
-    axs[row].legend(loc='upper right', frameon=False)
+    axs[row].legend(loc='upper right', frameon=False,
+                    handletextpad=0 if row > 0 else None)
 
     axs[row].grid(ls='solid', alpha=0.05, color='k')
 
 
 axs[-1].set_xlabel("$r$ (pixels)")
-axs[3].set_ylabel('$z$')
+axs[4].set_ylabel('$z$')
 
 for ax, letter in zip(axs, 'abcdefghi'):
     ax.grid(ls='solid', alpha=0.05, color='k')
     if case == 'gaussian':
-        ax.set_ylim(-0.25, 1.25)
+        ax.set_ylim(-0.2, 1.3)
+        ax.set_yticks([0, 0.5, 1])
         ax.set_xlim(0, n*0.74)
     else:
         ax.set_xlim(0, 1)
 
-    ax.annotate(letter + ')', xy=(0.02, 0.88), xytext=(0, 0),
+    ax.annotate(letter + ')', xy=(0.02, 0.86), xytext=(0, 0),
                 textcoords='offset points', xycoords='axes fraction',
                 weight='bold')
 
diff --git a/doc/transform_methods/daun-basis.py b/doc/transform_methods/daun-basis.py
new file mode 100644
index 00000000..52ccda6b
--- /dev/null
+++ b/doc/transform_methods/daun-basis.py
@@ -0,0 +1,90 @@
+# -*- coding: UTF-8 -*-
+from __future__ import print_function, division
+
+import numpy as np
+from scipy.linalg import solve_banded
+
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
+
+from abel.tools.polynomial import PiecewisePolynomial as PP
+
+n = 10
+r0 = 2.3
+
+# radial array with 20 point per pixel and half-integer point split into
+# ±ε pairs for sharp transitions for degree = 0
+r = np.linspace(0, n - 1, (n - 1) * 20 + 1) + 1e-5
+r = np.sort(np.concatenate((r, 0.5 + np.arange(0, n - 1) - 1e-5)))
+
+
+def f(x):
+    return np.exp(-(x - r0)**2 / 2) + np.exp(-(-x - r0)**2 / 2) + \
+           np.exp(-(x - 7)**2 * 3)
+
+
+def p(degree, j):
+    if degree == 0:
+        return PP(r, [(j - 1/2, j + 1/2, [1], j)])
+    elif degree == 1:
+        return PP(r, [(j - 1, j, [1,  1], j),
+                      (j, j + 1, [1, -1], j)])
+    elif degree == 2:
+        return PP(r, [(j - 1,   j - 1/2, [0, 0,  2], j - 1),
+                      (j - 1/2, j + 1/2, [1, 0, -2], j),
+                      (j + 1/2, j + 1,   [0, 0,  2], j + 1)])
+    else:  # degree == 3
+        return PP(r, [(j - 1, j, [1, 0, -3, -2], j),
+                      (j, j + 1, [1, 0, -3,  2], j)])
+
+
+def q(j):
+    return PP(r, [(j - 1, j, [0, 1,  2, 1], j),
+                  (j, j + 1, [0, 1, -2, 1], j)])
+
+
+plt.figure(figsize=(5, 7))
+
+colors = cm.rainbow(np.linspace(1, 0, n - 1))
+
+for degree in range(4):
+    plt.subplot(4, 1, 1 + degree)
+    plt.title('degree = ' + str(degree),
+              fontdict={'fontsize': plt.rcParams['axes.labelsize'],
+                        'fontweight': 'bold'},
+              loc='left')
+
+    F = f(np.arange(n))
+    if degree == 3:
+        D = np.concatenate(([0], 3 * (F[2:] - F[:-2]), [0]))
+        G = solve_banded((1, 1), ([0, 0] + [1] * (n - 2),
+                                  [4] * n,
+                                  [1] * (n - 2) + [0, 0]), D)
+
+    total = np.zeros_like(r)
+    for j, c in enumerate(colors):
+        p_j = F[j] * p(degree, j).func
+        plt.fill_between(r, 0, p_j, color=c, alpha=0.5)
+        total += p_j
+        if degree == 3:
+            q_j = G[j] * q(j).func
+            plt.plot(r, q_j, c=c)
+            total += q_j
+
+    plt.plot(r, f(r), c='k', label='test')
+    plt.plot(r, total, '--', c='k', lw=2, label='approx.')
+
+    plt.xlim((0, n - 1))
+    if degree == 3:
+        plt.xlabel('radius')
+
+    plt.ylim((-0.15, 1.05))
+    plt.ylabel('value')
+
+    if degree == 0:
+        plt.legend(loc='upper center', bbox_to_anchor=(0.55, 1))
+
+plt.tight_layout()
+
+#plt.savefig('daun-basis.svg')
+#plt.show()
diff --git a/doc/transform_methods/daun.rst b/doc/transform_methods/daun.rst
new file mode 100644
index 00000000..ea7a7f52
--- /dev/null
+++ b/doc/transform_methods/daun.rst
@@ -0,0 +1,340 @@
+.. _Daun:
+
+Daun
+====
+
+
+Introduction
+------------
+
+This suite of methods is based on the deconvolution procedure with Tikhonov
+regularization described by Daun at al. [1]_ and extends it with additional,
+smoother, approximations and regularization types.
+
+
+How it works
+------------
+
+The original method formulates the numerical Abel transform in a form
+equivalent to the “onion peeling” method, where the original distribution is
+approximated with a `step function
+<https://en.wikipedia.org/wiki/Step_function>`__ (piecewise constant function
+with 1-pixel-wide radial intervals). The forward transform thus can be
+described by a system on linear equations in a matrix form
+
+.. math::
+    \mathbf A_\text{OP} \mathbf x = \mathbf b,
+    :label: forward
+
+where the vector :math:`\mathbf x` consists of the original distribution values
+sampled at a uniform radial grid :math:`r_i = i \Delta r`, the vector
+:math:`\mathbf b` consists of the projection values sampled at the
+corresponding uniform grid :math:`y_i = i \Delta r`, and the matrix
+:math:`\mathbf A_\text{OP}` corresponds to the “onion peeling” forward Abel
+transform. Its elements are
+
+.. math::
+    A_{\text{OP}, ij} = \begin{cases}
+        0, & j < i, \\
+        2 \Delta r \big[(j + 1/2)^2 - i^2\big]^{1/2}, & j = i, \\
+        2 \Delta r \Big(\big[(j + 1/2)^2 - i^2\big]^{1/2} -
+                        \big[(j - 1/2)^2 - i^2\big]^{1/2}\Big), & j > i
+    \end{cases}
+
+and represent contributions of each distribution interval to each projection
+interval.
+
+However, instead of performing the inverse transform by using the inverse
+matrix directly:
+
+.. math::
+    \mathbf x = \mathbf A_\text{OP}^{-1} \mathbf b,
+
+as it is done in the onion peeling method, the equation :eq:`forward` is solved
+by applying `Tikhonov regularization
+<https://en.wikipedia.org/wiki/Tikhonov_regularization>`__:
+
+.. math::
+    \tilde{\mathbf x} = \operatorname{arg\,min}_{\mathbf x} \left(
+        \|\mathbf A_\text{OP} \mathbf x - \mathbf b\|^2 +
+        \alpha \|\mathbf L \mathbf x\|^2
+    \right),
+    :label: minimization
+
+where :math:`\alpha` is the regularization parameter, and the finite-difference
+matrix
+
+.. math::
+    \mathbf L = \begin{bmatrix}
+        -1     &  1     & 0      & \cdots & 0      \\
+        0      & -1     & 1      & \ddots & \vdots \\
+        \vdots & \ddots & \ddots & \ddots & \vdots \\
+        0      & \cdots & 0      & -1     & 1
+    \end{bmatrix}
+
+(approximation of the derivative operator) is used as the Tikhonov matrix. The
+idea is that the admixture of the derivative norm to the minimization problem
+leads to a smoother solution. The regularization parameter :math:`\alpha`
+controls how much attention is paid to the derivative: when :math:`\alpha = 0`,
+the exact solution to :eq:`forward` is obtained, even if very noisy; when
+:math:`\alpha \to \infty`, the solution becomes very smooth, even if
+reproducing the data poorly. A reasonably chosen value of :math:`\alpha` can
+result in a significant suppression of high-frequency noise without noticeably
+affecting the signal.
+
+The minimization problem :eq:`minimization` leads again to a linear matrix
+equation, and the regularized inverse transform is obtained by using the
+regularized matrix inverse:
+
+.. math::
+    \tilde{\mathbf x} = \mathbf A_\text{Tik}^{-1} \mathbf b_\text{Tik},
+    \quad
+    A_\text{Tik} = (\mathbf A^T \mathbf A +
+                    \alpha \mathbf A \mathbf L^T \mathbf L),
+    \quad
+    \mathbf b_\text{Tik} = A^T \mathbf b.
+    :label: tikhonov
+
+(Note: here :math:`\mathbf x` and :math:`\mathbf b` are column vectors, but in
+PyAbel they are row vectors corresponding to image rows, so all the equations
+are transposed; moreover, instead of processing row vectors separately, they
+are transformed as the image matrix at once.)
+
+
+PyAbel additions
+----------------
+
+Basis sets
+^^^^^^^^^^
+
+The step-function approximation used in the original method implies a basis set
+consisting of rectangular functions
+
+.. math::
+    f_i(r) = \begin{cases}
+        1, & r \in [i - 1/2, i + 1/2], \\
+        0  & \text{otherwise}.
+    \end{cases}
+
+This approximation can be considered rather coarse, so in addition to these
+zero-degree piecewise polynomials we also implement basis sets consisting of
+piecewise polynomials up to 3rd degree. An example of a test function composed
+of broad and narrow Gaussian peaks and its approximations of various degrees is
+shown below:
+
+.. plot:: transform_methods/daun-basis.py
+    :align: center
+
+Here the solid black line is the test function, and the dashed black line is
+its approximation of degree :math:`n`, equal to the sum of the colored basis
+functions.
+
+degree = 0:
+    Rectangular functions produce a stepwise approximation. This is the only
+    approach mentioned in the original article and corresponds to the usual
+    “onion peeling” transform.
+degree = 1:
+    `Triangular functions
+    <https://en.wikipedia.org/wiki/Triangular_function>`__ produce a continuous
+    piecewise linear approximation. Besides being continuous (although not
+    smooth), this also corresponds to how numerical data is usually plotted
+    (with points connected by straight lines), so such plots would faithfully
+    convey the underlying method assumptions.
+degree = 2:
+    Piecewise quadratic functions
+
+    .. math::
+        f_i(r) = \begin{cases}
+            2[r - (i - 1)]^2, & r \in [i - 1,   i - 1/2], \\
+            1 - 2[r - i]^2,   & r \in [i - 1/2, i + 1/2], \\
+            2[r - (i + 1)]^2, & r \in [i + 1/2, i + 1], \\
+            0                 & \text{otherwise}.
+        \end{cases}
+
+    produce a smooth piecewise quadratic approximation. While resembling
+    :ref:`BASEX basis functions <BASEXcomp>` in shape, these are localized
+    within ±1 pixel, sum to unity (although produce oscillations on slopes),
+    and their projections are much faster to compute.
+degree = 3:
+    Combinations of `cubic Hermite basis functions
+    <https://en.wikipedia.org/wiki/Cubic_Hermite_spline#Interpolation_on_a_single_interval>`__
+    produce a cubic-spline approximation (with endpoint derivatives clamped to
+    zero for 2D smoothness). Offers the most accurate representation for
+    sufficiently smooth distributions, but produces ringing artifacts around
+    sharp features, which can result in negative interpolated intensities even
+    for non-negative data points.
+
+(The projections of all these basis functions are calculated as described in
+:ref:`Polynomials`.)
+
+In practice, however, the choice of the basis set has negligible effect on the
+transform results, as can be seen from an example :ref:`below
+<example_degrees>`. Nevertheless, cubic splines might be useful for
+transforming smooth functions, in which case they yield very accurate results.
+
+
+Regularization methods
+^^^^^^^^^^^^^^^^^^^^^^
+
+:math:`L_2` norm
+""""""""""""""""
+
+In addition to the original derivative (difference) Tikhonov regularization,
+PyAbel also implements the usual :math:`L_2` regularization, as in
+:ref:`BASEX`, with the identity matrix :math:`\mathbf I` used instead of
+:math:`\mathbf L` in :eq:`tikhonov`. The results are practically identical to
+the BASEX method, especially with **degree** = 2, except that the basis set is
+computed much faster.
+
+Non-negativity
+""""""""""""""
+
+A more substantial addition is the implementation of the non-negativity
+regularization. Namely, instead of solving the unconstrained quadratic problem
+:eq:`minimization`, non-negativity constraints are imposed on the original
+problem:
+
+.. math::
+    \tilde{\mathbf x} = \operatorname{arg\,min}_{\mathbf x \geqslant 0}
+        \|\mathbf A \mathbf x - \mathbf b\|^2.
+
+This `non-negative least-squares
+<https://en.wikipedia.org/wiki/Non-negative_least_squares>`__ solution yields
+the distribution without negative intensities that reproduces the input data as
+well as possible. In situations where the distribution must be non-negative,
+this is the best physically meaningful solution.
+
+The noise-filtering properties of this method come from the fact that noise in
+the inverse Abel transform is strongly oscillating, so if negative-going spikes
+are forbidden in the solution, the positive-going spikes must also be reduced
+in order to preserve the overall intensity. Thus the method is most beneficial
+for very noisy images, for which linear methods produce a large amount of noise
+reaching negative values. For clean images of non-negative distributions, the
+constrained solution exactly matches the solution of the original problem
+:eq:`forward`. And unlike Tikhonov regularization, it does not blur legitimate
+sharp features in any case.
+
+Notice that constrained quadratic minimization remains a *non-linear* problem.
+This has several important implications. First, it is much more computationally
+challenging, so that transforming a megapixel image takes many seconds instead
+of several milliseconds (and depends on the image itself). Second, the average
+of transformed images is generally not equal to the transform of the averaged
+image. It is thus recommended to perform as much averaging (image
+symmetrization and summation of multiple images if applicable) as possible
+before applying the transform. In particular, using ``symmetry_axis=(0, 1)`` in
+:class:`abel.transform.Transform` would in fact require transforming only one
+quadrant, which is 4 times faster that transforming the whole image. Third, the
+method is only *asymptotically* `unbiased
+<https://en.wikipedia.org/wiki/Bias_of_an_estimator>`_, but for sufficiently
+noisy data can systematically shift or blur the intensities towards the
+symmetry axis. Thus while this regularization can be very helpful in revealing
+the structure in raw images, it is not recommended when further processing
+(like model fitting) is involved. In particular, for extraction of velocity and
+anisotropy distributions from velocity-map images, the :ref:`rBasex` method,
+possibly with its “positive” regularization, is more appropriate.
+
+
+When to use it
+--------------
+
+This method with default parameters (0th degree, 0 regularization parameter) is
+identical to the :doc:`“onion peeling” <onion_peeling>` method, but can also be
+used for the forward transform.
+
+The original (derivative/difference) Tikhonov regularization with non-zero
+regularization parameter helps to remove high-frequency oscillations from the
+transformed image. However, an excessively large regularization parameter can
+lead to oversmoothing and broadening of the useful signal and under/overshoots
+around sharp features. As recommended by Daun et al., by systematically
+adjusting the heuristic regularization parameter, the analyst can find a
+solution that represents an acceptable compromise between accuracy and
+regularity.
+
+The :math:`L_2` Tikhonov regularization approach is equivalent to that in the
+:ref:`BASEX` method and has the same use cases and [dis]advantages.
+
+The non-negativity regularization is recommended for visual inspection of very
+noisy images and images with sharp features without a broad background.
+However, due to its slowness, it cannot be used for real-time data processing.
+
+
+How to use it
+-------------
+
+The inverse Abel transform of a full image can be done with the
+:class:`abel.Transform <abel.transform.Transform>` class::
+
+    abel.Transform(myImage, method='daun').transform
+
+For the forward Abel transform, simply add :attr:`direction='forward'`::
+
+    abel.Transform(myImage, method='daun', direction='forward').transform
+
+Additional parameters can be passed through the :attr:`transform_options`
+parameter. For example, to use the original regularization method with the
+regularization parameter set to 100::
+
+    abel.Transform(myImage, method='daun',
+                   transform_options=dict{reg=100}).transform
+
+The :math:`L_2` regularization can be applied using ::
+
+    abel.Transform(myImage, method='daun',
+                   transform_options=dict{reg=('L2', 100)}).transform
+
+And the non-negative solution is obtained by ::
+
+    abel.Transform(myImage, method='daun',
+                   transform_options=dict{reg='nonneg'}).transform
+
+In this case, it is recommended to use symmetrization::
+
+    abel.Transform(myImage, method='daun',
+                   symmetry_axis=0,  # or symmetry_axis=(0, 1) if applicable
+                   transform_options=dict{reg='nonneg'}).transform
+
+unless independent inspection of all image parts is desired.
+
+The algorithm can be also accessed directly (to transform a right-side
+half-image or properly oriented quadrants) through the
+:func:`abel.daun.daun_transform()` function.
+
+.. note::
+    If you use any non-default options (degree, regularization), please cite
+    not only the article by Daun et al. and the PyAbel article, but also *this
+    PyAbel release* |zenodo|, because these capabilities are not present in the
+    original work by Daun et al. and were added to PyAbel after the RSI
+    publication.
+
+.. |zenodo| image:: https://zenodo.org/badge/30170345.svg
+    :target: https://zenodo.org/badge/latestdoi/30170345
+
+
+Examples
+--------
+
+Performance of various regularization methods for the Dribinski sample image
+with added Poissonian noise:
+
+.. plot:: ../examples/example_daun_reg.py
+
+:doc:`(source code) </example_daun_reg>`
+
+.. _example_degrees:
+
+The degree of basis-set polynomials has almost no effect on the results (shown
+here for :attr:`reg=0`):
+
+.. plot:: ../examples/example_daun_degree.py
+
+:doc:`(source code) </example_daun_degree>`
+
+
+Citation
+--------
+
+.. [1] \ K. J. Daun, K. A. Thomson, F. Liu, G. J. Smallwood,
+       “Deconvolution of axisymmetric flame properties using Tikhonov
+       regularization”,
+       `Appl. Opt. 45, 4638–4646 (2006)
+       <https://doi.org/10.1364/AO.45.004638>`_.
diff --git a/doc/transform_methods/rbasex-math.rst b/doc/transform_methods/rbasex-math.rst
index 18ae38e7..abdacc90 100644
--- a/doc/transform_methods/rbasex-math.rst
+++ b/doc/transform_methods/rbasex-math.rst
@@ -811,8 +811,7 @@ References
        CD\ :sub:`2`\ OH)”,
        Ph.D. dissertation, University of Southern California, 2012.
        (`ProQuest <https://search.proquest.com/docview/1289069738>`_,
-       `USC <http://digitallibrary.usc.edu/cdm/ref/collection/p15799coll3/id/
-       112619>`_).
+       `USC <https://digitallibrary.usc.edu/asset-management/2A3BF169XWB4>`_).
 
 .. [2] \ G. A. Garcia, L. Nahon, I. Powis,
        “Two-dimensional charged particle image inversion using a polar basis
diff --git a/doc/transform_methods/rbasex.rst b/doc/transform_methods/rbasex.rst
index 046a30ad..7b54549f 100644
--- a/doc/transform_methods/rbasex.rst
+++ b/doc/transform_methods/rbasex.rst
@@ -224,5 +224,4 @@ method from the PyAbel package”, using the current Zenodo DOI (see :ref:`READM
        CD\ :sub:`2`\ OH)”,
        Ph.D. dissertation, University of Southern California, 2012.
        (`ProQuest <https://search.proquest.com/docview/1289069738>`_,
-       `USC <https://digitallibrary.usc.edu/cdm/ref/collection/p15799coll3/id/
-       112619>`_).
+       `USC <https://digitallibrary.usc.edu/asset-management/2A3BF169XWB4>`_).
diff --git a/examples/example_GUI.py b/examples/example_GUI.py
index 02a936cf..7d8e39ed 100644
--- a/examples/example_GUI.py
+++ b/examples/example_GUI.py
@@ -24,8 +24,9 @@
 from six.moves import tkinter_tkfiledialog as filedialog
 
 
-Abel_methods = ['basex', 'direct', 'hansenlaw', 'linbasex', 'onion_bordas',
-                'onion_peeling', 'rbasex', 'three_point', 'two_point']
+Abel_methods = ['basex', 'daun(nonneg)', 'direct', 'hansenlaw', 'linbasex',
+                'onion_bordas', 'onion_peeling', 'rbasex', 'three_point',
+                'two_point']
 
 Abel_methods_forward = ['basex', 'direct', 'hansenlaw', 'rbasex']
 
@@ -392,6 +393,10 @@ def _transform(self):
                     tk.END,
                     "\nbasex: first time calculation of the basis"
                     " functions may take a while ...")
+            elif self.method == "daun(nonneg)":
+                self.text.insert(
+                    tk.END,
+                    "\ndaun(reg=nonneg): calculations will take a while...")
             elif self.method == "direct":
                 self.text.insert(
                    tk.END,
@@ -404,6 +409,10 @@ def _transform(self):
                     self.AIM = abel.Transform(
                         self.IM, method=self.method, direction=self.fi,
                         transform_options=dict(return_Beta=True))
+                elif self.method == 'daun(nonneg)':
+                    self.AIM = abel.Transform(
+                        self.IM, method='daun', direction=self.fi,
+                        transform_options=dict(reg='nonneg'))
                 else:
                     if self.method in Abel_methods_basis:
                         basis_dir = 'bases'
diff --git a/examples/example_daun_degree.py b/examples/example_daun_degree.py
new file mode 100644
index 00000000..7d7e3138
--- /dev/null
+++ b/examples/example_daun_degree.py
@@ -0,0 +1,78 @@
+# -*- coding: UTF-8 -*-
+from __future__ import division
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.interpolate import make_interp_spline, splev
+import abel
+
+# one quadrant of Dribinski sample image, its size and intensity distribution
+im = abel.tools.analytical.SampleImage().image
+q0 = abel.tools.symmetry.get_image_quadrants(im)[0]
+n = q0.shape[0]
+I0, _ = abel.tools.vmi.Ibeta(q0, origin='ll')
+
+# forward-transformed quadrant
+Q = abel.rbasex.rbasex_transform(q0, origin='ll', direction='forward')[0]
+# rescale intensities to 1000 max
+scale = 1000 / Q.max()
+q0 *= scale
+I0 *= scale
+Q *= scale
+# add Poissonian noise (comment out to compare the clean transform results)
+Q = np.random.RandomState(0).poisson(Q)
+
+# array for corresponding intensity distributions
+Is = []
+
+plt.figure(figsize=(7, 5))
+
+# transformed images
+for degree in range(4):
+    plt.subplot(3, 4, 1 + degree)
+    plt.axis('off')
+
+    q = abel.daun.daun_transform(Q, degree=degree)
+    plt.imshow(q, clim=(-3, 3), cmap='seismic')
+    plt.text(n / 2, 0, 'degree=' + str(degree), ha='center', va='top')
+
+    I, _ = abel.tools.vmi.Ibeta(q, origin='ll')
+    Is.append(I)
+
+# pixel subdivisions for smooth plots
+sub = 10
+rsub = np.linspace(0, n, sub * n + 1)
+
+# plots of intensity distributions
+for plot in [2, 3]:
+    plt.subplot(3, 1, plot)
+    plt.axhline(0, c='k', ls=':', lw=1)
+    plt.plot(I0, '--k', lw=1)
+
+    # degree = 0: plot with steps
+    plt.step(np.arange(n), Is[0], lw=1, label='0', where='mid')
+    # degree = 1: plot with lines
+    plt.plot(Is[1], lw=1, label='1')
+    # degree = 2: plot using parabolic segments
+    r1, r2 = np.arange(sub // 2) / sub, np.arange(sub // 2, 0, -1) / sub
+    b = np.concatenate((2 * r1**2, 1 - 2 * r2**2, 1 - 2 * r1**2, 2 * r2**2))
+    I = np.zeros_like(rsub)
+    for m in range(1, n):
+        i0 = sub * m
+        I[i0 - sub: i0 + sub] += Is[2][m] * b
+    plt.plot(rsub, I, lw=1, label='2')
+    # degree = 3: plot with cubic splines
+    spl = make_interp_spline(np.arange(n), Is[3], bc_type='clamped')
+    plt.plot(rsub, splev(rsub, spl), lw=1, label='3')
+
+    plt.xlim((0, n - 1))
+    plt.yticks([])
+    if plot == 2:  # full y range
+        plt.xticks([])
+        plt.legend(loc='upper center', bbox_to_anchor=(0.25, 1))
+    else:  # magnified
+        plt.ylim((-2000, 8000))
+
+plt.subplots_adjust(left=0.01, right=0.98,
+                    bottom=0.05, top=1,
+                    wspace=0, hspace=0.03)
+plt.show()
diff --git a/examples/example_daun_reg.py b/examples/example_daun_reg.py
new file mode 100644
index 00000000..3cb1856f
--- /dev/null
+++ b/examples/example_daun_reg.py
@@ -0,0 +1,66 @@
+# -*- coding: UTF-8 -*-
+from __future__ import division
+import numpy as np
+import matplotlib.pyplot as plt
+import abel
+
+# one quadrant of Dribinski sample image, its size and intensity distribution
+im = abel.tools.analytical.SampleImage().image
+q0 = abel.tools.symmetry.get_image_quadrants(im)[0]
+n = q0.shape[0]
+I0, _ = abel.tools.vmi.Ibeta(q0, origin='ll')
+
+# forward-transformed quadrant
+Q = abel.rbasex.rbasex_transform(q0, origin='ll', direction='forward')[0]
+# rescale intensities to 1000 max
+scale = 1000 / Q.max()
+q0 *= scale
+I0 *= scale
+Q *= scale
+# add Poissonian noise
+Q = np.random.RandomState(0).poisson(Q)
+
+# regularization parameters
+regs = [None,
+        ('diff', 100),  # RMS optimal is ~50
+        ('L2c', 50),    # RMS optimal is ~25
+        'nonneg']
+
+# array for corresponding intensity distributions
+Is = []
+
+plt.figure(figsize=(7, 5))
+
+# transformed images
+for i, reg in enumerate(regs):
+    plt.subplot(3, 4, 1 + i)
+    plt.axis('off')
+
+    q = abel.daun.daun_transform(Q, degree=1, reg=reg)
+    plt.imshow(q, clim=(-3, 3), cmap='seismic')
+    plt.text(n / 2, 0, 'reg=' + repr(reg), ha='center', va='top')
+
+    I, _ = abel.tools.vmi.Ibeta(q, origin='ll')
+    Is.append(I)
+
+# plots of intensity distributions
+for plot in [2, 3]:
+    plt.subplot(3, 1, plot)
+    plt.axhline(0, c='k', ls=':', lw=1)
+    plt.plot(I0, '--k', lw=1)
+
+    for i, reg in enumerate(regs):
+        plt.plot(Is[i], lw=1, label=repr(reg))
+
+    plt.xlim((0, n - 1))
+    plt.yticks([])
+    if plot == 2:  # full y range
+        plt.xticks([])
+        plt.legend(loc='upper center', bbox_to_anchor=(0.25, 1))
+    else:  # magnified
+        plt.ylim((-2000, 8000))
+
+plt.subplots_adjust(left=0.01, right=0.98,
+                    bottom=0.05, top=1,
+                    wspace=0, hspace=0.03)
+plt.show()


Test patch:
diff --git a/abel/tests/test_daun.py b/abel/tests/test_daun.py
new file mode 100644
index 00000000..b5842381
--- /dev/null
+++ b/abel/tests/test_daun.py
@@ -0,0 +1,253 @@
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+
+import os.path
+
+import numpy as np
+from numpy.testing import assert_allclose
+from scipy.linalg import inv, toeplitz
+
+import abel
+from abel.daun import _bs_daun, get_bs_cached, cache_cleanup, daun_transform
+from abel.tools.polynomial import PiecewisePolynomial as PP
+from abel.tools.analytical import GaussianAnalytical
+
+
+DATA_DIR = os.path.join(os.path.split(__file__)[0], 'data')
+
+
+def daun_bs(n, degree):
+    """Reference implementation of basis-set generation"""
+    r = np.arange(float(n))
+
+    A = np.empty((n, n))
+    for j in range(n):
+        if degree == 0:
+            p = PP(r, [(j - 1/2, j + 1/2, [1], j)])
+        elif degree == 1:
+            p = PP(r, [(j - 1, j, [1,  1], j),
+                       (j, j + 1, [1, -1], j)])
+        elif degree == 2:
+            p = PP(r, [(j - 1,   j - 1/2, [0, 0,  2], j - 1),
+                       (j - 1/2, j + 1/2, [1, 0, -2], j),
+                       (j + 1/2, j + 1,   [0, 0,  2], j + 1)])
+        else:  # degree == 3:
+            p = PP(r, [(j - 1, j, [1, 0, -3, -2], j),
+                       (j, j + 1, [1, 0, -3,  2], j)])
+        A[j] = p.abel
+
+    if degree == 3:
+        B = np.empty((n, n))
+        for j in range(n):
+            B[j] = PP(r, [(j - 1, j, [0, 1,  2, 1], j),
+                          (j, j + 1, [0, 1, -2, 1], j)]).abel
+        C = toeplitz([4, 1] + [0] * (n - 2))
+        C[1, 0] = C[-2, -1] = 0
+        D = toeplitz([0, 3] + [0] * (n - 2), [0, -3] + [0] * (n - 2))
+        D[1, 0] = D[-2, -1] = 0
+        A += D.dot(inv(C)).dot(B)
+
+    return A
+
+
+def test_daun_bs():
+    """Check basis-set generation for all degrees"""
+    n = 10
+    for degree in range(4):
+        A = _bs_daun(n, degree)
+        Aref = daun_bs(n, degree)
+        assert_allclose(A, Aref, err_msg='-> degree = ' + str(degree))
+
+
+def test_daun_bs_cache():
+    """Check basis-set cache handling"""
+
+    # changing size
+    n1, n2 = 20, 30
+    for degree in range(4):
+        # forward
+        cache_cleanup()
+        bs2_ref = get_bs_cached(n2, degree, direction='forward')
+        cache_cleanup()
+        bs1_ref = get_bs_cached(n1, degree, direction='forward')
+        bs2 = get_bs_cached(n2, degree, direction='forward')
+        bs1 = get_bs_cached(n1, degree, direction='forward')
+        assert_allclose(bs1, bs1_ref, atol=1e-15,
+                        err_msg='-> forward: n<, degree = {}'.format(degree))
+        assert_allclose(bs2, bs2_ref, atol=1e-15,
+                        err_msg='-> forward: n>, degree = {}'.format(degree))
+
+        # inverse
+        for reg_type in [None, 'diff', 'L2', 'nonneg']:
+            cache_cleanup()
+            bs2_ref = get_bs_cached(n2, degree, reg_type, 1)
+            cache_cleanup()
+            bs1_ref = get_bs_cached(n1, degree, reg_type, 1)
+            bs2 = get_bs_cached(n2, degree, reg_type, 1)
+            bs1 = get_bs_cached(n1, degree, reg_type, 1)
+            assert_allclose(bs1, bs1_ref, atol=1e-15,
+                            err_msg='-> inverse: n<, degree = {}, reg_type = {}'.
+                                    format(degree, reg_type))
+            assert_allclose(bs2, bs2_ref, atol=1e-15,
+                            err_msg='-> inverse: n>, degree = {}, reg_type = {}'.
+                                    format(degree, reg_type))
+
+    # changing degree
+    for direction in ['forward', 'inverse']:
+        cache_cleanup()
+        bs0_ref = get_bs_cached(n1, 0)
+        cache_cleanup()
+        bs1_ref = get_bs_cached(n1, 1)
+        bs0 = get_bs_cached(n1, 0)
+        bs1 = get_bs_cached(n1, 1)
+        assert_allclose(bs0, bs0_ref, atol=1e-15,
+                        err_msg='-> {}: degree>'.
+                                format(direction, degree, reg_type))
+        assert_allclose(bs1, bs1_ref, atol=1e-15,
+                        err_msg='-> {}: degree<'.
+                                format(direction, degree, reg_type))
+
+    # changing regularization strength
+    bs_ref = []
+    for reg in [0, 1, 2]:
+        cache_cleanup()
+        bs_ref.append(get_bs_cached(n1, strength=reg))
+    bs = []
+    for reg in [0, 1, 2]:
+        bs.append(get_bs_cached(n1, strength=reg))
+    assert_allclose(bs[0], bs_ref[0], atol=1e-15, err_msg='-> reg = 0 after 2')
+    assert_allclose(bs[1], bs_ref[1], atol=1e-15, err_msg='-> reg = 1 after 0')
+    assert_allclose(bs[2], bs_ref[2], atol=1e-15, err_msg='-> reg = 2 after 1')
+
+
+def get_basis_file_name(n, degree):
+    return os.path.join(DATA_DIR, 'daun_basis_{}_{}.npy'.format(n, degree))
+
+
+def test_daun_bs_disk_cache():
+    """Check basis-set disk caching"""
+    n1, n2 = 10, 20
+
+    f_n1_0 = get_basis_file_name(n1, 0)
+    f_n2_0 = get_basis_file_name(n2, 0)
+    f_n1_1 = get_basis_file_name(n1, 1)
+
+    for fname in [f_n1_0, f_n2_0, f_n1_1]:
+        if os.path.exists(fname):
+            os.remove(fname)
+
+    # saving
+    cache_cleanup()
+    bs_n1_0 = get_bs_cached(n1, 0, basis_dir=DATA_DIR)
+    assert os.path.exists(f_n1_0), \
+           'Basis set n = {}, degree = {} was not saved!'.format(n1, 0)
+
+    bs_n2_0 = get_bs_cached(n2, 0, basis_dir=DATA_DIR)
+    assert os.path.exists(f_n2_0), \
+           'Basis set n = {}, degree = {} was not saved!'.format(n2, 0)
+
+    bs_n1_1 = get_bs_cached(n1, 1, basis_dir=DATA_DIR)
+    assert os.path.exists(f_n1_1), \
+           'Basis set n = {}, degree = {} was not saved!'.format(n1, 1)
+
+    # loading
+    cache_cleanup()
+    bs = get_bs_cached(n1, 0, basis_dir=DATA_DIR)
+    assert_allclose(bs_n1_0, bs, err_msg='Loaded basis set ' +
+                    'n = {}, degree = {} differs from saved!'.format(n1, 0))
+
+    cache_cleanup()
+    bs = get_bs_cached(n2, 0, basis_dir=DATA_DIR)
+    assert_allclose(bs_n2_0, bs, err_msg='Loaded basis set ' +
+                    'n = {}, degree = {} differs from saved!'.format(n2, 0))
+
+    cache_cleanup()
+    bs = get_bs_cached(n1, 1, basis_dir=DATA_DIR)
+    assert_allclose(bs_n1_1, bs, err_msg='Loaded basis set ' +
+                    'n = {}, degree = {} differs from saved!'.format(n1, 1))
+
+    # crop-loading
+    os.remove(f_n1_0)
+    cache_cleanup()
+    bs = get_bs_cached(n1, 0, basis_dir=DATA_DIR)
+    assert_allclose(bs_n1_0, bs, err_msg='Loaded cropped basis set differs')
+
+    # clean-up
+    os.remove(f_n2_0)
+    os.remove(f_n1_1)
+
+
+def test_daun_shape():
+    n = 21
+    x = np.ones((n, n), dtype='float32')
+    recon = daun_transform(x, verbose=False)
+    assert recon.shape == (n, n)
+
+
+def test_daun_zeros():
+    n = 21
+    x = np.zeros((n, n), dtype='float32')
+    recon = abel.daun.daun_transform(x, verbose=False)
+    assert_allclose(recon, 0)
+
+
+def test_daun_gaussian():
+    """Check inverse Daun transform of a gaussian"""
+    n = 100
+    r_max = n - 1
+
+    ref = GaussianAnalytical(n, r_max, symmetric=False, sigma=30)
+    tr = np.tile(ref.abel[None, :], (n, 1))  # make a 2D array from 1D
+
+    for degree, tol in [(0, 2e-3), (1, 3e-4), (2, 9e-4), (3, 4e-5)]:
+        recon = daun_transform(tr, degree=degree, verbose=False)
+        recon = recon[n // 2 + n % 2]
+        assert_allclose(recon, ref.func, atol=tol,
+                        err_msg='-> degree = ' + str(degree))
+
+    for reg, tol in ([None, 2e-3],
+                     [1, 1e-3],
+                     [('diff', 1), 1e-3],
+                     [('L2', 1), 6e-4],  # for r > 0
+                     [('L2c', 1), 6e-4],
+                     ['nonneg', 2e-3]):
+        recon = daun_transform(tr, reg=reg, verbose=False)
+        recon = recon[n // 2 + n % 2]
+        skip = 1 if isinstance(reg, tuple) and reg[0] == 'L2' else 0
+        assert_allclose(recon[skip:], ref.func[skip:], atol=tol,
+                        err_msg='-> reg = ' + repr(reg))
+
+    # test dr
+    recon = daun_transform(tr, dr=0.5, verbose=False)[0]
+    assert_allclose(recon, 2 * ref.func, atol=3e-3, err_msg='-> dr = 0.5')
+
+
+def test_daun_forward_gaussian():
+    """Check forward Daun transform of a gaussian"""
+    n = 100
+    r_max = n - 1
+
+    ref = GaussianAnalytical(n, r_max, symmetric=False, sigma=30)
+    tr = np.tile(ref.func[None, :], (n, 1))  # make a 2D array from 1D
+
+    for degree, tol in [(0, 3e-2), (1, 7e-3), (2, 2e-2), (3, 7e-4)]:
+        proj = daun_transform(tr, degree=degree, direction='forward',
+                              verbose=False)
+        proj = proj[n // 2 + n % 2]
+        assert_allclose(proj, ref.abel, atol=tol,
+                        err_msg='-> degree = ' + str(degree))
+
+    # test dr
+    proj = daun_transform(tr, dr=0.5, direction='forward', verbose=False)[0]
+    assert_allclose(proj, ref.abel / 2, atol=2e-2, err_msg='-> dr = 0.5')
+
+
+if __name__ == '__main__':
+    test_daun_bs()
+    test_daun_bs_cache()
+    test_daun_bs_disk_cache()
+    test_daun_shape()
+    test_daun_zeros()
+    test_daun_gaussian()
+    test_daun_forward_gaussian()