NIH_Project_Testing_count_AI.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "645bbb21-2d17-43e6-acd7-fcda0a7f3314",
   "metadata": {},
   "outputs": [],
   "source": [
    "###SQUID PY###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e553503a-dda3-48ee-8bc4-73715ec66e31",
   "metadata": {},
   "outputs": [],
   "source": [
    "pip install git+https://github.com/scverse/squidpy@main"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ca02444-0a1a-4309-b9dc-5982b984f1cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "import squidpy as sq\n",
    "import cv2 \n",
    "\n",
    "import xarray as xr\n",
    "import squidpy as sq\n",
    "\n",
    "# # load the H&E stained tissue image and crop to a smaller segment\n",
    "# img = sq.datasets.visium_hne_image_crop()\n",
    "# crop = img.crop_corner(0, 0, size=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "3931ebcd-685c-4a99-b3de-9ad92246af37",
   "metadata": {},
   "outputs": [],
   "source": [
    "img = cv2.imread(r\"C:\\Users\\cheukjy\\WC25_LLEDGE2_test\\Test_2_H_tiles\\267_29_H_test_2.tif\")\n",
    "img = cv2.cvtColor(img,cv2.COLOR_BGR2RGB)\n",
    "img2 = cv2.cvtColor(img,cv2.COLOR_RGB2GRAY)\n",
    "\n",
    "cv2.imshow(\"test\",img2)\n",
    "cv2.waitKey(0)\n",
    "cv2.destroyAllWindows() \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b86ffe58-9f2d-466f-9a38-f51a4727adb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n",
    "crop.show(\"image_smooth\", cmap=\"gray\", ax=axes[0])\n",
    "axes[1].imshow(crop[\"image_smooth\"][:, :, 0, 0] < 90)\n",
    "_ = sns.histplot(np.array(crop[\"image_smooth\"]).flatten(), bins=50, ax=axes[2])\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94f92b5b-d15a-4bcb-93ca-f8ef67c1f663",
   "metadata": {},
   "outputs": [],
   "source": [
    "sq.im.segment(img=crop, layer=\"image_smooth\", method=\"watershed\", thresh=90, geq=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f134c7df-e1eb-4033-b54c-f664ec7351ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(crop)\n",
    "print(f\"Number of segments in crop: {len(np.unique(crop['segmented_watershed']))}\")\n",
    "\n",
    "fig, axes = plt.subplots(1, 2)\n",
    "crop.show(\"image\", channel=0, ax=axes[0])\n",
    "_ = axes[0].set_title(\"H&E\")\n",
    "crop.show(\"segmented_watershed\", cmap=\"jet\", interpolation=\"none\", ax=axes[1])\n",
    "_ = axes[1].set_title(\"segmentation\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8527542b-198a-48f1-8bb1-205e1880b699",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b71c8380-9f45-493b-b8d5-9695d8b9d229",
   "metadata": {},
   "outputs": [],
   "source": [
    "###TESTING MYSELF####"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b472590a-546f-4aef-80e7-4e6982fa0961",
   "metadata": {},
   "outputs": [],
   "source": [
    "# libraries below are from 16\n",
    "import cv2 \n",
    "from skimage.filters import sobel\n",
    "from skimage import io \n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from skimage import img_as_float\n",
    "# Libraries below are from 35\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy import ndimage\n",
    "from skimage import measure, color, io\n",
    "from skimage.segmentation import clear_border\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ce8f8802-5f20-4f68-aa3b-d7c6b19710b1",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\anaconda3\\Lib\\site-packages\\paramiko\\pkey.py:82: CryptographyDeprecationWarning: TripleDES has been moved to cryptography.hazmat.decrepit.ciphers.algorithms.TripleDES and will be removed from this module in 48.0.0.\n",
      "  \"cipher\": algorithms.TripleDES,\n",
      "C:\\ProgramData\\anaconda3\\Lib\\site-packages\\paramiko\\transport.py:219: CryptographyDeprecationWarning: Blowfish has been moved to cryptography.hazmat.decrepit.ciphers.algorithms.Blowfish and will be removed from this module in 45.0.0.\n",
      "  \"class\": algorithms.Blowfish,\n",
      "C:\\ProgramData\\anaconda3\\Lib\\site-packages\\paramiko\\transport.py:243: CryptographyDeprecationWarning: TripleDES has been moved to cryptography.hazmat.decrepit.ciphers.algorithms.TripleDES and will be removed from this module in 48.0.0.\n",
      "  \"class\": algorithms.TripleDES,\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "missing required 'shape' or 'dtype' argument",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m~\\AppData\\Local\\Temp\\ipykernel_17140\\826036945.py\u001b[0m in \u001b[0;36m?\u001b[1;34m()\u001b[0m\n\u001b[0;32m     75\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     76\u001b[0m \u001b[0moutput_file\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mclose\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     77\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     78\u001b[0m \u001b[1;31m# Display the results\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 79\u001b[1;33m \u001b[0mtiff\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mimwrite\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mimg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     80\u001b[0m \u001b[0mtiff\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mimwrite\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mgray\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     81\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfigure\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfigsize\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m10\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m10\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     82\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msubplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m121\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mimshow\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mimg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\anaconda3\\Lib\\site-packages\\tifffile\\tifffile.py\u001b[0m in \u001b[0;36m?\u001b[1;34m(file, data, bigtiff, byteorder, imagej, ome, shaped, append, shape, dtype, photometric, planarconfig, extrasamples, volumetric, tile, rowsperstrip, bitspersample, compression, compressionargs, predictor, subsampling, jpegtables, colormap, description, datetime, resolution, resolutionunit, subfiletype, software, metadata, extratags, contiguous, truncate, align, maxworkers, returnoffset)\u001b[0m\n\u001b[0;32m   1215\u001b[0m     \"\"\"\n\u001b[0;32m   1216\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0mdata\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1217\u001b[0m         \u001b[1;31m# write empty file\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1218\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mshape\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mNone\u001b[0m \u001b[1;32mor\u001b[0m \u001b[0mdtype\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1219\u001b[1;33m             \u001b[1;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"missing required 'shape' or 'dtype' argument\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1220\u001b[0m         \u001b[0mdtype\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnumpy\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdtype\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdtype\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1221\u001b[0m         \u001b[0mshape\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtuple\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1222\u001b[0m         \u001b[0mdatasize\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mproduct\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mitemsize\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mValueError\u001b[0m: missing required 'shape' or 'dtype' argument"
     ]
    }
   ],
   "source": [
    "import cv2\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy import ndimage\n",
    "from skimage import measure, color, io\n",
    "from skimage.segmentation import clear_border\n",
    "import tifffile as tiff\n",
    "\n",
    "# Read the image\n",
    "img = cv2.imread(r\"C:\\Users\\cheukjy\\WC25_LLEDGE2_test\\Test_2_H_tiles\\267_29_H_test_2.tif\")\n",
    "pixels_to_um = 0.100\n",
    "# Convert the image to grayscale\n",
    "gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n",
    "\n",
    "# Threshold the image to binary\n",
    "ret1, thresh = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY+cv2.THRESH_OTSU)\n",
    "\n",
    "# Morphological operations to remove small noise - opening\n",
    "kernel = np.ones((3,3),np.uint8)\n",
    "opening = cv2.morphologyEx(thresh, cv2.MORPH_OPEN, kernel, iterations = 2)\n",
    "opening = clear_border(opening)  # Remove edge-touching grains\n",
    "\n",
    "# Dilate to obtain sure background\n",
    "sure_bg = cv2.dilate(opening, kernel, iterations=10)\n",
    "\n",
    "# Distance transform\n",
    "dist_transform = cv2.distanceTransform(opening, cv2.DIST_L2, 5)\n",
    "\n",
    "# Threshold the distance transform to get sure foreground\n",
    "ret2, sure_fg = cv2.threshold(dist_transform, 0.5*dist_transform.max(), 255, 0)\n",
    "\n",
    "# Convert to uint8 and subtract to get unknown regions\n",
    "sure_fg = np.uint8(sure_fg)\n",
    "unknown = cv2.subtract(sure_bg, sure_fg)\n",
    "\n",
    "# Marker labeling\n",
    "ret3, markers = cv2.connectedComponents(sure_fg)\n",
    "\n",
    "# Add 10 to all labels to ensure the background is not zero\n",
    "markers = markers + 10\n",
    "\n",
    "# Mark the unknown region with zero\n",
    "markers[unknown == 255] = 0\n",
    "\n",
    "# Apply watershed\n",
    "markers = cv2.watershed(img, markers)\n",
    "img[markers == -1] = [0,255,255]  # Mark boundaries in yellow\n",
    "\n",
    "# Convert the labeled image to RGB\n",
    "img2 = color.label2rgb(markers, bg_label=0)\n",
    "\n",
    "# Analyze the regions\n",
    "regions = measure.regionprops(markers, intensity_image=gray)\n",
    "\n",
    "# Output properties to a CSV file\n",
    "propList = ['Area', 'equivalent_diameter', 'orientation', 'MajorAxisLength',\n",
    "            'MinorAxisLength', 'Perimeter', 'MinIntensity', 'MeanIntensity', 'MaxIntensity']\n",
    "\n",
    "output_file = open('cell_measurements.csv', 'w')\n",
    "output_file.write(',' + \",\".join(propList) + '\\n')\n",
    "\n",
    "for region_props in regions:\n",
    "    output_file.write(str(region_props.label))\n",
    "    for prop in propList:\n",
    "        if prop == 'Area':\n",
    "            to_print = region_props[prop] * pixels_to_um**2\n",
    "        elif prop == 'orientation':\n",
    "            to_print = region_props[prop] * 57.2958  # Convert to degrees\n",
    "        elif 'Intensity' not in prop:\n",
    "            to_print = region_props[prop] * pixels_to_um\n",
    "        else:\n",
    "            to_print = region_props[prop]\n",
    "        output_file.write(',' + str(to_print))\n",
    "    output_file.write('\\n')\n",
    "\n",
    "output_file.close()\n",
    "\n",
    "# Display the results\n",
    "tiff.imwrite(img)\n",
    "tiff.imwrite(gray)\n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.subplot(121), plt.imshow(img)\n",
    "plt.title('Overlay on original image'), plt.xticks([]), plt.yticks([])\n",
    "plt.subplot(122), plt.imshow(img2)\n",
    "plt.title('Colored Grains'), plt.xticks([]), plt.yticks([])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc26605a-7b3b-4d57-ad0a-41a17955e896",
   "metadata": {},
   "outputs": [],
   "source": [
    "img = cv2.imread(r\"C:\\Users\\cheukjy\\WC25_LLEDGE2_test\\Test_2_H_tiles\\267_29_H_test_2.tif\")\n",
    "\n",
    "gray = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)\n",
    "\n",
    "pixels_to_um = 0.100 # 1 pixel = 454 nm (got this from the metadata of original image)\n",
    "\n",
    "#Threshold image to binary using OTSU. ALl thresholded pixels will be set to 255\n",
    "# ret1, thresh = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY+cv2.THRESH_OTSU)\n",
    "\n",
    "plt.imshow(gray)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7d9b359e-6dc5-4636-8f0e-1cce3da7f9d1",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'np' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[1], line 3\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[38;5;66;03m# Morphological operations to remove small noise - opening\u001b[39;00m\n\u001b[0;32m      2\u001b[0m \u001b[38;5;66;03m#To remove holes we can use closing\u001b[39;00m\n\u001b[1;32m----> 3\u001b[0m kernel \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mones((\u001b[38;5;241m3\u001b[39m,\u001b[38;5;241m3\u001b[39m),np\u001b[38;5;241m.\u001b[39muint8)\n\u001b[0;32m      4\u001b[0m opening \u001b[38;5;241m=\u001b[39m cv2\u001b[38;5;241m.\u001b[39mmorphologyEx(thresh,cv2\u001b[38;5;241m.\u001b[39mMORPH_OPEN,kernel, iterations \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m2\u001b[39m)\n\u001b[0;32m      5\u001b[0m opening \u001b[38;5;241m=\u001b[39m clear_border(opening) \u001b[38;5;66;03m# to remove the edge-touching grains\u001b[39;00m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'np' is not defined"
     ]
    }
   ],
   "source": [
    "# Morphological operations to remove small noise - opening\n",
    "#To remove holes we can use closing\n",
    "kernel = np.ones((3,3),np.uint8)\n",
    "opening = cv2.morphologyEx(thresh,cv2.MORPH_OPEN,kernel, iterations = 2)\n",
    "opening = clear_border(opening) # to remove the edge-touching grains\n",
    "\n",
    "\n",
    "sure_bg = cv2.dilate(opening,kernel,iterations=10)\n",
    "\n",
    "\n",
    "dist_transform = cv2.distanceTransform(opening,cv2.DIST_L2,5)\n",
    "\n",
    "\n",
    "#Let us threshold the dist transform by starting at 1/2 its max value.\n",
    "#print(dist_transform.max()) gives about 21.9\n",
    "ret2, sure_fg = cv2.threshold(dist_transform,0.5*dist_transform.max(),255,0)\n",
    "\n",
    "sure_fg = np.uint8(sure_fg)\n",
    "unknown = cv2.subtract(sure_bg,sure_fg)\n",
    "\n",
    "#Now we create a marker and label the regions inside. \n",
    "\n",
    "#For markers let us use ConnectedComponents. \n",
    "ret3, markers = cv2.connectedComponents(sure_fg)\n",
    "\n",
    "\n",
    "markers = markers+10\n",
    "\n",
    "# Now, mark the region of unknown with zero\n",
    "markers[unknown==255] = 0\n",
    "\n",
    "markers = cv2.watershed(img,markers)\n",
    "\n",
    "#Let us color boundaries in yellow. \n",
    "img[markers == -1] = [0,255,255]  \n",
    "\n",
    "img2 = color.label2rgb(markers, bg_label=0)\n",
    "\n",
    "\n",
    "regions = measure.regionprops(markers, intensity_image=gray)\n",
    "\n",
    "#Can print various parameters for all objects\n",
    "#Best way is to output all properties to a csv file\n",
    "#Let us pick which ones we want to export. \n",
    "\n",
    "propList = ['Area',\n",
    "            'equivalent_diameter', #Added... verify if it works\n",
    "            'orientation', #Added, verify if it works. Angle btwn x-axis and major axis.\n",
    "            'MajorAxisLength',\n",
    "            'MinorAxisLength',\n",
    "            'Perimeter',\n",
    "            'MinIntensity',\n",
    "            'MeanIntensity',\n",
    "            'MaxIntensity']    \n",
    "    \n",
    "\n",
    "output_file = open('cell_measurements.csv', 'w')\n",
    "output_file.write(',' + \",\".join(propList) + '\\n') #join strings in array by commas, leave first cell blank\n",
    "#First cell blank to leave room for header (column names)\n",
    "\n",
    "for region_props in regions:\n",
    "    #output cluster properties to the excel file\n",
    "    output_file.write(str(region_props['Label']))\n",
    "    for i,prop in enumerate(propList):\n",
    "        if(prop == 'Area'): \n",
    "            to_print = region_props[prop]*pixels_to_um**2   #Convert pixel square to um square\n",
    "        elif(prop == 'orientation'): \n",
    "            to_print = region_props[prop]*57.2958  #Convert to degrees from radians\n",
    "        elif(prop.find('Intensity') < 0):          # Any prop without Intensity in its name\n",
    "            to_print = region_props[prop]*pixels_to_um\n",
    "        else: \n",
    "            to_print = region_props[prop]     #Reamining props, basically the ones with Intensity in its name\n",
    "        output_file.write(',' + str(to_print))\n",
    "    output_file.write('\\n')\n",
    "\n",
    "\n",
    "plt.output_file.close()\n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.subplot(121), plt.imshow(img)\n",
    "plt.title('Overlay on original image'), plt.xticks([]), plt.yticks([])\n",
    "plt.subplot(122), plt.imshow(img2)\n",
    "plt.title('Colored Grains'), plt.xticks([]), plt.yticks([])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae251318-8db9-4761-acd6-88b10b9c1cc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "img1= cv2.imread(r\"C:\\Users\\cheukjy\\WC25_LLEDGE2_test\\Test_2_H_tiles\\267_29_H_test_2.tif\")\n",
    " \n",
    "# img=img1[:,:,:]\n",
    "img = cv2.cvtColor(img1,cv2.COLOR_BGR2GRAY)\n",
    "# cv2.imshow(\"Black and white\",img)\n",
    "#imshow(\"name of tab opening\", file_that_will_be_opened)\n",
    "\n",
    "# pixels_to_um = 0,10 #1 pixel = 500 nm (got this from the metadata of the original image)\n",
    "\n",
    "#doing this only shows the blue in the image \n",
    "cv2.imread(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49385d2d-557a-4aca-a659-d7959aa27790",
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy import ndimage\n",
    "from skimage import measure, color, io\n",
    "\n",
    "img1 = cv2.imread(r\"C:\\Users\\cheukjy\\WC25_LLEDGE2_test\\Test_2_H_tiles\\267_29_H_test_2.tif\")\n",
    "\n",
    "img = cv2.cvtColor(img1,cv2.COLOR_BGR2GRAY)\n",
    "cv2.imshow(\"Black and white\",img)\n",
    "\n",
    "pixels_to_um = .100 \n",
    "\n",
    "ret1, thresh = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY+cv2.THRESH_OTSU)\n",
    "\n",
    "kernel = np.ones((3,3),np.uint8)\n",
    "opening = cv2.morphologyEx(thresh,cv2.MORPH_OPEN,kernel, iterations = 2)\n",
    "\n",
    "from skimage.segmentation import clear_border\n",
    "opening = clear_border(opening) #Remove edge touching grains\n",
    "\n",
    "sure_bg = cv2.dilate(opening,kernel,iterations=2)\n",
    "cv2.imshow(\"Sure Background\", sure_bg)\n",
    "cv2.waitKey(0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b76f2fb-e7ff-4e2a-be1e-53e1c86a33aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "img1 = cv2.imread(r\"C:\\Users\\cheukjy\\WC25_LLEDGE2_test\\Test_2_H_tiles\\267_29_H_test_2.tif\")\n",
    "#Extract only blue channel as DAPI / nuclear (blue) staining is the best\n",
    "#channel to perform cell count.\n",
    "img=img1[:, :, 0]#Blue channel. Image equivalent to grey image.\n",
    "\n",
    "\n",
    "pixels_to_um = 0.454 # 1 pixel = 454 nm (got this from the metadata of original image)\n",
    "\n",
    "#Threshold image to binary using OTSU. ALl thresholded pixels will be set to 255\n",
    "ret1, thresh = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY+cv2.THRESH_OTSU)\n",
    "\n",
    "# Morphological operations to remove small noise - opening\n",
    "#To remove holes we can use closing\n",
    "kernel = np.ones((3,3),np.uint8)\n",
    "opening = cv2.morphologyEx(thresh,cv2.MORPH_OPEN,kernel, iterations = 2)\n",
    "\n",
    "from skimage.segmentation import clear_border\n",
    "opening = clear_border(opening) #Remove edge touching grains\n",
    "#Check the total regions found before and after applying this. \n",
    "\n",
    "#Now we know that the regions at the center of cells is for sure cells\n",
    "#The region far away is background.\n",
    "#We need to extract sure regions. For that we can use erode. \n",
    "#But we have cells touching, so erode alone will not work. \n",
    "#To separate touching objects, the best approach would be distance transform and then thresholding.\n",
    "\n",
    "# let us start by identifying sure background area\n",
    "# dilating pixes a few times increases cell boundary to background. \n",
    "# This way whatever is remaining for sure will be background. \n",
    "#The area in between sure background and foreground is our ambiguous area. \n",
    "#Watershed should find this area for us.\n",
    "# - in other words, you're looking at the background. and seeing if the pixels \n",
    "#are too far awawy \n",
    "sure_bg = cv2.dilate(opening,kernel,iterations=2)\n",
    "\n",
    "# Finding sure foreground area using distance transform and thresholding\n",
    "#intensities of the points inside the foreground regions are changed to \n",
    "#distance their respective distances from the closest 0 value (boundary).\n",
    "#https://www.tutorialspoint.com/opencv/opencv_distance_transformation.htm\n",
    "## - what is the distance between each black area?\n",
    "dist_transform = cv2.distanceTransform(opening,cv2.DIST_L2,5)\n",
    "\n",
    "#Let us threshold the dist transform by starting at 1/2 its max value.\n",
    "#print(dist_transform.max()) gives about 21.9\n",
    "# Distance transform is when you are looking at how finely tuned you want to\n",
    "# filter only showing the white spaces on the image. The larger the number you put, the less white spots you see\n",
    "# this is helpful to only see cells \n",
    "ret2, sure_fg = cv2.threshold(dist_transform,0.5*dist_transform.max(),255,0)\n",
    "\n",
    "#Later you realize that 0.25* max value will not separate the cells well.\n",
    "#High value like 0.7 will not recognize some cells. 0.5 seems to be a good compromize\n",
    "\n",
    "# Unknown ambiguous region is nothing but bkground - foreground\n",
    "sure_fg = np.uint8(sure_fg)\n",
    "unknown = cv2.subtract(sure_bg,sure_fg)\n",
    "\n",
    "#Now we create a marker and label the regions inside. \n",
    "# For sure regions, both foreground and background will be labeled with positive numbers.\n",
    "# Unknown regions will be labeled 0. \n",
    "#For markers let us use ConnectedComponents. \n",
    "ret3, markers = cv2.connectedComponents(sure_fg)\n",
    "\n",
    "#One problem rightnow is that the entire background pixels is given value 0.\n",
    "#This means watershed considers this region as unknown.\n",
    "#So let us add 10 to all labels so that sure background is not 0, but 10\n",
    "markers = markers+10\n",
    "\n",
    "# Now, mark the region of unknown with zero\n",
    "markers[unknown==255] = 0\n",
    "#plt.imshow(markers, cmap='jet')   #Look at the 3 distinct regions.\n",
    "\n",
    "#Now we are ready for watershed filling. \n",
    "markers = cv2.watershed(img1,markers)\n",
    "#The boundary region will be marked -1\n",
    "#https://docs.opencv.org/3.3.1/d7/d1b/group__imgproc__misc.html#ga3267243e4d3f95165d55a618c65ac6e1\n",
    "\n",
    "#Let us color boundaries in yellow. \n",
    "img1[markers == -1] = [0,255,255]  \n",
    "\n",
    "img2 = color.label2rgb(markers, bg_label=0)\n",
    "\n",
    "cv2.imshow('Overlay on original image', img1)\n",
    "cv2.imshow('Colored Grains', img2)\n",
    "cv2.waitKey(0)\n",
    "\n",
    "#Now, time to extract properties of detected cells\n",
    "# regionprops function in skimage measure module calculates useful parameters for each object.\n",
    "regions = measure.regionprops(markers, intensity_image=img1)\n",
    "\n",
    "#Can print various parameters for all objects\n",
    "for prop in regions:\n",
    "    print('Label: {} Area: {}'.format(prop.label, prop.area))\n",
    "\n",
    "#Best way is to output all properties to a csv file\n",
    "#Let us pick which ones we want to export. \n",
    "\n",
    "propList = ['Area',\n",
    "            'equivalent_diameter', #Added... verify if it works\n",
    "            'orientation', #Added, verify if it works. Angle btwn x-axis and major axis.\n",
    "            'MajorAxisLength',\n",
    "            'MinorAxisLength',\n",
    "            'Perimeter',\n",
    "            'MinIntensity',\n",
    "            'MeanIntensity',\n",
    "            'MaxIntensity']    \n",
    "    \n",
    "\n",
    "output_file = open('cell_measurements.csv', 'w')\n",
    "output_file.write('Cell #,' + \",\".join(propList) + '\\n') #join strings in array by commas, leave first cell blank\n",
    "#First cell blank to leave room for header (column names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4b70653-ca3f-4ca8-a5be-8b7fe4fd0a55",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import cv2 as cv\n",
    "from matplotlib import pyplot as plt\n",
    " \n",
    "img = cv.imread(r'C:\\Users\\cheukjy\\WC25_LLEDGE2_test\\Test_2_H_tiles\\267_29_H_test_2.tif')\n",
    "assert img is not None, \"file could not be read, check with os.path.exists()\"\n",
    "gray = cv.cvtColor(img,cv.COLOR_BGR2GRAY)\n",
    "ret, thresh = cv.threshold(gray,0,255,cv.THRESH_BINARY_INV+cv.THRESH_OTSU)"
   ]
  }
 ],
 "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}