diff --git a/reference/reference_head_fbp.py b/reference/reference_head_fbp.py
index 4750a9f..4747de2 100644
--- a/reference/reference_head_fbp.py
+++ b/reference/reference_head_fbp.py
@@ -2,6 +2,7 @@
 
 import numpy as np
 import odl
+from odl.contrib import fom
 
 mu_water = 0.02
 photons_per_pixel = 10000.0
@@ -48,7 +49,7 @@
 # Reconstruct using FBP
 recon = pseudoinverse(-np.log(epsilon + data) / mu_water)
 
-print('psnr = {}'.format(odl.util.psnr(phantom, recon)))
+print('psnr = {}'.format(fom.psnr(phantom, recon)))
 
 # Display results
 phantom.show('phantom', clim=[0.8, 1.2])
diff --git a/reference/reference_head_tv.py b/reference/reference_head_tv.py
index 3f55498..67e8ef4 100644
--- a/reference/reference_head_tv.py
+++ b/reference/reference_head_tv.py
@@ -2,6 +2,7 @@
 
 import numpy as np
 import odl
+from odl.contrib import fom
 
 mu_water = 0.02
 photons_per_pixel = 10000.0
@@ -11,7 +12,7 @@
 # Create ODL data structures
 size = 512
 space = odl.uniform_discr([-128, -128], [128, 128], [size, size],
-                          dtype='float32', weighting='const')
+                          dtype='float32', weighting=1.0)
 
 # Make a fan beam geometry with flat detector
 # Angles: uniformly spaced, n = 1000, min = 0, max = 2 * pi
@@ -56,7 +57,7 @@
 op = odl.BroadcastOperator(nonlinear_operator, gradient)
 
 # Do not use the g functional, set it to zero.
-g = odl.solvers.ZeroFunctional(op.domain)
+f = odl.solvers.ZeroFunctional(op.domain)
 
 # Create functionals for the dual variable
 
@@ -67,7 +68,7 @@
 l1_norm = 0.0002 * odl.solvers.L1Norm(gradient.range)
 
 # Combine functionals, order must correspond to the operator K
-f = odl.solvers.SeparableSum(data_discr, l1_norm)
+g = odl.solvers.SeparableSum(data_discr, l1_norm)
 
 
 # --- Select solver parameters and solve using Chambolle-Pock --- #
@@ -84,14 +85,14 @@
 gamma = 0.01
 
 # Pass callback to the solver to display intermediate results
-callback = (odl.solvers.CallbackPrint(lambda x: odl.util.psnr(phantom, x)) &
+callback = (odl.solvers.CallbackPrint(lambda x: fom.psnr(x, phantom)) &
             odl.solvers.CallbackShow(clim=[0.8, 1.2]))
 
 odl.solvers.pdhg(
     x, f, g, op, tau=tau, sigma=sigma, niter=niter, gamma=gamma,
     callback=callback)
 
-print('psnr = {}'.format(odl.util.psnr(phantom, x)))
+print('psnr = {}'.format(fom.psnr(phantom, x)))
 
 # Display images
 x.show('head TV', clim=[0.8, 1.2])
diff --git a/reference/reference_shepp_fbp.py b/reference/reference_shepp_fbp.py
index 8cff868..5f2db17 100644
--- a/reference/reference_shepp_fbp.py
+++ b/reference/reference_shepp_fbp.py
@@ -2,6 +2,7 @@
 
 import numpy as np
 import odl
+from odl.contrib import fom
 
 
 # Create ODL data structures
@@ -31,7 +32,7 @@
 
 recon = pseudoinverse(data)
 
-print('psnr = {}'.format(odl.util.psnr(phantom, recon)))
+print('psnr = {}'.format(fom.psnr(phantom, recon)))
 
 # Display images
 data.show('Data')
diff --git a/reference/reference_shepp_huber.py b/reference/reference_shepp_huber.py
index 1d067a3..2f57a54 100644
--- a/reference/reference_shepp_huber.py
+++ b/reference/reference_shepp_huber.py
@@ -2,6 +2,7 @@
 
 import numpy as np
 import odl
+from odl.contrib import fom
 
 
 # Create ODL data structures
@@ -61,7 +62,7 @@
 
 
 # Optionally pass callback to the solver to display intermediate results
-callback = (odl.solvers.CallbackPrint(lambda x: odl.util.psnr(phantom, x)) &
+callback = (odl.solvers.CallbackPrint(lambda x: fom.psnr(phantom, x)) &
             odl.solvers.CallbackShow(clim=[0.1, 0.4]))
 
 # Choose a starting point
@@ -73,7 +74,7 @@
     hessinv_estimate=hessinv_estimate,
     callback=callback)
 
-print('psnr = {}'.format(odl.util.psnr(phantom, x)))
+print('psnr = {}'.format(fom.psnr(phantom, x)))
 
 # Display images
 x.show('Shepp-Logan TV windowed', clim=[0.1, 0.4])
diff --git a/reference/reference_shepp_tv.py b/reference/reference_shepp_tv.py
index b3fde25..dc46f3e 100644
--- a/reference/reference_shepp_tv.py
+++ b/reference/reference_shepp_tv.py
@@ -2,7 +2,7 @@
 
 import numpy as np
 import odl
-from odl.contrib.fom import psnr
+from odl.contrib import fom
 
 
 # Create ODL data structures
@@ -41,7 +41,7 @@
 op = odl.BroadcastOperator(operator, gradient)
 
 # Do not use the g functional, set it to zero.
-g = odl.solvers.ZeroFunctional(op.domain)
+f = odl.solvers.ZeroFunctional(op.domain)
 
 # Create functionals for the dual variable
 
@@ -52,22 +52,13 @@
 l1_norm = 0.3 * odl.solvers.L1Norm(gradient.range)
 
 # Combine functionals, order must correspond to the operator K
-f = odl.solvers.SeparableSum(l2_norm, l1_norm)
+g = odl.solvers.SeparableSum(l2_norm, l1_norm)
 
 
 # --- Select solver parameters and solve using Chambolle-Pock --- #
 
-
-# Estimated operator norm, add 10 percent to ensure ||K||_2^2 * sigma * tau < 1
-op_norm = 1.1 * odl.power_method_opnorm(op)
-
-niter = 1000  # Number of iterations
-tau = 0.1  # Step size for the primal variable
-sigma = 1.0 / (op_norm ** 2 * tau)  # Step size for the dual variable
-gamma = 0.1
-
 # Optionally pass callback to the solver to display intermediate results
-callback = (odl.solvers.CallbackPrint(lambda x: psnr(phantom, x)) &
+callback = (odl.solvers.CallbackPrint(lambda x: fom.psnr(x, phantom)) &
             odl.solvers.CallbackShow(clim=[0.1, 0.4]))
 
 # Choose a starting point
@@ -75,10 +66,10 @@
 
 # Run the algorithm
 odl.solvers.pdhg(
-    x, f, g, op, tau=tau, sigma=sigma, niter=niter, gamma=gamma,
-    callback=callback)
+    x, f, g, op, niter=1000, gamma=0.3,
+    callback=None)
 
-print('psnr = {}'.format(psnr(phantom, x)))
+print('psnr = {}'.format(fom.psnr(x, phantom)))
 
 # Display images
 x.show('Shepp-Logan TV windowed', clim=[0.1, 0.4])
diff --git a/reference/reference_shepp_tv_isotropic.py b/reference/reference_shepp_tv_isotropic.py
index 581a2fa..2030838 100644
--- a/reference/reference_shepp_tv_isotropic.py
+++ b/reference/reference_shepp_tv_isotropic.py
@@ -2,7 +2,7 @@
 
 import numpy as np
 import odl
-
+from odl.contrib import fom
 
 # Create ODL data structures
 size = 128
@@ -10,7 +10,7 @@
                           dtype='float32')
 
 # Creat parallel beam geometry
-geometry = odl.tomo.parallel_beam_geometry(space, angles=30)
+geometry = odl.tomo.parallel_beam_geometry(space, num_angles=30)
 
 # Create ray transform operator
 operator = odl.tomo.RayTransform(space, geometry)
@@ -41,7 +41,7 @@
 op = odl.BroadcastOperator(operator, gradient)
 
 # Do not use the g functional, set it to zero.
-g = odl.solvers.ZeroFunctional(op.domain)
+f = odl.solvers.ZeroFunctional(op.domain)
 
 # Create functionals for the dual variable
 
@@ -49,25 +49,17 @@
 l2_norm = odl.solvers.L2NormSquared(operator.range).translated(data)
 
 # Isotropic TV-regularization i.e. the l1-norm
-l1_norm = 0.3 * odl.solvers.L1Norm(gradient.range)
+l1_norm = 0.26 * odl.solvers.GroupL1Norm(gradient.range)
 
 # Combine functionals, order must correspond to the operator K
-f = odl.solvers.SeparableSum(l2_norm, l1_norm)
+g = odl.solvers.SeparableSum(l2_norm, l1_norm)
 
 
 # --- Select solver parameters and solve using Chambolle-Pock --- #
 
 
-# Estimated operator norm, add 10 percent to ensure ||K||_2^2 * sigma * tau < 1
-op_norm = 1.1 * odl.power_method_opnorm(op)
-
-niter = 1000  # Number of iterations
-tau = 0.1  # Step size for the primal variable
-sigma = 1.0 / (op_norm ** 2 * tau)  # Step size for the dual variable
-gamma = 0.1
-
 # Optionally pass callback to the solver to display intermediate results
-callback = (odl.solvers.CallbackPrint(lambda x: odl.util.psnr(phantom, x)) &
+callback = (odl.solvers.CallbackPrint(lambda x: fom.psnr(x, phantom)) &
             odl.solvers.CallbackShow(clim=[0.1, 0.4]))
 
 # Choose a starting point
@@ -75,10 +67,10 @@
 
 # Run the algorithm
 odl.solvers.pdhg(
-    x, f, g, op, tau=tau, sigma=sigma, niter=niter, gamma=gamma,
+    x, f, g, op, niter=1000, gamma=0.3,
     callback=None)
 
-print('psnr = {}'.format(odl.util.psnr(phantom, x)))
+print('psnr = {}'.format(fom.psnr(x, phantom)))
 
 # Display images
 x.show('Shepp-Logan TV windowed', clim=[0.1, 0.4])