diff --git a/misc/phase_retrieval.ipynb b/misc/phase_retrieval.ipynb
new file mode 100644
index 0000000..651fe07
--- /dev/null
+++ b/misc/phase_retrieval.ipynb
@@ -0,0 +1,1097 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### PaganinProcessor examples\n",
+    "The code in this notebook contains some examples and comparisons of using the `PaganinProcessor` phase retrieval methods in CIL "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This notebook requires CIL v24.1.0 or greater, check the version below"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import cil\n",
+    "print(cil.__version__)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Load some dependencies from CIL"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from cil.utilities import dataexample\n",
+    "from cil.processors import PaganinProcessor, Slicer, Binner, TransmissionAbsorptionConverter, Padder, RingRemover\n",
+    "from cil.utilities.display import show2D\n",
+    "from cil.recon import FDK, FBP\n",
+    "from cil.utilities.jupyter import islicer\n",
+    "from cil.io.utilities import HDF5_utilities\n",
+    "from cil.framework import AcquisitionGeometry, AcquisitionData"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This notebook also requires TomoPy, numpy and matplotlib"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "from tomopy.prep.phase import retrieve_phase\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "##### Parallel beam data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "First we test the PaganinProcessor with parallel beam data. In the following cell\n",
+    "- Get a test parallel beam dataset \n",
+    "- Perform a filtered back projection reconstruction using `FBP()`"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data = dataexample.SIMULATED_PARALLEL_BEAM_DATA.get()\n",
+    "data.reorder(order='tigre')\n",
+    "data.geometry.config.units = 'um'\n",
+    "data_abs = -1*data.log()\n",
+    "ig = data.geometry.get_ImageGeometry()\n",
+    "fbp =  FBP(data_abs, ig)\n",
+    "recon = fbp.run(verbose=0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next we repeat the above steps with the PaganinProcessor"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Choose the parameters to be used in the phase retrieval\n",
+    "- Delta and beta are the real and complex part of the material refractive index. These can be found for x-ray wavelengths at  https://henke.lbl.gov/optical_constants/getdb2.html. Here we just use some exagerated values to demonstrate the effect.\n",
+    "- The experiment peak energy in default units eV"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "delta = 1\n",
+    "beta = 0.002\n",
+    "energy = 40000"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Run the phase retrieval using the `PaganinProcessor` which is implemented based on [Paganin 2002](https://onlinelibrary.wiley.com/doi/10.1046/j.1365-2818.2002.01010.x). The processor returns the material retrieved thickness, removing the effect of phase in the image\n",
+    "\n",
+    "$$\n",
+    "T(x,y) = - \\frac{1}{\\mu}\\ln\\left (\\mathcal{F}^{-1}\\left \n",
+    "        (\\frac{\\mathcal{F}\\left ( M^2I_{norm}(x, y,z = \\Delta) \\right )}{1 + \n",
+    "          \\alpha\\left ( k_x^2 + k_y^2 \\right )}  \\right )\\right )\n",
+    "$$\n",
+    "\n",
+    "where\n",
+    "- $T$, is the sample thickness,\n",
+    "- $\\mu = \\frac{4\\pi\\beta}{\\lambda}$ is the material linear \n",
+    "attenuation coefficient where $\\beta$ is the complex part of the \n",
+    "material refractive index and $\\lambda=\\frac{hc}{E}$ is the probe \n",
+    "wavelength,\n",
+    "- $M$ is the magnification at the detector,\n",
+    "- $I_{norm}$ is the input image which is expected to be the \n",
+    "normalised transmission data, \n",
+    "- $\\Delta$ is the propagation distance,\n",
+    "- $\\alpha = \\frac{\\Delta\\delta}{\\mu}$ is a parameter determining \n",
+    "the strength of the filter to be applied in Fourier space where \n",
+    "$\\delta` is the real part of the deviation of the material \n",
+    "refractive index from 1 \n",
+    "- $k_x, k_y = \\left ( \\frac{2\\pi p}{N_xW}, \\frac{2\\pi q}{N_yW} \n",
+    "\\right )$ where $p$ and $q$ are co-ordinates in a Fourier \n",
+    "mesh in the range $-N_x/2$ to $N_x/2$ for an image with \n",
+    "size $N_x, N_y$ and pixel size $W$.\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In the following cell:\n",
+    "- Run the `PaganinProcessor` to retrieve $T(x,y)$ from the test dataset\n",
+    "- Reconstruct using FBP"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "processor = PaganinProcessor(delta=delta, beta=beta, energy=energy)\n",
+    "processor.set_input(data)\n",
+    "thickness = processor.get_output(override_geometry={'propagation_distance':10})\n",
+    "fbp =  FBP(thickness, ig)\n",
+    "recon_thickness = fbp.run(verbose=0)\n",
+    "\n",
+    "# calculate mu to get recon_attenuation with the same scaling as the original image\n",
+    "attenuation = thickness*processor.mu\n",
+    "fbp =  FBP(attenuation, ig)\n",
+    "recon_attenuation = fbp.run(verbose=0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next we test the phase retrieval using `PaganinProcessor(full_retrieval=False)`. In this implementation, the same filter is applied in Fourier space but the $-log()$ is not applied. \n",
+    "$$\n",
+    "I_{filt} = \\mathcal{F}^{-1}\\left (\\frac{\\mathcal{F}\\left ( \n",
+    "        I(x, y,z = \\Delta) \\right )}\n",
+    "        {1 - \\alpha\\left ( k_x^2 + k_y^2 \\right )}  \\right )\n",
+    "$$\n",
+    "This gives flexibility to apply a Paganin-like filter but doesn't require data that has already been converted from trans.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Run PaganinProcessor as a filter using full_retrieval=False on the absorption data and reconstruct\n",
+    "processor = PaganinProcessor(delta=delta, beta=beta, energy=energy, full_retrieval=False)\n",
+    "processor.set_input(data_abs)\n",
+    "filtered_image = processor.get_output(override_geometry={'propagation_distance':10})\n",
+    "fbp =  FBP(filtered_image, ig)\n",
+    "recon_filter = fbp.run(verbose=0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For comparison run Tomopy phase retreival with raw data, then convert to absorption and reconstruct"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "tomopy_alpha = (1/(delta/beta))/(4*np.pi**2)\n",
+    "data_tomopy = data.copy()\n",
+    "data_tmp = retrieve_phase(data.array, pixel_size=processor.pixel_size, dist=processor.propagation_distance, energy=energy/1000, alpha=tomopy_alpha)\n",
+    "data_tomopy.fill(data_tmp)\n",
+    "data_tomopy = -1*data_tomopy.log()\n",
+    "ig = data_tomopy.geometry.get_ImageGeometry()\n",
+    "fbp =  FBP(data_tomopy, ig)\n",
+    "recon_tomopy = fbp.run(verbose=0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Also run Tomopy phase retreival with absorption data and reconstruct"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "tomopy_alpha = (1/(delta/beta))/(4*np.pi**2)\n",
+    "data_tomopy_abs = data_abs.copy()\n",
+    "data_tmp = retrieve_phase(data_abs.array, pixel_size=processor.pixel_size, dist=processor.propagation_distance, energy=energy/1000, alpha=tomopy_alpha)\n",
+    "data_tomopy_abs.fill(data_tmp)\n",
+    "ig = data_tomopy_abs.geometry.get_ImageGeometry()\n",
+    "fbp =  FBP(data_tomopy_abs, ig)\n",
+    "recon_tomopy_abs = fbp.run(verbose=0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Compare the reconstructions"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vertical_slice = 67\n",
+    "x_range = slice(50,90)\n",
+    "y_range = slice(50,90)\n",
+    "\n",
+    "show2D([recon.array[vertical_slice,x_range,y_range], recon_thickness.array[vertical_slice,x_range,y_range], recon_attenuation.array[vertical_slice,x_range,y_range], recon_filter.array[vertical_slice,x_range,y_range], recon_tomopy.array[vertical_slice,x_range,y_range], recon_tomopy_abs.array[vertical_slice,x_range,y_range]],\n",
+    "title=['Original image', 'Phase retrieval - thickness', 'Phase retrieval - scaled by mu', 'Phase retrieval - full_retrieval=False', 'Tomopy phase retrieval transmission', 'Tomopy phase retrieval absorption'], \n",
+    "axis_labels = ('horizontal_y', 'horizontal_x'), num_cols=3)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Compare a cross-section through the reconstruction"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axs = plt.subplots(1,2, figsize=(12,5))\n",
+    "ax = axs[0]\n",
+    "vertical_slice = 67\n",
+    "y_slice = 70\n",
+    "x_range = range(50,90)\n",
+    "ax.plot(x_range, recon.array[vertical_slice,y_slice,x_range])\n",
+    "ax.plot(x_range, recon_attenuation.array[vertical_slice,y_slice,x_range])\n",
+    "ax.plot(x_range, recon_filter.array[vertical_slice,y_slice,x_range])\n",
+    "ax.plot(x_range, recon_tomopy.array[vertical_slice,y_slice,x_range],'--')\n",
+    "ax.plot(x_range, recon_tomopy_abs.array[vertical_slice,y_slice,x_range],'--')\n",
+    "\n",
+    "ax.set_xlabel('Horizontal x')\n",
+    "ax.set_ylabel('Intensity')\n",
+    "ax.set_title('Line Profile at horizontal_y=' + str(y_slice) + ', vertical slice=' + str(vertical_slice))\n",
+    "\n",
+    "ax = axs[1]\n",
+    "x_slice = 70\n",
+    "y_range = range(50,90)\n",
+    "ax.plot(y_range, recon.array[vertical_slice,y_range,x_slice])\n",
+    "ax.plot(y_range, recon_attenuation.array[vertical_slice,y_range,x_slice])\n",
+    "ax.plot(y_range, recon_filter.array[vertical_slice,y_range,x_slice])\n",
+    "ax.plot(y_range, recon_tomopy.array[vertical_slice,y_range,x_slice],'--')\n",
+    "ax.plot(y_range, recon_tomopy_abs.array[vertical_slice,y_range,x_slice],'--')\n",
+    "\n",
+    "ax.set_xlabel('Horizontal y')\n",
+    "ax.set_ylabel('Intensity')\n",
+    "ax.set_title('Line Profile at horizontal_x=' + str(x_slice) + ', vertical slice=' + str(vertical_slice))\n",
+    "ax.legend(['Original', 'Phase retrieved - scaled by mu', 'Phase retrieved - full_retrieval=False', 'TomoPy on transmission data', 'TomoPy on absorption data'])\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We see that all methods blur the result in comparison to the original reconstruction. The scaled phase retrieval in CIL matches the Tomopy method performed on transmission data and the filter in CIL matches the Tomopy method performed on absorption data\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can approximate the signal to noise of each reconstruction as the mean divided by the standard deviation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(\"Original reconstruction SNR = \" + str(recon.mean()/recon.array.std()))\n",
+    "print(\"Phase retrieved reconstruction SNR = \" + str(recon_attenuation.mean()/recon_attenuation.array.std()))\n",
+    "print(\"Phase retrieved (full_retrieval=False) reconstruction SNR = \" + str(recon_filter.mean()/recon_filter.array.std()))\n",
+    "print(\"TomoPy on transmission data reconstruction SNR = \" + str(recon_tomopy.mean()/recon_tomopy.array.std()))\n",
+    "print(\"TomoPy on absorption data reconstruction SNR = \" + str(recon_tomopy_abs.mean()/recon_tomopy_abs.array.std()))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In all cases, the phase retrieval improves the SNR"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "\n",
+    "##### Cone beam data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "delta = 1\n",
+    "beta = 0.0001\n",
+    "energy = 40000"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "With cone beam data, the magnification $M$ has an effect on the phase retrieval\\\n",
+    "$ T = -\\frac{1}{\\mu}\\ln(F^{-1}\\frac{F(M^2 I_{norm}(x,y,z=\\Delta))}{1+\\frac{\\Delta\\lambda\\delta}{4\\pi\\beta}(k_x^2+k_y^2)/M})$\\\n",
+    "The $M^2$ on top means sometimes we get a number larger than 1 inside the $\\ln$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Get some cone beam data and perform reconstruction without phase retrieval"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data = dataexample.SIMULATED_CONE_BEAM_DATA.get()\n",
+    "data.geometry.config.units = 'um'\n",
+    "print('Magnification = ' + str(data.geometry.magnification))\n",
+    "data.reorder(order='tigre')\n",
+    "data_abs = -1*data.log()\n",
+    "ig = data.geometry.get_ImageGeometry()\n",
+    "fdk =  FDK(data_abs, ig)\n",
+    "recon = fdk.run(verbose=0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Run phase retrieval on raw data and reconstruct"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "processor = PaganinProcessor(delta=delta, beta=beta, energy=energy)\n",
+    "processor.set_input(data)\n",
+    "thickness = processor.get_output()\n",
+    "recon_thickness = fdk.run(verbose=0)\n",
+    "# calculate mu to get recon_attenuation with the same scaling as the original image\n",
+    "attenuation = thickness*processor.mu\n",
+    "fdk =  FDK(attenuation, ig)\n",
+    "recon_attenuation = fdk.run(verbose=0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Run PaganinProcessor as a filter using `get_output(full_retrieval=False)` on the absorption data and reconstruct"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "processor = PaganinProcessor(delta=delta, beta=beta, energy=energy, full_retrieval=False)\n",
+    "processor.set_input(data_abs)\n",
+    "filtered_image = processor.get_output()\n",
+    "\n",
+    "fdk =  FDK(filtered_image, ig)\n",
+    "recon_filter = fdk.run(verbose=0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "TomoPy can only be used with parallel beam data so we do not use it for comparison here"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Compare the reconstructions"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vertical_slice = 67\n",
+    "x_range = slice(50,90)\n",
+    "y_range = slice(50,90)\n",
+    "show2D([recon.array[vertical_slice,x_range,y_range], recon_thickness.array[vertical_slice,x_range,y_range], recon_attenuation.array[vertical_slice,x_range,y_range], recon_filter.array[vertical_slice,x_range,y_range]],\n",
+    "title=['Original image', 'Phase retrieval - thickness', 'Phase retrieval - scaled by mu', 'Phase retrieval - full_retrieval=False'],\n",
+    "axis_labels = ('horizontal_y', 'horizontal_x'))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Compare the cross-sections through the reconstructions"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axs = plt.subplots(1,2, figsize=(12,5))\n",
+    "ax = axs[0]\n",
+    "vertical_slice = 67\n",
+    "y_slice = 70\n",
+    "x_range = range(50,90)\n",
+    "ax.plot(recon.array[vertical_slice, y_slice, x_range])\n",
+    "ax.plot(recon_attenuation.array[vertical_slice, y_slice, x_range])\n",
+    "ax.plot(recon_filter.array[vertical_slice, y_slice, x_range])\n",
+    "\n",
+    "ax.set_xlabel('Horizontal x')\n",
+    "ax.set_ylabel('Intensity')\n",
+    "ax.set_title('Line Profile at horizontal_y=' + str(y_slice) + ', vertical slice=' + str(vertical_slice))\n",
+    "\n",
+    "ax = axs[1]\n",
+    "x_slice = 70\n",
+    "y_range = range(50,90)\n",
+    "ax.plot(recon.array[vertical_slice,y_range,x_slice])\n",
+    "ax.plot(recon_attenuation.array[vertical_slice,y_range,x_slice])\n",
+    "ax.plot(recon_filter.array[vertical_slice,y_range,x_slice])\n",
+    "\n",
+    "ax.set_xlabel('Horizontal y')\n",
+    "ax.set_ylabel('Intensity')\n",
+    "ax.set_title('Line Profile at horizontal_x=' + str(x_slice) + ', vertical slice=' + str(vertical_slice))\n",
+    "ax.legend(['Original', 'Phase retrieval - scaled by mu', 'Phase retrieval - full_retrieval=False'])\n",
+    "\n",
+    "plt.tight_layout()\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can see that both methods blur the original image. The phase retrieval method becomes negative because of the values > 1 are passed into the negative log. This may be an indication that the full phase retrieval is not valid for this experimental setup, in which case using `full_retrieval = False` may be more useful as it just applies a Paganin-like filter to the data "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Compare the signal to noise ratio"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(\"Original reconstruction SNR = \" + str(recon.mean()/recon.array.std()))\n",
+    "print(\"Phase retrieved reconstruction SNR = \" + str(recon_attenuation.mean()/recon_attenuation.array.std()))\n",
+    "print(\"Phase retrieved (full_retrieval=False) reconstruction SNR = \" + str(recon_filter.mean()/recon_filter.array.std()))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The negative values in the phase retrival skew the result but with full_retrieval=False the SNR is improved"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "##### Generalised Paganin method"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The generalised Paganin method is implemented in CIL following the description in https://iopscience.iop.org/article/10.1088/2040-8986/abbab9 \\\n",
+    "When features in the image are close to the Nyquist frequency of the system, a more generalised form of the Pagnin filter can be used which preserves these high frequency features while still boosting SNR. This may have a similar effect to applying an unsharp mask after the normal Paganin phase retrieval. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$$\n",
+    "T(x,y) =  -\\frac{1}{\\mu}\\ln\\left (\\mathcal{F}^{-1}\\left (\\frac{\n",
+    "        \\mathcal{F}\\left ( M^2I_{norm}(x, y,z = \\Delta) \\right )}{1 - \\frac{2\n",
+    "        \\alpha}{W^2}\\left ( \\cos(Wk_x) + \\cos(Wk_y) -2 \\right )}  \\right )\n",
+    "        \\right )\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Choose delta and beta"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "delta = 1\n",
+    "beta = 0.001"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Get the simulated parallel data and perform the reconstruction without phase retrieval"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data = dataexample.SIMULATED_PARALLEL_BEAM_DATA.get()\n",
+    "data.geometry.config.units = 'um'\n",
+    "data.reorder(order='tigre')\n",
+    "data_abs = -1*data.log()\n",
+    "ig = data.geometry.get_ImageGeometry()\n",
+    "fbp =  FBP(data_abs, ig)\n",
+    "recon = fbp.run(verbose=0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Run phase retrival using the original Paganin method and reconstruct"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "processor = PaganinProcessor(delta=delta, beta=beta, energy=energy)\n",
+    "processor.set_input(data)\n",
+    "thickness = processor.get_output(override_geometry={'propagation_distance':10})\n",
+    "fbp =  FBP(thickness, ig)\n",
+    "recon_thickness = fbp.run(verbose=0)\n",
+    "\n",
+    "# calculate mu to get recon_attenuation with the same scaling as the original image\n",
+    "attenuation = thickness*processor.mu\n",
+    "fbp =  FBP(attenuation, ig)\n",
+    "recon_attenuation = fbp.run(verbose=0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Run phase retrieval on the data using the generalised Paganin method and reconstruct"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "processor = PaganinProcessor(delta=delta, beta=beta, energy=energy, filter_type='generalised_paganin_method')\n",
+    "processor.set_input(data)\n",
+    "thickness_GPM = processor.get_output(override_geometry={'propagation_distance':10})\n",
+    "fbp =  FBP(thickness_GPM, ig)\n",
+    "recon_thickness_GPM = fbp.run(verbose=0)\n",
+    "\n",
+    "# calculate mu to get recon_attenuation with the same scaling as the original image\n",
+    "attenuation_GPM = thickness_GPM*processor.mu\n",
+    "fbp =  FBP(attenuation_GPM, ig)\n",
+    "recon_attenuation_GPM = fbp.run(verbose=0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Plot cross-sections through the results"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axs = plt.subplots(1,2, figsize=(12,5))\n",
+    "ax = axs[0]\n",
+    "vertical_slice = 67\n",
+    "y_slice = 70\n",
+    "x_range = range(50,90)\n",
+    "ax.plot(recon.array[vertical_slice, y_slice, x_range])\n",
+    "ax.plot(recon_attenuation.array[vertical_slice, y_slice, x_range])\n",
+    "ax.plot(recon_attenuation_GPM.array[vertical_slice, y_slice, x_range])\n",
+    "\n",
+    "ax.set_xlabel('Horizontal x')\n",
+    "ax.set_ylabel('Intensity')\n",
+    "ax.set_title('Line Profile at horizontal_x=' + str(x_slice) + ', vertical slice=' + str(vertical_slice))\n",
+    "ax.legend(['Original', 'Phase retrieval', 'Phase retrieval - GPM'])\n",
+    "\n",
+    "ax = axs[1]\n",
+    "x_slice = 70\n",
+    "y_range = range(50,90)\n",
+    "ax.plot(recon.array[vertical_slice,y_range,x_slice])\n",
+    "ax.plot(recon_attenuation.array[vertical_slice,y_range,x_slice])\n",
+    "ax.plot(recon_attenuation_GPM.array[vertical_slice,y_range,x_slice])\n",
+    "\n",
+    "ax.set_xlabel('Horizontal y')\n",
+    "ax.set_ylabel('Intensity')\n",
+    "ax.set_title('Line Profile at horizontal_x=' + str(x_slice) + ', vertical slice=' + str(vertical_slice))\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Check the SNR"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(\"Original reconstruction SNR = \" + str(recon.mean()/recon.array.std()))\n",
+    "print(\"Phase retrieved reconstruction SNR = \" + str(recon_attenuation.mean()/recon_attenuation.array.std()))\n",
+    "print(\"Phase retrieved with GPM reconstruction SNR = \" + str(recon_attenuation_GPM.mean()/recon_attenuation_GPM.array.std()))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The GPM has slightly improved resolution of the sample features while maintaining the SNR boost"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### TomoBank example"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This example uses dataset tomo_00068 from TomoBank https://tomobank.readthedocs.io/en/latest/source/data/docs.data.phasecontrast.html#wet-sample which can be retrieved using:\n",
+    "`wget https://g-a0400.fd635.8443.data.globus.org/tomo_00068/tomo_00068.h5`\n",
+    "\n",
+    "The data were collected at Syrmep beamline of the Elettra synchotron. A description of the experiment is given in https://link.springer.com/chapter/10.1007/978-3-319-19387-8_70"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Load the file, you may need to change the filename to the path where you downloaded it"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "filename = 'tomo_00068.h5' \n",
+    "data = HDF5_utilities.read(filename=filename, dset_path='/exchange/data')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Construct a CIL AcquisitionData object using the parameters in https://tomobank.readthedocs.io/en/latest/source/data/docs.data.phasecontrast.html#wet-sample"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pixel_size = 0.0041 #mm\n",
+    "propagation_distance = 150 #mm\n",
+    "angles = HDF5_utilities.read(filename=filename, dset_path='/exchange/theta')\n",
+    "ag = AcquisitionGeometry.create_Parallel3D(detector_position=[0, propagation_distance, 0], units='mm').set_panel([np.shape(data)[2],np.shape(data)[1]], pixel_size=pixel_size).set_angles(angles)\n",
+    "data = AcquisitionData(data, deep_copy=False, geometry = ag)\n",
+    "data.reorder(order='tigre')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Apply a centre of rotation correction"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data.geometry.set_centre_of_rotation(1463.5-(data.shape[2]/2), distance_units='pixels')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Use `islicer` to view all slices so we can choose a region to crop"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "islicer(data)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Crop the data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "processor = Slicer(roi={'horizontal':(600,2500,1)})\n",
+    "processor.set_input(data)\n",
+    "data = processor.get_output()\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Check the cropped data looks sensible"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "show2D(data)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next we bin a few angles"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print('Data shape before: ' + str(data.shape))\n",
+    "processor = Binner(roi={'angle':(None, None, 3)})\n",
+    "processor.set_input(data)\n",
+    "data = processor.get_output()\n",
+    "print('Data shape after: ' + str(data.shape))\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Run phase retrieval on the raw data\n",
+    "\n",
+    "Parameters from https://tomobank.readthedocs.io/en/latest/source/data/docs.data.phasecontrast.html#wet-sample"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "delta = 1\n",
+    "beta = 1e-1\n",
+    "energy = 14000\n",
+    "\n",
+    "processor = PaganinProcessor(delta=delta, beta=beta, energy=energy)\n",
+    "processor.set_input(data)\n",
+    "thickness = processor.get_output()\n",
+    "\n",
+    "# calculate mu to get recon_attenuation with the same scaling as the original image\n",
+    "data_phase = thickness*processor.mu"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Run phase retrieval using the generalised Paganin method on the raw data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "processor = PaganinProcessor(delta=delta, beta=beta, energy=energy, \n",
+    "                             filter_type='generalised_paganin_method')\n",
+    "processor.set_input(data)\n",
+    "thickness = processor.get_output()\n",
+    "\n",
+    "# calculate mu to get recon_attenuation with the same scaling as the original image\n",
+    "data_phase_generalised = thickness*processor.mu"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Get a slice of each data set"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vertical_slice = 900\n",
+    "data_phase = data_phase.get_slice(vertical=vertical_slice)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vertical_slice = 900\n",
+    "data_phase_generalised = data_phase_generalised.get_slice(vertical=vertical_slice)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For comparison just run TransmissionAbsorptionConverter on the same slice of the original data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data_slice = data.get_slice(vertical=vertical_slice)\n",
+    "\n",
+    "processor = TransmissionAbsorptionConverter()\n",
+    "processor.set_input(data_slice)\n",
+    "processor.get_output(out=data_slice)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Perform some extra processing steps on both datasets"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Pad the data\n",
+    "ig = data_slice.geometry.get_ImageGeometry()\n",
+    "padsize = 1000\n",
+    "data_slice = Padder.linear_ramp(pad_width={'horizontal': padsize}, end_values=0)(data_slice)\n",
+    "data_phase = Padder.linear_ramp(pad_width={'horizontal': padsize}, end_values=0)(data_phase)\n",
+    "data_phase_generalised = Padder.linear_ramp(pad_width={'horizontal': padsize}, end_values=0)(data_phase_generalised)\n",
+    "\n",
+    "# Ring remover\n",
+    "N_decompositions = 20\n",
+    "wavelet_filter_name = 'db20'\n",
+    "sigma = 5.5\n",
+    "\n",
+    "processor = RingRemover(N_decompositions, wavelet_filter_name, sigma, info=False) \n",
+    "processor.set_input(data_slice)\n",
+    "data_slice = processor.get_output()\n",
+    "\n",
+    "processor = RingRemover(N_decompositions, wavelet_filter_name, sigma, info=False) \n",
+    "processor.set_input(data_phase)\n",
+    "data_phase = processor.get_output()\n",
+    "\n",
+    "processor = RingRemover(N_decompositions, wavelet_filter_name, sigma, info=False) \n",
+    "processor.set_input(data_phase_generalised)\n",
+    "data_phase_generalised = processor.get_output()\n",
+    "\n",
+    "# Reconstruct\n",
+    "fbp =  FBP(data_slice, ig)\n",
+    "recon = fbp.run(verbose=0)\n",
+    "\n",
+    "fbp =  FBP(data_phase, ig)\n",
+    "recon_phase = fbp.run(verbose=0)\n",
+    "\n",
+    "fbp =  FBP(data_phase_generalised, ig)\n",
+    "recon_phase_g = fbp.run(verbose=0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Compare the reconstructions"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "show2D([recon, recon_phase, recon_phase_g], ['Original','Paganin Method', 'Generalised Paganin Method'], num_cols=3)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Plot a difference map between the Paganin Method and Generalised Paganin Method reconstructions"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "show2D(recon_phase.array[700:950,700:950]-recon_phase_g.array[700:950,700:950], title='Paganin Method - Generalised Paganin Method', cmap='seismic')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We see some of the high resolution details (e.g. the edges) are preserved in the GPM"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Compare SNR"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(\"Original reconstruction SNR = \" + str(recon.mean()/recon.array.std()))\n",
+    "print(\"Phase retrieved reconstruction SNR = \" + str(recon_phase.mean()/recon_phase.array.std()))\n",
+    "print(\"Phase retrieved with GPM reconstruction SNR = \" + str(recon_phase_g.mean()/recon_phase_g.array.std()))"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "cil_tests",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}