SproutAngio_main.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a11fd5c1-3bcb-4ad8-bb21-9dbbbe9f81af",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1- Basic parameters:\n",
    "image_path = 'test_data/VEGFA_10ng.czi'\n",
    "median_filtering = (7,7,7)\n",
    "second_filtering = (5,5)\n",
    "nuclei_median_filtering = (3,3,3)\n",
    "# bead radius:\n",
    "rad_ = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "32162734-7fc0-4119-94fd-05621e830d8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2- Sprout segmentation:\n",
    "# Importing necessary libraries and defining functions for the analysis:\n",
    "import numpy as np\n",
    "from czifile import CziFile\n",
    "from skimage.filters import threshold_li\n",
    "import scipy\n",
    "import matplotlib.pyplot as plt\n",
    "from skimage.measure import label   \n",
    "from scipy import ndimage\n",
    "from fil_finder import FilFinder2D\n",
    "import astropy.units as u\n",
    "from skimage.morphology import skeletonize as skeletonize_sci\n",
    "\n",
    "def getLargestCC(segmentation):\n",
    "    labels = label(segmentation)\n",
    "    assert(labels.max() != 0)\n",
    "    largestCC = labels == np.argmax(np.bincount(labels.flat)[1:])+1\n",
    "    return largestCC\n",
    "\n",
    "def modified_max_inscribed_circle(bw, f):    \n",
    "    D = ndimage.distance_transform_edt(bw)\n",
    "    Rs = -np.sort(-D, axis=None)\n",
    "    R = Rs[0]\n",
    "    RInds = np.argsort(-D, axis=None)\n",
    "    RInds = RInds[Rs >= f*R]\n",
    "    [cy, cx] = np.unravel_index(RInds, D.shape)\n",
    "    return R, cx, cy\n",
    "\n",
    "def create_circular_mask(h, w, center=None, radius=None):\n",
    "    if center is None: \n",
    "        center = (int(w/2), int(h/2))\n",
    "    if radius is None: \n",
    "        radius = min(center[0], center[1], w-center[0], h-center[1])\n",
    "\n",
    "    Y, X = np.ogrid[:h, :w]\n",
    "    dist_from_center = np.sqrt((X - center[0])**2 + (Y-center[1])**2)\n",
    "\n",
    "    mask = dist_from_center <= radius\n",
    "    return mask\n",
    "\n",
    "with CziFile(image_path) as im:\n",
    "    image_array = im.asarray()\n",
    "\n",
    "# checking the shape of the array\n",
    "print(image_array.shape)\n",
    "\n",
    "# defining the correct channels into new arrays based on the above (3) shape\n",
    "im_sprout = image_array[0,0,0,0,:,:,:,0]\n",
    "im_nuclei = image_array[0,0,0,1,:,:,:,0]\n",
    "\n",
    "# Opening napari\n",
    "import napari\n",
    "viewer = napari.view_image(im_sprout, name='im_sprout')\n",
    "\n",
    "# median filtering the sprout channel\n",
    "med_sprout = im_sprout\n",
    "med_sprout = scipy.ndimage.median_filter(med_sprout, size=median_filtering)\n",
    "\n",
    "viewer.add_image(med_sprout, name=\"median_filtered_sprout\")    \n",
    "\n",
    "# taking the maximum projection of the sprout channel and checking if median filtering size is ok\n",
    "max_proj_0 = np.max(im_sprout, axis=0)\n",
    "viewer.add_image(max_proj_0, name=\"max_projection_sprout\") \n",
    "\n",
    "max_projection = np.max(med_sprout, axis=0)\n",
    "viewer.add_image(max_projection, name=\"med_filtered_max_proj\")  \n",
    "\n",
    "# thresholding the max projection and filling the holes\n",
    "from skimage.filters import threshold_otsu\n",
    "bw_mproj = max_projection >= threshold_otsu(max_projection)\n",
    "viewer.add_image(bw_mproj, name=\"threshold_otsu\")  \n",
    "\n",
    "largest_cc = getLargestCC(bw_mproj)\n",
    "viewer.add_image(largest_cc, name=\"largest_connected_component\")  \n",
    "\n",
    "bw = ndimage.binary_fill_holes(largest_cc) * 1\n",
    "viewer.add_image(bw, name=\"fill_holes\")  \n",
    "\n",
    "# median filtering to make the sprouts smoother. \n",
    "# If the sprouts are not smooth, the analysis might find the skeleton wrong.\n",
    "med_bw = scipy.ndimage.median_filter(bw, size=second_filtering)\n",
    "viewer.add_image(med_bw, name=\"median_filtering\")\n",
    "\n",
    "skeleton_sci = skeletonize_sci(med_bw)\n",
    "# Check bw, med_bw, skeleton images to see if you need to change the size of median filtering above, before passing this part. \n",
    "viewer.add_labels(skeleton_sci, name=\"skeletonized\")\n",
    "\n",
    "# Fil-finder analysis of skeletons\n",
    "# As long as it shows the branches, it is ok to continue with possible errors.\n",
    "fil = FilFinder2D(med_bw, ang_scale=0.1*u.deg, beamwidth=1*u.arcsec, mask=med_bw)\n",
    "fil.preprocess_image(flatten_percent=90)\n",
    "fil.medskel(verbose=True)\n",
    "fil.analyze_skeletons(skel_thresh=10*u.pix, prune_criteria=\"length\")\n",
    "branch_tables = fil.branch_tables()\n",
    "\n",
    "# defining the sprout area, bead and labeling the sprouts\n",
    "bb = fil.filament_extents\n",
    "num_skl = fil.branch_properties['number'].squeeze()\n",
    "pix_skl = fil.branch_properties['pixels'][0]\n",
    "pix_skl = [x+bb[0][0] for x in pix_skl]\n",
    "\n",
    "R, cx, cy = modified_max_inscribed_circle(bw, f=0.9)\n",
    "sz = bw.shape\n",
    "bead = np.zeros(sz, dtype=bool)\n",
    "if len(cx) > 200: sp = 10\n",
    "else: sp = 20 \n",
    "for i in range(0, len(cx), sp):\n",
    "    # You can change radius by changing rad_ value in the beginning to get the correct seperation of sprouts\n",
    "    circ_ = create_circular_mask(sz[0], sz[1], center=(cx[i],cy[i]), radius=R*rad_)\n",
    "    bead = np.logical_or(bead, circ_)\n",
    "viewer.add_image(bead, name=\"bead\")\n",
    "\n",
    "bw_rmv_bead = np.logical_and(med_bw, np.logical_not(bead))\n",
    "labels = label(bw_rmv_bead)\n",
    "lbls = np.unique(labels[(fil.skeleton).astype(np.bool_)])\n",
    "lbls_inx = np.nonzero(lbls)[0]\n",
    "\n",
    "L = 1\n",
    "sprouts = np.zeros(sz, dtype=np.uint8)\n",
    "for i in lbls_inx:\n",
    "    sprouts[labels==lbls[i]] = L\n",
    "    L += 1\n",
    "viewer.add_labels(sprouts, name=\"sprouts\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "71bdd35b-41e0-44b6-9616-b88a350a6610",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3- Nuclei segmentation\n",
    "# importing the necessary libraries, median filtering the nuclei array and creating a maximum projection image\n",
    "from skimage.transform import resize\n",
    "from skimage.segmentation import morphological_chan_vese, checkerboard_level_set\n",
    "from skimage.measure import label, regionprops\n",
    "from scipy.ndimage import gaussian_filter\n",
    "from skimage.morphology import local_maxima\n",
    "import skimage.graph\n",
    "\n",
    "def shortest_path(bw, st_p, end_p):         \n",
    "    costs = np.where(bw, 1, 1000)\n",
    "    path, cost = skimage.graph.route_through_array(\n",
    "        costs, start=st_p, end=end_p, fully_connected=True)\n",
    "    return path, cost\n",
    "\n",
    "viewer.add_image(im_nuclei, name=\"nuclei_channel\")\n",
    "med_nuclei = im_nuclei\n",
    "for _ in range(2):\n",
    "    med_nuclei = scipy.ndimage.median_filter(med_nuclei, size=nuclei_median_filtering)\n",
    "viewer.add_image(med_nuclei, name=\"median_filtered_nuclei\")\n",
    "\n",
    "max_projection_n = np.max(med_nuclei, axis=0)\n",
    "viewer.add_image(max_projection_n, name=\"max_proj_nuclei\")\n",
    "\n",
    "# Nuclei segmentation continues, omitting bead part (bead2) and other beads/migrated cells (bw_not)\n",
    "bead2 = np.zeros(sz, dtype=bool)\n",
    "bw_not = np.logical_not(bw)\n",
    "\n",
    "if len(cx) > 200: sp = 10\n",
    "else: sp = 20 \n",
    "for i in range(0, len(cx), sp):\n",
    "    # You can change radius to R*n in order to get the correct seperation of sprouts\n",
    "    circ_2 = create_circular_mask(sz[0], sz[1], center=(cx[i],cy[i]), radius=R*1)\n",
    "    bead2 = np.logical_or(bead2, circ_2)\n",
    "\n",
    "\n",
    "med_nuclei2 = med_nuclei.copy()\n",
    "for i in range(med_nuclei.shape[0]):\n",
    "    med_nuclei2[i,:,:][np.where(bead2)] = 0\n",
    "    med_nuclei2[i,:,:][np.where(bw_not)] = 0\n",
    "    #np.logical_and(med_nuclei2[i,:,:], np.logical_not(bead2))\n",
    "\n",
    "#viewer.add_image(med_nuclei2, name=\"med_nuclei2\")\n",
    "\n",
    "# Downsize the image if it is large, so that the segmentation will not take so much time depending on the computing power.\n",
    "if med_nuclei.shape[1] >= 2048:\n",
    "    md = 4\n",
    "elif med_nuclei.shape[1] == 1024:\n",
    "    md = 2\n",
    "else:\n",
    "    md = 1\n",
    "med_n_ds = (med_nuclei2[:, ::md, ::md]).astype(np.float64)\n",
    "\n",
    "# Chan-Vese segmentation:\n",
    "# increasing/decreasing iterations for more precise detection of nuclei.\n",
    "# Decrease iterations if not enough nuclei segmented. For details you can check scikit image chan-vese segmentation.\n",
    "max_proj_med_nuclei2 = np.max(med_nuclei2, axis=0)\n",
    "mask_ = np.zeros(med_nuclei2.shape)\n",
    "for i in range(med_nuclei2.shape[0]):\n",
    "    mask_[i,:,:] = max_proj_med_nuclei2\n",
    "mask = (mask_[:, ::md, ::md]).astype(np.float64)\n",
    "\n",
    "bw_n = morphological_chan_vese(image=med_n_ds, num_iter=10, init_level_set=mask, lambda1=1, lambda2=1) \n",
    "\n",
    "# Histogram of objects. To select the area threshold and correct object afterwards.\n",
    "bw_nn = bw_n.astype(np.bool_)\n",
    "bw_nn = (resize(bw_nn, output_shape=(med_nuclei.shape), order=0)).astype(np.bool_)\n",
    "viewer.add_labels(bw_nn, name=\"Chan-Vese_segmentation\")\n",
    "\n",
    "# if satisfied continue. \n",
    "# It would be good if segmentation covers all the nuclei and you observe a seperation between nuclei groups.\n",
    "# We will seperate the individual nucleus soon.\n",
    "\n",
    "lbl_n = label(bw_nn, connectivity=1)\n",
    "stats = regionprops(lbl_n)\n",
    "label_n, area_n, eqds_n, maj_axis_n, box_n, centroid_n, image_n = [], [], [], [], [], [], []\n",
    "for p in stats:\n",
    "    label_n.append(p.label)\n",
    "    area_n.append(p.area)\n",
    "    eqds_n.append(p.equivalent_diameter)\n",
    "    maj_axis_n.append(p.major_axis_length)\n",
    "    box_n.append(p.bbox)\n",
    "    centroid_n.append(p.centroid)\n",
    "    image_n.append(p.image)\n",
    "\n",
    "objects_library = {\"label_n\":label_n, \"image_n\":image_n, \"area_n\":area_n,\n",
    "                   \"eqds_n\":eqds_n, \"maj_Axis_n\":maj_axis_n, \"box_n\":box_n, \"centroid_n\":centroid_n}\n",
    "\n",
    "# You can change the range of the histogram to see the details\n",
    "plt.figure(), plt.hist(area_n, bins=200, range=(0,2000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "150b6e12-9d40-441f-ba52-2dc9dbace9bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4- Lumen analysis and results\n",
    "# Nuclei location detection:\n",
    "# Based on histogram, you can change area_threshold:\n",
    "area_threshold = [30, 400]\n",
    "correct_obj = (np.array(objects_library[\"area_n\"])>area_threshold[0]) & (np.array(objects_library[\"area_n\"])<area_threshold[1])\n",
    "in_range_eqds = np.array(objects_library[\"eqds_n\"])[correct_obj]\n",
    "cents = np.array(objects_library[\"centroid_n\"])[correct_obj]\n",
    "cents_f = cents.astype(int)\n",
    "rep = np.zeros(med_nuclei.shape, dtype=np.bool_)\n",
    "for i in cents_f:    \n",
    "    rep[tuple(i)]=True\n",
    "large_obj = np.array(objects_library[\"area_n\"])>=area_threshold[1]\n",
    "intensity_thr = 0.1\n",
    "theta_gd = np.mean(in_range_eqds)\n",
    "himg = med_nuclei / np.max(med_nuclei)\n",
    "\n",
    "for i in range(len(large_obj)):\n",
    "    if large_obj[i]: \n",
    "        #print(i)           \n",
    "        bb = objects_library[\"box_n\"][i]\n",
    "        D = himg[bb[0]:bb[3], bb[1]:bb[4], bb[2]:bb[5]]        \n",
    "        B = bw_nn[bb[0]:bb[3], bb[1]:bb[4], bb[2]:bb[5]]\n",
    "        Dfil = gaussian_filter(np.pad(D, [5]), sigma=0)\n",
    "        Dfil = Dfil[5:-5, 5:-5, 5:-5]               \n",
    "        reg_max = local_maxima(Dfil, connectivity=3, indices=False, allow_borders=False)\n",
    "        reg_max = np.logical_and(reg_max, B)        \n",
    "        reg_max_idx = np.argwhere(reg_max) \n",
    "        if len(reg_max_idx):      \n",
    "            reg_max_idx = [x for x in reg_max_idx if Dfil[tuple(x)]>=intensity_thr]\n",
    "            reg_max_value = [Dfil[tuple(x)] for x in reg_max_idx] \n",
    "            ranked = np.argsort(reg_max_value)[::-1]            \n",
    "            reg_max_idx = [reg_max_idx[x] for x in ranked]\n",
    "            flags = np.ones(len(reg_max_idx), dtype=bool)\n",
    "            for j in range(len(reg_max_idx)):\n",
    "                if flags[j]:\n",
    "                    for k in range(j+1, len(reg_max_idx)):                                                                                         \n",
    "                        path_seq, cost = shortest_path(B, reg_max_idx[j], reg_max_idx[k])\n",
    "                        if cost<theta_gd:\n",
    "                            flags[k] = False    \n",
    "            reg_max_idx = np.array(reg_max_idx)[flags]\n",
    "        for x in reg_max_idx:\n",
    "            x = tuple((x[0]+bb[0], x[1]+bb[1], x[2]+bb[2]))\n",
    "            rep[tuple(x)]=True\n",
    "vis_rep = np.max(rep, axis=0)\n",
    "\n",
    "#print(rep)\n",
    "\n",
    "[rr, cc] = np.where(vis_rep)\n",
    "#plt.figure(), plt.imshow(max_projection_n)\n",
    "#plt.scatter(cc, rr, s=5, c='yellow')\n",
    "\n",
    "dots = []\n",
    "for i in range(len(rr)):\n",
    "    dots.append((rr[i],cc[i]))\n",
    "\n",
    "max_proj_nuclei = np.max(med_nuclei2, axis=0)\n",
    "viewer.add_image(max_proj_nuclei, name=\"maximum_projection_nuclei\")\n",
    "viewer.add_points(dots, name=\"nuclei\", size=5, face_color='yellow')\n",
    "\n",
    "# Lumen space analysis: Regional width and paired nuclei distance analysis\n",
    "from scipy.ndimage import distance_transform_edt as distance\n",
    "from scipy.spatial.distance import cdist\n",
    "def curve_parametrization(L):\n",
    "    # not entirely correct\n",
    "    s = L[0]\n",
    "    dist = np.sum((L-s)**2, axis=1)\n",
    "    max_inx = np.argmax(dist)\n",
    "    s = L[max_inx]\n",
    "    dist = np.sum((L-s)**2, axis=1)\n",
    "    ranked = np.argsort(dist)[::-1]\n",
    "    return L[ranked]\n",
    "def point_in_lineside (points, curve):\n",
    "    \"\"\"d=(x−x1)(y2−y1)−(y−y1)(x2−x1), where [x,y] are points to assess,\n",
    "    and [x1 y1], [x2 y2] form the line segment\n",
    "    if d>0 or d<0 the point is in either side of the line\"\"\"\n",
    "    dist_to_curve = cdist(points, curve, 'euclidean')\n",
    "    proj = np.argmin(dist_to_curve, axis=1)\n",
    "    ll = len(L)-1\n",
    "    if ll in proj: proj[proj==ll] = ll-1\n",
    "    proj_ep = proj+1\n",
    "    d1 = L[proj_ep] - L[proj]\n",
    "    d2 = rep_coord - L[proj]\n",
    "    d = (d2[:,1] * d1[:,0]) - (d2[:,0] *  d1[:,1])\n",
    "    return d\n",
    "\n",
    "skeleton = pix_skl\n",
    "skeleton = [curve_parametrization(L) for L in skeleton]\n",
    "num_sp = sprouts.max()\n",
    "num_skl = len(skeleton)\n",
    "mproj = np.logical_or(bead, sprouts)\n",
    "D = distance(mproj)\n",
    "rep_2d = np.max(rep, axis=0)\n",
    "mproj_n = np.max(med_nuclei, axis=0)\n",
    "plt.imshow(np.max(med_nuclei, axis=0))\n",
    "i_list = []\n",
    "skel_len_list = []\n",
    "tip_width_list = []\n",
    "stalk_width_list = []\n",
    "root_width_list = []\n",
    "paired_dist_list = []\n",
    "\n",
    "for i in range(1, num_sp+1):\n",
    "    sprout_name = str(i)\n",
    "    prop = {\"sp_skeleton\":[]}\n",
    "    sp = sprouts==i    \n",
    "    skl_intsc = np.zeros(num_skl)\n",
    "    for j in range(len(skeleton)):\n",
    "        L = skeleton[j]\n",
    "        intsc = np.array([sp[tuple(x)] for x in L])\n",
    "        skl_intsc[j] = intsc.sum()\n",
    "    max_inx = np.argmax(skl_intsc)\n",
    "    L = skeleton[max_inx]\n",
    "    intsc = np.array([sprouts[tuple(x)] for x in L])\n",
    "    non_z = intsc!=0\n",
    "    L = L[non_z,:]\n",
    "    if non_z[0]==False:\n",
    "        L = np.flipud(L)\n",
    "    prop = {\"sp_skeleton\": L}\n",
    "    prop[\"sp_skel_len\"] = len(L) #>>chnge\n",
    "    thick = np.array([D[tuple(x)] for x in L])\n",
    "    prop[\"thick\"] = thick\n",
    "    print(sprout_name,\"skeleton length is:\",len(L))\n",
    "    skel_len_list.append(len(L))\n",
    "    \n",
    "    if len(L)>30:\n",
    "        # You can change the section of the sprout the width measurement is made. \n",
    "        # For example tip_sp = int(0.30*len(L))... would measure the thickness \n",
    "        # a bit closer to the middle part compared to 0.25 part of the sprout \n",
    "        tip_sp = int(0.25*len(L)); tip_sp = range(tip_sp-2,tip_sp+2)\n",
    "        mid_sp = int(0.50*len(L)); mid_sp = range(mid_sp-2,mid_sp+2)\n",
    "        end_sp = int(0.75*len(L)); end_sp = range(end_sp-2,end_sp+2)\n",
    "        prop[\"tip_width\"] = np.mean(thick[tip_sp])*2\n",
    "        prop[\"stalk_width\"] = np.mean(thick[mid_sp])*2\n",
    "        prop[\"root_width\"] = np.mean(thick[end_sp])*2\n",
    "        \n",
    "        tip_width_list.append(prop[\"tip_width\"])\n",
    "        stalk_width_list.append(prop[\"stalk_width\"])\n",
    "        root_width_list.append(prop[\"root_width\"])\n",
    "        \n",
    "    else:\n",
    "        mid_sp = int(0.50*len(L)); mid_sp = range(mid_sp-1,mid_sp+1)\n",
    "        prop[\"stalk_width\"] = np.mean(thick[mid_sp])*2\n",
    "        stalk_width_list.append(prop[\"stalk_width\"])\n",
    "\n",
    "    rep_sp = np.logical_and(sp, rep_2d)\n",
    "    rep_coord = np.array(np.where(rep_sp)).T\n",
    "    if len(rep_coord)==0:\n",
    "        prop[\"paired_n\"] = []\n",
    "        prop[\"paired_dist\"] = []\n",
    "        continue\n",
    "    d = point_in_lineside (rep_coord, L)\n",
    "    rep1 = rep_coord[d>0]\n",
    "    rep2 = rep_coord[d<0]\n",
    "    if len(rep1)>len(rep2):\n",
    "        rep1 = rep_coord[d<0]\n",
    "        rep2 = rep_coord[d>0]\n",
    "    dist_to_curve = cdist(rep1, L, 'euclidean')\n",
    "    proj1 = L[np.argmin(dist_to_curve, axis=1)]\n",
    "    dist_to_curve = cdist(rep2, L, 'euclidean')\n",
    "    proj2 = L[np.argmin(dist_to_curve, axis=1)]\n",
    "    pair_dist = cdist(proj1, proj2, 'euclidean') \n",
    "\n",
    "    if len(pair_dist) > 0:\n",
    "        mv = np.min(pair_dist, axis=1)\n",
    "        mi = np.argmin(pair_dist, axis=1)\n",
    "        p = np.vstack((range(len(mi)), mi)).T\n",
    "        a = np.unique(mi)\n",
    "        paired_n = np.zeros((len(a), 2), dtype=int)\n",
    "        paired_dist = np.zeros(len(a))\n",
    "        plt.plot(L[:,1], L[:,0], \"r\")\n",
    "        plt.plot(rep1[:,1], rep1[:,0], '.', markersize=1, color=\"yellow\")\n",
    "        plt.plot(rep2[:,1], rep2[:,0], '.', markersize=1, color=\"green\")\n",
    "        plt.plot(proj1[:,1], proj1[:,0], '.', markersize=1, color=\"blue\")\n",
    "        plt.plot(proj2[:,1], proj2[:,0], '.', markersize=1, color=\"magenta\")\n",
    "        plt.axis('equal')\n",
    "    \n",
    "        for ii in range(len(a)):\n",
    "            bb = np.where(mi==a[ii])[0]\n",
    "            if len(bb)==1:\n",
    "                paired_n[ii,:] = p[bb, :]\n",
    "            else:\n",
    "                cc = np.argmin(mv[bb])\n",
    "                paired_n[ii,:] = p[bb[cc], :]\n",
    "            pp = paired_n[ii, :]\n",
    "            x1 = np.expand_dims((rep1[pp[0]]), axis=1).T\n",
    "            x2 = np.expand_dims((rep2[pp[1]]), axis=1).T\n",
    "        \n",
    "        # Putting a threshold value for paired distance according to width results:\n",
    "            if cdist(x1, x2, 'euclidean') < (np.max(thick)*2):\n",
    "                paired_dist[ii] = cdist(x1, x2, 'euclidean')\n",
    "            else:\n",
    "                continue\n",
    "        \n",
    "            # Putting a threshold value for paired distance according to width results:\n",
    "            if paired_dist[ii] < (np.max(thick)*2):\n",
    "                x = np.append(x1,x2, axis=0)\n",
    "                plt.plot(x[:,1], x[:,0])\n",
    "            else:\n",
    "                continue\n",
    "\n",
    "        prop[\"paired_n\"] = paired_n\n",
    "        prop[\"paired_dist\"] = paired_dist\n",
    "    \n",
    "        paired_dist_in_a_sprout=[]\n",
    "        for i in paired_dist:\n",
    "            if i > 0:\n",
    "                paired_dist_in_a_sprout.append(i)\n",
    "        print(sprout_name,\"paired_dist_in_this_sprout:\",paired_dist_in_a_sprout)\n",
    "\n",
    "        if len(paired_dist_in_a_sprout)>0:\n",
    "            paired_dist_average = sum(paired_dist_in_a_sprout)/len(paired_dist_in_a_sprout)\n",
    "            print(sprout_name,\"paired_dist_average:\",paired_dist_average)\n",
    "            paired_dist_list.append(paired_dist_average)\n",
    "        else:\n",
    "            continue\n",
    "    else:\n",
    "        continue\n",
    "    \n",
    "    fig_xyz = plt.savefig('analyzed_image', dpi = 300)  \n",
    "    i_list.append(str(i))\n",
    "\n",
    "nuclei_count = np.count_nonzero(rep)\n",
    "\n",
    "# Printing the measured results\n",
    "print(\"skeleton_length_measurements:\",skel_len_list)\n",
    "print()\n",
    "print(\"tip_width_measurements:\", tip_width_list)\n",
    "print(\"stalk_width_measurements:\", stalk_width_list)\n",
    "print(\"root_width_measurements:\", root_width_list)\n",
    "print()\n",
    "print(\"paired_dist_measurements per sprout:\", paired_dist_list)\n",
    "print()\n",
    "print(\"sprout_number:\", len(skel_len_list))\n",
    "print(\"skeleton_length_avg:\",sum(skel_len_list)/len(skel_len_list))\n",
    "print(\"nuclei number:\",nuclei_count)\n",
    "print()\n",
    "print(\"tip_width_avg\",sum(tip_width_list)/len(tip_width_list))\n",
    "print(\"stalk_width_avg\",sum(stalk_width_list)/len(stalk_width_list))\n",
    "print(\"root_width_avg\",sum(root_width_list)/len(root_width_list))\n",
    "print()\n",
    "if len(paired_dist_list)>0:\n",
    "    print(\"paired_dist average per bead:\",sum(paired_dist_list)/len(paired_dist_list))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ac496787-6647-42ad-a2d8-052e2cf804ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 5- VOLUMETRIC ANALYSIS\n",
    "med_n_sp = (med_sprout[:, ::2, ::2]).astype(np.float64)\n",
    "\n",
    "# Creating a mask which does not include the bead part only includes the sprout part for segmentation\n",
    "sprouts_2D = (resize(sprouts, output_shape=(med_n_sp.shape[1],med_n_sp.shape[2]), order=0)).astype(np.bool_)\n",
    "med_n_sp_copy = med_n_sp.copy()\n",
    "for i in range(med_n_sp_copy.shape[0]):\n",
    "    med_n_sp_copy[i,:,:] = sprouts_2D\n",
    "\n",
    "# Chan-Vese segmentation part. You can increase iterations if there are extra parts segmented\n",
    "# or decrease iterations if there are some target parts not segmented\n",
    "bw_sp = morphological_chan_vese(image=med_n_sp, num_iter=8, init_level_set=med_n_sp_copy, lambda1=1, lambda2=1)\n",
    "sprout_3D = (resize(bw_sp, output_shape=(med_sprout.shape), order=0)).astype(np.bool_)\n",
    "viewer.add_labels(sprout_3D)\n",
    "\n",
    "volume = np.count_nonzero(sprout_3D)\n",
    "print(\"volume:\",volume)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d7ddc482-2b42-41d5-8a61-9a2ff0f970c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 6- If you want the results in excel file:\n",
    "import xlsxwriter\n",
    "\n",
    "# Creating the excel file, once.\n",
    "workbook = xlsxwriter.Workbook('analysis_results.xlsx')\n",
    "worksheet = workbook.add_worksheet()\n",
    "\n",
    "# writing the data into excel sheet:\n",
    "worksheet.write('B1', str(\"total skeleton length\"))\n",
    "worksheet.write('B2', sum(skel_len_list))\n",
    "\n",
    "worksheet.write('C1', str(\"average skeleton length\"))\n",
    "worksheet.write('C2', sum(skel_len_list)/len(skel_len_list))\n",
    "\n",
    "worksheet.write('D1', str(\"sprout number\"))\n",
    "worksheet.write('D2', len(skel_len_list))\n",
    "\n",
    "worksheet.write('E1', str(\"nuclei number\"))\n",
    "worksheet.write('E2', nuclei_count)\n",
    "\n",
    "worksheet.write('F1', str(\"tip width\"))\n",
    "worksheet.write('F2', sum(tip_width_list)/len(tip_width_list))\n",
    "\n",
    "worksheet.write('G1', str(\"stalk width\"))\n",
    "worksheet.write('G2', sum(stalk_width_list)/len(stalk_width_list))\n",
    "\n",
    "worksheet.write('H1', str(\"root width\"))\n",
    "worksheet.write('H2', sum(root_width_list)/len(root_width_list))\n",
    "\n",
    "worksheet.write('I1', str(\"sprout volume\"))\n",
    "worksheet.write('I2', volume)\n",
    "\n",
    "worksheet.write('J1', str(\"average paired nuclei distance\"))\n",
    "if len(paired_dist_list)>0:\n",
    "    worksheet.write('J2', sum(paired_dist_list)/len(paired_dist_list))\n",
    "\n",
    "workbook.close()\n",
    "print(\"analysis is done\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05bce4e0-b8cc-4cfe-898f-fcf291e0779f",
   "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.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}