example_004_Nikon_NSIM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 004. Uncalibrated modulators: Nikon N-SIM three color immunofluorescence 2D-SIM image\n",
    "  \n",
    "Load raw Nikon N-SIM 2D-SIM data, determine SIM pattern infromation, and run reconstruction.\n",
    "\n",
    "Samples prepared by and imaging data acquired by Dr. Christophe Leterrier. COS cells immunofluorescence labeled for microtubules (640 nm ex.), clathrin (561 nm ex.), and actin (488 nm ex.).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import fft\n",
    "from pathlib import Path\n",
    "import mcsim.analysis.sim_reconstruction as sim\n",
    "from mcsim.analysis import analysis_tools\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import PowerNorm\n",
    "import napari\n",
    "import nd2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define data path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nd2_path = Path(\"data\", \"example_004\", \"raw_data\", \"neurocyte_lab_NSIM_data.nd2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parse raw Nikon N-SIM data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Nikon 2D SIM phases and angles\n",
    "n_phases = 3\n",
    "n_angles = 3\n",
    "\n",
    "# open ND2File\n",
    "try:\n",
    "    f = nd2.ND2File(nd2_path)\n",
    "except:\n",
    "    raise Exception('Cannot open specified file.')\n",
    "    \n",
    "# extract total image size\n",
    "# in N-SIM data, all 9 images are in one giant tiled image.\n",
    "# rows are angles, columns are phases\n",
    "height_pixels = f.attributes.heightPx\n",
    "width_pixels = f.attributes.widthPx\n",
    "\n",
    "# determine individual image size\n",
    "n_ypixels = int(height_pixels/n_angles)\n",
    "n_xpixels = int(height_pixels/n_phases)\n",
    "\n",
    "# extract experiment metadata\n",
    "sizes = f.sizes\n",
    "try:\n",
    "    n_channels = int(sizes['C'])\n",
    "except:\n",
    "    n_channels = 1\n",
    "try:\n",
    "    n_z_steps = int(sizes['Z'])\n",
    "except:\n",
    "    n_z_steps = 1  \n",
    "\n",
    "# extract per channel metadata\n",
    "wvl_ch_um = np.zeros(n_channels)\n",
    "NA_ch = np.zeros(n_channels)\n",
    "RI_ch = np.zeros(n_channels)\n",
    "pixel_size_ch_um = np.zeros([n_channels,3])\n",
    "for i in range(n_channels):\n",
    "    wvl_ch_um[i] = f.metadata.channels[i].channel.emissionLambdaNm / 1e3\n",
    "    NA_ch[i] = f.metadata.channels[i].microscope.objectiveNumericalAperture\n",
    "    RI_ch[i] = f.metadata.channels[i].microscope.immersionRefractiveIndex\n",
    "    pixel_size_ch_um[i,:] = np.flip(f.metadata.channels[i].volume.axesCalibration)\n",
    "\n",
    "# load data into memory\n",
    "data = f.asarray()\n",
    "\n",
    "# close ND2File\n",
    "f.close()\n",
    "\n",
    "# cut up raw ND2 data into individual channels, z planes, angles, and phases\n",
    "sim_images = np.zeros((n_channels,n_angles,n_phases,n_z_steps,n_ypixels,n_xpixels),dtype=np.uint16)\n",
    "for channel_idx in range(n_channels):\n",
    "    if n_z_steps == 1:\n",
    "        for angle_idx in range(n_angles):\n",
    "            for phase_idx in range(n_phases):\n",
    "                x_min = 0+(phase_idx*n_xpixels)\n",
    "                x_max = x_min + n_xpixels\n",
    "                y_min = 0+(angle_idx*n_ypixels)\n",
    "                y_max = y_min + n_ypixels\n",
    "                sim_images[channel_idx,angle_idx,phase_idx,0,:,:]=data[channel_idx,y_min:y_max,x_min:x_max]\n",
    "    else:\n",
    "        for z_idx in range(n_z_steps):\n",
    "            for angle_idx in range(n_angles):\n",
    "                for phase_idx in range(n_phases):\n",
    "                    x_min = 0+(phase_idx*n_xpixels)\n",
    "                    x_max = x_min + n_xpixels\n",
    "                    y_min = 0+(angle_idx*n_ypixels)\n",
    "                    y_max = y_min + n_ypixels\n",
    "                    sim_images[channel_idx,angle_idx,phase_idx,z_idx,:,:]=data[channel_idx,z_idx,y_min:y_max,x_min:x_max]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize raw Nikon N-SIM data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add images to napari viewer with scale information\n",
    "viewer = napari.view_image(data,name='Raw N-SIM',scale=(1,pixel_size_ch_um[0,1],pixel_size_ch_um[0,2]))\n",
    "viewer.add_image(sim_images,name='Parsed raw N-SIM',scale=(1,1,1,pixel_size_ch_um[0,1],pixel_size_ch_um[0,2]))\n",
    "\n",
    "# activate scale bar in physical units\n",
    "viewer.scale_bar.unit = 'um'\n",
    "viewer.scale_bar.visible = True\n",
    "\n",
    "# label sliders in napari viewer\n",
    "viewer.dims.axis_labels = ['c_parsed','a','p','c_raw','y','x']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate FFT of raw data to find frequency guesses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "z_idx=0\n",
    "for ch_idx in range(n_channels):\n",
    "\n",
    "    image_set = sim_images[ch_idx,:,:,z_idx,:,:]\n",
    "\n",
    "    dx = np.round(pixel_size_ch_um[ch_idx,1],3)\n",
    "    nx = n_xpixels\n",
    "    ny = n_ypixels\n",
    "    fxs = analysis_tools.get_fft_frqs(nx, dx)\n",
    "    df = fxs[1] - fxs[0]\n",
    "    fys = analysis_tools.get_fft_frqs(ny, dx)\n",
    "    ff = np.sqrt(np.expand_dims(fxs, axis=0)**2 + np.expand_dims(fys, axis=1)**2)\n",
    "    for ii in range(image_set.shape[0]):\n",
    "        ft = fft.fftshift(fft.fft2(fft.ifftshift(np.squeeze(image_set[ii, 0, :]))))\n",
    "\n",
    "        figh = plt.figure()\n",
    "        plt.title('Channel='+str(ch_idx)+\", Angle=\"+str(ii))\n",
    "        plt.imshow(np.abs(ft), norm=PowerNorm(gamma=0.1,vmin=65,vmax=1.5e7),\n",
    "                    extent=[fxs[0] - 0.5 * df, fxs[-1] + 0.5 * df, fys[-1] + 0.5 * df, fys[0] - 0.5 * df])\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define correct guesses for frequencies and phases based on FFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# frequency guesses from FFTs\n",
    "# 2D-SIM pattern angles are commonly 60 degrees (pi/3) rotated. Frequency, order, and relative rotation angle are often unknown.\n",
    "frequencies_ch = []\n",
    "frequencies_ch.append([[-3.3,-0.92], [-2.4,2.4], [-.85,-3.3]]) # ch 0\n",
    "frequencies_ch.append([[-3.7,-1.1], [-2.76,2.7], [-.95,-3.7]]) # ch 1\n",
    "frequencies_ch.append([[-4.6,-1.3], [-3.45,3.3], [-1.2,-4.6]]) # ch 2\n",
    "\n",
    "# phases guesses\n",
    "# 2D-SIM pattern phases are 120 degrees (2*pi/3) apart. Order and relative phase shift (for each angle) are often unknown.\n",
    "phases_ch =[]\n",
    "phases_ch.append([[0, 4*np.pi/3, 2*np.pi/3],[0, 4*np.pi/3, 2*np.pi/3],[0, 4*np.pi/3, 2*np.pi/3]]) # ch 0\n",
    "phases_ch.append([[0, 4*np.pi/3, 2*np.pi/3],[0, 4*np.pi/3, 2*np.pi/3],[0, 4*np.pi/3, 2*np.pi/3]]) # ch 1\n",
    "phases_ch.append([[0, 4*np.pi/3, 2*np.pi/3],[0, 4*np.pi/3, 2*np.pi/3],[0, 4*np.pi/3, 2*np.pi/3]]) # ch 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define SIM reconstruction parameters for all channels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define parameters for channels in order: [ch0,ch1,ch2]\n",
    "wiener_parameter=[0.4,0.4,.4]\n",
    "band0_exclusion_fraction=[0.5,0.5,0.5]\n",
    "backgrounds = [100,100,100]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Demonstrate SIM reconstruction output with incorrect parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# frequency guesses from FFTs -> use incorrect guess to demonstrate diagnostic output\n",
    "frequencies_ch_incorrect_guess = []\n",
    "frequencies_ch_incorrect_guess.append([[-4, 4], [-2,-4], [-3.5,4.5]]) # ch 0 intentionally incorrect guess\n",
    "frequencies_ch_incorrect_guess.append([[-3.7,-1.1], [-2.76,2.7], [-.95,-3.7]]) # ch 1\n",
    "frequencies_ch_incorrect_guess.append([[-4.6,-1.3], [-3.45,3.3], [-1.2,-4.6]]) # ch 2\n",
    "\n",
    "# phases guesses -> use incorrect guess to demonstrate diagnostic output\n",
    "phases_ch_incorrect_guess =[]\n",
    "phases_ch_incorrect_guess.append([[2*np.pi/3, 0, 4*np.pi/3],[2*np.pi/3, 0, 4*np.pi/3],[2*np.pi/3, 0, 4*np.pi/3]]) # ch 0 intentionally incorrect guess\n",
    "phases_ch_incorrect_guess.append([[0, 4*np.pi/3, 2*np.pi/3],[0, 4*np.pi/3, 2*np.pi/3],[0, 4*np.pi/3, 2*np.pi/3]]) # ch 1\n",
    "phases_ch_incorrect_guess.append([[0, 4*np.pi/3, 2*np.pi/3],[0, 4*np.pi/3, 2*np.pi/3],[0, 4*np.pi/3, 2*np.pi/3]]) # ch 2\n",
    "\n",
    "%matplotlib inline\n",
    "# path to save data\n",
    "save_path = Path(\"data\",\"example_004\",\"reconstruction_incorrect_guess\")\n",
    "\n",
    "# loop over all channels and z planes\n",
    "n_channels_limited = 1\n",
    "for ch_idx in range(n_channels_limited):\n",
    "    for z_idx in range(n_z_steps):\n",
    "\n",
    "        # define dictionary with physical parameters of system\n",
    "        physical_params = {'pixel_size': np.round(pixel_size_ch_um[ch_idx,1],3),\n",
    "                                        'na': NA_ch[ch_idx],\n",
    "                                        'wavelength': wvl_ch_um[ch_idx]}\n",
    "\n",
    "        # extract all angles/phases for current channel and z plane\n",
    "        sim_individual_channel = sim_images[ch_idx,:,:,z_idx,:,:]\n",
    "\n",
    "        # create mcSIM reconstruction object. See docstring for details on each parameter.\n",
    "        imgset = sim.SimImageSet(physical_params,\n",
    "                                sim_individual_channel,\n",
    "                                frq_guess=frequencies_ch_incorrect_guess[ch_idx],\n",
    "                                phases_guess=phases_ch_incorrect_guess[ch_idx],\n",
    "                                determine_amplitudes=True,\n",
    "                                wiener_parameter=wiener_parameter[ch_idx],\n",
    "                                fmax_exclude_band0=0*band0_exclusion_fraction[ch_idx],\n",
    "                                phase_estimation_mode=\"real-space\",\n",
    "                                gain=1,\n",
    "                                background=backgrounds[ch_idx],\n",
    "                                min_p2nr=0.8,\n",
    "                                max_phase_err=10*np.pi/180,\n",
    "                                save_dir=save_path,\n",
    "                                save_suffix='_z'+str(z_idx).zfill(3)+'_ch'+str(ch_idx).zfill(3))\n",
    "\n",
    "            \n",
    "        # perform reconstruction, plot figures, save, and clean up log file\n",
    "        imgset.reconstruct()\n",
    "        imgset.plot_figs()\n",
    "        imgset.save_imgs()\n",
    "        imgset.save_result()\n",
    "\n",
    "        # create variables to hold widefield and SIM SR images\n",
    "        if ch_idx == 0 and z_idx == 0:\n",
    "            wf_images_incorrect_guess = np.zeros((n_channels,n_z_steps,imgset.widefield.shape[0],imgset.widefield.shape[1]),dtype=np.float32)\n",
    "            SR_images_incorrect_guess = np.zeros((n_channels,n_z_steps,imgset.sim_sr.shape[0],imgset.sim_sr.shape[1]),dtype=np.float32)\n",
    "\n",
    "        # store widefield and SIM SR images for display\n",
    "        wf_images_incorrect_guess[ch_idx,z_idx,:]=imgset.widefield\n",
    "        SR_images_incorrect_guess[ch_idx,z_idx,:]=imgset.sim_sr\n",
    "\n",
    "        # clean up mcSIM reconstruction object\n",
    "        del sim_individual_channel, imgset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perform SIM reconstruction with correct parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "# path to save data\n",
    "save_path = Path(\"data\",\"example_004\",\"reconstruction\")\n",
    "\n",
    "# loop over all channels and z planes\n",
    "for ch_idx in range(n_channels):\n",
    "    for z_idx in range(n_z_steps):\n",
    "\n",
    "        # define dictionary with physical parameters of system\n",
    "        physical_params = {'pixel_size': np.round(pixel_size_ch_um[ch_idx,1],3),\n",
    "                                        'na': NA_ch[ch_idx],\n",
    "                                        'wavelength': wvl_ch_um[ch_idx]}\n",
    "\n",
    "        # extract all angles/phases for current channel and z plane\n",
    "        sim_individual_channel = sim_images[ch_idx,:,:,z_idx,:,:]\n",
    "\n",
    "        # create mcSIM reconstruction object. See docstring for details on each parameter.\n",
    "        imgset = sim.SimImageSet(physical_params,\n",
    "                                sim_individual_channel,\n",
    "                                frq_guess=frequencies_ch[ch_idx],\n",
    "                                phases_guess=phases_ch[ch_idx],\n",
    "                                determine_amplitudes=True,\n",
    "                                wiener_parameter=wiener_parameter[ch_idx],\n",
    "                                fmax_exclude_band0=band0_exclusion_fraction[ch_idx],\n",
    "                                phase_estimation_mode=\"real-space\",\n",
    "                                gain=1,\n",
    "                                background=backgrounds[ch_idx],\n",
    "                                min_p2nr=0.8,\n",
    "                                max_phase_err=10*np.pi/180,\n",
    "                                save_dir=save_path,\n",
    "                                save_suffix='_z'+str(z_idx).zfill(3)+'_ch'+str(ch_idx).zfill(3))\n",
    "\n",
    "            \n",
    "        # perform reconstruction, plot figures, save, and clean up log file\n",
    "        imgset.reconstruct()\n",
    "        imgset.plot_figs()\n",
    "        imgset.save_imgs()\n",
    "        imgset.save_result()\n",
    "\n",
    "        # create variables to hold widefield and SIM SR images\n",
    "        if ch_idx == 0 and z_idx == 0:\n",
    "            wf_images = np.zeros((n_channels,n_z_steps,imgset.widefield.shape[0],imgset.widefield.shape[1]),dtype=np.float32)\n",
    "            SR_images = np.zeros((n_channels,n_z_steps,imgset.sim_sr.shape[0],imgset.sim_sr.shape[1]),dtype=np.float32)\n",
    "\n",
    "        # store widefield and SIM SR images for display\n",
    "        wf_images[ch_idx,z_idx,:]=imgset.widefield\n",
    "        SR_images[ch_idx,z_idx,:]=imgset.sim_sr\n",
    "\n",
    "        # clean up mcSIM reconstruction object\n",
    "        del sim_individual_channel, imgset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Display results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define colormaps\n",
    "colormaps = ['bop purple', 'bop blue', 'bop orange']\n",
    "\n",
    "# add images to napari viewer with scale information, colormaps, and additive blending\n",
    "for ch_idx in range(n_channels):\n",
    "    if ch_idx == 0:\n",
    "        viewer = napari.view_image(wf_images[ch_idx,:],name='Widefield CH '+str(ch_idx).zfill(2),scale=(pixel_size_ch_um[ch_idx,1],pixel_size_ch_um[ch_idx,2]),colormap = colormaps[ch_idx],blending='additive',contrast_limits=[0,2**16-1])\n",
    "        viewer.add_image(SR_images[ch_idx,:],name='SR CH '+str(ch_idx).zfill(2),scale=(pixel_size_ch_um[ch_idx,1]/2,pixel_size_ch_um[ch_idx,2]/2),colormap = colormaps[ch_idx],blending='additive',contrast_limits=[0,2**16-1])\n",
    "    else:\n",
    "        viewer.add_image(wf_images[ch_idx,:],name='Widefield CH '+str(ch_idx).zfill(2),scale=(pixel_size_ch_um[ch_idx,1],pixel_size_ch_um[ch_idx,2]),colormap = colormaps[ch_idx],blending='additive',contrast_limits=[0,2**16-1])\n",
    "        viewer.add_image(SR_images[ch_idx,:],name='SR CH '+str(ch_idx).zfill(2),scale=(pixel_size_ch_um[ch_idx,1]/2,pixel_size_ch_um[ch_idx,2]/2),colormap = colormaps[ch_idx],blending='additive',contrast_limits=[0,2**16-1])\n",
    "\n",
    "# activate scale bar in physical units\n",
    "viewer.scale_bar.unit = 'um'\n",
    "viewer.scale_bar.visible = True"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "e966e8575e80f21eefd148db0808d44eea9284358157c0ee2fd32e18adfd524e"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}