From d8ab26a21d2204841eb66097bdd15be621ae47a4 Mon Sep 17 00:00:00 2001
From: Gregory Lee <grlee77@gmail.com>
Date: Fri, 29 Jan 2021 19:45:17 -0500
Subject: [PATCH 01/12] ENH: add float32 support to wiener,
 unsupervised_wiener, rolling_ball

---
 TODO.txt                                      |  3 ++
 skimage/restoration/_denoise.py               |  5 ++-
 skimage/restoration/deconvolution.py          | 42 +++++++++++++++----
 skimage/restoration/rolling_ball.py           |  5 ++-
 skimage/restoration/tests/test_denoise.py     | 10 +++--
 skimage/restoration/tests/test_inpaint.py     |  6 ++-
 skimage/restoration/tests/test_j_invariant.py |  7 +++-
 skimage/restoration/tests/test_restoration.py | 31 ++++++++++----
 .../restoration/tests/test_rolling_ball.py    | 10 +++++
 skimage/restoration/uft.py                    | 11 +++--
 10 files changed, 102 insertions(+), 28 deletions(-)

diff --git a/TODO.txt b/TODO.txt
index 17522c00a6a..08d30811e3b 100644
--- a/TODO.txt
+++ b/TODO.txt
@@ -102,6 +102,9 @@ Other (2022)
 ------------
 * Remove conditional import of ``scipy.fft`` in ``skimage._shared.fft`` once
   the minimum supported version of ``scipy`` reaches 1.4.
+* Remove unneeded `deconv_type` and astype call in `weiner` and
+  `unsupervised_wiener` in `skimage.restoration.deconvolution` once the minimum
+  supported version of ``scipy`` reaches 1.4.
 * Remove pillow version related warning for CVE when pillow > 7.1.0 in
   `skimage/io/_plugins/pil_plugin.py` and `skimage/io/collection.py`.
 
diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index 360d2d332d6..7502b59f599 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -3,7 +3,7 @@
 from math import ceil
 from .. import img_as_float
 from ._denoise_cy import _denoise_bilateral, _denoise_tv_bregman
-from .._shared.utils import warn
+from .._shared.utils import _float_type, warn
 import pywt
 import skimage.color as color
 from skimage.color.colorconv import ycbcr_from_rgb
@@ -645,6 +645,7 @@ def _scale_sigma_and_image_consistently(image, sigma, multichannel,
     """If the ``image`` is rescaled, also rescale ``sigma`` consistently.
 
     Images that are not floating point will be rescaled via ``img_as_float``.
+    Half-precision images will be promoted to single precision.
     """
     if multichannel:
         if isinstance(sigma, numbers.Number) or sigma is None:
@@ -666,6 +667,8 @@ def _scale_sigma_and_image_consistently(image, sigma, multichannel,
                          for s in sigma]
             elif sigma is not None:
                 sigma *= scale_factor
+    elif image.dtype == np.float16:
+        image = image.astype(np.float32)
     return image, sigma
 
 
diff --git a/skimage/restoration/deconvolution.py b/skimage/restoration/deconvolution.py
index 121609d2f77..21f158f7395 100644
--- a/skimage/restoration/deconvolution.py
+++ b/skimage/restoration/deconvolution.py
@@ -6,6 +6,7 @@
 from scipy.signal import convolve
 
 from . import uft
+from .._shared.utils import _float_type
 
 __keywords__ = "restoration, image, deconvolution"
 
@@ -116,6 +117,10 @@ def wiener(image, psf, balance, reg=None, is_real=True, clip=True):
         reg, _ = uft.laplacian(image.ndim, image.shape, is_real=is_real)
     if not np.iscomplexobj(reg):
         reg = uft.ir2tf(reg, image.shape, is_real=is_real)
+    float_type = _float_type(image)
+    image = image.astype(float_type, copy=False)
+    psf = psf.real.astype(float_type, copy=False)
+    reg = reg.real.astype(float_type, copy=False)
 
     if psf.shape != reg.shape:
         trans_func = uft.ir2tf(psf, image.shape, is_real=is_real)
@@ -130,6 +135,13 @@ def wiener(image, psf, balance, reg=None, is_real=True, clip=True):
     else:
         deconv = uft.uifft2(wiener_filter * uft.ufft2(image))
 
+    # TODO: can remove astype call below once minimum SciPy >= 1.4
+    if deconv.dtype.kind == 'c':
+        deconv_type = np.promote_types(float_type, np.complex64)
+    else:
+        deconv_type = float_type
+    deconv = deconv.astype(deconv_type, copy=False)
+
     if clip:
         deconv[deconv > 1] = 1
         deconv[deconv < -1] = -1
@@ -238,6 +250,10 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True,
         reg, _ = uft.laplacian(image.ndim, image.shape, is_real=is_real)
     if not np.iscomplexobj(reg):
         reg = uft.ir2tf(reg, image.shape, is_real=is_real)
+    float_type = _float_type(image)
+    image = image.astype(float_type, copy=False)
+    psf = psf.real.astype(float_type, copy=False)
+    reg = reg.real.astype(float_type, copy=False)
 
     if psf.shape != reg.shape:
         trans_fct = uft.ir2tf(psf, image.shape,  is_real=is_real)
@@ -245,9 +261,9 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True,
         trans_fct = psf
 
     # The mean of the object
-    x_postmean = np.zeros(trans_fct.shape)
+    x_postmean = np.zeros(trans_fct.shape, dtype=float_type)
     # The previous computed mean in the iterative loop
-    prev_x_postmean = np.zeros(trans_fct.shape)
+    prev_x_postmean = np.zeros(trans_fct.shape, dtype=float_type)
 
     # Difference between two successive mean
     delta = np.NAN
@@ -263,9 +279,9 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True,
     # The Fourier transform may change the image.size attribute, so we
     # store it.
     if is_real:
-        data_spectrum = uft.urfft2(image.astype(float))
+        data_spectrum = uft.urfft2(image)
     else:
-        data_spectrum = uft.ufft2(image.astype(float))
+        data_spectrum = uft.ufft2(image)
 
     # Gibbs sampling
     for iteration in range(params['max_iter']):
@@ -273,9 +289,11 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True,
 
         # weighting (correlation in direct space)
         precision = gn_chain[-1] * atf2 + gx_chain[-1] * areg2  # Eq. 29
+        _rand1 = np.random.standard_normal(data_spectrum.shape)
+        _rand2 = np.random.standard_normal(data_spectrum.shape)
         excursion = np.sqrt(0.5) / np.sqrt(precision) * (
-            np.random.standard_normal(data_spectrum.shape) +
-            1j * np.random.standard_normal(data_spectrum.shape))
+            _rand1.astype(float_type, copy=False) +
+            1j * _rand2.astype(float_type, copy=False))
 
         # mean Eq. 30 (RLS for fixed gn, gamma0 and gamma1 ...)
         wiener_filter = gn_chain[-1] * np.conj(trans_fct) / precision
@@ -319,6 +337,13 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True,
     else:
         x_postmean = uft.uifft2(x_postmean)
 
+    # TODO: remove astype call below once minimum SciPy >= 1.4
+    if x_postmean.dtype.kind == 'c':
+        deconv_type = np.promote_types(float_type, np.complex64)
+    else:
+        deconv_type = float_type
+    x_postmean = x_postmean.astype(deconv_type, copy=False)
+
     if clip:
         x_postmean[x_postmean > 1] = 1
         x_postmean[x_postmean < -1] = -1
@@ -364,7 +389,7 @@ def richardson_lucy(image, psf, iterations=50, clip=True, filter_epsilon=None):
     ----------
     .. [1] https://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution
     """
-    float_type = np.promote_types(image.dtype, np.float32)
+    float_type = _float_type(image)
     image = image.astype(float_type, copy=False)
     psf = psf.astype(float_type, copy=False)
     im_deconv = np.full(image.shape, 0.5, dtype=float_type)
@@ -373,7 +398,8 @@ def richardson_lucy(image, psf, iterations=50, clip=True, filter_epsilon=None):
     for _ in range(iterations):
         conv = convolve(im_deconv, psf, mode='same')
         if filter_epsilon:
-            relative_blur = np.where(conv < filter_epsilon, 0, image / conv)
+            relative_blur = np.where(conv < filter_epsilon, 0,
+                                     image / conv)
         else:
             relative_blur = image / conv
         im_deconv *= convolve(relative_blur, psf_mirror, mode='same')
diff --git a/skimage/restoration/rolling_ball.py b/skimage/restoration/rolling_ball.py
index 85a03ee8d70..5142ac0e08c 100644
--- a/skimage/restoration/rolling_ball.py
+++ b/skimage/restoration/rolling_ball.py
@@ -1,5 +1,6 @@
 import numpy as np
 
+from .._shared.utils import _float_type
 from ._rolling_ball_cy import apply_kernel, apply_kernel_nan
 
 
@@ -78,7 +79,8 @@ def rolling_ball(image, *, radius=100, kernel=None,
     """
 
     image = np.asarray(image)
-    img = image.astype(float)
+    float_type = _float_type(image)
+    img = image.astype(float_type, copy=False)
 
     if num_threads is None:
         num_threads = 0
@@ -86,6 +88,7 @@ def rolling_ball(image, *, radius=100, kernel=None,
     if kernel is None:
         kernel = ball_kernel(radius, image.ndim)
 
+    kernel = kernel.astype(float_type)
     kernel_shape = np.asarray(kernel.shape)
     kernel_center = (kernel_shape // 2)
     center_intensity = kernel[tuple(kernel_center)]
diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py
index 09b2f113061..936018c0530 100644
--- a/skimage/restoration/tests/test_denoise.py
+++ b/skimage/restoration/tests/test_denoise.py
@@ -11,6 +11,7 @@
 from skimage._shared import testing
 from skimage._shared.testing import (assert_equal, assert_almost_equal,
                                      assert_warns, assert_)
+from skimage._shared.utils import _float_type
 from skimage._shared._warnings import expected_warnings
 from distutils.version import LooseVersion as Version
 
@@ -35,22 +36,24 @@
 astro_odd = astro[:, :-1]
 
 
-def test_denoise_tv_chambolle_2d():
+@testing.parametrize('dtype', [np.float32, np.float64])
+def test_denoise_tv_chambolle_2d(dtype):
     # astronaut image
-    img = astro_gray.copy()
+    img = astro_gray.astype(dtype, copy=True)
     # add noise to astronaut
     img += 0.5 * img.std() * np.random.rand(*img.shape)
     # clip noise so that it does not exceed allowed range for float images.
     img = np.clip(img, 0, 1)
     # denoise
     denoised_astro = restoration.denoise_tv_chambolle(img, weight=0.1)
+    assert denoised_astro.dtype == _float_type(img)
     # which dtype?
     assert_(denoised_astro.dtype in [float, np.float32, np.float64])
     from scipy import ndimage as ndi
     grad = ndi.morphological_gradient(img, size=((3, 3)))
     grad_denoised = ndi.morphological_gradient(denoised_astro, size=((3, 3)))
     # test if the total variation has decreased
-    assert_(grad_denoised.dtype == float)
+    assert_(grad_denoised.dtype == _float_type(img))
     assert_(np.sqrt((grad_denoised**2).sum()) < np.sqrt((grad**2).sum()))
 
 
@@ -554,6 +557,7 @@ def test_wavelet_denoising_scaling(case, dtype, convert2ycbcr,
                                            multichannel=multichannel,
                                            convert2ycbcr=convert2ycbcr,
                                            rescale_sigma=True)
+    assert denoised.dtype == _float_type(noisy)
 
     data_range = x.max() - x.min()
     psnr_noisy = peak_signal_noise_ratio(x, noisy, data_range=data_range)
diff --git a/skimage/restoration/tests/test_inpaint.py b/skimage/restoration/tests/test_inpaint.py
index 03033689e57..0019343152b 100644
--- a/skimage/restoration/tests/test_inpaint.py
+++ b/skimage/restoration/tests/test_inpaint.py
@@ -6,14 +6,16 @@
 from skimage._shared.testing import assert_allclose
 
 
-def test_inpaint_biharmonic_2d():
-    img = np.tile(np.square(np.linspace(0, 1, 5)), (5, 1))
+@testing.parametrize('dtype', [np.float32, np.float64])
+def test_inpaint_biharmonic_2d(dtype):
+    img = np.tile(np.square(np.linspace(0, 1, 5, dtype=dtype)), (5, 1))
     mask = np.zeros_like(img)
     mask[2, 2:] = 1
     mask[1, 3:] = 1
     mask[0, 4:] = 1
     img[np.where(mask)] = 0
     out = inpaint.inpaint_biharmonic(img, mask)
+    assert out.dtype == img.dtype
     ref = np.array(
         [[0., 0.0625, 0.25000000, 0.5625000, 0.73925058],
          [0., 0.0625, 0.25000000, 0.5478048, 0.76557821],
diff --git a/skimage/restoration/tests/test_j_invariant.py b/skimage/restoration/tests/test_j_invariant.py
index 5630eb82c94..6d8a313c7f9 100644
--- a/skimage/restoration/tests/test_j_invariant.py
+++ b/skimage/restoration/tests/test_j_invariant.py
@@ -1,6 +1,7 @@
 import functools
 import numpy as np
 
+from skimage._shared import testing
 from skimage._shared.testing import assert_
 from skimage.data import binary_blobs
 from skimage.data import camera, chelsea
@@ -27,10 +28,12 @@ def test_invariant_denoise():
     assert_(denoised_mse < original_mse)
 
 
-def test_invariant_denoise_color():
+@testing.parametrize('dtype', [np.float32, np.float64])
+def test_invariant_denoise_color(dtype):
     denoised_img_color = _invariant_denoise(
-            noisy_img_color, _denoise_wavelet,
+            noisy_img_color.astype(dtype), _denoise_wavelet,
             denoiser_kwargs=dict(multichannel=True))
+    assert denoised_img_color.dtype == dtype
 
     denoised_mse = mse(denoised_img_color, test_img_color)
     original_mse = mse(noisy_img_color, test_img_color)
diff --git a/skimage/restoration/tests/test_restoration.py b/skimage/restoration/tests/test_restoration.py
index cc2781bbe53..6313798a2ef 100644
--- a/skimage/restoration/tests/test_restoration.py
+++ b/skimage/restoration/tests/test_restoration.py
@@ -2,6 +2,7 @@
 import pytest
 from scipy.signal import convolve2d
 from scipy import ndimage as ndi
+from skimage._shared import testing
 from skimage._shared.testing import fetch
 from skimage.color import rgb2gray
 from skimage.data import astronaut, camera
@@ -11,46 +12,60 @@
 test_img = util.img_as_float(camera())
 
 
-def test_wiener():
-    psf = np.ones((5, 5)) / 25
+@testing.parametrize('dtype', [np.float32, np.float64])
+def test_wiener(dtype):
+    psf = np.ones((5, 5), dtype=dtype) / 25
     data = convolve2d(test_img, psf, 'same')
     np.random.seed(0)
     data += 0.1 * data.std() * np.random.standard_normal(data.shape)
+    data = data.astype(dtype, copy=False)
     deconvolved = restoration.wiener(data, psf, 0.05)
+    assert deconvolved.dtype == dtype
 
     path = fetch('restoration/tests/camera_wiener.npy')
-    np.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-3)
+    atol = 1e-5 if dtype == np.float32 else 0
+    np.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-3,
+                               atol=atol)
 
     _, laplacian = uft.laplacian(2, data.shape)
     otf = uft.ir2tf(psf, data.shape, is_real=False)
+    assert otf.real.dtype == dtype
     deconvolved = restoration.wiener(data, otf, 0.05,
                                      reg=laplacian,
                                      is_real=False)
+    assert deconvolved.real.dtype == dtype
     np.testing.assert_allclose(np.real(deconvolved),
                                np.load(path),
-                               rtol=1e-3)
+                               rtol=1e-3, atol=atol)
 
 
-def test_unsupervised_wiener():
-    psf = np.ones((5, 5)) / 25
+@testing.parametrize('dtype', [np.float32, np.float64])
+def test_unsupervised_wiener(dtype):
+    psf = np.ones((5, 5), dtype=dtype) / 25
     data = convolve2d(test_img, psf, 'same')
     np.random.seed(0)
     data += 0.1 * data.std() * np.random.standard_normal(data.shape)
+    data = data.astype(dtype, copy=False)
     deconvolved, _ = restoration.unsupervised_wiener(data, psf)
+    assert deconvolved.dtype == dtype
 
     path = fetch('restoration/tests/camera_unsup.npy')
-    np.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-3)
+    atol = 1e-5 if dtype == np.float32 else 0
+    np.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-3,
+                               atol=atol)
 
     _, laplacian = uft.laplacian(2, data.shape)
     otf = uft.ir2tf(psf, data.shape, is_real=False)
+    assert otf.real.dtype == dtype
     np.random.seed(0)
     deconvolved = restoration.unsupervised_wiener(
         data, otf, reg=laplacian, is_real=False,
         user_params={"callback": lambda x: None})[0]
+    assert deconvolved.real.dtype == dtype
     path = fetch('restoration/tests/camera_unsup2.npy')
     np.testing.assert_allclose(np.real(deconvolved),
                                np.load(path),
-                               rtol=1e-3)
+                               rtol=1e-3, atol=atol)
 
 
 def test_image_shape():
diff --git a/skimage/restoration/tests/test_rolling_ball.py b/skimage/restoration/tests/test_rolling_ball.py
index b23eca923fc..a8dd87a5fcf 100644
--- a/skimage/restoration/tests/test_rolling_ball.py
+++ b/skimage/restoration/tests/test_rolling_ball.py
@@ -7,6 +7,7 @@
 import pytest
 
 from skimage import data, io
+from skimage._shared import testing
 from skimage.restoration import rolling_ball
 from skimage.restoration.rolling_ball import ellipsoid_kernel
 
@@ -18,6 +19,15 @@ def test_ellipsoid_const():
     assert np.allclose(img - background, np.zeros_like(img))
 
 
+@testing.parametrize('dtype', [np.float32, np.float64])
+def test_rolling_ball_dtype(dtype):
+    img = 155 * np.ones((100, 100), dtype=dtype)
+    kernel = ellipsoid_kernel((25, 53), 50).astype(dtype, copy=False)
+    background = rolling_ball(img, kernel=kernel)
+    assert background.dtype == img.dtype
+    assert np.allclose(img - background, np.zeros_like(img))
+
+
 def test_nan_const():
     img = 123 * np.ones((100, 100), dtype=float)
     img[20, 20] = np.nan
diff --git a/skimage/restoration/uft.py b/skimage/restoration/uft.py
index 69029638f55..631df9ad1f7 100644
--- a/skimage/restoration/uft.py
+++ b/skimage/restoration/uft.py
@@ -390,7 +390,8 @@ def ir2tf(imp_resp, shape, dim=None, is_real=True):
     if not dim:
         dim = imp_resp.ndim
     # Zero padding and fill
-    irpadded = np.zeros(shape)
+    irpadded_dtype = imp_resp.dtype if imp_resp.dtype.kind == 'f' else float
+    irpadded = np.zeros(shape, dtype=irpadded_dtype)
     irpadded[tuple([slice(0, s) for s in imp_resp.shape])] = imp_resp
     # Roll for zero convention of the fft to avoid the phase
     # problem. Work with odd and even size.
@@ -400,9 +401,13 @@ def ir2tf(imp_resp, shape, dim=None, is_real=True):
                                shift=-int(np.floor(axis_size / 2)),
                                axis=axis)
     if is_real:
-        return fft.rfftn(irpadded, axes=range(-dim, 0))
+        out = fft.rfftn(irpadded, axes=range(-dim, 0))
     else:
-        return fft.fftn(irpadded, axes=range(-dim, 0))
+        out = fft.fftn(irpadded, axes=range(-dim, 0))
+
+    # TODO: remove .astype call onces SciPy >= 1.4 is required
+    cplx_dtype = np.promote_types(irpadded_dtype, np.complex64)
+    return out.astype(cplx_dtype, copy=False)
 
 
 def laplacian(ndim, shape, is_real=True):

From ff0ef76d5326b9e98c1ea4f655a3eadeb59c8dd7 Mon Sep 17 00:00:00 2001
From: Gregory Lee <grlee77@gmail.com>
Date: Thu, 4 Feb 2021 17:34:49 -0500
Subject: [PATCH 02/12] pep8 fixes

---
 skimage/restoration/deconvolution.py          | 4 ++--
 skimage/restoration/tests/test_j_invariant.py | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/skimage/restoration/deconvolution.py b/skimage/restoration/deconvolution.py
index 21f158f7395..2ae7067fd87 100644
--- a/skimage/restoration/deconvolution.py
+++ b/skimage/restoration/deconvolution.py
@@ -292,8 +292,8 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True,
         _rand1 = np.random.standard_normal(data_spectrum.shape)
         _rand2 = np.random.standard_normal(data_spectrum.shape)
         excursion = np.sqrt(0.5) / np.sqrt(precision) * (
-            _rand1.astype(float_type, copy=False) +
-            1j * _rand2.astype(float_type, copy=False))
+            _rand1.astype(float_type, copy=False)
+            + 1j * _rand2.astype(float_type, copy=False))
 
         # mean Eq. 30 (RLS for fixed gn, gamma0 and gamma1 ...)
         wiener_filter = gn_chain[-1] * np.conj(trans_fct) / precision
diff --git a/skimage/restoration/tests/test_j_invariant.py b/skimage/restoration/tests/test_j_invariant.py
index 6d8a313c7f9..7c7c617f9a7 100644
--- a/skimage/restoration/tests/test_j_invariant.py
+++ b/skimage/restoration/tests/test_j_invariant.py
@@ -31,8 +31,8 @@ def test_invariant_denoise():
 @testing.parametrize('dtype', [np.float32, np.float64])
 def test_invariant_denoise_color(dtype):
     denoised_img_color = _invariant_denoise(
-            noisy_img_color.astype(dtype), _denoise_wavelet,
-            denoiser_kwargs=dict(multichannel=True))
+        noisy_img_color.astype(dtype), _denoise_wavelet,
+        denoiser_kwargs=dict(multichannel=True))
     assert denoised_img_color.dtype == dtype
 
     denoised_mse = mse(denoised_img_color, test_img_color)

From b4d905ab9679be27d4d3a1765e7c0b36cb5f2cfc Mon Sep 17 00:00:00 2001
From: Gregory Lee <grlee77@gmail.com>
Date: Tue, 13 Apr 2021 17:18:25 -0400
Subject: [PATCH 03/12] incorporate reviewer feedback

_float_type -> _supported_float_type

test with additional dtypes (float16, float128)

suppress RuntimeWarnings in richardson_lucy
---
 skimage/restoration/_denoise.py               |  6 ++-
 skimage/restoration/deconvolution.py          | 15 +++---
 skimage/restoration/inpaint.py                |  6 +++
 skimage/restoration/rolling_ball.py           |  4 +-
 skimage/restoration/tests/test_denoise.py     | 21 +++++----
 skimage/restoration/tests/test_inpaint.py     |  6 ++-
 skimage/restoration/tests/test_restoration.py | 47 ++++++++++++-------
 .../restoration/tests/test_rolling_ball.py    | 15 ++----
 skimage/restoration/uft.py                    | 10 ++--
 9 files changed, 79 insertions(+), 51 deletions(-)

diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index 7502b59f599..9e50682916d 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -3,7 +3,7 @@
 from math import ceil
 from .. import img_as_float
 from ._denoise_cy import _denoise_bilateral, _denoise_tv_bregman
-from .._shared.utils import _float_type, warn
+from .._shared.utils import _supported_float_type, warn
 import pywt
 import skimage.color as color
 from skimage.color.colorconv import ycbcr_from_rgb
@@ -469,6 +469,10 @@ def denoise_tv_chambolle(image, weight=0.1, eps=2.e-4, n_iter_max=200,
     if not im_type.kind == 'f':
         image = img_as_float(image)
 
+    # enforce float16->float32 and float128->float64
+    float_dtype = _supported_float_type(image.dtype)
+    image = image.astype(float_dtype, copy=False)
+
     if multichannel:
         out = np.zeros_like(image)
         for c in range(image.shape[-1]):
diff --git a/skimage/restoration/deconvolution.py b/skimage/restoration/deconvolution.py
index 2ae7067fd87..fbf615ef782 100644
--- a/skimage/restoration/deconvolution.py
+++ b/skimage/restoration/deconvolution.py
@@ -6,7 +6,7 @@
 from scipy.signal import convolve
 
 from . import uft
-from .._shared.utils import _float_type
+from .._shared.utils import _supported_float_type
 
 __keywords__ = "restoration, image, deconvolution"
 
@@ -117,7 +117,7 @@ def wiener(image, psf, balance, reg=None, is_real=True, clip=True):
         reg, _ = uft.laplacian(image.ndim, image.shape, is_real=is_real)
     if not np.iscomplexobj(reg):
         reg = uft.ir2tf(reg, image.shape, is_real=is_real)
-    float_type = _float_type(image)
+    float_type = _supported_float_type(image.dtype)
     image = image.astype(float_type, copy=False)
     psf = psf.real.astype(float_type, copy=False)
     reg = reg.real.astype(float_type, copy=False)
@@ -250,13 +250,13 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True,
         reg, _ = uft.laplacian(image.ndim, image.shape, is_real=is_real)
     if not np.iscomplexobj(reg):
         reg = uft.ir2tf(reg, image.shape, is_real=is_real)
-    float_type = _float_type(image)
+    float_type = _supported_float_type(image.dtype)
     image = image.astype(float_type, copy=False)
     psf = psf.real.astype(float_type, copy=False)
     reg = reg.real.astype(float_type, copy=False)
 
     if psf.shape != reg.shape:
-        trans_fct = uft.ir2tf(psf, image.shape,  is_real=is_real)
+        trans_fct = uft.ir2tf(psf, image.shape, is_real=is_real)
     else:
         trans_fct = psf
 
@@ -389,7 +389,7 @@ def richardson_lucy(image, psf, iterations=50, clip=True, filter_epsilon=None):
     ----------
     .. [1] https://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution
     """
-    float_type = _float_type(image)
+    float_type = _supported_float_type(image.dtype)
     image = image.astype(float_type, copy=False)
     psf = psf.astype(float_type, copy=False)
     im_deconv = np.full(image.shape, 0.5, dtype=float_type)
@@ -398,8 +398,9 @@ def richardson_lucy(image, psf, iterations=50, clip=True, filter_epsilon=None):
     for _ in range(iterations):
         conv = convolve(im_deconv, psf, mode='same')
         if filter_epsilon:
-            relative_blur = np.where(conv < filter_epsilon, 0,
-                                     image / conv)
+            with np.errstate(invalid='ignore'):
+                relative_blur = np.where(conv < filter_epsilon, 0,
+                                         image / conv)
         else:
             relative_blur = image / conv
         im_deconv *= convolve(relative_blur, psf_mirror, mode='same')
diff --git a/skimage/restoration/inpaint.py b/skimage/restoration/inpaint.py
index ddcce13f63b..e4e02abb883 100644
--- a/skimage/restoration/inpaint.py
+++ b/skimage/restoration/inpaint.py
@@ -5,6 +5,7 @@
 import scipy.ndimage as ndi
 from scipy.ndimage import laplace
 import skimage
+from .._shared.utils import _supported_float_type
 from ..measure import label
 
 
@@ -124,6 +125,11 @@ def inpaint_biharmonic(image, mask, multichannel=False):
         raise TypeError('Masked arrays are not supported')
 
     image = skimage.img_as_float(image)
+
+    # float16->float32 and float128->float64
+    float_dtype = _supported_float_type(image.dtype)
+    image = image.astype(float_dtype, copy=False)
+
     mask = mask.astype(bool)
 
     # Split inpainting mask into independent regions
diff --git a/skimage/restoration/rolling_ball.py b/skimage/restoration/rolling_ball.py
index 5142ac0e08c..1149fad0d92 100644
--- a/skimage/restoration/rolling_ball.py
+++ b/skimage/restoration/rolling_ball.py
@@ -1,6 +1,6 @@
 import numpy as np
 
-from .._shared.utils import _float_type
+from .._shared.utils import _supported_float_type
 from ._rolling_ball_cy import apply_kernel, apply_kernel_nan
 
 
@@ -79,7 +79,7 @@ def rolling_ball(image, *, radius=100, kernel=None,
     """
 
     image = np.asarray(image)
-    float_type = _float_type(image)
+    float_type = _supported_float_type(image.dtype)
     img = image.astype(float_type, copy=False)
 
     if num_threads is None:
diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py
index 936018c0530..7388509f3ee 100644
--- a/skimage/restoration/tests/test_denoise.py
+++ b/skimage/restoration/tests/test_denoise.py
@@ -11,7 +11,7 @@
 from skimage._shared import testing
 from skimage._shared.testing import (assert_equal, assert_almost_equal,
                                      assert_warns, assert_)
-from skimage._shared.utils import _float_type
+from skimage._shared.utils import _supported_float_type
 from skimage._shared._warnings import expected_warnings
 from distutils.version import LooseVersion as Version
 
@@ -36,7 +36,8 @@
 astro_odd = astro[:, :-1]
 
 
-@testing.parametrize('dtype', [np.float32, np.float64])
+@testing.parametrize('dtype',
+                     [np.float16, np.float32, np.float64, np.float128])
 def test_denoise_tv_chambolle_2d(dtype):
     # astronaut image
     img = astro_gray.astype(dtype, copy=True)
@@ -46,15 +47,19 @@ def test_denoise_tv_chambolle_2d(dtype):
     img = np.clip(img, 0, 1)
     # denoise
     denoised_astro = restoration.denoise_tv_chambolle(img, weight=0.1)
-    assert denoised_astro.dtype == _float_type(img)
-    # which dtype?
-    assert_(denoised_astro.dtype in [float, np.float32, np.float64])
+    assert denoised_astro.dtype == _supported_float_type(img)
+
     from scipy import ndimage as ndi
+
+    # Convert to a floating point type supported by scipy.ndimage
+    float_dtype = _supported_float_type(img)
+    img = img.astype(float_dtype, copy=False)
+
     grad = ndi.morphological_gradient(img, size=((3, 3)))
     grad_denoised = ndi.morphological_gradient(denoised_astro, size=((3, 3)))
     # test if the total variation has decreased
-    assert_(grad_denoised.dtype == _float_type(img))
-    assert_(np.sqrt((grad_denoised**2).sum()) < np.sqrt((grad**2).sum()))
+    assert grad_denoised.dtype == float_dtype
+    assert np.sqrt((grad_denoised**2).sum()) < np.sqrt((grad**2).sum())
 
 
 def test_denoise_tv_chambolle_multichannel():
@@ -557,7 +562,7 @@ def test_wavelet_denoising_scaling(case, dtype, convert2ycbcr,
                                            multichannel=multichannel,
                                            convert2ycbcr=convert2ycbcr,
                                            rescale_sigma=True)
-    assert denoised.dtype == _float_type(noisy)
+    assert denoised.dtype == _supported_float_type(noisy)
 
     data_range = x.max() - x.min()
     psnr_noisy = peak_signal_noise_ratio(x, noisy, data_range=data_range)
diff --git a/skimage/restoration/tests/test_inpaint.py b/skimage/restoration/tests/test_inpaint.py
index 0019343152b..60845ccf4f6 100644
--- a/skimage/restoration/tests/test_inpaint.py
+++ b/skimage/restoration/tests/test_inpaint.py
@@ -4,9 +4,11 @@
 
 from skimage._shared import testing
 from skimage._shared.testing import assert_allclose
+from skimage._shared.utils import _supported_float_type
 
 
-@testing.parametrize('dtype', [np.float32, np.float64])
+@testing.parametrize('dtype',
+                     [np.float16, np.float32, np.float64, np.float128])
 def test_inpaint_biharmonic_2d(dtype):
     img = np.tile(np.square(np.linspace(0, 1, 5, dtype=dtype)), (5, 1))
     mask = np.zeros_like(img)
@@ -15,7 +17,7 @@ def test_inpaint_biharmonic_2d(dtype):
     mask[0, 4:] = 1
     img[np.where(mask)] = 0
     out = inpaint.inpaint_biharmonic(img, mask)
-    assert out.dtype == img.dtype
+    assert out.dtype == _supported_float_type(img)
     ref = np.array(
         [[0., 0.0625, 0.25000000, 0.5625000, 0.73925058],
          [0., 0.0625, 0.25000000, 0.5478048, 0.76557821],
diff --git a/skimage/restoration/tests/test_restoration.py b/skimage/restoration/tests/test_restoration.py
index 6313798a2ef..be0741084f2 100644
--- a/skimage/restoration/tests/test_restoration.py
+++ b/skimage/restoration/tests/test_restoration.py
@@ -4,6 +4,7 @@
 from scipy import ndimage as ndi
 from skimage._shared import testing
 from skimage._shared.testing import fetch
+from skimage._shared.utils import _supported_float_type
 from skimage.color import rgb2gray
 from skimage.data import astronaut, camera
 from skimage import restoration, util
@@ -12,7 +13,19 @@
 test_img = util.img_as_float(camera())
 
 
-@testing.parametrize('dtype', [np.float32, np.float64])
+def _get_rtol_atol(dtype):
+    rtol = 1e-3
+    atol = 0
+    if dtype == np.float16:
+        rtol = 1e-2
+        atol = 1e-3
+    elif dtype == np.float32:
+        atol = 1e-5
+    return rtol, atol
+
+
+@testing.parametrize('dtype',
+                     [np.float16, np.float32, np.float64, np.float128])
 def test_wiener(dtype):
     psf = np.ones((5, 5), dtype=dtype) / 25
     data = convolve2d(test_img, psf, 'same')
@@ -20,26 +33,27 @@ def test_wiener(dtype):
     data += 0.1 * data.std() * np.random.standard_normal(data.shape)
     data = data.astype(dtype, copy=False)
     deconvolved = restoration.wiener(data, psf, 0.05)
-    assert deconvolved.dtype == dtype
+    assert deconvolved.dtype == _supported_float_type(dtype)
 
+    rtol, atol = _get_rtol_atol(dtype)
     path = fetch('restoration/tests/camera_wiener.npy')
-    atol = 1e-5 if dtype == np.float32 else 0
-    np.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-3,
+    np.testing.assert_allclose(deconvolved, np.load(path), rtol=rtol,
                                atol=atol)
 
     _, laplacian = uft.laplacian(2, data.shape)
     otf = uft.ir2tf(psf, data.shape, is_real=False)
-    assert otf.real.dtype == dtype
+    assert otf.real.dtype == _supported_float_type(dtype)
     deconvolved = restoration.wiener(data, otf, 0.05,
                                      reg=laplacian,
                                      is_real=False)
-    assert deconvolved.real.dtype == dtype
+    assert deconvolved.real.dtype == _supported_float_type(dtype)
     np.testing.assert_allclose(np.real(deconvolved),
                                np.load(path),
-                               rtol=1e-3, atol=atol)
+                               rtol=rtol, atol=atol)
 
 
-@testing.parametrize('dtype', [np.float32, np.float64])
+@testing.parametrize('dtype',
+                     [np.float16, np.float32, np.float64, np.float128])
 def test_unsupervised_wiener(dtype):
     psf = np.ones((5, 5), dtype=dtype) / 25
     data = convolve2d(test_img, psf, 'same')
@@ -47,25 +61,25 @@ def test_unsupervised_wiener(dtype):
     data += 0.1 * data.std() * np.random.standard_normal(data.shape)
     data = data.astype(dtype, copy=False)
     deconvolved, _ = restoration.unsupervised_wiener(data, psf)
-    assert deconvolved.dtype == dtype
+    assert deconvolved.dtype == _supported_float_type(dtype)
 
+    rtol, atol = _get_rtol_atol(dtype)
     path = fetch('restoration/tests/camera_unsup.npy')
-    atol = 1e-5 if dtype == np.float32 else 0
-    np.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-3,
+    np.testing.assert_allclose(deconvolved, np.load(path), rtol=rtol,
                                atol=atol)
 
     _, laplacian = uft.laplacian(2, data.shape)
     otf = uft.ir2tf(psf, data.shape, is_real=False)
-    assert otf.real.dtype == dtype
+    assert otf.real.dtype == _supported_float_type(dtype)
     np.random.seed(0)
     deconvolved = restoration.unsupervised_wiener(
         data, otf, reg=laplacian, is_real=False,
         user_params={"callback": lambda x: None})[0]
-    assert deconvolved.real.dtype == dtype
+    assert deconvolved.real.dtype == _supported_float_type(dtype)
     path = fetch('restoration/tests/camera_unsup2.npy')
     np.testing.assert_allclose(np.real(deconvolved),
                                np.load(path),
-                               rtol=1e-3, atol=atol)
+                               rtol=rtol, atol=atol)
 
 
 def test_image_shape():
@@ -102,7 +116,8 @@ def test_richardson_lucy():
     np.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-3)
 
 
-@pytest.mark.parametrize('dtype_image', [np.float32, np.float64])
+@pytest.mark.parametrize('dtype_image', [np.float16, np.float32, np.float64,
+                                         np.float128])
 @pytest.mark.parametrize('dtype_psf', [np.float32, np.float64])
 def test_richardson_lucy_filtered(dtype_image, dtype_psf):
     if dtype_image == np.float64:
@@ -117,7 +132,7 @@ def test_richardson_lucy_filtered(dtype_image, dtype_psf):
 
     deconvolved = restoration.richardson_lucy(data, psf, 5,
                                               filter_epsilon=1e-6)
-    assert deconvolved.dtype == data.dtype
+    assert deconvolved.dtype == _supported_float_type(data.dtype)
 
     path = fetch('restoration/tests/astronaut_rl.npy')
     np.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-3,
diff --git a/skimage/restoration/tests/test_rolling_ball.py b/skimage/restoration/tests/test_rolling_ball.py
index a8dd87a5fcf..a3a51f902b7 100644
--- a/skimage/restoration/tests/test_rolling_ball.py
+++ b/skimage/restoration/tests/test_rolling_ball.py
@@ -12,20 +12,15 @@
 from skimage.restoration.rolling_ball import ellipsoid_kernel
 
 
-def test_ellipsoid_const():
-    img = 155 * np.ones((100, 100), dtype=np.uint8)
+@testing.parametrize('dtype',
+                     [np.uint8, np.int32, np.float16, np.float32, np.float64,
+                      np.float128])
+def test_ellipsoid_const(dtype):
+    img = 155 * np.ones((100, 100), dtype=dtype)
     kernel = ellipsoid_kernel((25, 53), 50)
     background = rolling_ball(img, kernel=kernel)
     assert np.allclose(img - background, np.zeros_like(img))
-
-
-@testing.parametrize('dtype', [np.float32, np.float64])
-def test_rolling_ball_dtype(dtype):
-    img = 155 * np.ones((100, 100), dtype=dtype)
-    kernel = ellipsoid_kernel((25, 53), 50).astype(dtype, copy=False)
-    background = rolling_ball(img, kernel=kernel)
     assert background.dtype == img.dtype
-    assert np.allclose(img - background, np.zeros_like(img))
 
 
 def test_nan_const():
diff --git a/skimage/restoration/uft.py b/skimage/restoration/uft.py
index 631df9ad1f7..59e19f03084 100644
--- a/skimage/restoration/uft.py
+++ b/skimage/restoration/uft.py
@@ -22,6 +22,7 @@
 
 import numpy as np
 from .._shared.fft import fftmodule as fft
+from .._shared.utils import _supported_float_type
 
 __keywords__ = "fft, Fourier Transform, orthonormal, unitary"
 
@@ -390,7 +391,7 @@ def ir2tf(imp_resp, shape, dim=None, is_real=True):
     if not dim:
         dim = imp_resp.ndim
     # Zero padding and fill
-    irpadded_dtype = imp_resp.dtype if imp_resp.dtype.kind == 'f' else float
+    irpadded_dtype = _supported_float_type(imp_resp)
     irpadded = np.zeros(shape, dtype=irpadded_dtype)
     irpadded[tuple([slice(0, s) for s in imp_resp.shape])] = imp_resp
     # Roll for zero convention of the fft to avoid the phase
@@ -400,10 +401,9 @@ def ir2tf(imp_resp, shape, dim=None, is_real=True):
             irpadded = np.roll(irpadded,
                                shift=-int(np.floor(axis_size / 2)),
                                axis=axis)
-    if is_real:
-        out = fft.rfftn(irpadded, axes=range(-dim, 0))
-    else:
-        out = fft.fftn(irpadded, axes=range(-dim, 0))
+
+    func = fft.rfftn if is_real else fft.fftn
+    out = func(irpadded, axes=(range(-dim, 0)))
 
     # TODO: remove .astype call onces SciPy >= 1.4 is required
     cplx_dtype = np.promote_types(irpadded_dtype, np.complex64)

From fb83eed72aeeebcd5e7c4f6336b8337ecc1db85d Mon Sep 17 00:00:00 2001
From: Gregory Lee <grlee77@gmail.com>
Date: Tue, 13 Apr 2021 17:47:03 -0400
Subject: [PATCH 04/12] float128 unavailable on some platforms

---
 skimage/restoration/tests/test_denoise.py      | 10 ++++++++--
 skimage/restoration/tests/test_inpaint.py      |  2 +-
 skimage/restoration/tests/test_restoration.py  |  2 +-
 skimage/restoration/tests/test_rolling_ball.py |  3 +--
 4 files changed, 11 insertions(+), 6 deletions(-)

diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py
index 7388509f3ee..2c5fa05d4de 100644
--- a/skimage/restoration/tests/test_denoise.py
+++ b/skimage/restoration/tests/test_denoise.py
@@ -36,8 +36,14 @@
 astro_odd = astro[:, :-1]
 
 
-@testing.parametrize('dtype',
-                     [np.float16, np.float32, np.float64, np.float128])
+float_dtypes = [np.float16, np.float32, np.float64]
+try:
+    float_dtypes += [np.float128]
+except AttributeError:
+    pass
+
+
+@testing.parametrize('dtype',float_dtypes)
 def test_denoise_tv_chambolle_2d(dtype):
     # astronaut image
     img = astro_gray.astype(dtype, copy=True)
diff --git a/skimage/restoration/tests/test_inpaint.py b/skimage/restoration/tests/test_inpaint.py
index 60845ccf4f6..16f79c2c510 100644
--- a/skimage/restoration/tests/test_inpaint.py
+++ b/skimage/restoration/tests/test_inpaint.py
@@ -8,7 +8,7 @@
 
 
 @testing.parametrize('dtype',
-                     [np.float16, np.float32, np.float64, np.float128])
+                     [np.float16, np.float32, np.float64])
 def test_inpaint_biharmonic_2d(dtype):
     img = np.tile(np.square(np.linspace(0, 1, 5, dtype=dtype)), (5, 1))
     mask = np.zeros_like(img)
diff --git a/skimage/restoration/tests/test_restoration.py b/skimage/restoration/tests/test_restoration.py
index be0741084f2..b8fb1e55607 100644
--- a/skimage/restoration/tests/test_restoration.py
+++ b/skimage/restoration/tests/test_restoration.py
@@ -25,7 +25,7 @@ def _get_rtol_atol(dtype):
 
 
 @testing.parametrize('dtype',
-                     [np.float16, np.float32, np.float64, np.float128])
+                     [np.float16, np.float32, np.float64])
 def test_wiener(dtype):
     psf = np.ones((5, 5), dtype=dtype) / 25
     data = convolve2d(test_img, psf, 'same')
diff --git a/skimage/restoration/tests/test_rolling_ball.py b/skimage/restoration/tests/test_rolling_ball.py
index a3a51f902b7..520ca9c451d 100644
--- a/skimage/restoration/tests/test_rolling_ball.py
+++ b/skimage/restoration/tests/test_rolling_ball.py
@@ -13,8 +13,7 @@
 
 
 @testing.parametrize('dtype',
-                     [np.uint8, np.int32, np.float16, np.float32, np.float64,
-                      np.float128])
+                     [np.uint8, np.int32, np.float16, np.float32, np.float64])
 def test_ellipsoid_const(dtype):
     img = 155 * np.ones((100, 100), dtype=dtype)
     kernel = ellipsoid_kernel((25, 53), 50)

From f1cb2d6655bbe231a570059a68f968fdb0924c65 Mon Sep 17 00:00:00 2001
From: Gregory Lee <grlee77@gmail.com>
Date: Tue, 13 Apr 2021 17:58:39 -0400
Subject: [PATCH 05/12] remove additional float128 case

---
 skimage/restoration/tests/test_restoration.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/skimage/restoration/tests/test_restoration.py b/skimage/restoration/tests/test_restoration.py
index b8fb1e55607..71a24ba8568 100644
--- a/skimage/restoration/tests/test_restoration.py
+++ b/skimage/restoration/tests/test_restoration.py
@@ -53,7 +53,7 @@ def test_wiener(dtype):
 
 
 @testing.parametrize('dtype',
-                     [np.float16, np.float32, np.float64, np.float128])
+                     [np.float16, np.float32, np.float64])
 def test_unsupervised_wiener(dtype):
     psf = np.ones((5, 5), dtype=dtype) / 25
     data = convolve2d(test_img, psf, 'same')
@@ -116,8 +116,7 @@ def test_richardson_lucy():
     np.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-3)
 
 
-@pytest.mark.parametrize('dtype_image', [np.float16, np.float32, np.float64,
-                                         np.float128])
+@pytest.mark.parametrize('dtype_image', [np.float16, np.float32, np.float64])
 @pytest.mark.parametrize('dtype_psf', [np.float32, np.float64])
 def test_richardson_lucy_filtered(dtype_image, dtype_psf):
     if dtype_image == np.float64:

From 9561041f87d22c0b58c933ae931aadfcff42b1a5 Mon Sep 17 00:00:00 2001
From: Gregory Lee <grlee77@gmail.com>
Date: Wed, 12 May 2021 17:58:18 -0400
Subject: [PATCH 06/12] apply style suggestion

---
 skimage/restoration/deconvolution.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/skimage/restoration/deconvolution.py b/skimage/restoration/deconvolution.py
index fbf615ef782..d9f3bf23a61 100644
--- a/skimage/restoration/deconvolution.py
+++ b/skimage/restoration/deconvolution.py
@@ -290,10 +290,11 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True,
         # weighting (correlation in direct space)
         precision = gn_chain[-1] * atf2 + gx_chain[-1] * areg2  # Eq. 29
         _rand1 = np.random.standard_normal(data_spectrum.shape)
+        _rand1 = _rand1.astype(float_type, copy=False)
         _rand2 = np.random.standard_normal(data_spectrum.shape)
+        _rand2 = _rand2.astype(float_type, copy=False)
         excursion = np.sqrt(0.5) / np.sqrt(precision) * (
-            _rand1.astype(float_type, copy=False)
-            + 1j * _rand2.astype(float_type, copy=False))
+            _rank1 + 1j * _rand2)
 
         # mean Eq. 30 (RLS for fixed gn, gamma0 and gamma1 ...)
         wiener_filter = gn_chain[-1] * np.conj(trans_fct) / precision

From cba69bbf12ed3a8dd12e4377f6a0820ab97c8d68 Mon Sep 17 00:00:00 2001
From: Gregory Lee <grlee77@gmail.com>
Date: Wed, 12 May 2021 18:02:45 -0400
Subject: [PATCH 07/12] fix _invariant_denoise for float16 inputs

---
 skimage/restoration/j_invariant.py            | 8 ++++++++
 skimage/restoration/tests/test_j_invariant.py | 5 +++--
 2 files changed, 11 insertions(+), 2 deletions(-)

diff --git a/skimage/restoration/j_invariant.py b/skimage/restoration/j_invariant.py
index 4883d0a513d..7f5f212bd71 100644
--- a/skimage/restoration/j_invariant.py
+++ b/skimage/restoration/j_invariant.py
@@ -4,10 +4,12 @@
 import numpy as np
 from scipy import ndimage as ndi
 
+from .._shared.utils import _supported_float_type
 from ..metrics import mean_squared_error
 from ..util import img_as_float
 
 
+
 def _interpolate_image(image, *, multichannel=False):
     """Replacing each pixel in ``image`` with the average of its neighbors.
 
@@ -112,9 +114,15 @@ def _invariant_denoise(image, denoise_function, *, stride=4,
         Denoised image, of same shape as `image`.
     """
     image = img_as_float(image)
+
+    # promote float16->float32 if needed
+    float_dtype = _supported_float_type(image.dtype)
+    image = image.astype(float_dtype, copy=False)
+
     if denoiser_kwargs is None:
         denoiser_kwargs = {}
 
+
     if 'multichannel' in denoiser_kwargs:
         multichannel = denoiser_kwargs['multichannel']
     else:
diff --git a/skimage/restoration/tests/test_j_invariant.py b/skimage/restoration/tests/test_j_invariant.py
index ea362c8c23b..dbe37aa97fc 100644
--- a/skimage/restoration/tests/test_j_invariant.py
+++ b/skimage/restoration/tests/test_j_invariant.py
@@ -3,6 +3,7 @@
 import pytest
 
 from skimage._shared.testing import assert_, expected_warnings
+from skimage._shared.utils import _supported_float_type
 from skimage.data import binary_blobs
 from skimage.data import camera, chelsea
 from skimage.metrics import mean_squared_error as mse
@@ -29,7 +30,7 @@ def test_invariant_denoise():
     assert_(denoised_mse < original_mse)
 
 
-@pytest.mark.parametrize('dtype', [np.float32, np.float64])
+@pytest.mark.parametrize('dtype', [np.float16, np.float32, np.float64])
 def test_invariant_denoise_color(dtype):
     denoised_img_color = _invariant_denoise(
         noisy_img_color.astype(dtype), _denoise_wavelet,
@@ -37,7 +38,7 @@ def test_invariant_denoise_color(dtype):
     denoised_mse = mse(denoised_img_color, test_img_color)
     original_mse = mse(noisy_img_color, test_img_color)
     assert denoised_mse < original_mse
-    assert denoised_img_color.dtype == dtype
+    assert denoised_img_color.dtype == _supported_float_type(dtype)
 
 
 def test_invariant_denoise_color_deprecated():

From b713a87a46d58a8613abd296e47fcda81447d733 Mon Sep 17 00:00:00 2001
From: Gregory Lee <grlee77@gmail.com>
Date: Wed, 12 May 2021 18:29:47 -0400
Subject: [PATCH 08/12] fix typo

---
 skimage/restoration/deconvolution.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/restoration/deconvolution.py b/skimage/restoration/deconvolution.py
index d9f3bf23a61..c3a7bc8f97f 100644
--- a/skimage/restoration/deconvolution.py
+++ b/skimage/restoration/deconvolution.py
@@ -294,7 +294,7 @@ def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True,
         _rand2 = np.random.standard_normal(data_spectrum.shape)
         _rand2 = _rand2.astype(float_type, copy=False)
         excursion = np.sqrt(0.5) / np.sqrt(precision) * (
-            _rank1 + 1j * _rand2)
+            _rand1 + 1j * _rand2)
 
         # mean Eq. 30 (RLS for fixed gn, gamma0 and gamma1 ...)
         wiener_filter = gn_chain[-1] * np.conj(trans_fct) / precision

From d32f0731f282e1497fe6fa1af6cbfe6935a859be Mon Sep 17 00:00:00 2001
From: Gregory Lee <grlee77@gmail.com>
Date: Wed, 12 May 2021 19:15:13 -0400
Subject: [PATCH 09/12] Update test_denoise_nl_means_multichannel

- use a test image with more structure
- fix the random seed
- tune denoising parameters a bit

With this updates the difference between 3D vs 2D is more pronounced.

test float16 dtype
---
 skimage/restoration/non_local_means.py    |  4 +--
 skimage/restoration/tests/test_denoise.py | 30 ++++++++++++++---------
 2 files changed, 21 insertions(+), 13 deletions(-)

diff --git a/skimage/restoration/non_local_means.py b/skimage/restoration/non_local_means.py
index cf574d18545..ee9ef2bf61a 100644
--- a/skimage/restoration/non_local_means.py
+++ b/skimage/restoration/non_local_means.py
@@ -161,10 +161,10 @@ def denoise_nl_means(image, patch_size=7, patch_distance=11, h=0.1,
         preserve_range = True
 
     image = convert_to_float(image, preserve_range)
-
-    kwargs = dict(s=patch_size, d=patch_distance, h=h, var=sigma * sigma)
     if not image.flags.c_contiguous:
         image = np.ascontiguousarray(image)
+
+    kwargs = dict(s=patch_size, d=patch_distance, h=h, var=sigma * sigma)
     if channel_axis is not None:  # 2-D images
         if fast_mode:
             return _fast_nl_means_denoising_2d(image, **kwargs)
diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py
index 85732dc382a..c1e40d4328f 100644
--- a/skimage/restoration/tests/test_denoise.py
+++ b/skimage/restoration/tests/test_denoise.py
@@ -480,28 +480,36 @@ def test_denoise_nl_means_3d(fast_mode, dtype):
 
 
 @pytest.mark.parametrize('fast_mode', [False, True])
-@pytest.mark.parametrize('dtype', ['float64', 'float32'])
+@pytest.mark.parametrize('dtype', ['float64', 'float32', 'float16'])
 @pytest.mark.parametrize('channel_axis', [0, -1])
 def test_denoise_nl_means_multichannel(fast_mode, dtype, channel_axis):
     # for true 3D data, 3D denoising is better than denoising as 2D+channels
-    img = np.zeros((13, 10, 8), dtype=dtype)
-    img[6, 4:6, 2:-2] = 1.
-    sigma = 0.3
-    imgn = img + sigma * np.random.randn(*img.shape)
+    dtype = np.float64
+    rstate = np.random.RandomState(5)
+
+    # synthetic 3d volume
+    img = data.binary_blobs(length=32, n_dim=3, seed=5)
+    img = img[:, :24, :16].astype(dtype, copy=False)
+
+    sigma = 0.2
+    imgn = img + sigma * rstate.randn(*img.shape)
     imgn = imgn.astype(dtype)
+
+    # test 3D denoising (channel_axis = None)
+    denoised_ok_multichannel = restoration.denoise_nl_means(
+        imgn, 3, 2, h=0.6 * sigma, sigma=sigma, fast_mode=fast_mode,
+        channel_axis=None)
+
+    # set a channel axis: one dimension is (incorrectly) considered "channels"
     imgn = np.moveaxis(imgn, -1, channel_axis)
     denoised_wrong_multichannel = restoration.denoise_nl_means(
-        imgn, 3, 4, 0.6 * sigma, fast_mode=fast_mode,
+        imgn, 3, 2, h=0.6 * sigma, sigma=sigma, fast_mode=fast_mode,
         channel_axis=channel_axis
     )
-    denoised_ok_multichannel = restoration.denoise_nl_means(
-        imgn, 3, 4, 0.6 * sigma, fast_mode=fast_mode, channel_axis=None)
     denoised_wrong_multichannel = np.moveaxis(
         denoised_wrong_multichannel, channel_axis, -1
     )
-    denoised_ok_multichannel = np.moveaxis(
-        denoised_ok_multichannel, channel_axis, -1
-    )
+
     psnr_wrong = peak_signal_noise_ratio(img, denoised_wrong_multichannel)
     psnr_ok = peak_signal_noise_ratio(img, denoised_ok_multichannel)
     assert_(psnr_ok > psnr_wrong)

From bdb5adfb41a949c7b1972003e09b30dc175988eb Mon Sep 17 00:00:00 2001
From: "Gregory R. Lee" <grlee77@gmail.com>
Date: Fri, 21 May 2021 10:32:04 -0400
Subject: [PATCH 10/12] Apply suggestions from code review

Co-authored-by: Riadh Fezzani <rfezzani@gmail.com>
---
 skimage/restoration/_denoise.py      | 4 ++--
 skimage/restoration/deconvolution.py | 1 -
 skimage/restoration/uft.py           | 1 -
 3 files changed, 2 insertions(+), 4 deletions(-)

diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index e6143d0ee01..5bb4fd66e6d 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -10,8 +10,8 @@
 from .._shared import utils
 from .._shared.utils import _supported_float_type, warn
 from ._denoise_cy import _denoise_bilateral, _denoise_tv_bregman
-from skimage import color
-from skimage.color.colorconv import ycbcr_from_rgb
+from .. import color
+from ..color.colorconv import ycbcr_from_rgb
 
 
 def _gaussian_weight(array, sigma_squared, *, dtype=float):
diff --git a/skimage/restoration/deconvolution.py b/skimage/restoration/deconvolution.py
index c3a7bc8f97f..40859fb141c 100644
--- a/skimage/restoration/deconvolution.py
+++ b/skimage/restoration/deconvolution.py
@@ -8,7 +8,6 @@
 from . import uft
 from .._shared.utils import _supported_float_type
 
-__keywords__ = "restoration, image, deconvolution"
 
 
 def wiener(image, psf, balance, reg=None, is_real=True, clip=True):
diff --git a/skimage/restoration/uft.py b/skimage/restoration/uft.py
index 59e19f03084..1cc4d91b712 100644
--- a/skimage/restoration/uft.py
+++ b/skimage/restoration/uft.py
@@ -24,7 +24,6 @@
 from .._shared.fft import fftmodule as fft
 from .._shared.utils import _supported_float_type
 
-__keywords__ = "fft, Fourier Transform, orthonormal, unitary"
 
 
 def ufftn(inarray, dim=None):

From 67efd1cc3484f03127d73b9941199e66c4594ca3 Mon Sep 17 00:00:00 2001
From: Gregory Lee <grlee77@gmail.com>
Date: Fri, 21 May 2021 10:34:44 -0400
Subject: [PATCH 11/12] prefer pytest.mark.parameterize

---
 skimage/restoration/tests/test_restoration.py  | 7 ++-----
 skimage/restoration/tests/test_rolling_ball.py | 7 ++++---
 2 files changed, 6 insertions(+), 8 deletions(-)

diff --git a/skimage/restoration/tests/test_restoration.py b/skimage/restoration/tests/test_restoration.py
index 71a24ba8568..3d1e4056eaf 100644
--- a/skimage/restoration/tests/test_restoration.py
+++ b/skimage/restoration/tests/test_restoration.py
@@ -2,7 +2,6 @@
 import pytest
 from scipy.signal import convolve2d
 from scipy import ndimage as ndi
-from skimage._shared import testing
 from skimage._shared.testing import fetch
 from skimage._shared.utils import _supported_float_type
 from skimage.color import rgb2gray
@@ -24,8 +23,7 @@ def _get_rtol_atol(dtype):
     return rtol, atol
 
 
-@testing.parametrize('dtype',
-                     [np.float16, np.float32, np.float64])
+@pytest.mark.parametrize('dtype', [np.float16, np.float32, np.float64])
 def test_wiener(dtype):
     psf = np.ones((5, 5), dtype=dtype) / 25
     data = convolve2d(test_img, psf, 'same')
@@ -52,8 +50,7 @@ def test_wiener(dtype):
                                rtol=rtol, atol=atol)
 
 
-@testing.parametrize('dtype',
-                     [np.float16, np.float32, np.float64])
+@pytest.mark.parametrize('dtype', [np.float16, np.float32, np.float64])
 def test_unsupervised_wiener(dtype):
     psf = np.ones((5, 5), dtype=dtype) / 25
     data = convolve2d(test_img, psf, 'same')
diff --git a/skimage/restoration/tests/test_rolling_ball.py b/skimage/restoration/tests/test_rolling_ball.py
index 520ca9c451d..650332ee607 100644
--- a/skimage/restoration/tests/test_rolling_ball.py
+++ b/skimage/restoration/tests/test_rolling_ball.py
@@ -7,13 +7,14 @@
 import pytest
 
 from skimage import data, io
-from skimage._shared import testing
 from skimage.restoration import rolling_ball
 from skimage.restoration.rolling_ball import ellipsoid_kernel
 
 
-@testing.parametrize('dtype',
-                     [np.uint8, np.int32, np.float16, np.float32, np.float64])
+@pytest.mark.parametrize(
+    'dtype',
+    [np.uint8, np.int32, np.float16, np.float32, np.float64]
+)
 def test_ellipsoid_const(dtype):
     img = 155 * np.ones((100, 100), dtype=dtype)
     kernel = ellipsoid_kernel((25, 53), 50)

From 2eae9be733f33d2e1ec3d050d84b9dc09b65d5f0 Mon Sep 17 00:00:00 2001
From: Gregory Lee <grlee77@gmail.com>
Date: Fri, 21 May 2021 10:35:34 -0400
Subject: [PATCH 12/12] prefer pytest.mark.parameterize

---
 skimage/restoration/tests/test_restoration.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/skimage/restoration/tests/test_restoration.py b/skimage/restoration/tests/test_restoration.py
index 3d1e4056eaf..12fd6ed8131 100644
--- a/skimage/restoration/tests/test_restoration.py
+++ b/skimage/restoration/tests/test_restoration.py
@@ -1,12 +1,13 @@
 import numpy as np
 import pytest
-from scipy.signal import convolve2d
 from scipy import ndimage as ndi
+from scipy.signal import convolve2d
+
+from skimage import restoration, util
 from skimage._shared.testing import fetch
 from skimage._shared.utils import _supported_float_type
 from skimage.color import rgb2gray
 from skimage.data import astronaut, camera
-from skimage import restoration, util
 from skimage.restoration import uft
 
 test_img = util.img_as_float(camera())