Updated example and added OMP options

examples/PGURE-SVT HyperSpy Demo.ipynb

@@ -4,12 +4,20 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# PGURE-SVT HyperSpy Demonstration\n",
+    "# PGURE-SVT Demonstration\n",
     "\n",
     "### Tom Furnival ([[email protected]](mailto:[email protected]))\n",
     "\n",
     "PGURE-SVT is an algorithm designed to denoise image sequences acquired in microscopy. It exploits the correlations between consecutive frames to form low-rank matrices, which are then recovered using a technique known as nuclear norm minimization. An unbiased risk estimator for mixed Poisson-Gaussian noise is used to automate the selection of the regularization parameter, while robust noise and motion estimation maintain broad applicability to many different types of microscopy.\n",
     "\n",
+    "You can read more about the algorithm and its applications in:\n",
+    "\n",
+    "> T Furnival, RK Leary, PA Midgley _\"Denoising time-resolved microscopy sequences with singular value thresholding.\"_ (Article in press). DOI:[10.1016/j.ultramic.2016.05.005](http://dx.doi.org/10.1016/j.ultramic.2016.05.005)\n",
+    "\n",
+    "The source code is available to download from [http://tjof2.github.io/pgure-svt/](http://tjof2.github.io/pgure-svt/).\n",
+    "\n",
+    "### Use with HyperSpy\n",
+    "\n",
     "This example notebook shows how PGURE-SVT can be combined with [HyperSpy](http://hyperspy.org), which is an open-source Python library that makes signal handling and processing straightforward in Python, with a friendly API. While you can use `pguresvt.pguresvt.SVT` to denoise a numpy array directly, `pguresvt.hspysvt.HSPYSVT` can directly denoise a HyperSpy signal."
    ]
   },
@@ -22,33 +30,15 @@
    "outputs": [],
    "source": [
     "# Configures the plotting backend\n",
-    "#%matplotlib inline\n",
-    "%matplotlib qt"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
-   "outputs": [],
-   "source": [
+    "%matplotlib qt\n",
+    "\n",
+    "# Import NumPy and the HyperSpy API\n",
     "import numpy as np\n",
-    "# Import the HyperSpy API\n",
-    "import hyperspy.api as hs"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
-   "outputs": [],
-   "source": [
+    "import hyperspy.api as hs\n",
+    "\n",
     "# Import the HyperSpy wrapper for PGURE-SVT\n",
-    "from pguresvt import hspysvt"
+    "from pguresvt import hspysvt\n",
+    "from pguresvt.pguresvt import PoissonGaussianNoiseGenerator"
    ]
   },
   {
@@ -80,7 +70,14 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Now we corrupt it with a mixture of Poisson and Gaussian noise."
+    "Now we corrupt the movie with using PGURE-SVT's built-in noise generator for mixed Poisson-Gaussian noise, according to the equation:\n",
+    "<img src=\"mixed_noise_figure.png\">\n",
+    "where the parameters are defined as:\n",
+    "```\n",
+    "alpha = detector gain\n",
+    "mu = detector offset\n",
+    "sigma = detector noise\n",
+    "```"
    ]
   },
   {
@@ -91,25 +88,19 @@
    },
    "outputs": [],
    "source": [
-    "# First extract the data and rescale to [0,1] range\n",
-    "clean = movie._data_aligned_with_axes\n",
-    "clean = clean / np.amax(clean)\n",
-    "\n",
     "# Detector gain\n",
-    "gain = 0.1\n",
+    "detector_gain = 0.1\n",
     "# Detector offset\n",
-    "offset = 0.1\n",
-    "# Detector variance\n",
-    "sigma = 0.1\n",
-    "\n",
-    "def addnoise(x):\n",
-    "    return gain * np.random.poisson(x / gain) + offset + sigma * np.random.randn()\n",
-    "addnoise = np.vectorize(addnoise, otypes=[np.float])\n",
-    "\n",
-    "noisy = addnoise(clean)\n",
-    "\n",
+    "detector_offset = 0.1\n",
+    "# Detector noise\n",
+    "detector_sigma = 0.1\n",
+    "\n",
+    "noisy = PoissonGaussianNoiseGenerator(movie._data_aligned_with_axes,\n",
+    "                                      alpha=detector_gain,\n",
+    "                                      mu=detector_offset,\n",
+    "                                      sigma=detector_sigma)\n",
     "noisy_movie = hs.signals.Image(noisy)\n",
-    "noisy_movie.plot()"
+    "noisy_movie.plot(navigator='slider')"
    ]
   },
   {
@@ -131,10 +122,13 @@
     "                arpssize=7, \n",
     "                tol=1e-7,\n",
     "                median=5,\n",
-    "                hotpixelthreshold=10)\n",
+    "                hotpixelthreshold=10,\n",
+    "                numthreads=1)\n",
     "```\n",
     "\n",
-    "In this example we do not use the noise estimation procedure, and instead provide the known parameters to the algorithm directly. This information is used by the PGURE optimizer to calculate the threshold."
+    "In this example we do not use the noise estimation procedure, and instead provide the known parameters to the algorithm directly. This information is used by the PGURE optimizer to calculate the threshold.\n",
+    "\n",
+    "*Note:* If you have a multicore machine, you can set `numthreads > 1` to speed up the calculation."
    ]
   },
   {
@@ -149,10 +143,11 @@
     "# Initialize with suggested parameters\n",
     "svt = hspysvt.HSPYSVT(patchsize=4,\n",
     "                      estimatenoise=False,\n",
-    "                      alpha=gain,\n",
-    "                      mu=offset,\n",
-    "                      sigma=sigma,\n",
-    "                      tol=1e-6)"
+    "                      alpha=detector_gain,\n",
+    "                      mu=detector_offset,\n",
+    "                      sigma=detector_sigma,\n",
+    "                      tol=1e-5,\n",
+    "                      numthreads=2)"
    ]
   },
   {
@@ -174,16 +169,18 @@
     "denoised_movie = svt.denoise(noisy_movie)\n",
     "\n",
     "# Plot denoised data\n",
-    "denoised_movie.plot()"
+    "denoised_movie.plot(navigator='slider')"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### 2. Time-resolved ADF-STEM image sequence\n",
+    "## 2. Time-resolved ADF-STEM image sequence\n",
     "\n",
-    "In this example we apply PGURE-SVT to an experimental dataset of a nanoparticle acquired using ADF-STEM. Here the noise levels are not known, so a noise estimation procedure is used."
+    "In this example we apply PGURE-SVT to an experimental dataset of a nanoparticle acquired using ADF-STEM. For larger images, you can use the `patchoverlap` parameter to control the trade-off between speed and accuracy of the denoising procedure.\n",
+    "\n",
+    "Here only the detector offset (`mu`) was known beforehand, so a noise estimation procedure is used for the other values."
    ]
   },
   {
@@ -194,19 +191,21 @@
    },
    "outputs": [],
    "source": [
-    "# Load example dataset\n",
+    "# Load example dataset and plot\n",
     "expt_movie = hs.load(\"experimentalsequence.tif\")\n",
+    "expt_movie.plot(navigator='slider')\n",
     "\n",
     "# Initialize with suggested parameters\n",
     "expt_svt = hspysvt.HSPYSVT(patchsize=4,\n",
-    "                           patchoverlap=2,                     \n",
-    "                           tol=1e-6)\n",
+    "                           patchoverlap=2,\n",
+    "                           mu=0.1,\n",
+    "                           numthreads=2)\n",
     "\n",
     "# Run the denoising\n",
     "denoised_movie = expt_svt.denoise(expt_movie)\n",
     "\n",
     "# Plot denoised data\n",
-    "denoised_movie.plot()"
+    "denoised_movie.plot(navigator='slider')"
    ]
   }
  ],