rgb2ihc.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2 as cv\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy import *\n",
    "from skimage.color import rgb2hed, hed2rgb,separate_stains, combine_stains\n",
    "from skimage.exposure import rescale_intensity\n",
    "from matplotlib.colors import LinearSegmentedColormap\n",
    "from skimage import data\n",
    "import numpy as np\n",
    "from mpl_toolkits.mplot3d.axes3d import Axes3D #--- For 3D plot\n",
    "#https://github.com/scikit-image/scikit-image/blob/0e8cff47a533e240508e2485032dac78220ac6c9/skimage/color/colorconv.py#L1375"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 定义工具函数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def norm_by_row(M):\n",
    "    for k in range(M.shape[1]):\n",
    "        M[k,:] /= np.sqrt(np.sum(np.power(M[k,:],2)))\n",
    "    return M\n",
    "\n",
    "def showbychan(im_ihc):\n",
    "    for k in range(3):\n",
    "        plt.figure()\n",
    "        plt.imshow(im_ihc[:, :, k], cmap=\"gray\")\n",
    "        \n",
    "def rgbdeconv(rgb, conv_matrix, C=0):\n",
    "    rgb = rgb.copy().astype(float)\n",
    "    rgb += C\n",
    "    print(rgb.shape)\n",
    "    print(np.reshape(-np.log10(rgb), (-1, 3)).shape)\n",
    "    stains = np.reshape(-np.log10(rgb), (-1, 3)) @ conv_matrix\n",
    "    return np.reshape(stains, rgb.shape)\n",
    "\n",
    "def hecconv(stains, conv_matrix, C=0):\n",
    "#     from skimage.exposure import rescale_intensity\n",
    "    stains = stains.astype(float)\n",
    "    logrgb2 = -np.reshape(stains, (-1, 3)) @ conv_matrix\n",
    "    rgb2 = np.power(10, logrgb2)\n",
    "    return np.reshape(rgb2 - C, stains.shape)\n",
    "\n",
    "def surf(matIn, name=\"fig\", div = (50, 50), SIZE = (8, 6)):\n",
    "    x = np.arange(0, matIn.shape[0])\n",
    "    y = np.arange(0, matIn.shape[1])\n",
    "    x, y = np.meshgrid(y, x)\n",
    "    fig = plt.figure(figsize = SIZE)\n",
    "    ax = Axes3D(fig)\n",
    "    ax.plot_surface(x, y, matIn, rstride=div[0], cstride=div[1], cmap='jet')\n",
    "    plt.title(name)\n",
    "    plt.show()\n",
    "    \n",
    "def in_range(d):\n",
    "    return (0, np.max(cv.GaussianBlur(d.copy(), (3,3), 0)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 读取文件/常量定义"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "H_DAB = array([\n",
    "    [0.65,0.70,0.29],\n",
    "    [0.07, 0.99, 0.11],\n",
    "    [0.27,0.57,0.78]\n",
    "])\n",
    "\n",
    "H_Mou = H_DAB.copy()\n",
    "H_Mou[2,:] = np.cross(H_DAB[0,:], H_DAB[1,:])\n",
    "\n",
    "H_ki67 = H_DAB.copy()\n",
    "H_ki67[1,:] = np.cross(H_DAB[0,:], H_DAB[2,:])\n",
    "\n",
    "im_hex = cv.imread(\"/Users/cunyuan/DATA/ji1024_orig/4d/val1024/chips/13hepatches.tif\", cv.CV_32F)\n",
    "im_ki67 = cv.imread(\"/Users/cunyuan/DATA/ji1024_orig/4d/val1024/dab/ki67patches13.tif\", cv.CV_32F)\n",
    "\n",
    "im_hex = cv.imread(\"/Users/cunyuan/DATA/Kimura/qupath-proj/tiles/01_17-7930_Ki67_HE_1024.0_1.0/53_01_17-7930_Ki67_HE (d=1.0, x=25360.0, y=88432.0, w=1024.0, h=1024.0, z=53).tif\", cv.CV_32F)\n",
    "im_ki67 = cv.imread(\"/Users/cunyuan/DATA/Kimura/qupath-proj/tiles/01_17-7930_Ki67_IHC_1024.0_1.0/53_01_17-7930_Ki67_IHC (d=1.0, x=25360.0, y=88432.0, w=1024.0, h=1024.0, z= 53).tif\", cv.CV_32F)\n",
    "\n",
    "im_hex = cv.imread(\"/Users/cunyuan/DATA/Kimura/qupath-proj/tiles/0.36/results/200/2502/he_seq_0c.tif\")\n",
    "im_ki67 = cv.imread(\"/Users/cunyuan/DATA/Kimura/qupath-proj/tiles/0.36/results/200/2502/test_seq_0c.tif\")\n",
    "\n",
    "im_mask = cv.cvtColor(cv.imread(\"/Users/cunyuan/DATA/ji1024_orig/4d/val1024/masks/13hepatches_mask.tif\",\n",
    "                                cv.CV_32F), cv.COLOR_BGR2GRAY)\n",
    "\n",
    "im_hex = cv.cvtColor(im_hex, cv.COLOR_BGR2RGB)/255.\n",
    "im_ki67 = cv.cvtColor(im_ki67, cv.COLOR_BGR2RGB)/255.\n",
    "\n",
    "cmap_hema = LinearSegmentedColormap.from_list('mycmap', ['white', 'navy'])\n",
    "cmap_eosin = LinearSegmentedColormap.from_list('mycmap', ['white','darkviolet'])\n",
    "cmap_dab = LinearSegmentedColormap.from_list('mycmap', ['white','saddlebrown'])\n",
    "cmap_res = LinearSegmentedColormap.from_list('mycmap', ['white','green'])\n",
    "print(\"Trans. H.E\", H_Mou)\n",
    "print(\"Trans. Ki67\", H_ki67)\n",
    "plt.imshow(im_ki67)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 实验：染色解组"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "H = H_Mou\n",
    "Hinv = linalg.inv(norm_by_row(H))\n",
    "\n",
    "img = im_hex\n",
    "img[img==0] = 1E-6\n",
    "im_sepa_hex=abs(rgbdeconv(img, Hinv))\n",
    "h = im_sepa_hex[:,:,0];e = im_sepa_hex[:,:,1];d_r = im_sepa_hex[:,:,2];\n",
    "\n",
    "fig = plt.figure(figsize=(10,10));\n",
    "plt.subplot(221);plt.imshow(img);plt.title(\"Input\");plt.axis('off')\n",
    "axis=plt.subplot(222);plt.imshow(rescale_intensity(h, in_range=in_range(h)), cmap=cmap_hema);plt.title(\"Hema.\"),plt.axis('off')\n",
    "plt.subplot(223);plt.imshow(rescale_intensity(e, in_range=in_range(e)), cmap=cmap_eosin);plt.title(\"Eosin\");plt.axis('off')\n",
    "plt.subplot(224);plt.imshow(rescale_intensity(d_r, in_range=in_range(d_r)), cmap=cmap_res);plt.title(\"Residual\");plt.axis('off')\n",
    "fig.tight_layout()\n",
    "# print(im_sepa_hex)\n",
    "\n",
    "surf(d_r, \"residual\", div=(50,50))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h1 = rescale_intensity(h, in_range=in_range(h))\n",
    "plt.figure(figsize=(10,10));plt.imshow(h1, cmap=cmap_hema)\n",
    "surf(h1, \"Rescaled Hematoxylin\")\n",
    "plt.figure(figsize=(10,10));plt.imshow(h1>0.3, cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img_r = img[:,:,0]\n",
    "plt.figure(figsize=(10,10));plt.imshow(img_r, cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from cv2 import ximgproc\n",
    "def gfLabel(img_r, mode=-1, hard_thresh=0.5, geps=1000, grad=5,\n",
    "           MOP_SIZE=3, MCL_SIZE=3, MORPH_ITER=2):\n",
    "    img_r = rescale_intensity(img_r, in_range=in_range(img_r))\n",
    "    d1 = img_r\n",
    "    d2 = (d1*255).astype(uint8)\n",
    "    \n",
    "    d2 = cv.morphologyEx(d2, op=cv.MORPH_OPEN, kernel=np.ones((MOP_SIZE,MOP_SIZE), uint8), iterations = MORPH_ITER)\n",
    "#     d2 = cv.morphologyEx(d2, op=cv.MORPH_CLOSE, kernel=np.ones((MCL_SIZE, MCL_SIZE), uint8), iterations=MORPH_ITER)\n",
    "\n",
    "    guidedDAB = ximgproc.guidedFilter(guide=d2, src=d2, radius=grad, eps=geps, dDepth=-1)\n",
    "\n",
    "    if mode==-1:\n",
    "        d1 *= (guidedDAB < (hard_thresh*255))\n",
    "        # _, _d1 = cv.threshold(guidedDAB,25,255,cv.THRESH_BINARY+cv.THRESH_OTSU)\n",
    "        # d1 = (guidedDAB/255)*(_d1/255)\n",
    "        d2 = (d1*255).astype(uint8)\n",
    "    else:\n",
    "        d1 *= (guidedDAB > (hard_thresh*255))\n",
    "        d2 = (d1*255).astype(uint8)\n",
    "    guide = d2\n",
    "\n",
    "    guidedDAB = ximgproc.guidedFilter(guide=guide, src=d2, radius=grad, eps=geps, dDepth=-1)\n",
    "\n",
    "    _,guidedDAB = cv.threshold((guidedDAB).astype(uint8),0,255,cv.THRESH_BINARY+cv.THRESH_OTSU)\n",
    "\n",
    "    guidedDAB = ((guidedDAB) > 0).astype(uint8)\n",
    "    plt.figure(figsize=(20,10));plt.subplot(121);plt.imshow(guidedDAB, cmap = \"gray\");plt.axis('off')\n",
    "\n",
    "    return guidedDAB\n",
    "\n",
    "guidedDAB = gfLabel(img_r)\n",
    "plt.figure(figsize=(10,10));plt.imshow(guidedDAB, cmap = \"gray\");plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "im_mask =gfLabel(h, mode=1, hard_thresh = 0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 实验：细胞分割"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx= 10\n",
    "im_he = im_hex.copy()\n",
    "im_ihc = im_ki67.copy()\n",
    "\n",
    "thresh = im_mask\n",
    "\n",
    "# noise removal\n",
    "kernel = np.ones((3,3),np.uint8)\n",
    "# opening = cv.morphologyEx(thresh,cv.MORPH_OPEN,kernel, iterations = 10)\n",
    "opening = thresh\n",
    "# sure background area\n",
    "sure_bg = np.uint8(cv.dilate(opening,kernel,iterations=5))\n",
    "\n",
    "# Finding sure foreground area\n",
    "dist_transform = cv.distanceTransform(opening,cv.DIST_L2,3)\n",
    "ret, sure_fg = cv.threshold(dist_transform,0.3*dist_transform.max(),255,0)\n",
    "# ret, sure_fg = cv.threshold((dist_transform).astype(uint8),0,dist_transform.max(),cv.THRESH_BINARY+cv.THRESH_OTSU)\n",
    "# sure_fg  = cv.adaptiveThreshold((dist_transform).astype(uint8), dist_transform.max(), cv.ADAPTIVE_THRESH_GAUSSIAN_C, cv.THRESH_BINARY, 3,5)\n",
    "sure_fg = np.uint8(sure_fg)\n",
    "unknown = cv.subtract(sure_bg,sure_fg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(unknown)\n",
    "unknown.max()\n",
    "plt.figure();plt.imshow(dist_transform)\n",
    "plt.figure();plt.imshow(sure_bg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(30,10))\n",
    "plt.subplot(131);plt.imshow(im_ihc)\n",
    "plt.subplot(132);plt.imshow(im_he)\n",
    "plt.subplot(133);plt.imshow(opening)\n",
    "plt.figure();plt.imshow(im_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Marker labelling\n",
    "ret, markers = cv.connectedComponents(sure_fg)\n",
    "# Add one to all labels so that sure background is not 0, but 1\n",
    "markers = markers+1\n",
    "# Now, mark the region of unknown with zero\n",
    "markers[unknown==1] = 0\n",
    "\n",
    "plt.figure(figsize=(4,4), dpi=300);plt.imshow(markers, cmap='jet');plt.axis(\"off\")\n",
    "surf(markers>1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "markers = cv.watershed(np.uint8(im_ihc*255),markers)\n",
    "im_he[markers == -1] = [0,255,0]\n",
    "im_ihc[markers == -1] = [0,255,0]\n",
    "plt.figure(figsize=(8,8), dpi=300);plt.imshow(markers, cmap='jet');plt.axis(\"off\")\n",
    "# plt.figure(figsize=(8,4), dpi=300)\n",
    "# plt.subplot(121);plt.imshow(im_he);plt.axis('off')\n",
    "# plt.subplot(122);plt.imshow(im_ihc);plt.axis('off')\n",
    "np.unique(markers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ret, thresh = cv.threshold(im_mask, 127, 255, 0)\n",
    "contours, hierarchy = cv.findContours(im_mask, cv.RETR_TREE, cv.CHAIN_APPROX_SIMPLE)\n",
    "print(len(contours))\n",
    "plt.figure(figsize=(12,6), dpi=300)\n",
    "plt.subplot(121);plt.imshow(im_ihc);plt.axis('off');plt.tight_layout()\n",
    "plt.subplot(122);plt.imshow(cv.drawContours(im_ki67.copy(), contours, -1, (0,255,0), 1));plt.axis('off');plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "imci=abs(hecconv(im_sepa_hex, H))\n",
    "plt.figure(figsize = (10,10));\n",
    "plt.imshow(imci)\n",
    "\n",
    "plt.figure(figsize = (10,10));\n",
    "hmod = (abs(im_sepa_hex[:,:,0]))*im_mask/255\n",
    "plt.imshow(hmod, cmap=cmap_dab); plt.title(\"Pseudo DAB: Mask + Hema. Texture\")\n",
    "\n",
    "surf(hmod, \"Pseudo DAB surf.\", div=(50,50))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 实验：ki67染色解组及重组"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "tags": []
   },
   "outputs": [],
   "source": [
    "H = H_ki67\n",
    "Hinv = linalg.inv(norm_by_row(H))\n",
    "plt.figure();plt.imshow(im_ki67)\n",
    "img = np.float32(im_ki67)\n",
    "img[img==0] = 1E-6\n",
    "im_sepa_ki67=abs(rgbdeconv(img, Hinv))\n",
    "\n",
    "h = im_sepa_ki67[:,:,0];e_r = im_sepa_ki67[:,:,1];d = im_sepa_ki67[:,:,2];\n",
    "\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(10,10));\n",
    "plt.subplot(221);plt.imshow(img);plt.title(\"Input\");plt.axis('off')\n",
    "plt.subplot(222);plt.imshow(rescale_intensity(h, in_range=in_range(h)), cmap=cmap_hema);plt.title(\"Hema.\");plt.axis('off')\n",
    "plt.subplot(223);plt.imshow(rescale_intensity(e_r, in_range=in_range(e_r)), cmap=cmap_res);plt.title(\"Residual\");plt.axis('off')\n",
    "plt.subplot(224);plt.imshow(rescale_intensity(d, in_range=in_range(d)), cmap=cmap_dab);plt.title(\"DAB\");plt.axis('off')\n",
    "fig.tight_layout()\n",
    "\n",
    "surf(e_r, \"Residual (Eosin)\")\n",
    "surf(rescale_intensity(d, in_range=in_range(d)), \"DAB\")\n",
    "d.max()\n",
    "plt.imshow(rescale_intensity(d, in_range=in_range(d)), cmap=\"gray\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "h1 = rescale_intensity(h, in_range=in_range(h))\n",
    "d1 = rescale_intensity(d, in_range=in_range(d))\n",
    "plt.figure(figsize=(10,10));plt.imshow(d1, cmap=cmap_dab);plt.axis('off')\n",
    "plt.figure(figsize=(10,10));plt.imshow(d1>0.7, cmap=\"gray\")\n",
    "plt.figure(figsize=(10,10));plt.imshow(d1>0.5, cmap=\"gray\")\n",
    "plt.figure(figsize=(10,10));plt.imshow(d1>0.1, cmap=\"gray\")\n",
    "plt.figure(figsize=(10,10));plt.imshow(d1>0.05, cmap=\"gray\")\n",
    "plt.figure(figsize=(10,10));plt.imshow(h1>0.1, cmap=\"gray\")\n",
    "dmask = d1>0.05"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 引导滤波"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from cv2 import ximgproc\n",
    "\n",
    "d2 = (d1*255).astype(uint8)\n",
    "\n",
    "guide = d2\n",
    "\n",
    "geps = 1000; grad=5\n",
    "guidedDAB = ximgproc.guidedFilter(guide=guide, src=d2, radius=grad, eps=geps, dDepth=-1)\n",
    "gd = ximgproc.guidedFilter(guide=guide, src=d2, radius=grad, eps=geps, dDepth=-1)\n",
    "gd = rescale_intensity(gd, in_range=in_range(gd))\n",
    "\n",
    "fig = plt.figure(figsize=(20,20))\n",
    "axis = plt.subplot(1, 2, 2)\n",
    "axis.imshow(gd, cmap=cmap_dab)\n",
    "axis.set_title(\"Guided Filter Applied. Radius=%d, Eps=%d\"%(grad, geps))\n",
    "axis.axis('off')\n",
    "\n",
    "axis = plt.subplot(1,2,1)\n",
    "axis.imshow(d1, cmap=cmap_dab)\n",
    "axis.set_title(\"Orignal DAB\")\n",
    "axis.axis('off')\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "guidedDAB = d2\n",
    "plt.figure(figsize=(10,10));surf(gd)\n",
    "_,dmask = cv.threshold((guidedDAB*255).astype(uint8),0,255,cv.THRESH_BINARY+cv.THRESH_OTSU)\n",
    "\n",
    "guidedDAB = ((255-dmask) > 0).astype(uint8)\n",
    "plt.figure(figsize=(20,10));plt.subplot(121);plt.imshow(guidedDAB, cmap = \"gray\");plt.axis('off');plt.title(\"Thresh\")\n",
    "\n",
    "MOP_SIZE, MCL_SIZE = 7, 9\n",
    "guidedDAB = cv.morphologyEx(guidedDAB, op=cv.MORPH_OPEN, kernel=np.ones((MOP_SIZE,MOP_SIZE), uint8))\n",
    "# guidedDAB = cv.morphologyEx(guidedDAB, op=cv.MORPH_CLOSE, kernel=np.ones((MCL_SIZE, MCL_SIZE), uint8))\n",
    "\n",
    "plt.figure(figsize=(10,10));plt.imshow(guidedDAB, cmap = \"gray\");plt.axis('off');plt.title(\"MORPH proc\")\n",
    "\n",
    "im_mask = guidedDAB\n",
    "plt.figure(figsize = (10,10));\n",
    "hmod = (abs(im_sepa_hex[:,:,0]))*im_mask\n",
    "plt.imshow(hmod, cmap=cmap_dab); plt.title(\"Pseudo DAB: Mask + Hema. Texture\")\n",
    "surf(hmod, \"Pseudo DAB surf.\", div=(50,50))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "imci=hecconv(im_sepa_ki67, H)\n",
    "plt.figure(figsize=(10,10));plt.imshow(imci)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thresh = np.expand_dims(im_mask.astype(np.uint8), -1)\n",
    "\n",
    "im_he = cv.cvtColor(im_hex, cv.COLOR_BGR2RGB).astype(np.uint8)\n",
    "im_ihc = cv.cvtColor(im_ki67, cv.COLOR_BGR2RGB).astype(np.uint8)\n",
    "\n",
    "# noise removal\n",
    "kernel = np.ones((3,3),np.uint8)\n",
    "opening = cv.morphologyEx(thresh,cv.MORPH_OPEN,kernel, iterations = 5)\n",
    "\n",
    "# sure background area\n",
    "sure_bg = cv.dilate(opening,kernel,iterations=10)\n",
    "\n",
    "# Finding sure foreground area\n",
    "dist_transform = cv.distanceTransform(opening,cv.DIST_L2,3)\n",
    "ret, sure_fg = cv.threshold(dist_transform,0.5*dist_transform.max(),255,0)\n",
    "\n",
    "# Finding unknown region\n",
    "sure_fg = np.uint8(sure_fg)\n",
    "unknown = cv.subtract(sure_bg,sure_fg)\n",
    "\n",
    "\n",
    "plt.figure(figsize=(30,10))\n",
    "plt.subplot(131);plt.imshow(im_ihc)\n",
    "plt.subplot(132);plt.imshow(im_he)\n",
    "plt.subplot(133);plt.imshow(opening)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(dist_transform, cmap=\"gray\");plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Marker labelling\n",
    "ret, markers = cv.connectedComponents(sure_fg)\n",
    "# Add one to all labels so that sure background is not 0, but 1\n",
    "markers = markers+1\n",
    "# Now, mark the region of unknown with zero\n",
    "markers[unknown==255] = 0\n",
    "\n",
    "plt.figure(figsize=(10,10));plt.imshow(markers, cmap='jet')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "markers = cv.watershed(im_he,markers)\n",
    "im_he[markers == -1] = [0,255,0]\n",
    "im_ihc[markers == -1] = [0,255,0]\n",
    "plt.imshow(markers>1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(40,20))\n",
    "plt.subplot(121);plt.imshow(im_he);plt.axis('off')\n",
    "plt.subplot(122);plt.imshow(im_ihc);plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,10));plt.imshow(markers);plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 实验：染色浓度调整"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from skimage.exposure import rescale_intensity\n",
    "\n",
    "# Rescale hematoxylin and DAB signals and give them a fluorescence look\n",
    "\n",
    "\n",
    "stain_scale_factor = np.array([0.8, 0.05, 0.5]).T # to adjust the instensity of dab in pseudo image\n",
    "\n",
    "im_sepa_hex[np.where(im_sepa_hex==inf)] = 0\n",
    "h = im_sepa_hex[:, :, 0] # hema.\n",
    "r = im_sepa_hex[:, :, 1] # residual. (no eosin applied)\n",
    "\n",
    "zdh = np.dstack((h, r, hmod))\n",
    "\n",
    "zdh = (zdh.reshape(-1,3)*stain_scale_factor).reshape(im_sepa_hex.shape)\n",
    "\n",
    "fig = plt.figure(figsize=(20,20))\n",
    "axis = plt.subplot(1, 2, 1)\n",
    "axis.imshow(hecconv(zdh, H))\n",
    "axis.set_title(\"Pseudo Hematoxylin-DAB\")\n",
    "axis.axis('off')\n",
    "\n",
    "axis = plt.subplot(1,2,2)\n",
    "axis.imshow(img)\n",
    "axis.set_title(\"True Hematoxylin-DAB\")\n",
    "axis.axis('off')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 实验：染色浓度配准"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def reject_outliers(data, m=3):\n",
    "    dt = data.copy()\n",
    "    dt[(abs(data - np.mean(data)) > m * np.std(data))] = np.mean(data)\n",
    "    return dt\n",
    "\n",
    "h3 = reject_outliers(h)\n",
    "d3 = reject_outliers(d)\n",
    "\n",
    "h1 = rescale_intensity(reject_outliers(h), \n",
    "                       in_range=in_range(reject_outliers(h)))\n",
    "d1 = rescale_intensity(d, in_range(reject_outliers(d)))\n",
    "plt.imshow(h3, cmap=cmap_hema)\n",
    "plt.figure();plt.imshow(d1, cmap=cmap_dab)\n",
    "plt.figure();plt.imshow(d3, cmap=cmap_dab)\n",
    "surf(h1, \"Rescl. H\")\n",
    "surf(guidedDAB*d1, \"Rescl. D\")\n",
    "plt.figure(); plt.hist(d1.reshape(-1,), 50)\n",
    "plt.figure(); plt.hist((d1[d1>0.1]).reshape(-1,), 50)\n",
    "plt.figure(); plt.hist(d3.reshape(-1,), 50)\n",
    "plt.figure(); plt.hist((gd/255).reshape(-1,), 50)\n",
    "\n",
    "plt.figure(); plt.hist(h.reshape(-1,), 50)\n",
    "plt.figure(); plt.hist(h3.reshape(-1,), 50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,10))\n",
    "lbl = np.dstack([h1,  guidedDAB*d1, zeros_like(h1),])\n",
    "lbl = np.dstack([h1,  e_r, guidedDAB*d1])\n",
    "lbl_toshow = lbl\n",
    "plt.imshow(lbl_toshow)\n",
    "plt.axis('off')\n",
    "h.max()\n",
    "plt.figure();plt.hist(h.reshape(-1,), 50)\n",
    "plt.figure();plt.hist(d.reshape(-1,), 50)\n",
    "plt.imsave(\"lb.png\", lbl_toshow)\n",
    "lbl = np.dstack([h,  e_r, d])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 蒙板阈值测试"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from PIL import Image\n",
    "im = Image.open(\"./lb.png\");im"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "recover_ki67 = hecconv(lbl, H_ki67)\n",
    "plt.imshow(recover_ki67);plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "surf(h1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "surf(d1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "goupi = cv.cvtColor(cv.imread(\"/Users/cunyuan/DATA/Kimura/qupath-proj/tiles/Label_256/1506_01_17-7930_Ki67_IHC (d=1.0, x=29968.0, y=86384.0, w=256.0, h=256.0, z= 1506).png\", cv.CV_32F), cv.COLOR_BGR2GRAY)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "goupi.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "surf(goupi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.unique(goupi.reshape(-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g3 = [1018,\n",
    "984,\n",
    "931,\n",
    "841,\n",
    "692,\n",
    "517,\n",
    "311,\n",
    "114,\n",
    "15]\n",
    "_4d = [1057,\n",
    "921,\n",
    "734,\n",
    "202,\n",
    "1,\n",
    "0,\n",
    "0,\n",
    "0,\n",
    "0]\n",
    "th = [0.01,\n",
    "0.02,\n",
    "0.04,\n",
    "0.09,\n",
    "0.16,\n",
    "0.25,\n",
    "0.36,\n",
    "0.49,\n",
    "0.64\n",
    "]\n",
    "_4d = np."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,8))\n",
    "plt.plot(th, g3, th, _4d)\n",
    "plt.grid()\n",
    "plt.xlabel(\"Threshold\")\n",
    "plt.ylabel(\"# of tiles\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.1-final"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}