2D_analysis.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "2D-analysis.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ncFHiRF9REx6"
      },
      "outputs": [],
      "source": [
        "!pip install SimpleITK\n",
        "!pip install skan \n",
        "!pip install tifffile\n",
        "!pip install plotly==5.3.1"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "from google.colab import drive\n",
        "drive.mount('/content/drive')"
      ],
      "metadata": {
        "id": "rtdDGBuyRL1l"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import SimpleITK as sitk\n",
        "from skimage.morphology import skeletonize, thin, medial_axis, skeletonize_3d\n",
        "from scipy import ndimage\n",
        "from skan import skeleton_to_csgraph, Skeleton\n",
        "from skan import summarize\n",
        "from skan import draw\n",
        "import tifffile as tiff\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "from mpl_toolkits.mplot3d import Axes3D\n",
        "import numpy as np\n",
        "import cv2\n",
        "import plotly.express as px\n",
        "import plotly.graph_objects as go\n",
        "import pandas as pd\n",
        "import glob\n",
        "from tqdm.notebook import tnrange"
      ],
      "metadata": {
        "id": "Ybn_et51ROYz"
      },
      "execution_count": 18,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def load_img_from_tiff(path2img):\n",
        "    \"\"\"\n",
        "    Parameters\n",
        "    ----------\n",
        "    path2img: str \n",
        "        path to image Tiff file\n",
        "    \n",
        "    Returns\n",
        "    -------\n",
        "    img_array: np.array\n",
        "        image data in numpy format\n",
        "    \"\"\"\n",
        "    img = sitk.ReadImage(path2img)\n",
        "    img_array = sitk.GetArrayFromImage(img)\n",
        "    return(img_array)"
      ],
      "metadata": {
        "id": "yhknJLeORQmu"
      },
      "execution_count": 19,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "single_thresh = load_img_from_tiff(\"/content/drive/MyDrive/mydata/sea_urchin_data/Test n.1\")\n",
        "# single_thresh = cv2.resize(single_thresh, dsize=(256,256), interpolation=cv2.INTER_CUBIC)\n",
        "\n",
        "single_thresh = single_thresh.astype('float64')\n",
        "single_thresh = (single_thresh - np.min(single_thresh))/np.ptp(single_thresh)\n",
        "skimage_skeleton, dist = medial_axis(single_thresh, return_distance=True)\n",
        "skan_skeleton = Skeleton(dist*skimage_skeleton)\n",
        "\n",
        "lee94_skimage_skeleton = skeletonize(single_thresh, method='lee')\n",
        "lee94_skimage_skeleton=lee94_skimage_skeleton/255\n",
        "lee_skan_skeleton=Skeleton(dist*lee94_skimage_skeleton)\n",
        "df=summarize(lee_skan_skeleton)"
      ],
      "metadata": {
        "id": "KabW62fqRTal"
      },
      "execution_count": 20,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#rename columns:\n",
        "df = df.rename(columns={'image-coord-src-0': 'src-y', \n",
        "                        'image-coord-src-1': 'src-x',\n",
        "                        'image-coord-dst-0': 'dst-y',\n",
        "                        'image-coord-dst-1': 'dst-x',\n",
        "                        'mean-pixel-value': 'thickness'})\n",
        "\n",
        "#drop columns\n",
        "df.drop([#'mean-pixel-value',\n",
        "         'stdev-pixel-value',\n",
        "         'coord-src-0',\n",
        "         'coord-src-1',\n",
        "         'coord-dst-0',\n",
        "         'coord-dst-1'], axis=1, inplace=True)"
      ],
      "metadata": {
        "id": "Ytr6c5vqRU4G"
      },
      "execution_count": 21,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def get_slope(x1, y1, x2, y2):\n",
        "    if x1==x2:\n",
        "        return \"\"\n",
        "    m = (y2-y1)/(x2-x1)\n",
        "    return m\n",
        "\n",
        "#initiate blank columns\n",
        "df[\"slope\"] = 0\n",
        "df[\"tortuosity\"] = 0\n",
        "\n",
        "#compute and assign values of slope and tortuosity\n",
        "for index, row in df.iterrows():\n",
        "    df.loc[index, 'slope']=get_slope(row[\"src-x\"], row[\"src-y\"], row[\"dst-x\"], row[\"dst-y\"])\n",
        "    df.loc[index, 'tortuosity']=row[\"branch-distance\"]/row[\"euclidean-distance\"]\n",
        "\n",
        "df.to_csv(\"single_summary.csv\")"
      ],
      "metadata": {
        "id": "T0M4grM2RZN8"
      },
      "execution_count": 22,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "branches = pd.DataFrame(df['thickness'])\n",
        "branches['source_node_id']=df['node-id-src']\n",
        "branches['destination_node_id']=df['node-id-dst']\n",
        "branches.to_csv(\"branches_lee.csv\", index=False)\n",
        "branches.tail()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 205
        },
        "id": "sOUf4Lt4ReKJ",
        "outputId": "0db09860-9718-4b56-eea1-063eab00b285"
      },
      "execution_count": 23,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "     thickness  source_node_id  destination_node_id\n",
              "726   3.552071            7704                 7757\n",
              "727   2.118389            7704                 7783\n",
              "728   2.404508            7711                 7881\n",
              "729   2.033724            7757                 7880\n",
              "730   2.342973            7773                 7800"
            ],
            "text/html": [
              "\n",
              "  <div id=\"df-81701f75-861e-4d9f-82ce-71a9ac11d7a9\">\n",
              "    <div class=\"colab-df-container\">\n",
              "      <div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>thickness</th>\n",
              "      <th>source_node_id</th>\n",
              "      <th>destination_node_id</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>726</th>\n",
              "      <td>3.552071</td>\n",
              "      <td>7704</td>\n",
              "      <td>7757</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>727</th>\n",
              "      <td>2.118389</td>\n",
              "      <td>7704</td>\n",
              "      <td>7783</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>728</th>\n",
              "      <td>2.404508</td>\n",
              "      <td>7711</td>\n",
              "      <td>7881</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>729</th>\n",
              "      <td>2.033724</td>\n",
              "      <td>7757</td>\n",
              "      <td>7880</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>730</th>\n",
              "      <td>2.342973</td>\n",
              "      <td>7773</td>\n",
              "      <td>7800</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>\n",
              "      <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-81701f75-861e-4d9f-82ce-71a9ac11d7a9')\"\n",
              "              title=\"Convert this dataframe to an interactive table.\"\n",
              "              style=\"display:none;\">\n",
              "        \n",
              "  <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n",
              "       width=\"24px\">\n",
              "    <path d=\"M0 0h24v24H0V0z\" fill=\"none\"/>\n",
              "    <path d=\"M18.56 5.44l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94zm-11 1L8.5 8.5l.94-2.06 2.06-.94-2.06-.94L8.5 2.5l-.94 2.06-2.06.94zm10 10l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94z\"/><path d=\"M17.41 7.96l-1.37-1.37c-.4-.4-.92-.59-1.43-.59-.52 0-1.04.2-1.43.59L10.3 9.45l-7.72 7.72c-.78.78-.78 2.05 0 2.83L4 21.41c.39.39.9.59 1.41.59.51 0 1.02-.2 1.41-.59l7.78-7.78 2.81-2.81c.8-.78.8-2.07 0-2.86zM5.41 20L4 18.59l7.72-7.72 1.47 1.35L5.41 20z\"/>\n",
              "  </svg>\n",
              "      </button>\n",
              "      \n",
              "  <style>\n",
              "    .colab-df-container {\n",
              "      display:flex;\n",
              "      flex-wrap:wrap;\n",
              "      gap: 12px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert {\n",
              "      background-color: #E8F0FE;\n",
              "      border: none;\n",
              "      border-radius: 50%;\n",
              "      cursor: pointer;\n",
              "      display: none;\n",
              "      fill: #1967D2;\n",
              "      height: 32px;\n",
              "      padding: 0 0 0 0;\n",
              "      width: 32px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert:hover {\n",
              "      background-color: #E2EBFA;\n",
              "      box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "      fill: #174EA6;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert {\n",
              "      background-color: #3B4455;\n",
              "      fill: #D2E3FC;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert:hover {\n",
              "      background-color: #434B5C;\n",
              "      box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
              "      filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
              "      fill: #FFFFFF;\n",
              "    }\n",
              "  </style>\n",
              "\n",
              "      <script>\n",
              "        const buttonEl =\n",
              "          document.querySelector('#df-81701f75-861e-4d9f-82ce-71a9ac11d7a9 button.colab-df-convert');\n",
              "        buttonEl.style.display =\n",
              "          google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "\n",
              "        async function convertToInteractive(key) {\n",
              "          const element = document.querySelector('#df-81701f75-861e-4d9f-82ce-71a9ac11d7a9');\n",
              "          const dataTable =\n",
              "            await google.colab.kernel.invokeFunction('convertToInteractive',\n",
              "                                                     [key], {});\n",
              "          if (!dataTable) return;\n",
              "\n",
              "          const docLinkHtml = 'Like what you see? Visit the ' +\n",
              "            '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n",
              "            + ' to learn more about interactive tables.';\n",
              "          element.innerHTML = '';\n",
              "          dataTable['output_type'] = 'display_data';\n",
              "          await google.colab.output.renderOutput(dataTable, element);\n",
              "          const docLink = document.createElement('div');\n",
              "          docLink.innerHTML = docLinkHtml;\n",
              "          element.appendChild(docLink);\n",
              "        }\n",
              "      </script>\n",
              "    </div>\n",
              "  </div>\n",
              "  "
            ]
          },
          "metadata": {},
          "execution_count": 23
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "nodes_1 = pd.DataFrame(df[\"node-id-src\"])\n",
        "nodes_2=pd.DataFrame(df[\"node-id-dst\"])\n",
        "\n",
        "nodes_1 = nodes_1.rename(columns={'node-id-src': 'node_id'})\n",
        "nodes_2 = nodes_2.rename(columns={'node-id-dst': 'node_id'})\n",
        "\n",
        "nodes_1[\"node_coordinate_x\"]=df['src-x']\n",
        "nodes_1[\"node_coordinate_y\"]=df['src-y']\n",
        "\n",
        "nodes_2[\"node_coordinate_x\"]=df['dst-x']\n",
        "nodes_2[\"node_coordinate_y\"]=df['dst-y']\n",
        "\n",
        "final_nodes=nodes_1.append(nodes_2)\n",
        "\n",
        "final_nodes=final_nodes.drop_duplicates(subset=['node_id',], keep='last').reset_index()\n",
        "final_nodes = final_nodes.drop(['index'], axis=1)\n",
        "final_nodes\n",
        "\n",
        "final_nodes.to_csv(\"nodes_lee.csv\", index=False)"
      ],
      "metadata": {
        "id": "5FdH0n0OSXEe"
      },
      "execution_count": 24,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# generate_node_image (includes dst nodes)\n",
        "xs_=np.array(final_nodes['node_coordinate_x'])\n",
        "ys_=np.array(final_nodes['node_coordinate_y'])\n",
        "node_image = np.zeros([305,305], dtype=np.uint8)\n",
        "# print(node_image.shape)\n",
        "node_image[ys_.astype(np.uint16), xs_.astype(np.uint16)] = 4."
      ],
      "metadata": {
        "id": "fXTG6FAQSZ7X"
      },
      "execution_count": 25,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#compute porosity\n",
        "def compute_porosity_2d(img):\n",
        "    number_of_white_pix = np.sum(img == 1)  \n",
        "    number_of_black_pix = np.sum(img == 0)\n",
        "    porosity=number_of_black_pix/(number_of_white_pix+number_of_black_pix)\n",
        "    return(porosity)\n",
        "\n",
        "compute_porosity_2d(single_thresh)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "3izh6GSoSbrN",
        "outputId": "9e972054-7011-48d4-d894-5e89e02d639a"
      },
      "execution_count": 26,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0.5397688793335125"
            ]
          },
          "metadata": {},
          "execution_count": 26
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "plt.rcParams[\"figure.figsize\"] = (15,15)\n",
        "f, axarr = plt.subplots(2,2)\n",
        "\n",
        "axarr[0,0].imshow(single_thresh+lee94_skimage_skeleton, cmap='jet')\n",
        "axarr[0,0].set_title(\"Skeleton\")\n",
        "axarr[0,1].imshow(node_image+lee94_skimage_skeleton+single_thresh, cmap='inferno')\n",
        "axarr[0,1].set_title(\"Branch Nodes\")\n",
        "axarr[1,0].imshow(dist)\n",
        "axarr[1,0].set_title(\"Distance from nearest void\")\n",
        "axarr[1,1].imshow(dist*lee94_skimage_skeleton+single_thresh)\n",
        "axarr[1,1].set_title(\"Distance from nearest void (skeleton)\")"
      ],
      "metadata": {
        "id": "WZjoUPB-Sdlb"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#branch path coords and thicknesses -> csv\n",
        "branch_path_coords=[]\n",
        "\n",
        "branch_thics=skan_skeleton.path_means()\n",
        "# refer https://jni.github.io/skan/api/skan.csr.html\n",
        "\n",
        "for i in range(len(branch_thics)):\n",
        "    path_arr=skan_skeleton.path_coordinates(i).tolist()\n",
        "    # print(path_arr.tolist())\n",
        "    #replace newlines**** w commas\n",
        "    branch_path_coords.append(path_arr)\n",
        "    # break\n",
        "\n",
        "branch_path_coords_and_thickness=pd.DataFrame(branch_thics, columns=['thickness'])\n",
        "branch_path_coords_and_thickness['path_coordinates']=branch_path_coords\n",
        "# branch_path_coords_and_thickness['thickness']=branch_thics\n",
        "\n",
        "branch_path_coords_and_thickness\n",
        "branch_path_coords_and_thickness.to_csv(\"branch_path_coords_and_thickness.csv\", index=False)"
      ],
      "metadata": {
        "id": "bpKMP8aqSfGe"
      },
      "execution_count": 28,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Thickness, src_x, src_y, dst_x, dst_y -> csv\n",
        "\n",
        "thickness_and_coordinates = pd.DataFrame(df['thickness'])\n",
        "thickness_and_coordinates['start_x']=df['src-x']\n",
        "thickness_and_coordinates['start_y']=df['src-y']\n",
        "thickness_and_coordinates['end_x']=df['dst-x']\n",
        "thickness_and_coordinates['end_y']=df['dst-y']\n",
        "\n",
        "thickness_and_coordinates.tail()\n",
        "thickness_and_coordinates.to_csv(\"thickness_and_endpoint_coordinates.csv\", index=False )"
      ],
      "metadata": {
        "id": "9qdhtUciSk4v"
      },
      "execution_count": 29,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        ""
      ],
      "metadata": {
        "id": "6LkPRCZAVuPo"
      },
      "execution_count": 29,
      "outputs": []
    }
  ]
}