diff --git a/odl/solvers/functional/default_functionals.py b/odl/solvers/functional/default_functionals.py
index d23b43a387d..e9049d2c9a8 100644
--- a/odl/solvers/functional/default_functionals.py
+++ b/odl/solvers/functional/default_functionals.py
@@ -22,7 +22,7 @@
     proximal_l1, proximal_convex_conj_l1, proximal_l2, proximal_convex_conj_l2,
     proximal_l2_squared, proximal_const_func, proximal_box_constraint,
     proximal_convex_conj, proximal_convex_conj_kl,
-    proximal_convex_conj_kl_cross_entropy,
+    proximal_convex_conj_kl_cross_entropy, proximal_huber_norm,
     combine_proximals)
 from odl.util import conj_exponent
 
@@ -35,7 +35,7 @@
            'QuadraticForm',
            'NuclearNorm', 'IndicatorNuclearNormUnitBall',
            'ScalingFunctional', 'IdentityFunctional',
-           'MoreauEnvelope')
+           'MoreauEnvelope', 'HuberNorm')
 
 
 class LpNorm(Functional):
@@ -2292,6 +2292,134 @@ def gradient(self):
                 (1 / self.sigma) * self.functional.proximal(self.sigma))
 
 
+class HuberNorm(Functional):
+
+    """The Huber norm functional.
+
+    Notes
+    -----
+    The functional :math:`F` with smoothing :math:`\\epsilon>0` is given by
+
+    .. math::
+        F(x)
+        =
+        \\begin{cases}
+            \\frac{1}{2 \\epsilon} x^2 & \\text{if } |x| \leq \\epsilon
+            \\\\
+            |x| - \\frac{\\epsilon}{2} & \\text{if } |x| > \\epsilon,
+        \\end{cases}
+
+    and the Huber norm is the integral of this functional over the domain,
+    i.e., the Huber norm is given by
+
+    .. math::
+        \\int_\Omega F(x) dx.
+
+    In the discrete case, this becomes
+
+    .. math::
+        \\sum_{i=1}^n F(x_i).
+
+    """
+
+    def __init__(self, space, epsilon):
+        """Initialize a new instance.
+
+        Parameters
+        ----------
+        space : `DiscreteLp` or `FnBase`
+            Domain of the functional.
+        epsilon : float
+            The parameter of the Huber norm functional.
+
+        Examples
+        --------
+        Example of initializing the Huber norm functional
+
+        >>> space = odl.uniform_discr(0, 1, 14)
+        >>> epsilon = 0.1
+        >>> huber_norm = odl.solvers.HuberNorm(space, epsilon)
+
+        Check that if all elements are > epsilon we get the L1-norm modified
+        epsilon/2 * the mass of the space
+        >>> l1_norm = odl.solvers.L1Norm(space)
+        >>> element = space.one()
+        >>> constant = epsilon/2 * element.inner(element)
+        >>> (huber_norm(element) - (l1_norm(element) - constant)) < 1e-5
+        True
+
+        Check that if all elements are < epsilon we get the L2-norm modified
+        with weight 1/(2*epsilon)
+        >>> l2_norm = odl.solvers.L2Norm(space)
+        >>> element = (epsilon/2) * space.one()
+        >>> (huber_norm(element) - 1/(2*epsilon) * l2_norm(element)) < 1e-5
+        True
+        """
+        self.__epsilon = float(epsilon)
+        super(HuberNorm, self).__init__(space=space, linear=False,
+                                        grad_lipschitz=2)
+
+    @property
+    def epsilon(self):
+        """The smoothing parameter of the Huber norm functional."""
+        return self.__epsilon
+
+    def _call(self, x):
+        """Return ``self(x)``."""
+        indices = x.ufuncs.absolute().asarray() < self.epsilon
+        indices = np.float32(indices)
+
+        tmp = ((x * indices)**2 / (2.0 * self.epsilon) +
+               (x.ufuncs.absolute() - self.epsilon / 2.0) * (1-indices))
+
+        return tmp.inner(self.domain.one())
+
+    @property
+    def gradient(self):
+        """Gradient operator of the functional."""
+        functional = self
+
+        class HuberNormGradient(Operator):
+
+            """The gradient operator of this functional."""
+
+            def __init__(self):
+                """Initialize a new instance."""
+                super(HuberNormGradient, self).__init__(functional.domain,
+                                                        functional.domain,
+                                                        linear=False)
+
+            # TODO: Update this call. Might not work for PorductSpaces
+            def _call(self, x):
+                """Apply the gradient operator to the given point."""
+                indices = x.ufuncs.absolute().asarray() < functional.epsilon
+                indices = np.float32(indices)
+
+                to_return = ((x * indices) / (functional.epsilon) +
+                             (x).ufuncs.sign() * (1-indices))
+
+                return to_return
+
+        return HuberNormGradient()
+
+    @property
+    def proximal(self):
+        """Return the ``proximal factory`` of the functional.
+
+        See Also
+        --------
+        odl.solvers.proximal_huber_norm : `proximal factory` for the Huber
+            norm.
+        """
+        return proximal_huber_norm(space=self.domain, epsilon=self.epsilon)
+
+    def __repr__(self):
+        """Return ``repr(self)``."""
+        return '{}({!r}, {})'.format(self.__class__.__name__,
+                                     self.domain,
+                                     self.epsilon)
+
+
 if __name__ == '__main__':
     # pylint: disable=wrong-import-position
     from odl.util.testutils import run_doctests
diff --git a/odl/solvers/nonsmooth/proximal_operators.py b/odl/solvers/nonsmooth/proximal_operators.py
index 0e1c6258950..1e58b22d8ba 100644
--- a/odl/solvers/nonsmooth/proximal_operators.py
+++ b/odl/solvers/nonsmooth/proximal_operators.py
@@ -40,7 +40,8 @@
            'proximal_l1', 'proximal_convex_conj_l1',
            'proximal_l2', 'proximal_convex_conj_l2',
            'proximal_l2_squared', 'proximal_convex_conj_l2_squared',
-           'proximal_convex_conj_kl', 'proximal_convex_conj_kl_cross_entropy')
+           'proximal_convex_conj_kl', 'proximal_convex_conj_kl_cross_entropy',
+           'proximal_huber_norm')
 
 
 def combine_proximals(*factory_list):
@@ -1406,6 +1407,67 @@ def _call(self, x, out):
     return ProximalConvexConjKLCrossEntropy
 
 
+def proximal_huber_norm(space, epsilon):
+    """Proximal factory of the Huber norm.
+
+    Parameters
+    ----------
+    space : `FnBase`
+        Space X which is the domain of the functional F
+    epsilon : float
+        The parameter of the Huber norm functional.
+
+    Returns
+    -------
+    prox_factory : function
+        Factory for the proximal operator to be initialized.
+
+    See Also
+    --------
+    odl.solvers.HuberNorm : the Huber norm functional
+
+    Notes
+    -----
+    The proximal operator is given by given by the proximal operator of
+    1/(2 * epsilon) * L2 norm in points that are <= epsilon, and by the
+    proximal operator of the l1 norm in points that are > epsilon.
+    """
+
+    epsilon = float(epsilon)
+
+    class ProximalHuberNorm(Operator):
+
+        """Proximal operator of conjugate of cross entropy KL divergence."""
+
+        def __init__(self, sigma):
+            """Initialize a new instance.
+
+            Parameters
+            ----------
+            sigma : positive float
+            """
+            self.sigma = float(sigma)
+            super(ProximalHuberNorm, self).__init__(
+                domain=space, range=space, linear=False)
+
+        def _call(self, x, out):
+            """Apply the operator to ``x`` and stores the result in ``out``."""
+            l2_indices = x.ufuncs.absolute().asarray() < epsilon + self.sigma
+            l2_indices = np.float32(l2_indices)
+
+            pos_indices = np.float32(x.asarray() > 0)
+
+            tmp = (epsilon / (epsilon + self.sigma) * x * l2_indices +
+                   (x - self.sigma) * (1 - l2_indices) * pos_indices +
+                   (x + self.sigma) * (1 - l2_indices) * (1 - pos_indices))
+
+            out.assign(tmp)
+
+            return out
+
+    return ProximalHuberNorm
+
+
 if __name__ == '__main__':
     # pylint: disable=wrong-import-position
     from odl.util.testutils import run_doctests
diff --git a/odl/test/largescale/solvers/nonsmooth/default_functionals_slow_test.py b/odl/test/largescale/solvers/nonsmooth/default_functionals_slow_test.py
index c3fea29152e..fbd8c5d80cd 100644
--- a/odl/test/largescale/solvers/nonsmooth/default_functionals_slow_test.py
+++ b/odl/test/largescale/solvers/nonsmooth/default_functionals_slow_test.py
@@ -32,7 +32,7 @@
 func_params = ['l1', 'l2', 'l2^2', 'kl', 'kl_cross_ent', 'const',
                'groupl1-1', 'groupl1-2',
                'nuclearnorm-1-1', 'nuclearnorm-1-2', 'nuclearnorm-1-inf',
-               'quadratic', 'linear']
+               'quadratic', 'linear', 'huber']
 
 func_ids = [" functional='{}' ".format(p) for p in func_params]
 
@@ -73,6 +73,8 @@ def functional(request, linear_offset, quadratic_offset, dual):
                                          vector=space.one(), constant=0.623)
     elif name == 'linear':
         func = odl.solvers.QuadraticForm(vector=space.one(), constant=0.623)
+    elif name == 'huber':
+        func = odl.solvers.HuberNorm(space, epsilon=0.162)
     else:
         assert False
 
@@ -126,6 +128,14 @@ def test_proximal_defintion(functional, stepsize):
         pytest.skip('functional has no call method')
         return
 
+    # Special case not covered above. Cconj not implemetned is not seen through
+    # the quadratic perturbation
+    if (isinstance(functional, odl.solvers.FunctionalQuadraticPerturb) and
+            isinstance(functional.functional,
+                       FunctionalDefaultConvexConjugate)):
+        pytest.skip('functional has no call method')
+        return
+
     # No implementation of the proximal for convex conj of
     # FunctionalQuadraticPerturb unless the quadratic term is 0.
     if (isinstance(functional, odl.solvers.FunctionalQuadraticPerturb) and
@@ -189,6 +199,23 @@ def test_convex_conj_defintion(functional):
         pytest.skip('functional has no call')
         return
 
+    # Special case not covered above. Cconj not implemetned is not seen through
+    # the quadratic perturbation
+    if (isinstance(functional, odl.solvers.FunctionalQuadraticPerturb) and
+            isinstance(functional.functional,
+                       FunctionalDefaultConvexConjugate)):
+        pytest.skip('functional has no call')
+        return
+
+    # Same but with translation on the functional
+    if (isinstance(functional, odl.solvers.FunctionalTranslation) and
+            (isinstance(functional.functional,
+                       FunctionalDefaultConvexConjugate) or
+             isinstance(functional.functional.convex_conj,
+                       FunctionalDefaultConvexConjugate))):
+        pytest.skip('functional has no call')
+        return
+
     f_convex_conj = functional.convex_conj
     if isinstance(f_convex_conj, FunctionalDefaultConvexConjugate):
         pytest.skip('functional has no convex conjugate')
diff --git a/odl/test/solvers/functional/functional_test.py b/odl/test/solvers/functional/functional_test.py
index 9512589d09c..ba505dccac6 100644
--- a/odl/test/solvers/functional/functional_test.py
+++ b/odl/test/solvers/functional/functional_test.py
@@ -49,7 +49,7 @@ def space(request, fn_impl):
 func_params = ['l1 ', 'l2', 'l2^2', 'constant', 'zero', 'ind_unit_ball_1',
                'ind_unit_ball_2', 'ind_unit_ball_pi', 'ind_unit_ball_inf',
                'product', 'quotient', 'kl', 'kl_cc', 'kl_cross_ent',
-               'kl_cc_cross_ent']
+               'kl_cc_cross_ent', 'huber']
 func_ids = [" functional='{}' ".format(p) for p in func_params]
 
 
@@ -91,6 +91,8 @@ def functional(request, space):
         func = odl.solvers.functional.KullbackLeiblerCrossEntropy(space)
     elif name == 'kl_cc_cross_ent':
         func = odl.solvers.KullbackLeiblerCrossEntropy(space).convex_conj
+    elif name == 'huber':
+        func = odl.solvers.HuberNorm(space, epsilon=0.1)
     else:
         assert False