Image Windowing and Matrix Factorization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6bcce3c2",
   "metadata": {},
   "source": [
    "# Image Windowing and Matrix Factorization\n",
    "\n",
    "### Rama Vasudevan\n",
    "\n",
    "### Px Credits: Gerd Duscher, Suhas Somnath, S. Mani Valleti, Maxim Ziatdinov, Rama Vasudevan\n",
    "\n",
    "Let's look at some image analytics here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "233273c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "!pip install -U pyNSID sidpy SciFiReaders nanonispy gwyfile pycroscopy\n",
    "#!pip install pip install threadpoolctl==3.1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b992e8e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import pyNSID\n",
    "import sidpy as sid\n",
    "import SciFiReaders as sr\n",
    "from sidpy.proc.fitter import SidFitter"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73507c41",
   "metadata": {},
   "source": [
    "# Second example - Image analysis\n",
    "\n",
    "The pycroscopy package has some tools for generic image analysis, as well as wrappers around common machine learning methods. These include matrix and tensor factorization techniques. Let us explore one example.\n",
    "\n",
    "First we will import a microscopy image, and then we will perform image windowing\n",
    "\n",
    "We will then use matrix factorization to analyze the spatial distribution of different phases\n",
    "\n",
    "This is explained in <a href = \"https://pubs.acs.org/doi/full/10.1021/acs.nanolett.6b02130\">this article</a>.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "baaf0522",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pyNSID\n",
    "\n",
    "import sidpy as sid\n",
    "import SciFiReaders as sr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b5d824c",
   "metadata": {},
   "outputs": [],
   "source": [
    "dm3_file = r'data/bto_atomic.dm3'\n",
    "png_file = r'data/re48.tiff'\n",
    "\n",
    "dm3_reader = sr.DM3Reader(dm3_file)\n",
    "img_reader = sr.ImageReader(png_file)\n",
    "\n",
    "data = dm3_reader.read()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6946005d",
   "metadata": {},
   "outputs": [],
   "source": [
    "data.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65a6c0fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let's crop the image and try again\n",
    "# Make a sidpy dataset\n",
    "\n",
    "x_dim = data.dim_0.values\n",
    "y_dim = data.dim_1.values\n",
    "\n",
    "data_cropped = data[:][200:-200,10:550]\n",
    "\n",
    "data_set = sid.Dataset.from_array(data_cropped, title='BTO_STEM')\n",
    "\n",
    "# Set the data type\n",
    "data_set.data_type = sid.DataType.IMAGE\n",
    "\n",
    "# Add quantity and units\n",
    "data_set.units = 'counts'\n",
    "data_set.quantity = 'Intensity'\n",
    "\n",
    "y_dim = y_dim[10:550]\n",
    "x_dim = x_dim[200:-200]\n",
    "\n",
    "# Add dimension info\n",
    "data_set.set_dimension(0, sid.Dimension(x_dim,\n",
    "                                        name='x',\n",
    "                                        units='nm', quantity='x',\n",
    "                                        dimension_type='spatial'))\n",
    "data_set.set_dimension(1, sid.Dimension(y_dim,\n",
    "                                        name='y',\n",
    "                                        units='nm', quantity='y',\n",
    "                                        dimension_type='spatial'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b360dc78",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_set.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "805bce73",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from pycroscopy.image import ImageWindowing\n",
    "\n",
    "parms_dict = {}\n",
    "parms_dict['window_step_x'] = 16\n",
    "parms_dict['window_step_y'] = 16\n",
    "parms_dict['window_size_x'] = 64\n",
    "parms_dict['window_size_y'] = 64\n",
    "parms_dict['mode'] = 'fft'\n",
    "parms_dict['filter'] = 'hamming'\n",
    "parms_dict['zoom_factor'] = 2\n",
    "parms_dict['interpol_factor'] = 1\n",
    "iw = ImageWindowing(parms_dict)\n",
    "windows = iw.MakeWindows(data_set)\n",
    "windows = np.abs(np.log(np.abs(windows)))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b896200c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "windows.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f9ad82f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pycroscopy.learn.ml.matrix_factor import MatrixFactor\n",
    "mfactor = MatrixFactor(np.abs(windows), method = 'pca',n_components = 3 )\n",
    "output = mfactor.do_fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9880227",
   "metadata": {},
   "outputs": [],
   "source": [
    "abundances = output[0]\n",
    "components = output[1]\n",
    "abund = np.array(abundances)\n",
    "comps = np.array(components)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61908793",
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "\n",
    "fig, axes = plt.subplots(nrows=1, ncols=mfactor.ncomp, figsize = (10,3))\n",
    "for ind, ax in enumerate(axes.flat):\n",
    "    im1 = ax.imshow(comps[ind,:,:])\n",
    "    ax.set_title('Component #' + str(ind))\n",
    "    divider = make_axes_locatable(ax)\n",
    "    cax = divider.append_axes('right', size='5%', pad=0.05)\n",
    "    fig.colorbar(im1, cax=cax, orientation='vertical')\n",
    "    ax.axis('off')\n",
    "fig.tight_layout()\n",
    "\n",
    "fig, axes = plt.subplots(nrows=1, ncols=mfactor.ncomp, figsize = (10,3))\n",
    "for ind, ax in enumerate(axes.flat):\n",
    "    im1 = ax.imshow(abund[:,:,ind])\n",
    "    divider = make_axes_locatable(ax)\n",
    "    cax = divider.append_axes('right', size='5%', pad=0.05)\n",
    "    fig.colorbar(im1, cax=cax, orientation='vertical')\n",
    "    ax.axis('off')\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "005f6d1c",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "1. Play with the number of components present and see how it affects the result. Similarly, play with the different types of factorization ('nmf', 'pca', 'nfindr'). How do the results differ?\n",
    "\n",
    "2. Think about how we might want to determine how many components we should use. HINT: Look at reconstruction errors, or similarity of the components.\n",
    "\n",
    "[RKV to do: write the code for structural similarity of the components, if any two have a similarity beyond a certain level then we can cut off there]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26e5eb78",
   "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
}