python/tutorials/tutorial_workspace_manipulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Tutorial: Workspace Manipulations\n",
    "In this tutorial we demonstrate how to create a workspace and perform basic operations on it, including cropping, rotation, thresholding."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we import puma. Note that in order to run pumapy, the pumapy conda environment must first be activated, by executing \"conda activate puma\" in a terminal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pumapy as puma\n",
    "import os\n",
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A workspace is the datastructure at the basis of both PuMA and pumapy and it is basically a container for the material sample that you want to analyze. A workspace is made of little cubes, or voxels (i.e. 3D pixels), holding a value. This simple element definition (formally called Cartesian grid) allows for very fast operations. Inside a workspace object, two different arrays are defined: one called \"matrix\" and the other called \"orientation\". Both of these are nothing but a 3D Numpy array for the matrix (X,Y,Z dimensions of the domain) and a 4D Numpy array for the orientation (dimensions of X,Y,Z,3 for vectors throughout the domain). \n",
    "\n",
    "Next we show the different ways we have implemented to define a workspace class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape of workspace 1: (10, 11, 12)\n",
      "\n",
      "Shape of workspace 2: (20, 31, 212)\n",
      "\n",
      "Matrix shape of workspace 3: (5, 6, 2)\n",
      "Orientation shape of workspace 3: (5, 6, 2, 3)\n",
      "Display Workspace 3 matrix\n",
      "\n",
      "3D Workspace:\n",
      "  o---> y\n",
      "  |\n",
      "x v\n",
      "[(:,:,0)\n",
      "[[1 1 1 1 1 1]\n",
      "[1 1 1 1 1 1]\n",
      "[1 1 1 1 1 1]\n",
      "[1 1 1 1 1 1]\n",
      "[1 1 1 1 1 1]\n",
      "\n",
      "(:,:,1)\n",
      "[[1 1 1 1 1 1]\n",
      "[1 1 1 1 1 1]\n",
      "[1 1 1 1 1 1]\n",
      "[1 1 1 1 1 1]\n",
      "[1 1 1 1 1 1]]\n",
      "\n",
      " Display Workspace 3 orientation\n",
      "\n",
      "3D Orientation:\n",
      "  o---> y\n",
      "  |\n",
      "x v\n",
      "[(:,:,0)\n",
      "[[(0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0)]\n",
      "[(0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0)]\n",
      "[(0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0)]\n",
      "[(0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0)]\n",
      "[(0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0)]\n",
      "\n",
      "(:,:,1)\n",
      "[[(0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0)]\n",
      "[(0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0)]\n",
      "[(0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0)]\n",
      "[(0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0)]\n",
      "[(0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0) (0.4, 2.0, 5.0)]]\n",
      "\n",
      " Display Workspace 5\n",
      "\n",
      "3D Workspace:\n",
      "  o---> y\n",
      "  |\n",
      "x v\n",
      "[(:,:,0)\n",
      "[[3 4 4 0 3 3 1 1 3 1]\n",
      "[1 2 1 1 2 4 2 3 1 4]\n",
      "[0 0 3 0 2 0 2 1 4 0]\n",
      "[0 2 3 2 0 4 4 1 1 0]\n",
      "[0 0 2 1 4 1 3 1 1 2]\n",
      "[2 3 3 2 0 4 0 0 2 2]\n",
      "[3 0 4 4 2 3 1 1 3 0]\n",
      "[0 1 3 4 2 0 1 4 2 1]\n",
      "[0 3 0 0 1 3 4 4 3 3]\n",
      "[4 2 0 3 2 0 1 1 4 0]\n",
      "\n",
      "(:,:,1)\n",
      "[[4 2 0 2 3 1 2 4 4 2]\n",
      "[2 0 4 4 0 1 0 4 4 3]\n",
      "[0 0 2 3 1 3 4 0 3 4]\n",
      "[3 0 2 1 2 0 2 4 0 2]\n",
      "[0 0 4 4 1 0 0 0 2 0]\n",
      "[0 0 0 2 0 0 2 0 0 1]\n",
      "[2 2 2 3 4 3 2 4 1 3]\n",
      "[0 3 4 3 0 4 0 3 3 2]\n",
      "[1 0 2 2 1 0 3 2 0 3]\n",
      "[3 3 2 0 3 1 1 3 0 0]\n",
      "\n",
      "(:,:,2)\n",
      "[[0 1 4 2 2 4 2 1 1 0]\n",
      "[4 2 0 2 3 0 1 4 0 3]\n",
      "[1 3 4 3 3 2 4 1 4 0]\n",
      "[0 4 2 3 3 1 1 3 2 3]\n",
      "[4 4 0 1 2 0 1 3 3 1]\n",
      "[3 3 2 2 2 2 2 3 0 4]\n",
      "[4 4 2 4 0 2 1 3 0 3]\n",
      "[2 4 1 4 3 3 4 3 1 4]\n",
      "[1 4 4 1 1 4 0 3 2 1]\n",
      "[4 3 0 1 2 0 1 2 0 4]]\n"
     ]
    }
   ],
   "source": [
    "# defines a workspace full of zeros of shape 10x11x12\n",
    "ws1 = puma.Workspace.from_shape((10, 11, 12))\n",
    "print(\"Shape of workspace 1: {}\\n\".format(ws1.matrix.shape))\n",
    "\n",
    "# defines a workspace of shape 10x11x12, full of a custom value (in this case ones)\n",
    "ws2 = puma.Workspace.from_shape_value((20, 31, 212), 1)\n",
    "print(\"Shape of workspace 2: {}\\n\".format(ws2.matrix.shape))\n",
    "\n",
    "# defines a workspace of shape 10x11x12, full of a custom value (in this case ones)\n",
    "ws3 = puma.Workspace.from_shape_value_vector((5, 6, 2), 1, (0.4, 2, 5))\n",
    "print(\"Matrix shape of workspace 3: {}\".format(ws3.matrix.shape))\n",
    "print(\"Orientation shape of workspace 3: {}\".format(ws3.orientation.shape))\n",
    "print(\"Display Workspace 3 matrix\")\n",
    "ws3.show_matrix()\n",
    "print(\"\\n Display Workspace 3 orientation\")\n",
    "ws3.show_orientation()\n",
    "\n",
    "# we can also convert a Numpy array into a Workspace as follows:\n",
    "array = np.random.randint(5, size=(10, 10, 10))\n",
    "ws4 = puma.Workspace.from_array(array)\n",
    "\n",
    "# finally, we can also create an empty workspace object and assign its matrix directly as:\n",
    "ws5 = puma.Workspace()\n",
    "ws5.matrix = np.random.randint(5, size=(10, 10, 3))\n",
    "print(\"\\n Display Workspace 5\")\n",
    "ws5.show_matrix()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is important to keep the first three dimensions (X,Y,Z) of the matrix and orientation class variables the same. This is automatically enforced by using the class methods, but it is not when assigning them directly as in the last two examples. \n",
    "\n",
    "We can import a tomography image directly into a workspace: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Importing /Users/fsemerar/Documents/PuMA_playground/puma-dev/python/tests/testdata/200_FiberForm.tif ... Done\n"
     ]
    }
   ],
   "source": [
    "ws_raw = puma.import_3Dtiff(\"python/tests/testdata/200_FiberForm.tif\", 1.3e-6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The voxel length of the workspace can either be set during import of a 3D tiff, or manually afterwards, as shown below: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "ws_raw.voxel_length = 1.3e-6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visualize its slices by running the command below. By scrolling on top of the plot, you can slice through the material along the z axis. You can also use the left/right arrows on the keyboard to skip +/-10 slices or the up/down arrows to skip +/-100 slices. In addition, on the bottom of the plot, the (x,y) coordinates are shown along with the corresponding grayscale value. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "85ddd74acb754f559be3c8f0379279c6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<pumapy.visualization.slicer.PlotSlicer at 0x7fc7b1281990>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "puma.plot_slices(ws_raw)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we show how to manipulate the domain, e.g. crop, rescale, resize and rotate it. \n",
    "\n",
    "An approach to crop a domain is the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape of original workspace: (200, 200, 200)\n",
      "Shape of cropped workspace: (30, 165, 20)\n"
     ]
    }
   ],
   "source": [
    "ws_copy = ws_raw.copy()\n",
    "ws_copy.matrix = ws_copy.matrix[10:40, 35:, -20:]\n",
    "print(\"Shape of original workspace: {}\".format(ws_raw.get_shape()))\n",
    "print(\"Shape of cropped workspace: {}\".format(ws_copy.get_shape()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, it is important to not fall in the trap of referencing the same Numpy array. Here is an example of how YOU SHOULDN'T perform cropping:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "85\n",
      "85\n"
     ]
    }
   ],
   "source": [
    "ws_bad = puma.Workspace()\n",
    "ws_bad.matrix = ws_raw[10:40, 35:, -20:]\n",
    "ws_bad[0, 0, 0] = np.random.randint(0, 255)\n",
    "print(ws_raw.matrix[10, 35, -20])\n",
    "print(ws_bad.matrix[0, 0, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see from the output, now both the original Workspace and the newly created one share the same Numpy array for the matrix class variable (the second one is only a section of it). This way, when one is changed, the other one is changed as well. It is important to make a copy of a domain if the original workspace needs to be kept.\n",
    "\n",
    "Next, we show how we can rescale a domain by a factor or resize it to a specified size. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rescaled workspace size: (100, 100, 100)\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9b915c5afcba422d82e69b358ed00836",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<pumapy.visualization.slicer.CompareSlicer at 0x7fc7a082ab50>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ws_copy = ws_raw.copy()\n",
    "ws_copy.rescale(scale=0.5, segmented=False)\n",
    "\n",
    "# Notice that now the axis have different limits\n",
    "puma.compare_slices(ws_raw, ws_copy)"
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}