deconvolution-subset-correctness-investigations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook looks at the impact of flux-normalising H, and also verifies that the deconvolution is stable if I fix certain planes in the reconstruction, subtract their effect from the input image, and then repeat the deconvolution for the remaining planes. That is something I intend to rely on to speed up my PIV analysis of Nils' huge dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import multiprocessing\n",
    "import tifffile\n",
    "import sys, time, os, csv\n",
    "import cProfile, pstats\n",
    "from jutils import tqdm_alias as tqdm\n",
    "\n",
    "import psfmatrix, lfimage\n",
    "import projector, lfdeconv\n",
    "import special_fftconvolve as special\n",
    "import jutils as util\n",
    "import py_light_field as plf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "projectorClass = projector.Projector_allC\n",
    "maxiter = 8\n",
    "inputImage = lfimage.LoadLightFieldTiff('/Users/jonny/Movies/Nils files/Rectified/Left/Cam_Left_40_X1_N19.tif')\n",
    "hMatrix = psfmatrix.LoadMatrix('PSFmatrix/reducedPSFMatrix_M22.2NA0.5MLPitch125fml3125from-156to156zspacing4Nnum19lambda520n1.33.mat')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stability investigations\n",
    "Now that I have flux-normalised H, does the deconvolution converge to a stable image (importantly: with a stable absolute amplitude) as we run increasing numbers of iterations? The maximum values do seem to rise a little, but the sum does at least seem to remain approximately the same. **Conclusion**: I am reassured that the sum and maximum remain almost unchanged.\n",
    "\n",
    "Can we do an initial reconstruction, subtract the contribution from most of the planes to the camera image, and then perform subsequent reconstructions on just a few planes of interest in the dataset? (This is what I want to do with my flow analysis, where I don't want to be having to reconstruct the whole volume during my iterative optimization, when I am only interested in the artefact-corrupted planes). **Conclusion**: yes, that does seem to work - I get almost exactly the same reconstruction as in the original, full reconstruction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hMatrixSmall = psfmatrix.LoadMatrix('PSFmatrix/fdnormPSFMatrix_M22.2NA0.5MLPitch125fml3125from-12to0zspacing4Nnum19lambda520n1.33.mat')\n",
    "hMatrixSmall1 = psfmatrix.LoadMatrix('PSFmatrix/fdnormPSFMatrix_M22.2NA0.5MLPitch125fml3125from-12to-8zspacing4Nnum19lambda520n1.33.mat')\n",
    "hMatrixSmall2 = psfmatrix.LoadMatrix('PSFmatrix/fdnormPSFMatrix_M22.2NA0.5MLPitch125fml3125from-4to0zspacing4Nnum19lambda520n1.33.mat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate a synthetic four-plane image\n",
    "obj = tifffile.imread('Data/03_Reconstructed/exampleData/definitive_worm_crop_X15_backproject.tif')\n",
    "obj = obj[7:11,np.newaxis].astype(np.float32)\n",
    "im = lfdeconv.ForwardProjectACC(hMatrixSmall, obj, progress=tqdm, projector=projectorClass())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Deconvolve the full four planes\n",
    "Htf = lfdeconv.BackwardProjectACC(hMatrixSmall, im, progress=tqdm, projector=projectorClass())\n",
    "obj8 = lfdeconv.DeconvRL(hMatrixSmall, Htf, 8, Htf.copy(), logPrint=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sanity check that intensities remain approx the same as we iterate the deconvolution.\n",
    "if True:\n",
    "    obj10 = lfdeconv.DeconvRL(hMatrixSmall, Htf, 2, obj8, logPrint=False)\n",
    "    obj20 = lfdeconv.DeconvRL(hMatrixSmall, Htf, 10, obj10, logPrint=False)\n",
    "    \n",
    "    print(np.max(obj8), np.sum(obj8))\n",
    "    print(np.max(obj10), np.sum(obj10))\n",
    "    print(np.max(obj20), np.sum(obj20))\n",
    "\n",
    "    plt.figure()\n",
    "    plt.imshow(obj8[0,0])\n",
    "    plt.show()\n",
    "    plt.figure()\n",
    "    plt.imshow(obj10[0,0])\n",
    "    plt.show()    \n",
    "    plt.figure()\n",
    "    plt.imshow(obj20[0,0])\n",
    "    plt.show()        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Forward-project just the first, and second, halves of that reconstructed object\n",
    "imFirstHalf = lfdeconv.ForwardProjectACC(hMatrixSmall1, obj8[0:2], progress=tqdm, projector=projectorClass())\n",
    "imSecondHalf = lfdeconv.ForwardProjectACC(hMatrixSmall2, obj8[2:4], progress=tqdm, projector=projectorClass())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Subtract the forward-projection of the first two planes from our original camera image,\n",
    "# and try reconstructing the second two planes in isolation.\n",
    "imPartial = im - imFirstHalf\n",
    "obj8SecondHalf = lfdeconv.DeconvRL(hMatrixSmall2, None, 8, None, im=imPartial, logPrint=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare the reconstructions\n",
    "plt.imshow(obj8[2,0]); plt.colorbar(); plt.show()\n",
    "plt.imshow(obj8SecondHalf[0,0]); plt.colorbar(); plt.show()\n",
    "plt.imshow(obj8[2,0] - obj8SecondHalf[0,0]); plt.colorbar(); plt.show()\n",
    "\n",
    "plt.imshow(obj8[3,0]); plt.colorbar(); plt.show()\n",
    "plt.imshow(obj8SecondHalf[1,0]); plt.colorbar(); plt.show()\n",
    "plt.imshow(obj8[3,0] - obj8SecondHalf[1,0]); plt.colorbar(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare the deconvolutions as obtained from various different PSFs\n",
    "It's clear that the max bug does increase the artefacts in the native focal plane, although it's not catastrophic. It's hard to quantify how much difference the normalisation makes, but it is not a huge visual impact, and the relative differences seem to be on the 10% level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theirs = tifffile.imread('/Users/jonny/Movies/Nils files/Single-view-reconstructions/Cam_Left_40_uncompressed.tif')\n",
    "mineNorm = tifffile.imread('Data/03_Reconstructed/Left/Cam_Left_40_X1_N19.tif')\n",
    "mineUnnorm = tifffile.imread('Data/03_Reconstructed/Left/Cam_Left_40_X1_N19_unnorm.tif')\n",
    "mineFDnorm = tifffile.imread('Data/03_Reconstructed/Left/Cam_Left_40_X1_N19_fdnorm.tif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "for cc in [20, 38, 39, 40, 55]:\n",
    "    plt.figure(figsize=(10,10)); plt.title('theirs %d'%cc)\n",
    "    plt.imshow(theirs[cc])\n",
    "    plt.colorbar(); plt.show()\n",
    "    if False:\n",
    "        plt.figure(figsize=(10,10)); plt.title('mineNorm %d'%cc)\n",
    "        plt.imshow(mineNorm[cc])\n",
    "        plt.colorbar(); plt.show()\n",
    "        plt.figure(figsize=(10,10)); plt.title('mineUnnorm %d'%cc)\n",
    "        plt.imshow(mineUnnorm[cc])\n",
    "        plt.colorbar(); plt.show()\n",
    "    plt.figure(figsize=(10,10)); plt.title('mineFDnorm %d'%cc)\n",
    "    plt.imshow(mineFDnorm[cc])\n",
    "    plt.colorbar(); plt.show()\n",
    "    frac = mineUnnorm[cc]/mineFDnorm[cc]\n",
    "    plt.figure(figsize=(10,10)); plt.title('mine unnorm/fdnorm %d (%f)'%(cc, (np.max(frac)-np.min(frac))/np.min(frac)))\n",
    "    plt.imshow(frac)\n",
    "    plt.colorbar(); plt.show()\n",
    "    frac = mineUnnorm[cc]/theirs[cc]\n",
    "    plt.figure(figsize=(10,10)); plt.title('mine unnorm/theirs %d (%f)'%(cc, (np.max(frac)-np.min(frac))/np.min(frac)))\n",
    "    plt.imshow(frac)\n",
    "    plt.colorbar(); plt.show()    \n",
    "    \n",
    "plt.figure(figsize=(10,10)); plt.title('theirs MIP')\n",
    "plt.imshow(np.max(theirs, axis=0))\n",
    "plt.colorbar(); plt.show()\n",
    "plt.figure(figsize=(10,10)); plt.title('mine unnorm MIP')\n",
    "plt.imshow(np.max(mineUnnorm, axis=0))\n",
    "plt.colorbar(); plt.show()\n",
    "plt.figure(figsize=(10,10)); plt.title('mine FDnorm MIP')\n",
    "plt.imshow(np.max(mineFDnorm, axis=0))\n",
    "plt.colorbar(); plt.show()\n",
    "\n",
    "frac = np.max(mineUnnorm,axis=0)/np.max(mineFDnorm,axis=0)\n",
    "plt.figure(figsize=(10,10)); plt.title('MIP mine unnorm/fdnorm MIP %f)'%((np.max(frac)-np.min(frac))/np.min(frac)))\n",
    "plt.imshow(frac)\n",
    "plt.colorbar(); plt.show()\n",
    "frac = np.max(mineUnnorm,axis=0)/np.max(theirs,axis=0)\n",
    "plt.figure(figsize=(10,10)); plt.title('MIP mine unnorm/theirs MIP (%f)'%((np.max(frac)-np.min(frac))/np.min(frac)))\n",
    "plt.imshow(frac)\n",
    "plt.colorbar(); plt.show()    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.5.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}