miPolScope_fig1_u2os_fluorescence.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f8dbb2c7-38e8-42d8-9bf9-71bc42046cf7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import tifffile\n",
    "import napari\n",
    "\n",
    "from scipy.io import loadmat\n",
    "from skimage import transform\n",
    "\n",
    "from waveorder import waveorder_microscopy, fluorescence_microscopy, wavelet_softThreshold\n",
    "from waveorder.io import WaveorderReader"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "affec9b1-e0e6-4ee8-a360-e00841608961",
   "metadata": {},
   "source": [
    "# Table of contents\n",
    "- View dataset\n",
    "- Load data\n",
    "    - Load raw images\n",
    "    - Load calibration data\n",
    "    - Load background images\n",
    "- Recostruct fluorescence anisotropy channels\n",
    "    - Register images\n",
    "        - Crop edges\n",
    "        - View registered images\n",
    "    - Reconstruct Stokes images\n",
    "        - Initialize Stokes reconstructor\n",
    "        - Initialize fluorescence reconstructor\n",
    "        - Denoise intensity images\n",
    "            - View denoised intensity images\n",
    "        - Reconstruct Raw Stokes images\n",
    "        - Normalize Stokes images\n",
    "            - View normalized Stokes images and background images\n",
    "        - Correct background\n",
    "            - View background corrected Stokes images\n",
    "    - Compute anisotropy and orientation\n",
    "    - Compute deconvolved fluorescence intensity\n",
    "- View results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d4b5f02-4c85-4b35-b3aa-184b25c79e86",
   "metadata": {
    "tags": []
   },
   "source": [
    "# View dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "021d904b-6ea4-487c-bd06-6847ff1f9b84",
   "metadata": {},
   "outputs": [],
   "source": [
    "# assume data is in ~/Downloads folder\n",
    "data_dir = os.path.join(os.path.expanduser('~'), 'Downloads')\n",
    "data_path = os.path.join(data_dir, 'miPolScope_fig1_u2os_labelfree.zarr')\n",
    "\n",
    "# Check that data path exists\n",
    "if not os.path.exists(data_path):\n",
    "    raise ValueError('Data path does not exist.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6771ca36-e5d8-41ba-aa73-bfa28e068273",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "version mismatch: detected:FormatV01, requested:FormatV03\n",
      "version mismatch: detected:FormatV03, requested:FormatV01\n"
     ]
    }
   ],
   "source": [
    "viewer = napari.Viewer()\n",
    "layers = viewer.open(os.path.join(data_path, 'Row_0/Col_0/Pos_000'), plugin='napari-ome-zarr')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f6c38ae-7941-4eff-9e32-692487098d4a",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Load data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c2188fd-835b-45e6-99e4-8827a9765f34",
   "metadata": {},
   "source": [
    "## Load raw images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "70935830-0bfa-4743-9156-1851225d5b51",
   "metadata": {},
   "outputs": [],
   "source": [
    "wo_data = WaveorderReader(data_path, data_type='zarr')\n",
    "I = wo_data.get_zarr(0)\n",
    "n_timepoints, n_channels, n_slices, *img_size = I.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7c4d670d-d349-421b-9fdd-b6b7c3e43799",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of time points: 41\n",
      "Number of channels: 8\n",
      "Number of slices: 40\n",
      "Image size: [640, 640]\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    f'Number of time points: {n_timepoints}\\n'\n",
    "    f'Number of channels: {n_channels}\\n'\n",
    "    f'Number of slices: {n_slices}\\n'\n",
    "    f'Image size: {img_size}'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9c3d415-7f81-4361-80b0-78d874aeb9de",
   "metadata": {},
   "source": [
    "As demonstration, we will analyze only the first 5 timepoints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4ee336a2-e996-4573-b46b-0ea8dda4d11a",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_timepoints = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "bab182fd-5f1d-46c6-aa07-57dee44aaa9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data into memory\n",
    "I = np.array(I[:n_timepoints, 4:]) # fluorescence data is in the last four channels"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bb3354a-0780-4403-8567-85e3040db1de",
   "metadata": {},
   "source": [
    "## Load calibration data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c9ecf380",
   "metadata": {},
   "outputs": [],
   "source": [
    "cal_data = loadmat(os.path.join(data_path,'calibration_fluor.mat'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "63f3b269",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.transpose(cal_data['A'].astype('float'), (2, 3, 0, 1)) # A has shape (size_Y, size_X, N_channels, N_Stokes)\n",
    "black_level = cal_data['black_level'][0][0].astype('uint16')\n",
    "\n",
    "tform0 = transform.AffineTransform(cal_data['tform0'].T)\n",
    "tform45 = transform.AffineTransform(cal_data['tform45'].T)\n",
    "tform90 = transform.AffineTransform(cal_data['tform90'].T)\n",
    "tform135 = transform.AffineTransform(cal_data['tform135'].T)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a69c92e-859b-4ffe-a047-8409baa730df",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Load background images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "6b12c342",
   "metadata": {},
   "outputs": [],
   "source": [
    "S1_bg = cal_data['S1_bg'][0][0]\n",
    "S2_bg = cal_data['S2_bg'][0][0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6beb0c51-9b1e-4b71-bf1e-8d843689239f",
   "metadata": {},
   "source": [
    "# Reconstruct fluorescence anisotropy channels"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "695705ff-50d5-43cc-a346-dcb32730a0ff",
   "metadata": {},
   "source": [
    "## Register images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ace91f65",
   "metadata": {},
   "outputs": [],
   "source": [
    "I_registered = np.zeros((n_timepoints, 4, n_slices, *img_size), dtype='float')\n",
    "\n",
    "for t in range(n_timepoints):\n",
    "    for c, tform in enumerate((tform0, tform45, tform90, tform135)):\n",
    "        for z in range(n_slices):\n",
    "            I_registered[t,c,z] = transform.warp(I[t,c,z], tform.inverse, preserve_range=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fedd631-c1a2-4db2-b06e-0e4628674bcd",
   "metadata": {},
   "source": [
    "### Crop edges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "40a1e736",
   "metadata": {},
   "outputs": [],
   "source": [
    "I_registered = I_registered[..., 20:-20, 20:-20]\n",
    "img_size = I_registered.shape[-2:]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56083440-2890-48c6-bb9d-7343ba0d1412",
   "metadata": {},
   "source": [
    "### View registered images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "8c06383e",
   "metadata": {},
   "outputs": [],
   "source": [
    "viewer = napari.view_image(I_registered, contrast_limits=(99, 150))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abad6661-70c1-45e7-bd93-919f79657880",
   "metadata": {},
   "source": [
    "## Reconstruct Stokes images"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0444fd3",
   "metadata": {},
   "source": [
    "### Initialize Stokes reconstructor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "6cbcce68",
   "metadata": {},
   "outputs": [],
   "source": [
    "# z projection parameters\n",
    "z_chunk_size = 4\n",
    "n_slices = n_slices//z_chunk_size\n",
    "\n",
    "wavelength = 670 # in nm\n",
    "NA_obj = 1.2 # Numerical Aperture of Objective\n",
    "NA_illu = 0.4 # Numerical Aperture of Condenser\n",
    "n_objective_media = 1.33 # refractive index of objective immersion media\n",
    "mag = 30 # effective magnification\n",
    "n_slices = n_slices # number of slices in z-stack\n",
    "z_step_um = 0.25 * z_chunk_size # z-step size in um\n",
    "pad_z = 5 # slices to pad for phase reconstruction boundary artifacts\n",
    "pixel_size_um = 3.45 # camera pixel size in um\n",
    "bg_correction = 'None' # BG correction method: 'None', 'local_fit', 'global'\n",
    "mode = '3D' # phase reconstruction mode, '2D' or '3D'\n",
    "use_gpu = False\n",
    "gpu_id = 0 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "8ed88f82",
   "metadata": {},
   "outputs": [],
   "source": [
    "z_defocus = -(np.r_[:n_slices] - n_slices // 2) * z_step_um # assumes stack starts from the bottom\n",
    "swing = 0\n",
    "ps = pixel_size_um / mag\n",
    "\n",
    "reconstructor = waveorder_microscopy(img_dim=img_size,\n",
    "                                     lambda_illu=wavelength/1000,\n",
    "                                     ps=ps,\n",
    "                                     NA_obj=NA_obj,\n",
    "                                     NA_illu=NA_illu,\n",
    "                                     z_defocus=z_defocus,\n",
    "                                     chi=swing,\n",
    "                                     n_media=n_objective_media,\n",
    "                                     cali=True,\n",
    "                                     bg_option=bg_correction,\n",
    "                                     A_matrix=A,\n",
    "                                     QLIPP_birefringence_only=True,\n",
    "                                     pad_z=pad_z,\n",
    "                                     phase_deconv=mode,\n",
    "                                     illu_mode='BF',\n",
    "                                     use_gpu=use_gpu,\n",
    "                                     gpu_id=gpu_id)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5bb01c7",
   "metadata": {},
   "source": [
    "### Initialize fluorescence reconstructor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "40602c1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "lambda_emiss = [0.670]  # emission wavelength of the fluorescence channel (list, in um)\n",
    "\n",
    "fluor_reconstructor = fluorescence_microscopy((*img_size, n_slices), lambda_emiss, ps, z_step_um, NA_obj, \n",
    "                                                n_media=n_objective_media,\n",
    "                                                deconv_mode='3D-WF', \n",
    "                                                pad_z=3, \n",
    "                                                use_gpu=use_gpu, \n",
    "                                                gpu_id=gpu_id)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2fd08a9-5bba-4a14-a7ca-185865d03106",
   "metadata": {},
   "source": [
    "### Denoise intensity images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "acf30ca1",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\ivan.ivanov\\AppData\\Local\\Continuum\\miniconda3\\envs\\recorder\\lib\\site-packages\\pywt\\_multilevel.py:43: UserWarning: Level value of 2 is too high: all coefficients will experience boundary effects.\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "# z projection\n",
    "I_denoised = np.reshape(I_registered, (n_timepoints, 4, n_slices, z_chunk_size, *img_size)).mean(axis=3)\n",
    "\n",
    "# wavelet denoising\n",
    "for c in range(4):\n",
    "    I_denoised[:,c] = wavelet_softThreshold(I_denoised[:,c], 'db8', 1, level=2, axes=(1,2,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46f53900-3b85-4f83-985e-098da6fb42e7",
   "metadata": {},
   "source": [
    "#### View denoised intensity images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "8a15e8aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "viewer = napari.view_image(I_denoised, contrast_limits=(99, 150))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46d1242e-9a36-4283-86da-1cedbd4e414c",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Reconstruct Raw Stokes images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "7f26bc4f-03d2-4b71-9553-f0b82d0176ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "S_raw = np.zeros((n_timepoints, 3, n_slices, *img_size), dtype='float')\n",
    "\n",
    "for t in range(n_timepoints):\n",
    "    S_raw_ = reconstructor.Stokes_recon(np.moveaxis(I_denoised[t], 1, -1) - black_level)\n",
    "    S_raw[t] = np.moveaxis(S_raw_, -1, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3bed893-797c-412d-bcaf-4339fe07ed92",
   "metadata": {},
   "source": [
    "### Normalize Stokes images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "b32c3516",
   "metadata": {},
   "outputs": [],
   "source": [
    "S_norm = np.zeros_like(S_raw)\n",
    "\n",
    "for t in range(n_timepoints):\n",
    "    for z in range(n_slices):\n",
    "        S_norm[t,:,z] = reconstructor.Stokes_transform(S_raw[t,:,z])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef42793c-6080-4bae-bffe-d8897a8eda71",
   "metadata": {},
   "source": [
    "#### View normalized Stokes images and background images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "2f510cd3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Image layer 'S2_bg' at 0x156c1ef0d00>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "viewer = napari.Viewer()\n",
    "viewer.add_image(S_norm[:,0], name='S0_norm', colormap='gray', contrast_limits=(0, 120))\n",
    "viewer.add_image(S_norm[:,1], name='S1_norm', colormap='RdBu', visible=False, contrast_limits=(-0.25, 0.25))\n",
    "viewer.add_image(S_norm[:,2], name='S2_norm', colormap='RdBu', visible=False, contrast_limits=(-0.25, 0.25))\n",
    "\n",
    "viewer.add_image(S1_bg*np.ones(img_size), name='S1_bg', colormap='RdBu', visible=False, contrast_limits=(-0.25, 0.25))\n",
    "viewer.add_image(S2_bg*np.ones(img_size), name='S2_bg', colormap='RdBu', visible=False, contrast_limits=(-0.25, 0.25))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81b99e90-6512-488f-8e97-5c2311d82f19",
   "metadata": {},
   "source": [
    "### Correct background"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "416f75b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "S_corr = np.zeros_like(S_norm)\n",
    "\n",
    "S_corr[:,0] = S_norm[:,0]\n",
    "S_corr[:,1] = S_norm[:,1] - S1_bg\n",
    "S_corr[:,2] = S_norm[:,2] - S2_bg"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c263a48-7b96-47ce-833d-d0ffbf1dde42",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### View background corrected Stokes images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "91217848",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Image layer 'S2' at 0x156c2dfbe50>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "viewer = napari.Viewer()\n",
    "viewer.add_image(S_corr[:,0], name='S0', colormap='gray', contrast_limits=(0, 120))\n",
    "viewer.add_image(S_corr[:,1], name='S1', colormap='RdBu', visible=False, contrast_limits=(-0.25, 0.25))\n",
    "viewer.add_image(S_corr[:,2], name='S2', colormap='RdBu', visible=False, contrast_limits=(-0.25, 0.25))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2804b7df-7b10-44ab-8eed-b0e465fce2d7",
   "metadata": {},
   "source": [
    "## Compute anisotropy and orientation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "15bf477d",
   "metadata": {},
   "outputs": [],
   "source": [
    "anisotropy, orientation = fluor_reconstructor.Fluor_anisotropy_recon(S_corr[:,1], S_corr[:,2])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1c9de3c-a7c9-472a-8970-b88e09957b7d",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Compute deconvolved fluorescence intensity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "bf5cdeff-bd2a-4e29-b9fa-3a361ae46ab6",
   "metadata": {},
   "outputs": [],
   "source": [
    "blackground_level = [8]\n",
    "I_deconvolved = np.zeros((n_timepoints, n_slices, *img_size), dtype='float')\n",
    "\n",
    "for t in range(n_timepoints):\n",
    "    I_ = fluor_reconstructor.deconvolve_fluor_3D(np.moveaxis(S_raw[t,0], 0, -1), blackground_level, reg=[1e-2], verbose=False)\n",
    "    I_deconvolved[t] = np.moveaxis(I_, -1, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7beade26-39b9-43ca-a1ad-03d4e165c390",
   "metadata": {},
   "source": [
    "# View results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "d0e77e33",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Image layer 'orientation' at 0x156fe71bd00>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "viewer = napari.Viewer()\n",
    "viewer.add_image(I_deconvolved, name='deconvolved intensity', colormap='gray', contrast_limits=(0, 200))\n",
    "viewer.add_image(anisotropy, name='anisotropy', colormap='gray', visible=False, contrast_limits=(0, 0.2))\n",
    "viewer.add_image(orientation, name='orientation', colormap='gray', visible=False, contrast_limits=(0, np.pi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be78cd52-a32a-426b-aee0-a5f16e021aff",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}