diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index 86463b69..92d0d10e 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -173,3 +173,6 @@
 
 - Gregor Thalhammer
   Phase unwrapping integration
+
+- François Orieux
+  Image deconvolution http://research.orieux.fr
diff --git a/doc/examples/plot_restoration.py b/doc/examples/plot_restoration.py
new file mode 100644
index 00000000..fdffd953
--- /dev/null
+++ b/doc/examples/plot_restoration.py
@@ -0,0 +1,60 @@
+# -*- coding: utf-8 -*-
+"""
+=====================
+Deconvolution of Lena
+=====================
+
+In this example, we deconvolve a noisy version of Lena using Wiener
+and unsupervised Wiener algorithms. This algorithms are based on
+linear models that can't restore sharp edge as much as non-linear
+methods (like TV restoration) but are much faster.
+
+Wiener filter
+-------------
+The inverse filter based on the PSF (Point Spread Function),
+the prior regularisation (penalisation of high frequency) and the
+tradeoff between the data and prior adequacy. The regularization
+parameter must be hand tuned.
+
+Unsupervised Wiener
+-------------------
+This algorithm has a self-tuned regularisation parameters based on
+data learning. This is not common and based on the following
+publication. The algorithm is based on a iterative Gibbs sampler that
+draw alternatively samples of posterior conditionnal law of the image,
+the noise power and the image frequency power.
+
+.. [1] François Orieux, Jean-François Giovannelli, and Thomas
+       Rodet, "Bayesian estimation of regularization and point
+       spread function parameters for Wiener-Hunt deconvolution",
+       J. Opt. Soc. Am. A 27, 1593-1607 (2010)
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+
+from skimage import color, data, restoration
+
+lena = color.rgb2gray(data.lena())
+from scipy.signal import convolve2d as conv2
+psf = np.ones((5, 5)) / 25
+lena = conv2(lena, psf, 'same')
+lena += 0.1 * lena.std() * np.random.standard_normal(lena.shape)
+
+deconvolved, _ = restoration.unsupervised_wiener(lena, psf)
+
+fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8, 5))
+
+plt.gray()
+
+ax[0].imshow(lena, vmin=deconvolved.min(), vmax=deconvolved.max())
+ax[0].axis('off')
+ax[0].set_title('Data')
+
+ax[1].imshow(deconvolved)
+ax[1].axis('off')
+ax[1].set_title('Self tuned restoration')
+
+fig.subplots_adjust(wspace=0.02, hspace=0.2,
+                    top=0.9, bottom=0.05, left=0, right=1)
+
+plt.show()
diff --git a/skimage/__init__.py b/skimage/__init__.py
index caec81fe..6e44515d 100644
--- a/skimage/__init__.py
+++ b/skimage/__init__.py
@@ -29,6 +29,8 @@ measure
     Measurement of image properties, e.g., similarity and contours.
 morphology
     Morphological operations, e.g. opening or skeletonization.
+restoration
+    Deconvolution algorithms.
 segmentation
     Splitting an image into self-similar regions.
 transform
diff --git a/skimage/restoration/__init__.py b/skimage/restoration/__init__.py
new file mode 100644
index 00000000..324b0932
--- /dev/null
+++ b/skimage/restoration/__init__.py
@@ -0,0 +1,25 @@
+# -*- coding: utf-8 -*-
+"""Image restoration module.
+
+References
+----------
+.. [1] François Orieux, Jean-François Giovannelli, and Thomas
+       Rodet, "Bayesian estimation of regularization and point
+       spread function parameters for Wiener-Hunt deconvolution",
+       J. Opt. Soc. Am. A 27, 1593-1607 (2010)
+
+       http://www.opticsinfobase.org/josaa/abstract.cfm?URI=josaa-27-7-1593
+
+.. [2] Richardson, William Hadley, "Bayesian-Based Iterative Method of
+       Image Restoration". JOSA 62 (1): 55–59. doi:10.1364/JOSA.62.000055, 1972
+
+.. [3] B. R. Hunt "A matrix theory proof of the discrete
+       convolution theorem", IEEE Trans. on Audio and
+       Electroacoustics, vol. au-19, no. 4, pp. 285-288, dec. 1971
+"""
+
+from .deconvolution import wiener, unsupervised_wiener, richardson_lucy
+
+__all__ = ['wiener',
+           "unsupervised_wiener",
+           "richardson_lucy"]
diff --git a/skimage/restoration/deconvolution.py b/skimage/restoration/deconvolution.py
new file mode 100644
index 00000000..51c2c705
--- /dev/null
+++ b/skimage/restoration/deconvolution.py
@@ -0,0 +1,382 @@
+# -*- coding: utf-8 -*-
+# deconvolution.py --- Image deconvolution
+
+"""Implementations restoration functions"""
+
+from __future__ import division
+
+import numpy as np
+import numpy.random as npr
+from scipy.signal import convolve2d
+
+from . import uft
+
+__keywords__ = "restoration, image, deconvolution"
+
+
+def wiener(image, psf, balance, reg=None, is_real=True, clip=True):
+    """Wiener-Hunt deconvolution
+
+    Return the deconvolution with a Wiener-Hunt approach (i.e. with
+    Fourier diagonalisation).
+
+    Parameters
+    ----------
+    image : (M, N) ndarray
+       Input degraded image
+    psf : ndarray
+       Point Spread Function. This is assumed to be the impulse
+       response (input image space) if the data-type is real, or the
+       transfer function (Fourier space) if the data-type is
+       complex. There is no constraints on the shape of the impulse
+       response. The transfer function must be of shape `(M, N)` if
+       `is_real is True`, `(M, N // 2 + 1)` otherwise (see
+       `np.fft.rfftn`).
+    balance : float
+       The regularisation parameter value that tunes the balance
+       between the data adequacy that improve frequency restoration
+       and the prior adequacy that reduce frequency restoration (to
+       avoid noise artifact).
+    reg : ndarray, optional
+       The regularisation operator. The Laplacian by default. It can
+       be an impulse response or a transfer function, as for the
+       psf. Shape constraint is the same than for the `psf` parameter.
+    is_real : boolean, optional
+       True by default. Specify if ``psf`` and ``reg`` are provided
+       with hermitian hypothesis, that is only half of the frequency
+       plane is provided (due to the redundancy of Fourier transform
+       of real signal). It's apply only if ``psf`` and/or ``reg`` are
+       provided as transfer function.  For the hermitian property see
+       ``uft`` module or ``np.fft.rfftn``.
+    clip : boolean, optional
+       True by default. If true, pixel value of the result above 1 or
+       under -1 are thresholded for skimage pipeline
+       compatibility.
+
+    Returns
+    -------
+    im_deconv : (M, N) ndarray
+       The deconvolved image
+
+    Examples
+    --------
+    >>> from skimage import color, data, restoration
+    >>> lena = color.rgb2gray(data.lena())
+    >>> from scipy.signal import convolve2d
+    >>> psf = np.ones((5, 5)) / 25
+    >>> lena = convolve2d(lena, psf, 'same')
+    >>> lena += 0.1 * lena.std() * np.random.standard_normal(lena.shape)
+    >>> deconvolved_lena = restoration.wiener(lena, psf, 1100)
+
+    Notes
+    -----
+    This function applies the Wiener filter to a noisy and degraded
+    image by an impulse response (or PSF). If the data model is
+
+    .. math:: y = Hx + n
+
+    where :math:`n` is noise, :math:`H` the PSF and :math:`x` the
+    unknown original image, the Wiener filter is
+
+    .. math::
+       \hat x = F^\dag (|\Lambda_H|^2 + \lambda |\Lambda_D|^2)
+       \Lambda_H^\dag F y
+
+    where :math:`F` and :math:`F^\dag` are the Fourier and inverse
+    Fourier transfroms respectively, :math:`\Lambda_H` the transfer
+    function (or the Fourier transfrom of the PSF, see [Hunt] below)
+    and :math:`\Lambda_D` the filter to penalize the restored image
+    frequencies (Laplacian by default, that is penalization of high
+    frequency). The parameter :math:`\lambda` tunes the balance
+    between the data (that tends to increase high frequency, even
+    those coming from noise), and the regularization.
+
+    These methods are then specific to a prior model. Consequently,
+    the application or the true image nature must corresponds to the
+    prior model. By default, the prior model (Laplacian) introduce
+    image smoothness or pixel correlation. It can also be interpreted
+    as high-frequency penalization to compensate the instability of
+    the solution wrt. data (sometimes called noise amplification or
+    "explosive" solution).
+
+    Finally, the use of Fourier space implies a circulant property of
+    :math:`H`, see [Hunt].
+
+    References
+    ----------
+    .. [1] François Orieux, Jean-François Giovannelli, and Thomas
+           Rodet, "Bayesian estimation of regularization and point
+           spread function parameters for Wiener-Hunt deconvolution",
+           J. Opt. Soc. Am. A 27, 1593-1607 (2010)
+
+           http://www.opticsinfobase.org/josaa/abstract.cfm?URI=josaa-27-7-1593
+
+           http://research.orieux.fr/files/papers/OGR-JOSA10.pdf
+
+    .. [2] B. R. Hunt "A matrix theory proof of the discrete
+           convolution theorem", IEEE Trans. on Audio and
+           Electroacoustics, vol. au-19, no. 4, pp. 285-288, dec. 1971
+    """
+    if reg is None:
+        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)
+
+    if psf.shape != reg.shape:
+        trans_func = uft.ir2tf(psf, image.shape, is_real=is_real)
+    else:
+        trans_func = psf
+
+    wiener_filter = np.conj(trans_func) / (np.abs(trans_func)**2 +
+                                           balance * np.abs(reg)**2)
+    if is_real:
+        deconv = uft.uirfft2(wiener_filter * uft.urfft2(image))
+    else:
+        deconv = uft.uifft2(wiener_filter * uft.ufft2(image))
+
+    if clip:
+        deconv[deconv > 1] = 1
+        deconv[deconv < -1] = -1
+
+    return deconv
+
+
+def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True,
+                        clip=True):
+    """Unsupervised Wiener-Hunt deconvolution
+
+    Return the deconvolution with a Wiener-Hunt approach, where the
+    hyperparameters are automatically estimated. The algorithm is a
+    stochastic iterative process (Gibbs sampler) described in the
+    reference below. See also ``wiener`` function.
+
+    Parameters
+    ----------
+    image : (M, N) ndarray
+       The input degraded image
+    psf : ndarray
+       The impulse response (input image's space) or the transfer
+       function (Fourier space). Both are accepted. The transfer
+       function is recognize as being complex
+       (``np.iscomplexobj(psf)``).
+    reg : ndarray, optional
+       The regularisation operator. The Laplacian by default. It can
+       be an impulse response or a transfer function, as for the psf.
+    user_params : dict
+       dictionary of gibbs parameters. See below.
+    clip : boolean, optional
+       True by default. If true, pixel value of the result above 1 or
+       under -1 are thresholded for skimage pipeline
+       compatibility.
+
+    Returns
+    -------
+    x_postmean : (M, N) ndarray
+       The deconvolved image (the posterior mean).
+    chains : dict
+       The keys ``noise`` and ``prior`` contain the chain list of
+       noise and prior precision respectively.
+
+    Other parameters
+    ----------------
+    The keys of ``user_params`` are:
+
+    threshold : float
+       The stopping criterion: the norm of the difference between to
+       successive approximated solution (empirical mean of object
+       samples, see Notes section). 1e-4 by default.
+    burnin : int
+       The number of sample to ignore to start computation of the
+       mean. 100 by default.
+    min_iter : int
+       The minimum number of iterations. 30 by default.
+    max_iter : int
+       The maximum number of iterations if ``threshold`` is not
+       satisfied. 150 by default.
+    callback : callable (None by default)
+       A user provided callable to which is passed, if the function
+       exists, the current image sample for whatever purpose. The user
+       can store the sample, or compute other moments than the
+       mean. It has no influence on the algorithm execution and is
+       only for inspection.
+
+    Examples
+    --------
+    >>> from skimage import color, data, restoration
+    >>> lena = color.rgb2gray(data.lena())
+    >>> from scipy.signal import convolve2d
+    >>> psf = np.ones((5, 5)) / 25
+    >>> lena = convolve2d(lena, psf, 'same')
+    >>> lena += 0.1 * lena.std() * np.random.standard_normal(lena.shape)
+    >>> deconvolved_lena = restoration.unsupervised_wiener(lena, psf)
+
+    Notes
+    -----
+    The estimated image is design as the posterior mean of a
+    probability law (from a Bayesian analysis). The mean is defined as
+    a sum over all the possible images weighted by their respective
+    probability. Given the size of the problem, the exact sum is not
+    tractable. This algorithm use of MCMC to draw image under the
+    posterior law. The practical idea is to only draw high probable
+    image since they have the biggest contribution to the mean. At the
+    opposite, the lowest probable image are draw less often since
+    their contribution are low. Finally the empirical mean of these
+    samples give us an estimation of the mean, and an exact
+    computation with an infinite sample set.
+
+    References
+    ----------
+    .. [1] François Orieux, Jean-François Giovannelli, and Thomas
+           Rodet, "Bayesian estimation of regularization and point
+           spread function parameters for Wiener-Hunt deconvolution",
+           J. Opt. Soc. Am. A 27, 1593-1607 (2010)
+
+           http://www.opticsinfobase.org/josaa/abstract.cfm?URI=josaa-27-7-1593
+
+           http://research.orieux.fr/files/papers/OGR-JOSA10.pdf
+    """
+    params = {'threshold': 1e-4, 'max_iter': 200,
+              'min_iter': 30, 'burnin': 15, 'callback': None}
+    params.update(user_params or {})
+
+    if reg is None:
+        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)
+
+    if psf.shape != reg.shape:
+        trans_fct = uft.ir2tf(psf, image.shape,  is_real=is_real)
+    else:
+        trans_fct = psf
+
+    # The mean of the object
+    x_postmean = np.zeros(trans_fct.shape)
+    # The previous computed mean in the iterative loop
+    prev_x_postmean = np.zeros(trans_fct.shape)
+
+    # Difference between two successive mean
+    delta = np.NAN
+
+    # Initial state of the chain
+    gn_chain, gx_chain = [1], [1]
+
+    # The correlation of the object in Fourier space (if size is big,
+    # this can reduce computation time in the loop)
+    areg2 = np.abs(reg)**2
+    atf2 = np.abs(trans_fct)**2
+
+    # The Fourier transfrom may change the image.size attribut, so we
+    # store it.
+    if is_real:
+        data_spectrum = uft.urfft2(image.astype(np.float))
+    else:
+        data_spectrum = uft.ufft2(image.astype(np.float))
+
+    # Gibbs sampling
+    for iteration in range(params['max_iter']):
+        # Sample of Eq. 27 p(circX^k | gn^k-1, gx^k-1, y).
+
+        # weighting (correlation in direct space)
+        precision = gn_chain[-1] * atf2 + gx_chain[-1] * areg2  # Eq. 29
+        excursion = np.sqrt(0.5) / np.sqrt(precision) * (
+            np.random.standard_normal(data_spectrum.shape) +
+            1j * np.random.standard_normal(data_spectrum.shape))
+
+        # mean Eq. 30 (RLS for fixed gn, gamma0 and gamma1 ...)
+        wiener_filter = gn_chain[-1] * np.conj(trans_fct) / precision
+
+        # sample of X in Fourier space
+        x_sample = wiener_filter * data_spectrum + excursion
+        if params['callback']:
+            params['callback'](x_sample)
+
+        # sample of Eq. 31 p(gn | x^k, gx^k, y)
+        gn_chain.append(npr.gamma(image.size / 2,
+                                  2 / uft.image_quad_norm(data_spectrum -
+                                                          x_sample *
+                                                          trans_fct)))
+
+        # sample of Eq. 31 p(gx | x^k, gn^k-1, y)
+        gx_chain.append(npr.gamma((image.size - 1) / 2,
+                                  2 / uft.image_quad_norm(x_sample * reg)))
+
+        # current empirical average
+        if iteration > params['burnin']:
+            x_postmean = prev_x_postmean + x_sample
+
+        if iteration > (params['burnin'] + 1):
+            current = x_postmean / (iteration - params['burnin'])
+            previous = prev_x_postmean / (iteration - params['burnin'] - 1)
+
+            delta = np.sum(np.abs(current - previous)) / \
+                np.sum(np.abs(x_postmean)) / (iteration - params['burnin'])
+
+        prev_x_postmean = x_postmean
+
+        # stop of the algorithm
+        if (iteration > params['min_iter']) and (delta < params['threshold']):
+            break
+
+    # Empirical average \approx POSTMEAN Eq. 44
+    x_postmean = x_postmean / (iteration - params['burnin'])
+    if is_real:
+        x_postmean = uft.uirfft2(x_postmean)
+    else:
+        x_postmean = uft.uifft2(x_postmean)
+
+    if clip:
+        x_postmean[x_postmean > 1] = 1
+        x_postmean[x_postmean < -1] = -1
+
+    return (x_postmean, {'noise': gn_chain, 'prior': gx_chain})
+
+
+def richardson_lucy(image, psf, iterations=50, clip=True):
+    """Richardson-Lucy deconvolution.
+
+    Parameters
+    ----------
+    image : ndarray
+       Input degraded image
+    psf : ndarray
+       The point spread function
+    iterations : int
+       Number of iterations. This parameter play to role of
+       regularisation.
+    clip : boolean, optional
+       True by default. If true, pixel value of the result above 1 or
+       under -1 are thresholded for skimage pipeline
+       compatibility.
+
+    Returns
+    -------
+    im_deconv : ndarray
+       The deconvolved image
+
+    Examples
+    --------
+    >>> from skimage import color, data, restoration
+    >>> camera = color.rgb2gray(data.camera())
+    >>> from scipy.signal import convolve2d
+    >>> psf = np.ones((5, 5)) / 25
+    >>> camera = convolve2d(camera, psf, 'same')
+    >>> camera += 0.1 * camera.std() * np.random.standard_normal(camera.shape)
+    >>> deconvolved = restoration.richardson_lucy(camera, psf, 5)
+
+    References
+    ----------
+    .. [2] http://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution
+    """
+    image = image.astype(np.float)
+    psf = psf.astype(np.float)
+    im_deconv = 0.5 * np.ones(image.shape)
+    psf_mirror = psf[::-1, ::-1]
+    for _ in range(iterations):
+        relative_blur = image / convolve2d(im_deconv, psf, 'same')
+        im_deconv *= convolve2d(relative_blur, psf_mirror, 'same')
+
+    if clip:
+        im_deconv[im_deconv > 1] = 1
+        im_deconv[im_deconv < -1] = -1
+
+    return im_deconv
diff --git a/skimage/restoration/tests/camera_rl.npy b/skimage/restoration/tests/camera_rl.npy
new file mode 100644
index 00000000..10995f85
Binary files /dev/null and b/skimage/restoration/tests/camera_rl.npy differ
diff --git a/skimage/restoration/tests/camera_unsup.npy b/skimage/restoration/tests/camera_unsup.npy
new file mode 100644
index 00000000..12d2e8b8
Binary files /dev/null and b/skimage/restoration/tests/camera_unsup.npy differ
diff --git a/skimage/restoration/tests/camera_unsup2.npy b/skimage/restoration/tests/camera_unsup2.npy
new file mode 100644
index 00000000..5bb974db
Binary files /dev/null and b/skimage/restoration/tests/camera_unsup2.npy differ
diff --git a/skimage/restoration/tests/camera_wiener.npy b/skimage/restoration/tests/camera_wiener.npy
new file mode 100644
index 00000000..831b2be8
Binary files /dev/null and b/skimage/restoration/tests/camera_wiener.npy differ
diff --git a/skimage/restoration/tests/test_restoration.py b/skimage/restoration/tests/test_restoration.py
new file mode 100644
index 00000000..915640af
--- /dev/null
+++ b/skimage/restoration/tests/test_restoration.py
@@ -0,0 +1,64 @@
+from os.path import abspath, dirname, join as pjoin
+
+import numpy as np
+from scipy.signal import convolve2d
+
+import skimage
+from skimage.data import camera
+from skimage import restoration
+from skimage.restoration import uft
+
+test_img = skimage.img_as_float(camera())
+
+
+def test_wiener():
+    psf = np.ones((5, 5)) / 25
+    data = convolve2d(test_img, psf, 'same')
+    np.random.seed(0)
+    data += 0.1 * data.std() * np.random.standard_normal(data.shape)
+    deconvolved = restoration.wiener(data, psf, 0.05)
+
+    path = pjoin(dirname(abspath(__file__)), 'camera_wiener.npy')
+    np.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-3)
+
+    _, laplacian = uft.laplacian(2, data.shape)
+    otf = uft.ir2tf(psf, data.shape, is_real=False)
+    deconvolved = restoration.wiener(data, otf, 0.05,
+                                     reg=laplacian,
+                                     is_real=False)
+    np.testing.assert_allclose(np.real(deconvolved),
+                               np.load(path),
+                               rtol=1e-3)
+
+
+def test_unsupervised_wiener():
+    psf = np.ones((5, 5)) / 25
+    data = convolve2d(test_img, psf, 'same')
+    np.random.seed(0)
+    data += 0.1 * data.std() * np.random.standard_normal(data.shape)
+    deconvolved, _ = restoration.unsupervised_wiener(data, psf)
+
+    path = pjoin(dirname(abspath(__file__)), 'camera_unsup.npy')
+    np.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-3)
+
+    _, laplacian = uft.laplacian(2, data.shape)
+    otf = uft.ir2tf(psf, data.shape, is_real=False)
+    np.random.seed(0)
+    deconvolved = restoration.unsupervised_wiener(
+        data, otf, reg=laplacian, is_real=False,
+        user_params={"callback": lambda x: None})[0]
+    path = pjoin(dirname(abspath(__file__)), 'camera_unsup2.npy')
+    np.testing.assert_allclose(np.real(deconvolved),
+                               np.load(path),
+                               rtol=1e-3)
+
+
+def test_richardson_lucy():
+    psf = np.ones((5, 5)) / 25
+    data = convolve2d(test_img, psf, 'same')
+    np.random.seed(0)
+    data += 0.1 * data.std() * np.random.standard_normal(data.shape)
+    deconvolved = restoration.richardson_lucy(data, psf, 5)
+
+    path = pjoin(dirname(abspath(__file__)), 'camera_rl.npy')
+    np.testing.assert_allclose(deconvolved, np.load(path), rtol=1e-3)
diff --git a/skimage/restoration/uft.py b/skimage/restoration/uft.py
new file mode 100644
index 00000000..fedcc6f6
--- /dev/null
+++ b/skimage/restoration/uft.py
@@ -0,0 +1,445 @@
+# -*- coding: utf-8 -*-
+# uft.py --- Unitary fourier transform
+
+"""Function of unitary fourier transform and utilities
+
+This module implement unitary fourier transform, that is ortho-normal
+transform. They are especially and useful for convolution [1]: they
+respect the Parseval equality, the value of the null frequency is
+equal to
+
+.. math::  \frac{1}{\sqrt{n}} \sum_i x_i
+
+or the Fourier tranform have the same energy than the original image
+(see ``image_quad_norm`` function). The transform is applied from the
+last axes for performance reason (c order array). You may use directly
+the numpy.fft module for more sophisticated purpose.
+
+References
+----------
+.. [1] B. R. Hunt "A matrix theory proof of the discrete convolution
+       theorem", IEEE Trans. on Audio and Electroacoustics,
+       vol. au-19, no. 4, pp. 285-288, dec. 1971
+
+"""
+
+from __future__ import division, print_function
+
+import numpy as np
+
+__keywords__ = "fft, Fourier Transform, orthonormal, unitary"
+
+
+def ufftn(inarray, dim=None):
+    """N-dim unitary Fourier transform
+
+    Parameters
+    ----------
+    inarray : ndarray
+        The array to transform.
+    dim : int, optional
+        The ``dim`` last axis along wich to compute the transform. All
+        axes by default.
+
+    Returns
+    -------
+    outarray : ndarray (same shape than inarray)
+        The unitary N-D Fourier transform of ``inarray``.
+
+    Examples
+    --------
+    >>> input = np.ones((3, 3, 3))
+    >>> output = ufftn(input)
+    >>> np.allclose(np.sum(input) / np.sqrt(input.size), output[0, 0, 0])
+    True
+    >>> output.shape
+    (3, 3, 3)
+    """
+    if dim is None:
+        dim = inarray.ndim
+    outarray = np.fft.fftn(inarray, axes=range(-dim, 0))
+    return outarray / np.sqrt(np.prod(inarray.shape[-dim:]))
+
+
+def uifftn(inarray, dim=None):
+    """N-dim unitary inverse Fourier transform
+
+    Parameters
+    ----------
+    inarray : ndarray
+        The array to transform.
+    dim : int, optional
+        The ``dim`` last axis along wich to compute the transform. All
+        axes by default.
+
+    Returns
+    -------
+    outarray : ndarray (same shape than inarray)
+        The unitary inverse N-D Fourier transform of ``inarray``.
+
+    Examples
+    --------
+    >>> input = np.ones((3, 3, 3))
+    >>> output = uifftn(input)
+    >>> np.allclose(np.sum(input) / np.sqrt(input.size), output[0, 0, 0])
+    True
+    >>> output.shape
+    (3, 3, 3)
+    """
+    if dim is None:
+        dim = inarray.ndim
+    outarray = np.fft.ifftn(inarray, axes=range(-dim, 0))
+    return outarray * np.sqrt(np.prod(inarray.shape[-dim:]))
+
+
+def urfftn(inarray, dim=None):
+    """N-dim real unitary Fourier transform
+
+    This transform consider the Hermitian property of the transform on
+    real input
+
+    Parameters
+    ----------
+    inarray : ndarray
+        The array to transform.
+    dim : int, optional
+        The ``dim`` last axis along wich to compute the transform. All
+        axes by default.
+
+    Returns
+    -------
+    outarray : ndarray (the last dim as  N / 2 + 1 lenght)
+        The unitary N-D real Fourier transform of ``inarray``.
+
+    Notes
+    -----
+    The ``urfft`` functions assume an input array of real
+    values. Consequently, the output have an Hermitian property and
+    redondant values are not computed and returned.
+
+    Examples
+    --------
+    >>> input = np.ones((5, 5, 5))
+    >>> output = urfftn(input)
+    >>> np.allclose(np.sum(input) / np.sqrt(input.size), output[0, 0, 0])
+    True
+    >>> output.shape
+    (5, 5, 3)
+    """
+    if dim is None:
+        dim = inarray.ndim
+    outarray = np.fft.rfftn(inarray, axes=range(-dim, 0))
+    return outarray / np.sqrt(np.prod(inarray.shape[-dim:]))
+
+
+def uirfftn(inarray, dim=None, shape=None):
+    """N-dim real unitary Fourier transform
+
+    This transform consider the Hermitian property of the transform
+    from complex to real real input.
+
+    Parameters
+    ----------
+    inarray : ndarray
+        The array to transform.
+    dim : int, optional
+        The ``dim`` last axis along wich to compute the transform. All
+        axes by default.
+    shape : tuple of int
+        The shape of the output. The shape of ``rfft`` is ambiguous in
+        case of odd shape. In this case, the parameter must be
+        used. see ``np.fft.irfftn``.
+
+    Returns
+    -------
+    outarray : ndarray
+        The unitary N-D inverse real Fourier transform of ``inarray``.
+
+    Notes
+    -----
+    The ``uirfft`` function assume that output array is of real
+    values. Consequently, the input is assumed of having an Hermitian
+    property and redondant values are implicit.
+
+    Examples
+    --------
+    >>> input = np.ones((5, 5, 5))
+    >>> output = uirfftn(urfftn(input), shape=input.shape)
+    >>> np.allclose(input, output)
+    True
+    >>> output.shape
+    (5, 5, 5)
+    """
+    if dim is None:
+        dim = inarray.ndim
+    outarray = np.fft.irfftn(inarray, shape, axes=range(-dim, 0))
+    return outarray * np.sqrt(np.prod(outarray.shape[-dim:]))
+
+
+def ufft2(inarray):
+    """2-dim unitary Fourier transform
+
+    Compute the Fourier transform on the last 2 axes.
+
+    Parameters
+    ----------
+    inarray : ndarray
+        The array to transform.
+
+    Returns
+    -------
+    outarray : ndarray (same shape than inarray)
+        The unitary 2-D Fourier transform of ``inarray``.
+
+    See Also
+    --------
+    uifft2, ufftn, urfftn
+
+    Examples
+    --------
+    >>> input = np.ones((10, 128, 128))
+    >>> output = ufft2(input)
+    >>> np.allclose(np.sum(input[1, ...]) / np.sqrt(input[1, ...].size), output[1, 0, 0])
+    True
+    >>> output.shape
+    (10, 128, 128)
+    """
+    return ufftn(inarray, 2)
+
+
+def uifft2(inarray):
+    """2-dim inverse unitary Fourier transform
+
+    Compute the inverse Fourier transform on the last 2 axes.
+
+    Parameters
+    ----------
+    inarray : ndarray
+        The array to transform.
+
+    Returns
+    -------
+    outarray : ndarray (same shape than inarray)
+        The unitary 2-D inverse Fourier transform of ``inarray``.
+
+    See Also
+    --------
+    uifft2, uifftn, uirfftn
+
+    Examples
+    --------
+    >>> input = np.ones((10, 128, 128))
+    >>> output = uifft2(input)
+    >>> np.allclose(np.sum(input[1, ...]) / np.sqrt(input[1, ...].size), output[0, 0, 0])
+    True
+    >>> output.shape
+    (10, 128, 128)
+    """
+    return uifftn(inarray, 2)
+
+
+def urfft2(inarray):
+    """2-dim real unitary Fourier transform
+
+    Compute the real Fourier transform on the last 2 axes. This
+    transform consider the Hermitian property of the transform from
+    complex to real real input.
+
+    Parameters
+    ----------
+    inarray : ndarray
+        The array to transform.
+
+    Returns
+    -------
+    outarray : ndarray (the last dim as (N - 1) *2 lenght)
+        The unitary 2-D real Fourier transform of ``inarray``.
+
+    See Also
+    --------
+    ufft2, ufftn, urfftn
+
+    Examples
+    --------
+    >>> input = np.ones((10, 128, 128))
+    >>> output = urfft2(input)
+    >>> np.allclose(np.sum(input[1,...]) / np.sqrt(input[1,...].size), output[1, 0, 0])
+    True
+    >>> output.shape
+    (10, 128, 65)
+    """
+    return urfftn(inarray, 2)
+
+
+def uirfft2(inarray, shape=None):
+    """2-dim real unitary Fourier transform
+
+    Compute the real inverse Fourier transform on the last 2 axes.
+    This transform consider the Hermitian property of the transform
+    from complex to real real input.
+
+    Parameters
+    ----------
+    inarray : ndarray
+        The array to transform.
+
+    Returns
+    -------
+    outarray : ndarray (the last dim as (N - 1) *2 lenght)
+        The unitary 2-D inverse real Fourier transform of ``inarray``.
+
+    See Also
+    --------
+    urfft2, uifftn, uirfftn
+
+    Examples
+    --------
+    >>> input = np.ones((10, 128, 128))
+    >>> output = uirfftn(urfftn(input), shape=input.shape)
+    >>> np.allclose(input, output)
+    True
+    >>> output.shape
+    (10, 128, 128)
+    """
+    return uirfftn(inarray, 2, shape=shape)
+
+
+def image_quad_norm(inarray):
+    """Return quadratic norm of images in Fourier space
+
+    This function detect if the image suppose the hermitian property.
+
+    Parameters
+    ----------
+    inarray : ndarray
+        The images are supposed to be in the last two axes
+
+    Returns
+    -------
+    norm : float
+        The quadratic norm of ``inarray``.
+
+    Examples
+    --------
+    >>> input = np.ones((5, 5))
+    >>> image_quad_norm(ufft2(input)) == np.sum(np.abs(input)**2)
+    True
+    >>> image_quad_norm(ufft2(input)) == image_quad_norm(urfft2(input))
+    True
+    """
+    # If there is an hermitian symmetry
+    if inarray.shape[-1] != inarray.shape[-2]:
+        return 2 * np.sum(np.sum(np.abs(inarray)**2, axis=-1), axis=-1) - \
+            np.sum(np.abs(inarray[..., 0])**2, axis=-1)
+    else:
+        return np.sum(np.sum(np.abs(inarray)**2, axis=-1), axis=-1)
+
+
+def ir2tf(imp_resp, shape, dim=None, is_real=True):
+    """Compute the transfer function of IR
+
+    This function make the necessary correct zero-padding, zero
+    convention, correct fft2 etc... to compute the transfer function
+    of IR. To use with unitary Fourier transform for the signal (ufftn
+    or equivalent).
+
+    Parameters
+    ----------
+    imp_resp : ndarray
+       The impulsionnal responses.
+    shape : tuple of int
+       A tuple of integer corresponding to the target shape of the
+       tranfert function.
+    dim : int, optional
+        The ``dim`` last axis along wich to compute the transform. All
+        axes by default.
+    is_real : boolean (optionnal, default True)
+       If True, imp_resp is supposed real and the hermissian property
+       is used with rfftn Fourier transform.
+
+    Returns
+    -------
+    y : complex ndarray
+       The tranfert function of shape ``shape``.
+
+    See Also
+    --------
+    ufftn, uifftn, urfftn, uirfftn
+
+    Examples
+    --------
+    >>> np.all(np.array([[4, 0], [0, 0]]) == ir2tf(np.ones((2, 2)), (2, 2)))
+    True
+    >>> ir2tf(np.ones((2, 2)), (512, 512)).shape == (512, 257)
+    True
+    >>> ir2tf(np.ones((2, 2)), (512, 512), is_real=False).shape == (512, 512)
+    True
+
+    Notes
+    -----
+    The input array can be composed of multiple dimentionnal IR with
+    an arbitraru number of IR. The individual IR must be accesed
+    through first axes. The last ``dim`` axes of space definition. The
+    ``dim`` parameter must be specified to compute the transform only
+    along these last axes.
+    """
+    if not dim:
+        dim = imp_resp.ndim
+    # Zero padding and fill
+    irpadded = np.zeros(shape)
+    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.
+    for axis, axis_size in enumerate(imp_resp.shape):
+        if axis >= imp_resp.ndim - dim:
+            irpadded = np.roll(irpadded,
+                               shift=-int(np.floor(axis_size / 2)),
+                               axis=axis)
+    if is_real:
+        return np.fft.rfftn(irpadded, axes=range(-dim, 0))
+    else:
+        return np.fft.fftn(irpadded, axes=range(-dim, 0))
+
+
+def laplacian(ndim, shape, is_real=True):
+    """Return the transfer function of the Laplacian
+
+    Laplacian is the second order difference, on line and column.
+
+    Parameters
+    ----------
+    ndim : int
+        The dimension of the Laplacian
+    shape : tuple, shape
+        The support on which to compute the transfer function
+    is_real : boolean (optionnal, default True)
+       If True, imp_resp is supposed real and the hermissian property
+       is used with rfftn Fourier transform to return the transfer
+       function.
+
+    Returns
+    -------
+    tf : array_like, complex
+        The transfer function
+    impr : array_like, real
+        The Laplacian
+
+    Examples
+    --------
+    >>> tf, ir = laplacian(2, (32, 32))
+    >>> np.all(ir == np.array([[0, -1, 0], [-1, 4, -1], [0, -1, 0]]))
+    True
+    >>> np.all(tf == ir2tf(ir, (32, 32)))
+    True
+    """
+    impr = np.zeros([3] * ndim)
+    for dim in range(ndim):
+        idx = tuple([slice(1, 2)] * dim +
+                    [slice(None)] +
+                    [slice(1, 2)] * (ndim - dim - 1))
+        impr[idx] = np.array([-1.0,
+                              0.0,
+                              -1.0]).reshape([-1 if i == dim else 1
+                                              for i in range(ndim)])
+    impr[([slice(1, 2)] * ndim)] = 2.0 * ndim
+    return ir2tf(impr, shape, is_real=is_real), impr
diff --git a/skimage/setup.py b/skimage/setup.py
index 14f66897..962adb97 100644
--- a/skimage/setup.py
+++ b/skimage/setup.py
@@ -12,6 +12,7 @@ def configuration(parent_package='', top_path=None):
     config.add_subpackage('draw')
     config.add_subpackage('exposure')
     config.add_subpackage('feature')
+    config.add_subpackage('restoration')
     config.add_subpackage('filter')
     config.add_subpackage('graph')
     config.add_subpackage('io')