diff --git a/nifty/operators/__init__.py b/nifty/operators/__init__.py
index 44d0c545552e9a9b0f1a95efc526df2d00e65a5a..ed25e1c63cd1eee1b2b698d64b5abd1348b0843f 100644
--- a/nifty/operators/__init__.py
+++ b/nifty/operators/__init__.py
@@ -3,6 +3,7 @@ from .linear_operator import LinearOperator
from .diagonal_operator import DiagonalOperator
from .endomorphic_operator import EndomorphicOperator
from .fft_smoothing_operator import FFTSmoothingOperator
+from .direct_smoothing_operator import DirectSmoothingOperator
from .fft_operator import FFTOperator
from .inversion_enabler import InversionEnabler
from .composed_operator import ComposedOperator
diff --git a/nifty/operators/direct_smoothing_operator.py b/nifty/operators/direct_smoothing_operator.py
new file mode 100644
index 0000000000000000000000000000000000000000..73ae7c4c67b7a884e7ff144f9c2faca562316764
--- /dev/null
+++ b/nifty/operators/direct_smoothing_operator.py
@@ -0,0 +1,108 @@
+from __future__ import division
+from builtins import range
+import numpy as np
+
+from .endomorphic_operator import EndomorphicOperator
+from .. import dobj
+from ..spaces import PowerSpace
+from .. import utilities
+from .. import Field, DomainTuple
+
+
+class DirectSmoothingOperator(EndomorphicOperator):
+ def __init__(self, domain, sigma, log_distances=False,
+ space=None):
+ super(DirectSmoothingOperator, self).__init__()
+
+ self._domain = DomainTuple.make(domain)
+ self._space = utilities.infer_space(self._domain, space)
+ if not isinstance(self._domain[self._space], PowerSpace):
+ raise TypeError("PowerSpace needed")
+
+ self._sigma = float(sigma)
+ self._log_distances = log_distances
+ self._effective_smoothing_width = 3.01
+
+ self._axis = self._domain.axes[self._space][0]
+ if self._sigma != 0:
+ distances = self._domain[self._space].k_lengths
+ if self._log_distances:
+ distances = np.log(np.maximum(distances, 1e-15))
+
+ self._ibegin, self._nval, self._wgt = self._precompute(distances)
+
+ def _times(self, x):
+ if self._sigma == 0:
+ return x.copy()
+
+ return self._smooth(x)
+
+ @property
+ def domain(self):
+ return self._domain
+
+ @property
+ def self_adjoint(self):
+ return True
+
+ @property
+ def unitary(self):
+ return False
+
+ def _precompute(self, x):
+ """ Does precomputations for Gaussian smoothing on a 1D irregular grid.
+
+ Parameters
+ ----------
+ x: 1D floating point array or list containing the individual grid
+ positions. Points must be given in ascending order.
+
+
+ Returns
+ -------
+ ibegin: integer array of the same size as x
+ ibegin[i] is the minimum grid index to consider when computing the
+ smoothed value at grid index i
+ nval: integer array of the same size as x
+ nval[i] is the number of indices to consider when computing the
+ smoothed value at grid index i.
+ wgt: list with the same number of entries as x
+ wgt[i] is an array with nval[i] entries containing the
+ normalized smoothing weights.
+ """
+ dxmax = self._effective_smoothing_width*self._sigma
+
+ x = np.asarray(x)
+
+ ibegin = np.searchsorted(x, x-dxmax)
+ nval = np.searchsorted(x, x+dxmax) - ibegin
+
+ wgt = []
+ expfac = 1. / (2. * self._sigma*self._sigma)
+ for i in range(x.size):
+ if nval[i] > 0:
+ t = x[ibegin[i]:ibegin[i]+nval[i]]-x[i]
+ t = np.exp(-t*t*expfac)
+ t *= 1./np.sum(t)
+ wgt.append(t)
+ else:
+ wgt.append(np.array([]))
+
+ return ibegin, nval, wgt
+
+ def _smooth(self, x):
+ if dobj.distaxis(x.val) == self._axis:
+ val = dobj.redistribute(x.val, nodist=(self._axis,))
+ else:
+ val = x.val.copy()
+ lval = dobj.local_data(val)
+
+ for sl in utilities.get_slice_list(lval.shape, (self._axis,)):
+ inp = lval[sl]
+ out = np.zeros(inp.shape[0], dtype=inp.dtype)
+ for i in range(inp.shape[0]):
+ out[self._ibegin[i]:self._ibegin[i]+self._nval[i]] += \
+ inp[i] * self._wgt[i][:]
+ lval[sl] = out
+ val = dobj.redistribute(val, dobj.distaxis(x.val))
+ return Field(x.domain, val)
diff --git a/test/test_operators/test_smoothing_operator.py b/test/test_operators/test_smoothing_operator.py
index 79ad0c311aa096d0f5994c60e0c1231c667a5d5c..cbae4195bf470b00642da7cb400f25e71a2e20bb 100644
--- a/test/test_operators/test_smoothing_operator.py
+++ b/test/test_operators/test_smoothing_operator.py
@@ -83,3 +83,29 @@ class SmoothingOperator_Tests(unittest.TestCase):
mean=4, dtype=tp)
out = smo(inp)
assert_allclose(inp.sum(), out.sum(), rtol=tol, atol=tol)
+
+ @expand(product([100, 200], [False, True], [0., 1., 3.7],
+ [np.float64, np.complex128]))
+ def test_smooth_irregular1(self, sz, log, sigma, tp):
+ tol = _get_rtol(tp)
+ sp = ift.RGSpace(sz, harmonic=True)
+ bb = ift.PowerSpace.useful_binbounds(sp, logarithmic=log)
+ ps = ift.PowerSpace(sp, binbounds=bb)
+ smo = ift.DirectSmoothingOperator(ps, sigma=sigma)
+ inp = ift.Field.from_random(domain=ps, random_type='normal', std=1,
+ mean=4, dtype=tp)
+ out = smo(inp)
+ assert_allclose(inp.sum(), out.sum(), rtol=tol, atol=tol)
+
+ @expand(product([10, 15], [7, 10], [False, True], [0., 1., 3.7],
+ [np.float64, np.complex128]))
+ def test_smooth_irregular2(self, sz1, sz2, log, sigma, tp):
+ tol = _get_rtol(tp)
+ sp = ift.RGSpace([sz1, sz2], harmonic=True)
+ bb = ift.PowerSpace.useful_binbounds(sp, logarithmic=log)
+ ps = ift.PowerSpace(sp, binbounds=bb)
+ smo = ift.DirectSmoothingOperator(ps, sigma=sigma)
+ inp = ift.Field.from_random(domain=ps, random_type='normal', std=1,
+ mean=4, dtype=tp)
+ out = smo(inp)
+ assert_allclose(inp.sum(), out.sum(), rtol=tol, atol=tol)