3D_analysis.ipynb

{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "eJLAo-lZWV98"
      },
      "outputs": [],
      "source": [
        "!pip install SimpleITK\n",
        "!pip install skan \n",
        "!pip install tifffile\n",
        "!pip install imjoy"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Nniw7y-2WedD"
      },
      "outputs": [],
      "source": [
        "from google.colab import drive\n",
        "drive.mount('/content/drive')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "e8aZmFuxWcg2"
      },
      "outputs": [],
      "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",
        "from imjoy import api\n",
        "import asyncio\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\n",
        "import shutil\n",
        "from google.colab import files"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "mgUX7gmJXGlp"
      },
      "outputs": [],
      "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)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "EdihqFCbqHSz"
      },
      "outputs": [],
      "source": [
        "folder_id=\"Perforated Model 1 Threshold\"\n",
        "!cp -r \"/content/drive/MyDrive/echinoids-official-data/Perforato/Perforato Modello 1/Perforated Model 1 Threshold\" \"/content/\"\n",
        "\n",
        "## for zip files\n",
        "# filename=\"Galleried suture 1.zip\"\n",
        "# src=\"/content/drive/MyDrive/mydata/sea_urchin_data/\"+filename\n",
        "# dst=\"/content/\"+filename\n",
        "# shutil.copyfile(src, dst)\n",
        "# shutil.unpack_archive(\"/content/\"+filename, \"/content/\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "SUqb8KHtXMPq"
      },
      "outputs": [],
      "source": [
        "def tiffs_to_3d_numpy(folderpath):\n",
        "    '''\n",
        "    Converts tiff files to numpy in bulk, also sorts based on filename\n",
        "\n",
        "    args:\n",
        "        folderpath (str): path to folder containing tiff files\n",
        "\n",
        "    Returns:\n",
        "        np.array: array containing data from all tiff files in folderpath\n",
        "    '''\n",
        "\n",
        "    filelist = glob.glob(folderpath+'*.tiff') \n",
        "    filename_prefix = filelist[0][:-8]\n",
        "    filename_suffix = filelist[0][-5:]\n",
        "    number_ids=[]\n",
        "\n",
        "    for i in range(len(filelist)):\n",
        "        number_ids.append(int(filelist[i][-8:-5]))\n",
        "        filelist[i] = filelist[i][-8:-5]\n",
        "\n",
        "    filelist.sort()\n",
        "    sorted_filelist=[]\n",
        "\n",
        "    for j in range(len(filelist)):\n",
        "        sorted_filelist.append(filename_prefix+filelist[j]+filename_suffix)\n",
        "\n",
        "    val_seg = np.array([plt.imread(fname) for fname in sorted_filelist])\n",
        "\n",
        "    return(val_seg)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "mjxgcnx_XQe-"
      },
      "outputs": [],
      "source": [
        "# show segmentation data\n",
        "segmentation_data_np = tiffs_to_3d_numpy(\"/content/\"+folder_id+\"/\")\n",
        "images_to_show=np.array([segmentation_data_np])\n",
        "fig = px.imshow(images_to_show, \n",
        "                facet_col=0,\n",
        "                animation_frame=1,\n",
        "                color_continuous_scale =\"gray\")\n",
        "fig.show()"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "def compute_porosity(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",
        "porousity=compute_porosity(segmentation_data_np)\n",
        "print(\"porousity = \"+str(porousity*100)+\"%\")"
      ],
      "metadata": {
        "id": "NDll4iCVzrRr"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "p0dLhtueZJOx"
      },
      "outputs": [],
      "source": [
        "#compute 3d distance transfrom\n",
        "dist_trans_3d=ndimage.distance_transform_edt(segmentation_data_np, return_distances=True)\n",
        "\n",
        "#skeletonize\n",
        "lee_3d_skimage_skeleton=skeletonize(segmentation_data_np, method='lee')\n",
        "lee_3d_skan_skeleton=Skeleton(lee_3d_skimage_skeleton*dist_trans_3d)#*dist_trans_3d)\n",
        "df_3d=summarize(lee_3d_skan_skeleton)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "-XFRt0O_Zck1"
      },
      "outputs": [],
      "source": [
        "df_3d = df_3d.rename(columns={'image-coord-src-0': 'src-z', \n",
        "                        'image-coord-src-1': 'src-y',\n",
        "                        'image-coord-src-2': 'src-x',\n",
        "                        'image-coord-dst-0': 'dst-z',\n",
        "                        'image-coord-dst-1': 'dst-y',\n",
        "                        'image-coord-dst-2': 'dst-x',\n",
        "                        'mean-pixel-value': 'thickness'})\n",
        "\n",
        "#drop columns\n",
        "df_3d.drop([#'mean-pixel-value',\n",
        "         'stdev-pixel-value',\n",
        "         'coord-src-0',\n",
        "         'coord-src-1',\n",
        "         'coord-src-2',\n",
        "         'coord-dst-0',\n",
        "         'coord-dst-2',\n",
        "         'coord-dst-1'], axis=1, inplace=True)"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# path_coordinates=lee_3d_skan_skeleton.path_coordinates(1305)\n",
        "# path_coordinates\n",
        "dummy_=lee_3d_skan_skeleton.path_means()\n",
        "branch_path_coords=[]\n",
        "for i in tnrange(len(dummy_)):\n",
        "    path_arr=lee_3d_skan_skeleton.path_coordinates(i).tolist()\n",
        "    # print(path_arr.tolist())\n",
        "    #replace newlines**** w commas\n",
        "    branch_path_coords.append(path_arr)\n",
        "\n",
        "df_3d['path-coords']=branch_path_coords\n",
        "df_3d\n"
      ],
      "metadata": {
        "id": "2J0ZIcP1U4Py"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "o1bbRHBzZffR"
      },
      "outputs": [],
      "source": [
        "# branch avg thickness and endpoint node ids as csv\n",
        "segments = pd.DataFrame(df_3d['thickness'])\n",
        "segments['source_node_id']=df_3d['node-id-src']\n",
        "segments['destination_node_id']=df_3d['node-id-dst']\n",
        "segments['branch-distance']=df_3d['branch-distance']\n",
        "segments['euclidean-distance']=df_3d['euclidean-distance']\n",
        "segments['path-coords']=df_3d['path-coords']\n",
        "\n",
        "#compute tortuosity\n",
        "segments[\"tortuosity\"] = 0\n",
        "\n",
        "#compute and assign values of slope and tortuosity\n",
        "for index, row in segments.iterrows():\n",
        "    if row['path-coords'][-1] == row['path-coords'][0]: #set tortuosity as -1 for loop branch\n",
        "        segments.loc[index, 'tortuosity']= -1.\n",
        "    else:\n",
        "        segments.loc[index, 'tortuosity']=row[\"branch-distance\"]/row[\"euclidean-distance\"]\n",
        "    \n",
        "segments.to_csv(\"segments.csv\", index=False)\n",
        "segments.tail()"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "def get_coordination_num(img3d, coordinates):\n",
        "    x,y,z = coordinates\n",
        "    padded_img3d=np.pad(img3d, pad_width=1)\n",
        "    if img3d[z, y, x]==0:\n",
        "        print(coordinates, \"is not in the skeleton\")\n",
        "    else:\n",
        "        small_cube = padded_img3d[z:z+3, y:y+3, x:x+3]\n",
        "        number_of_white_pix = np.sum(small_cube == 1)  \n",
        "        # print(img3d[z, y, x])\n",
        "        return(number_of_white_pix-1)\n",
        "\n",
        "get_coordination_num(lee_3d_skimage_skeleton, (52, 98, 0))"
      ],
      "metadata": {
        "id": "4g8-pXGU7-uE"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "UhJ3TLPdZgSP"
      },
      "outputs": [],
      "source": [
        "# csv with every branch endpoints\n",
        "nodes_1_3d = pd.DataFrame(df_3d[\"node-id-src\"])\n",
        "nodes_2_3d=pd.DataFrame(df_3d[\"node-id-dst\"])\n",
        "\n",
        "nodes_1_3d = nodes_1_3d.rename(columns={'node-id-src': 'node_id'})\n",
        "nodes_2_3d = nodes_2_3d.rename(columns={'node-id-dst': 'node_id'})\n",
        "\n",
        "nodes_1_3d[\"node_coordinate_x\"]=df_3d['src-x']\n",
        "nodes_1_3d[\"node_coordinate_y\"]=df_3d['src-y']\n",
        "nodes_1_3d[\"node_coordinate_z\"]=df_3d['src-z']\n",
        "\n",
        "nodes_2_3d[\"node_coordinate_x\"]=df_3d['dst-x']\n",
        "nodes_2_3d[\"node_coordinate_y\"]=df_3d['dst-y']\n",
        "nodes_2_3d[\"node_coordinate_z\"]=df_3d['dst-z']\n",
        "\n",
        "final_nodes_3d=nodes_1_3d.append(nodes_2_3d)\n",
        "\n",
        "final_nodes_3d=final_nodes_3d.drop_duplicates(subset=['node_id',], keep='last').reset_index()\n",
        "final_nodes_3d = final_nodes_3d.drop(['index'], axis=1)\n",
        "final_nodes_3d"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "coordination_numbers=np.zeros(len(final_nodes_3d))\n",
        "for i in range(len(final_nodes_3d)):\n",
        "    x = int(final_nodes_3d['node_coordinate_x'][i])\n",
        "    y = int(final_nodes_3d['node_coordinate_y'][i])\n",
        "    z = int(final_nodes_3d['node_coordinate_z'][i])\n",
        "    coordination_numbers[i] = get_coordination_num(lee_3d_skimage_skeleton, (x, y, z))\n",
        "final_nodes_3d['coordination-number']=coordination_numbers"
      ],
      "metadata": {
        "id": "v9NsfBxO0lKZ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "final_nodes_3d.to_csv(folder_id+\" nodes.csv\", index=False)\n",
        "final_nodes_3d.head()\n",
        "files.download(folder_id+\" nodes.csv\") "
      ],
      "metadata": {
        "id": "cybthkSc0e4x"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "FRAtOFfwZkZ3"
      },
      "outputs": [],
      "source": [
        "# extract close nodes\n",
        "close_node_pair_coordinates=[]\n",
        "\n",
        "x1s_pre=np.array(final_nodes_3d['node_coordinate_x'])\n",
        "y1s_pre=np.array(final_nodes_3d['node_coordinate_y'])\n",
        "z1s_pre=np.array(final_nodes_3d['node_coordinate_z'])\n",
        "\n",
        "x2s_pre=x1s_pre\n",
        "y2s_pre=y1s_pre\n",
        "z2s_pre=z1s_pre\n",
        "\n",
        "for i in tnrange(len(final_nodes_3d)):\n",
        "    for j in range(len(final_nodes_3d)):\n",
        "        x1, y1, z1, x2, y2, z2 = (x1s_pre[i],\n",
        "                                  y1s_pre[i],\n",
        "                                  z1s_pre[i],\n",
        "                                  x2s_pre[j],\n",
        "                                  y2s_pre[j],\n",
        "                                  z2s_pre[j])\n",
        "\n",
        "        # print(x1, y1, z1, x2, y2, z2)\n",
        "\n",
        "        p1 = np.array([x1, y1, z1])\n",
        "        p2 = np.array([x2, y2, z2])\n",
        "\n",
        "        #compute distance between nodes\n",
        "        squared_dist = np.sum((p1-p2)**2, axis=0)\n",
        "        dist = np.sqrt(squared_dist)\n",
        "\n",
        "        if dist<2.0 and dist>0.0:\n",
        "            thickness_2=dist_trans_3d[int(z2)][int(y2)][int(x2)]\n",
        "            thickness_1=dist_trans_3d[int(z1)][int(y1)][int(x1)]\n",
        "            avg_thickness=np.mean([thickness_1, thickness_2])\n",
        "            close_node_pair_coordinates.append([(x1, y1, z1), (x2, y2, z2), avg_thickness])\n",
        "            #close_node_pair_coordinates -> [[(src coords), (dst coords), dst], [(),(),t], ..]\n",
        "\n",
        "            if thickness_1==0 or thickness_2==0:\n",
        "                print(\"check\", thickness_1, thickness_2)\n",
        "                print(\"coords \", x1, y1, z1, x2, y2, z2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "OImv4RGQZmUi"
      },
      "outputs": [],
      "source": [
        "# close nodes as pandas df\n",
        "fine_branches=pd.DataFrame(close_node_pair_coordinates, columns=[\"src\", \"dst\", \"thickness\"])\n",
        "\n",
        "fine_branches[['x1', 'y1', 'z1']] = pd.DataFrame(fine_branches['src'].tolist(), \n",
        "                                                 index=fine_branches.index)\n",
        "\n",
        "fine_branches[['x2', 'y2', 'z2']] = pd.DataFrame(fine_branches['dst'].tolist(), \n",
        "                                                 index=fine_branches.index)\n",
        "\n",
        "fine_branches.drop('src', axis=1, inplace=True)\n",
        "fine_branches.drop('dst', axis=1, inplace=True)\n",
        "\n",
        "fine_branches\n",
        "#now add intermediate nodes from branches"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Dx5MYWfxZphP"
      },
      "outputs": [],
      "source": [
        "# treat each pair of adjacent skeleton point as a branch\n",
        "dummy_=lee_3d_skan_skeleton.path_means()\n",
        "fine_branches_2=fine_branches\n",
        "for i in tnrange(len(dummy_)):\n",
        "    path_coordinates=lee_3d_skan_skeleton.path_coordinates(i)\n",
        "    # print(path_coordinates)\n",
        "\n",
        "    for c in range(len(path_coordinates)-1):\n",
        "        x1, y1, z1 = path_coordinates[c]\n",
        "        x2, y2, z2 = path_coordinates[c+1]\n",
        "        # print(\"src\",x1, y1, z1, \"dst\",x2, y2, z2)\n",
        "\n",
        "        thickness_2=dist_trans_3d[int(x2)][int(y2)][int(z2)]\n",
        "        thickness_1=dist_trans_3d[int(x1)][int(y1)][int(z1)]\n",
        "        avg_thickness=np.mean([thickness_1, thickness_2])\n",
        "\n",
        "        temp_dict={'thickness':avg_thickness,\n",
        "                   'x1':x1,\n",
        "                   'y1':y1,\n",
        "                   'z1':z1,\n",
        "                   'x2':x2,\n",
        "                   'y2':y2,\n",
        "                   'z2':z2}\n",
        "\n",
        "        fine_branches_2=fine_branches_2.append(temp_dict, ignore_index=True)\n",
        "\n",
        "fine_branches_2.to_csv(folder_id+\" finer branches.csv\", index=False)"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "files.download(folder_id+\" finer branches.csv\") "
      ],
      "metadata": {
        "id": "cHzzAQ4HX6rW"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "MIiE45bVZrge"
      },
      "outputs": [],
      "source": [
        "#compute 3d node image\n",
        "xs_3d=np.array(final_nodes_3d['node_coordinate_x'])\n",
        "ys_3d=np.array(final_nodes_3d['node_coordinate_y'])\n",
        "zs_3d=np.array(final_nodes_3d['node_coordinate_z'])\n",
        "node_image_3d = np.zeros([100,100,100], dtype=np.uint8)\n",
        "node_image_3d[zs_3d.astype(np.uint16), ys_3d.astype(np.uint16), xs_3d.astype(np.uint16)] = 4."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "c3GbAU9bZUgX"
      },
      "outputs": [],
      "source": [
        "images_to_show=np.array([dist_trans_3d, segmentation_data_np])\n",
        "fig = px.imshow(images_to_show, \n",
        "                facet_col=0,\n",
        "                animation_frame=1)\n",
        "# for i, label in enumerate(['2d distance transform', '3d distance transform', \"Difference\"]):\n",
        "#     fig.layout.annotations[i]['text'] = label\n",
        "fig.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "lxus08XUamKm"
      },
      "outputs": [],
      "source": [
        "# show thresh + nodes + skel as 2d slices\n",
        "images_to_show=np.array([lee_3d_skimage_skeleton+node_image_3d+segmentation_data_np, dist_trans_3d])\n",
        "fig = px.imshow(images_to_show, \n",
        "                facet_col=0,\n",
        "                animation_frame=1,\n",
        "                color_continuous_scale =\"inferno\")\n",
        "fig.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "UUXrFLTyamIP"
      },
      "outputs": [],
      "source": [
        "# matplotlib static 2d projection of 3d skel\n",
        "plt.rcParams[\"figure.figsize\"] = (10,10)\n",
        "fig = plt.figure()\n",
        "ax = fig.add_subplot(111, projection='3d')\n",
        "\n",
        "size = 21\n",
        "m = lee_3d_skimage_skeleton[:100]\n",
        "\n",
        "pos = np.where(m==1)\n",
        "ax.scatter(pos[0], pos[1], pos[2], c='green', alpha=0.3)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "-sp0nX98amFf"
      },
      "outputs": [],
      "source": [
        "# 3d skel as plotly scatter\n",
        "fig = go.Figure(data=[go.Scatter3d(\n",
        "    x=pos[0],\n",
        "    y=pos[1],\n",
        "    z=pos[2],\n",
        "    mode='markers',\n",
        "    marker=dict(\n",
        "        size=2,\n",
        "        color=np.ones_like(pos[0]),         \n",
        "        colorscale='Viridis',  \n",
        "        opacity=0.6\n",
        "    )\n",
        ")])\n",
        "\n",
        "# tight layout\n",
        "fig.update_layout(margin=dict(l=0, r=0, b=0, t=0))\n",
        "fig.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "gAiz7--tjZPf"
      },
      "outputs": [],
      "source": [
        "class ImJoyPlugin():\n",
        "    async def setup(self):\n",
        "        pass\n",
        "\n",
        "    async def run(self, ctx):\n",
        "        viewer = await api.createWindow(src=\"https://kaibu.org/#/app\")\n",
        "        volume = m\n",
        "        # view image\n",
        "        await viewer.view_image(volume, type=\"itk-vtk\", name=\"3D volume\")\n",
        "\n",
        "api.export(ImJoyPlugin())"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        ""
      ],
      "metadata": {
        "id": "KLK-cuGeuM6q"
      },
      "execution_count": null,
      "outputs": []
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [],
      "name": "alpha-3D-analysis.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}