sample_script.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import numpy as np\n",
    "import scipy.io as sio\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "import contexttimer\n",
    "%matplotlib notebook\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "import arrayfire as af\n",
    "af.set_device(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specify the path of your repo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from opticaltomography.opticsutil import compare3DStack, show3DStack\n",
    "from opticaltomography.opticsalg import PhaseObject3D, TomographySolver, AlgorithmConfigs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Specify parameters & load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Units in microns\n",
    "wavelength = 0.514\n",
    "n_measure = 1.0\n",
    "n_b = 1.0\n",
    "maginification = 80.\n",
    "dx = 6.5 / maginification\n",
    "dy = 6.5 / maginification\n",
    "dz = 3 * dx\n",
    "na = 0.65"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Make sure the path is correct\n",
    "#Illumination angle, change to [0.0] if only on-axis is needed:\n",
    "na_list = sio.loadmat(\"na_list_test.mat\")\n",
    "fx_illu_list = na_list[\"na_list\"][150:,0] / wavelength\n",
    "fy_illu_list = na_list[\"na_list\"][150:,1] / wavelength"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot object in z (y,x,z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "phantom = np.ones((400,400,10),dtype=\"complex64\") * n_b\n",
    "show3DStack(np.real(phantom), axis=2, clim=(np.min(np.real(phantom)), np.max(np.real(phantom))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fill in phantom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = np.meshgrid(np.linspace(-1,1,phantom.shape[0]), np.linspace(-1,1,phantom.shape[1]))\n",
    "r2 = x ** 2 + y ** 2\n",
    "phantom[...,4] += (r2 < 0.25 ** 2) * 0.1 / (2 * np.pi * dz / wavelength)\n",
    "show3DStack(np.real(phantom), axis=2, clim=(np.min(np.real(phantom)), np.max(np.real(phantom))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup solver objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "solver_params = dict(wavelength = wavelength, na = na, \\\n",
    "                     RI_measure = n_measure, sigma = 2 * np.pi * dz / wavelength,\\\n",
    "                     fx_illu_list = fx_illu_list, fy_illu_list = fy_illu_list,\\\n",
    "                     pad = True, pad_size = (25,25))\n",
    "phase_obj_3d = PhaseObject3D(shape=phantom.shape, voxel_size=(dy,dx,dz), RI=n_b, RI_obj=phantom)\n",
    "solver_obj   = TomographySolver(phase_obj_3d, **solver_params)\n",
    "# Forward simulation method\n",
    "# solver_obj.setScatteringMethod(model = \"MultiPhaseContrast\")\n",
    "solver_obj.setScatteringMethod(model = \"MultiBorn\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate forward prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with contexttimer.Timer() as timer:\n",
    "    forward_field_mb = solver_obj.forwardPredict(field=False)\n",
    "    print(timer.elapsed)   \n",
    "forward_field_mb = np.squeeze(forward_field_mb)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot\n",
    "%matplotlib notebook\n",
    "show3DStack(np.real(forward_field_mb), axis=2, clim=(np.min(np.real(forward_field_mb)), np.max(np.real(forward_field_mb))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving an inverse problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create a class for all inverse problem parameters\n",
    "configs            = AlgorithmConfigs()\n",
    "configs.batch_size = 1\n",
    "configs.method     = \"FISTA\"\n",
    "configs.restart    = True\n",
    "configs.max_iter   = 5\n",
    "# multislice stepsize\n",
    "# configs.stepsize   = 2e-4\n",
    "# multiborn stepsize\n",
    "configs.stepsize   = 10\n",
    "configs.error      = []\n",
    "configs.pure_real = True\n",
    "#total variation regularization\n",
    "configs.total_variation     = False\n",
    "configs.reg_tv              = 1.0 #lambda\n",
    "configs.max_iter_tv         = 15\n",
    "configs.order_tv            = 1\n",
    "configs.total_variation_gpu = True\n",
    "configs.total_variation_anisotropic = False\n",
    "\n",
    "# reconstruction method\n",
    "# solver_obj.setScatteringMethod(model = \"MultiPhaseContrast\")\n",
    "solver_obj.setScatteringMethod(model = \"MultiBorn\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "recon_obj_3d = solver_obj.solve(configs, forward_field_mb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "current_rec = recon_obj_3d\n",
    "cost = solver_obj.configs.error\n",
    "show3DStack(np.real(current_rec), axis=2, clim=(np.min(np.real(current_rec)), np.max(np.real(current_rec))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(np.log10(cost))"
   ]
  },
  {
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}