original.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2\n",
    "import matplotlib.pyplot as plt\n",
    "from PIL import Image\n",
    "import numpy as np\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_intensity_pictures(img,composite,outer_bool=False):  \n",
    "    _ , thresh = cv2.threshold(composite, 5, 10, cv2.THRESH_BINARY)\n",
    "    contours,_ = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)\n",
    "\n",
    "    #Find biggest contour\n",
    "    contours = [max(contours, key = cv2.contourArea)]\n",
    "\n",
    "    black = np.full(thresh.shape, 0, dtype=np.uint8)\n",
    "    if outer_bool:\n",
    "        return (black , composite, composite)\n",
    "\n",
    "    white_outer_blackbg = cv2.drawContours(black, contours, -1, 255, -1)\n",
    "    clean_composite = cv2.bitwise_and(composite, white_outer_blackbg)\n",
    "\n",
    "    left = (math.inf,0)\n",
    "    right = (0,0)\n",
    "    top = (0, math.inf)\n",
    "    bottom = (0,0)\n",
    "    threshold_intensity = 10\n",
    "    #Find left most edge\n",
    "    for y, row in enumerate(clean_composite):\n",
    "        for x, pixel in enumerate(row):\n",
    "            if pixel >= threshold_intensity and x < left[0]:\n",
    "                left = (x,y)\n",
    "            #Find right most edge\n",
    "            if pixel >= threshold_intensity and x > right[0]:\n",
    "                right = (x,y)\n",
    "            #Find top most edge\n",
    "            if pixel >= threshold_intensity and y < top[1]:\n",
    "                top = (x,y)\n",
    "            #Find bottom most edge\n",
    "            if pixel >= threshold_intensity and y > bottom[1]:\n",
    "                bottom = (x,y)\n",
    "\n",
    "    #Given four points find top left and bottom right of the rectangle\n",
    "    def find_top_left_bottom_right(p1,p2,p3,p4):\n",
    "        top_left = (min(p1[0],p2[0],p3[0],p4[0]),min(p1[1],p2[1],p3[1],p4[1]))\n",
    "        bottom_right = (max(p1[0],p2[0],p3[0],p4[0]),max(p1[1],p2[1],p3[1],p4[1]))\n",
    "        return top_left, bottom_right\n",
    "\n",
    "    def x_radius(p1,p2):\n",
    "        return (p2[0]-p1[0])/2\n",
    "    def y_radius(p1,p2):\n",
    "        return (p2[1]-p1[1])/2\n",
    "\n",
    "    top_left, bottom_right = find_top_left_bottom_right(left,right,top,bottom)\n",
    "\n",
    "    left_2 = (left[0] + x_radius(top_left,bottom_right)/2, left[1])\n",
    "    right_2 = (right[0] - x_radius(top_left,bottom_right)/2, right[1])\n",
    "    top_2 = (top[0], top[1] + y_radius(top_left,bottom_right)/2)\n",
    "    bottom_2 = (bottom[0], bottom[1] - y_radius(top_left,bottom_right)/2)\n",
    "    top_left_2, bottom_right_2 = find_top_left_bottom_right(left_2,right_2,top_2,bottom_2) \n",
    "    center = (int((top_left_2[0] + bottom_right_2[0])/2), int((top_left_2[1] + bottom_right_2[1])/2))\n",
    "\n",
    "    #Make a mask of the ellipse that sets everything outside the ellipse to 0\n",
    "    mask = np.zeros(composite.shape, dtype=np.uint8)\n",
    "    mask = cv2.ellipse(mask,\n",
    "                         center,\n",
    "                         (int(x_radius(top_left_2,bottom_right_2)), int(y_radius(top_left_2,bottom_right_2)))\n",
    "                         , 0, 0, 360, (255,255 ,255), -1)\n",
    "\n",
    "    #Apply the mask to the image\n",
    "    masked_image = cv2.bitwise_and(composite, mask)\n",
    "\n",
    "    #Find biggest contour\n",
    "    contours = [max(contours, key = cv2.contourArea)]\n",
    "\n",
    "    # Draw contours + hull results\n",
    "    mask_2 = np.zeros(thresh.shape, dtype=np.uint8)\n",
    "\n",
    "    cv2.drawContours(mask_2, contours, -1, (255,255 ,255), -1)\n",
    "    masked_image_2 = cv2.bitwise_and(composite, mask_2)\n",
    "\n",
    "    #Remove inner from outer\n",
    "    mask_1_1 = cv2.bitwise_not(mask)\n",
    "    outer_only = cv2.bitwise_and(masked_image_2, mask_1_1)\n",
    "\n",
    "    return (masked_image, outer_only, composite)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Get matrix of pixel intensities for red and blue for inner and outer\n",
    "def get_intensity_points(img,composite,outer_bool=False):  \n",
    "    _,thresh = cv2.threshold(composite, 5, 10, cv2.THRESH_BINARY)\n",
    "    contours,_ = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)\n",
    "    # composite = composite.astype(np.uint16)\n",
    "    #Find biggest contour\n",
    "    contours = [max(contours, key = cv2.contourArea)]\n",
    "\n",
    "    black = np.full(thresh.shape, 0, dtype=np.uint8)\n",
    "    white_outer_blackbg = cv2.drawContours(black, contours, -1, 255, -1)\n",
    "\n",
    "    clean_composite = cv2.bitwise_and(composite, white_outer_blackbg)\n",
    "\n",
    "    left = (math.inf,0)\n",
    "    right = (0,0)\n",
    "    top = (0, math.inf)\n",
    "    bottom = (0,0)\n",
    "    threshold_intensity = 10\n",
    "    #Find left most edge\n",
    "    for y, row in enumerate(clean_composite):\n",
    "        for x, pixel in enumerate(row):\n",
    "            if pixel >= threshold_intensity and x < left[0]:\n",
    "                left = (x,y)\n",
    "            #Find right most edge\n",
    "            if pixel >= threshold_intensity and x > right[0]:\n",
    "                right = (x,y)\n",
    "            #Find top most edge\n",
    "            if pixel >= threshold_intensity and y < top[1]:\n",
    "                top = (x,y)\n",
    "            #Find bottom most edge\n",
    "            if pixel >= threshold_intensity and y > bottom[1]:\n",
    "                bottom = (x,y)\n",
    "\n",
    "    #Given four points find top left and bottom right of the rectangle\n",
    "    def find_top_left_bottom_right(p1,p2,p3,p4):\n",
    "        top_left = (min(p1[0],p2[0],p3[0],p4[0]),min(p1[1],p2[1],p3[1],p4[1]))\n",
    "        bottom_right = (max(p1[0],p2[0],p3[0],p4[0]),max(p1[1],p2[1],p3[1],p4[1]))\n",
    "        return top_left, bottom_right\n",
    "\n",
    "    def x_radius(p1,p2):\n",
    "        return (p2[0]-p1[0])/2\n",
    "    def y_radius(p1,p2):\n",
    "        return (p2[1]-p1[1])/2\n",
    "\n",
    "    top_left, bottom_right = find_top_left_bottom_right(left,right,top,bottom)\n",
    "\n",
    "    color = (255, 0, 0)\n",
    "    left_2 = (left[0] + x_radius(top_left,bottom_right)/2, left[1])\n",
    "    right_2 = (right[0] - x_radius(top_left,bottom_right)/2, right[1])\n",
    "    top_2 = (top[0], top[1] + y_radius(top_left,bottom_right)/2)\n",
    "    bottom_2 = (bottom[0], bottom[1] - y_radius(top_left,bottom_right)/2)\n",
    "    top_left_2, bottom_right_2 = find_top_left_bottom_right(left_2,right_2,top_2,bottom_2) \n",
    "    center = (int((top_left_2[0] + bottom_right_2[0])/2), int((top_left_2[1] + bottom_right_2[1])/2))\n",
    "\n",
    "    ret, thresh = cv2.threshold(composite, 5, 10, cv2.THRESH_BINARY)\n",
    "    contours,_ = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)\n",
    "    #Find biggest contour\n",
    "    contours = [max(contours, key = cv2.contourArea)]\n",
    "\n",
    "    black = np.full(thresh.shape, 0, dtype=np.uint8)\n",
    "    white_outer_blackbg = cv2.drawContours(black, contours, -1, 255, -1)\n",
    "\n",
    "    white = np.full(img.shape, 255, dtype=np.uint8)\n",
    "    black_ellipse_whitebg= cv2.ellipse(white,\n",
    "                         center,\n",
    "                         (int(x_radius(top_left_2,bottom_right_2)), int(y_radius(top_left_2,bottom_right_2)))\n",
    "                         , 0, 0, 360, 0, -1)\n",
    "\n",
    "    outer_mask = cv2.bitwise_and(white_outer_blackbg, black_ellipse_whitebg)\n",
    "\n",
    "    #Make a mask of the ellipse that sets everything outside the ellipse to 0\n",
    "    inner_mask = np.full(img.shape, 0, dtype=np.uint8)\n",
    "    inner_mask = cv2.ellipse(inner_mask,\n",
    "                         center,\n",
    "                         (int(x_radius(top_left_2,bottom_right_2)), int(y_radius(top_left_2,bottom_right_2)))\n",
    "                         , 0, 0, 360, 255, -1)\n",
    "    \n",
    "    rows_ellipse, cols_ellipse  = np.where(inner_mask > 0)\n",
    "    pts_ellipse = (rows_ellipse, cols_ellipse)\n",
    "\n",
    "    rows_outer, cols_outer = np.where(outer_mask > 0)\n",
    "    pts_outer = (rows_outer, cols_outer)\n",
    "\n",
    "    #Get list of pixel intensities that are in the inner rectangle\n",
    "    def get_pixel_intensities(img, pts):\n",
    "        pixel_intensities = img[pts[0], pts[1]] \n",
    "        return pixel_intensities\n",
    "\n",
    "    outer = get_pixel_intensities(img, pts_outer)\n",
    "    inner = get_pixel_intensities(img, pts_ellipse)\n",
    "    \n",
    "    if outer_bool:\n",
    "        pts_outer = np.where(white_outer_blackbg > 0)\n",
    "        outer = get_pixel_intensities(img, pts_outer)\n",
    "        inner = [0 for i in range(len(outer))]\n",
    "\n",
    "    return (outer, inner)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_pictures_norm(img_red, img, composite, outer_bool = False):\n",
    "    _ , thresh = cv2.threshold(composite, 5, 10, cv2.THRESH_BINARY)\n",
    "    contours,_ = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)\n",
    "\n",
    "    #Find biggest contour\n",
    "    contours = [max(contours, key = cv2.contourArea)]\n",
    "\n",
    "    black = np.full(thresh.shape, 0, dtype=np.uint8)\n",
    "    white_outer_blackbg = cv2.drawContours(black, contours, -1, 255, -1)\n",
    "    clean_composite = cv2.bitwise_and(composite, white_outer_blackbg)\n",
    "\n",
    "    left = (math.inf,0)\n",
    "    right = (0,0)\n",
    "    top = (0, math.inf)\n",
    "    bottom = (0,0)\n",
    "    threshold_intensity = 10\n",
    "\n",
    "    #Find left most edge\n",
    "    for y, row in enumerate(clean_composite):\n",
    "        for x, pixel in enumerate(row):\n",
    "            if pixel >= threshold_intensity and x < left[0]:\n",
    "                left = (x,y)\n",
    "            #Find right most edge\n",
    "            if pixel >= threshold_intensity and x > right[0]:\n",
    "                right = (x,y)\n",
    "            #Find top most edge\n",
    "            if pixel >= threshold_intensity and y < top[1]:\n",
    "                top = (x,y)\n",
    "            #Find bottom most edge\n",
    "            if pixel >= threshold_intensity and y > bottom[1]:\n",
    "                bottom = (x,y)\n",
    "\n",
    "    #Given four points find top left and bottom right of the rectangle\n",
    "    def find_top_left_bottom_right(p1,p2,p3,p4):\n",
    "        top_left = (min(p1[0],p2[0],p3[0],p4[0]),min(p1[1],p2[1],p3[1],p4[1]))\n",
    "        bottom_right = (max(p1[0],p2[0],p3[0],p4[0]),max(p1[1],p2[1],p3[1],p4[1]))\n",
    "        return top_left, bottom_right\n",
    "\n",
    "    def x_radius(p1,p2):\n",
    "        return (p2[0]-p1[0])/2\n",
    "    def y_radius(p1,p2):\n",
    "        return (p2[1]-p1[1])/2\n",
    "\n",
    "    top_left, bottom_right = find_top_left_bottom_right(left,right,top,bottom)\n",
    "\n",
    "    left_2 = (left[0] + x_radius(top_left,bottom_right)/2, left[1])\n",
    "    right_2 = (right[0] - x_radius(top_left,bottom_right)/2, right[1])\n",
    "    top_2 = (top[0], top[1] + y_radius(top_left,bottom_right)/2)\n",
    "    bottom_2 = (bottom[0], bottom[1] - y_radius(top_left,bottom_right)/2)\n",
    "    top_left_2, bottom_right_2 = find_top_left_bottom_right(left_2,right_2,top_2,bottom_2) \n",
    "    center = (int((top_left_2[0] + bottom_right_2[0])/2), int((top_left_2[1] + bottom_right_2[1])/2))\n",
    "\n",
    "    black = np.full(thresh.shape, 0, dtype=np.uint8)\n",
    "    white_outer_blackbg = cv2.drawContours(black, contours, -1, 255, -1)\n",
    "\n",
    "    white = np.full(img.shape, 255, dtype=np.uint8)\n",
    "    black_ellipse_whitebg= cv2.ellipse(white,\n",
    "                        center,\n",
    "                        (int(x_radius(top_left_2,bottom_right_2)), int(y_radius(top_left_2,bottom_right_2)))\n",
    "                        , 0, 0, 360, 0, -1)\n",
    "\n",
    "    outer_mask = cv2.bitwise_and(white_outer_blackbg, black_ellipse_whitebg)\n",
    "\n",
    "    #Make a mask of the ellipse that sets everything outside the ellipse to 0\n",
    "    inner_mask = np.full(img.shape, 0, dtype=np.uint8)\n",
    "    inner_mask = cv2.ellipse(inner_mask,\n",
    "                        center,\n",
    "                        (int(x_radius(top_left_2,bottom_right_2)), int(y_radius(top_left_2,bottom_right_2)))\n",
    "                        , 0, 0, 360, 255, -1)\n",
    "\n",
    "    # Normalize on blue\n",
    "    # plt.imshow(\n",
    "    #     cv2.bitwise_and(outer_mask, cv2.normalize(img_red/img, None, 0, 255, cv2.NORM_MINMAX, dtype=cv2.CV_8U))\n",
    "    # )\n",
    "    # plt.show()\n",
    "\n",
    "    # plt.imshow(\n",
    "    #     cv2.bitwise_and(inner_mask, cv2.normalize(img_red/img, None, 0, 255, cv2.NORM_MINMAX, dtype=cv2.CV_8U))\n",
    "    # )\n",
    "    # plt.show()\n",
    "\n",
    "    # plt.imshow(\n",
    "    #     cv2.bitwise_and(white_outer_blackbg, cv2.normalize(img_red/img, None, 0, 255, cv2.NORM_MINMAX, dtype=cv2.CV_8U))\n",
    "    # )\n",
    "    # plt.show()\n",
    "\n",
    "    if outer_bool:\n",
    "        return (black , cv2.bitwise_and(white_outer_blackbg, cv2.normalize(img_red/img, None, 0, 255, cv2.NORM_MINMAX, dtype=cv2.CV_8U))\n",
    "                      , cv2.bitwise_and(white_outer_blackbg, cv2.normalize(img_red/img, None, 0, 255, cv2.NORM_MINMAX, dtype=cv2.CV_8U)))\n",
    "    return (cv2.bitwise_and(inner_mask, cv2.normalize(img_red/img, None, 0, 255, cv2.NORM_MINMAX, dtype=cv2.CV_8U))\n",
    "           ,cv2.bitwise_and(outer_mask, cv2.normalize(img_red/img, None, 0, 255, cv2.NORM_MINMAX, dtype=cv2.CV_8U))\n",
    "           ,cv2.bitwise_and(white_outer_blackbg, cv2.normalize(img_red/img, None, 0, 255, cv2.NORM_MINMAX, dtype=cv2.CV_8U)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "_, blue = cv2.imreadmulti('data/Blue-34a-dbl_0002.tiff', [], cv2.IMREAD_UNCHANGED)\n",
    "_, red = cv2.imreadmulti('data/Red-34a-dbl_0002.tiff', [], cv2.IMREAD_UNCHANGED)\n",
    "_, composite = cv2.imreadmulti('data/Composite-34a-dbl_0002.tiff', [], cv2.IMREAD_GRAYSCALE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "blue_composite = [(blue[i], composite[i]) for i in range(len(composite))]\n",
    "red_composite = [(red[i], composite[i]) for i in range(len(composite))]\n",
    "red_blue_composite = [(red[i], blue[i], composite[i]) for i in range(len(composite))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        ...,\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0]], dtype=uint8),\n",
       " array([[0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        ...,\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0]], dtype=uint8),\n",
       " array([[0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        ...,\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 0, 0, ..., 0, 0, 0]], dtype=uint8))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(red1, blue1, composite1) = red_blue_composite[10]\n",
    "get_pictures_norm(red1, blue1, composite1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Blue Red Composite\n",
    "images = []\n",
    "for idx, (rd, bl, cmp) in enumerate(red_blue_composite):\n",
    "    if idx < len(red_blue_composite)*0.05 or idx > len(red_blue_composite)*0.95:\n",
    "        point = get_pictures_norm(rd, bl, cmp, outer_bool=True)\n",
    "    else:\n",
    "        point = get_pictures_norm(rd, bl, cmp)\n",
    "    images.append(point)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'get_intensity_line' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m/Users/danielgeorge/Documents/work/bio/spheroid/clean/spheroid-analysis/original.ipynb Cell 9\u001b[0m in \u001b[0;36m<cell line: 3>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/danielgeorge/Documents/work/bio/spheroid/clean/spheroid-analysis/original.ipynb#X40sZmlsZQ%3D%3D?line=2'>3</a>\u001b[0m \u001b[39mfor\u001b[39;00m idx, (bl, cmp) \u001b[39min\u001b[39;00m \u001b[39menumerate\u001b[39m(blue_composite[\u001b[39m50\u001b[39m:\u001b[39m52\u001b[39m]):\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/danielgeorge/Documents/work/bio/spheroid/clean/spheroid-analysis/original.ipynb#X40sZmlsZQ%3D%3D?line=3'>4</a>\u001b[0m     \u001b[39mif\u001b[39;00m idx \u001b[39m<\u001b[39m \u001b[39mlen\u001b[39m(blue_composite)\u001b[39m*\u001b[39m\u001b[39m0.05\u001b[39m \u001b[39mor\u001b[39;00m idx \u001b[39m>\u001b[39m \u001b[39mlen\u001b[39m(blue_composite)\u001b[39m*\u001b[39m\u001b[39m0.95\u001b[39m:\n\u001b[0;32m----> <a href='vscode-notebook-cell:/Users/danielgeorge/Documents/work/bio/spheroid/clean/spheroid-analysis/original.ipynb#X40sZmlsZQ%3D%3D?line=4'>5</a>\u001b[0m         point \u001b[39m=\u001b[39m get_intensity_line(bl, cmp, outer_bool\u001b[39m=\u001b[39m\u001b[39mTrue\u001b[39;00m)\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/danielgeorge/Documents/work/bio/spheroid/clean/spheroid-analysis/original.ipynb#X40sZmlsZQ%3D%3D?line=5'>6</a>\u001b[0m     \u001b[39melse\u001b[39;00m:\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/danielgeorge/Documents/work/bio/spheroid/clean/spheroid-analysis/original.ipynb#X40sZmlsZQ%3D%3D?line=6'>7</a>\u001b[0m         point \u001b[39m=\u001b[39m get_intensity_line(bl, cmp)\n",
      "\u001b[0;31mNameError\u001b[0m: name 'get_intensity_line' is not defined"
     ]
    }
   ],
   "source": [
    "#Blue\n",
    "points_blue = []\n",
    "for idx, (bl, cmp) in enumerate(blue_composite[50:52]):\n",
    "    if idx < len(blue_composite)*0.05 or idx > len(blue_composite)*0.95:\n",
    "        point = get_intensity_line(bl, cmp, outer_bool=True)\n",
    "    else:\n",
    "        point = get_intensity_line(bl, cmp)\n",
    "    points_blue.append(point)\n",
    "\n",
    "# #Red\n",
    "# points_red = []\n",
    "# for idx, (rd, cmp) in enumerate(red_composite):\n",
    "#     if idx < len(red_composite)*0.05 or idx > len(red_composite)*0.95:\n",
    "#         point = get_intensity_points(rd, cmp, outer_bool=True)\n",
    "#     else:\n",
    "#         point = get_intensity_points(rd, cmp)\n",
    "#     points_red.append(point)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'points_blue' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m/Users/danielgeorge/Documents/work/bio/spheroid/clean/spheroid-analysis/original.ipynb Cell 9\u001b[0m in \u001b[0;36m<cell line: 8>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/danielgeorge/Documents/work/bio/spheroid/clean/spheroid-analysis/original.ipynb#X11sZmlsZQ%3D%3D?line=5'>6</a>\u001b[0m inner_mean \u001b[39m=\u001b[39m []\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/danielgeorge/Documents/work/bio/spheroid/clean/spheroid-analysis/original.ipynb#X11sZmlsZQ%3D%3D?line=6'>7</a>\u001b[0m outer_mean \u001b[39m=\u001b[39m []\n\u001b[0;32m----> <a href='vscode-notebook-cell:/Users/danielgeorge/Documents/work/bio/spheroid/clean/spheroid-analysis/original.ipynb#X11sZmlsZQ%3D%3D?line=7'>8</a>\u001b[0m \u001b[39mfor\u001b[39;00m i \u001b[39min\u001b[39;00m \u001b[39mrange\u001b[39m(\u001b[39mlen\u001b[39m(points_blue)):\n\u001b[1;32m      <a href='vscode-notebook-cell:/Users/danielgeorge/Documents/work/bio/spheroid/clean/spheroid-analysis/original.ipynb#X11sZmlsZQ%3D%3D?line=8'>9</a>\u001b[0m     \u001b[39mif\u001b[39;00m i \u001b[39m>\u001b[39m\u001b[39m=\u001b[39m \u001b[39mlen\u001b[39m(blue_composite)\u001b[39m*\u001b[39m\u001b[39m0.05\u001b[39m \u001b[39mand\u001b[39;00m i \u001b[39m<\u001b[39m\u001b[39m=\u001b[39m \u001b[39mlen\u001b[39m(blue_composite)\u001b[39m*\u001b[39m\u001b[39m0.95\u001b[39m:\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/danielgeorge/Documents/work/bio/spheroid/clean/spheroid-analysis/original.ipynb#X11sZmlsZQ%3D%3D?line=9'>10</a>\u001b[0m         \n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/danielgeorge/Documents/work/bio/spheroid/clean/spheroid-analysis/original.ipynb#X11sZmlsZQ%3D%3D?line=10'>11</a>\u001b[0m         \u001b[39m# plt.hist(points_red[i][0]/points_blue[i][0], bins=100, color='b', alpha=0.5, label='Outer',density=True)\u001b[39;00m\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/danielgeorge/Documents/work/bio/spheroid/clean/spheroid-analysis/original.ipynb#X11sZmlsZQ%3D%3D?line=11'>12</a>\u001b[0m         \u001b[39m# plt.hist(points_red[i][1]/points_blue[i][1], bins=100, color='r', alpha=0.5, label='Inner',density=True)\u001b[39;00m\n\u001b[1;32m     <a href='vscode-notebook-cell:/Users/danielgeorge/Documents/work/bio/spheroid/clean/spheroid-analysis/original.ipynb#X11sZmlsZQ%3D%3D?line=13'>14</a>\u001b[0m         plt\u001b[39m.\u001b[39mhist(np\u001b[39m.\u001b[39mlog(points_red[i][\u001b[39m0\u001b[39m]\u001b[39m/\u001b[39mpoints_blue[i][\u001b[39m0\u001b[39m]), bins\u001b[39m=\u001b[39m\u001b[39m100\u001b[39m, color\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39mb\u001b[39m\u001b[39m'\u001b[39m, alpha\u001b[39m=\u001b[39m\u001b[39m0.5\u001b[39m, label\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39mOuter\u001b[39m\u001b[39m'\u001b[39m,density\u001b[39m=\u001b[39m\u001b[39mTrue\u001b[39;00m)\n",
      "\u001b[0;31mNameError\u001b[0m: name 'points_blue' is not defined"
     ]
    }
   ],
   "source": [
    "#Blue vs red dist with images\n",
    "import imageio\n",
    "import os\n",
    "\n",
    "filenames = []\n",
    "inner_mean = []\n",
    "outer_mean = []\n",
    "for i in range(len(points_blue)):\n",
    "    if i >= len(blue_composite)*0.05 and i <= len(blue_composite)*0.95:\n",
    "        \n",
    "        # plt.hist(points_red[i][0]/points_blue[i][0], bins=100, color='b', alpha=0.5, label='Outer',density=True)\n",
    "        # plt.hist(points_red[i][1]/points_blue[i][1], bins=100, color='r', alpha=0.5, label='Inner',density=True)\n",
    "\n",
    "        plt.hist(np.log(points_red[i][0]/points_blue[i][0]), bins=100, color='b', alpha=0.5, label='Outer',density=True)\n",
    "        plt.axvline((np.log(points_red[i][0]/points_blue[i][0])).mean(), color='b', linestyle='dashed', linewidth=1)\n",
    "        outer_mean.append((np.log(points_red[i][0]/points_blue[i][0])).mean())\n",
    "        plt.hist(np.log(points_red[i][1]/points_blue[i][1]), bins=100, color='r', alpha=0.5, label='Inner',density=True)\n",
    "        plt.axvline((np.log(points_red[i][1]/points_blue[i][1])).mean(), color='r', linestyle='dashed', linewidth=1)\n",
    "        inner_mean.append((np.log(points_red[i][1]/points_blue[i][1])).mean())\n",
    "        plt.legend(loc='upper right')\n",
    "        plt.xlabel('Pixel Intensity')\n",
    "        plt.ylabel('Frequency')\n",
    "        plt.title(f'Frequency of Pixel Intensities at {i}')\n",
    "\n",
    "\n",
    "    # create file name and append it to a list\n",
    "    filename = f'{i}.png'\n",
    "    filenames.append(filename)\n",
    "    \n",
    "    # save frame\n",
    "    plt.savefig(filename)\n",
    "    plt.close()\n",
    "\n",
    "def pil_grid(images, max_horiz=np.iinfo(int).max):\n",
    "    n_images = len(images)\n",
    "    n_horiz = min(n_images, max_horiz)\n",
    "    h_sizes, v_sizes = [0] * n_horiz, [0] * (n_images // n_horiz)\n",
    "    for i, im in enumerate(images):\n",
    "        h, v = i % n_horiz, i // n_horiz\n",
    "        h_sizes[h] = max(h_sizes[h], im.size[0])\n",
    "        v_sizes[v] = max(v_sizes[v], im.size[1])\n",
    "    h_sizes, v_sizes = np.cumsum([0] + h_sizes), np.cumsum([0] + v_sizes)\n",
    "    im_grid = Image.new('RGB', (h_sizes[-1], v_sizes[-1]), color='white')\n",
    "    for i, im in enumerate(images):\n",
    "        im_grid.paste(im, (h_sizes[i % n_horiz], v_sizes[i // n_horiz]))\n",
    "    return im_grid\n",
    "\n",
    "# build gif\n",
    "with imageio.get_writer('bluevsreddistwimg.gif', mode='I') as writer:\n",
    "    for i in range(len(filenames)):\n",
    "        if i >= len(blue_composite)*0.05 and i <= len(blue_composite)*0.95:\n",
    "            imgs = [Image.fromarray(x) for x in images[i]] + [Image.open(filenames[i])]\n",
    "            pil_grid(imgs).save('temp.png')\n",
    "            image = imageio.imread('temp.png')\n",
    "            os.remove('temp.png')\n",
    "            writer.append_data(image)\n",
    "        \n",
    "# Remove files\n",
    "for filename in set(filenames):\n",
    "    os.remove(filename)\n",
    "#plt inner mean and outer mean\n",
    "plt.plot(inner_mean, label='Inner', color='r')\n",
    "plt.plot(outer_mean, label='Outer', color='b')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Density scatter plot\n",
    "import imageio\n",
    "import os\n",
    "\n",
    "filenames = []\n",
    "for i in range(len(points_blue)):\n",
    "    if i >= len(blue_composite)*0.05 and i <= len(blue_composite)*0.95:\n",
    "        plt.scatter(points_blue[i][0], points_red[i][0] , c='b', marker='o', label='Outer',alpha=0.01)\n",
    "        plt.scatter(points_blue[i][1], points_red[i][1], c='r', marker='o', label='Inner',alpha=0.01)\n",
    "        plt.legend(loc='upper left')\n",
    "        plt.xlabel('Blue')\n",
    "        plt.ylabel('Red')\n",
    "        ax = plt.gca()\n",
    "        #set log scale\n",
    "        ax.set_xscale('log')\n",
    "        ax.set_yscale('log')\n",
    "        plt.title(f'Blue vs Red Intensity at {i}')\n",
    "\n",
    "\n",
    "    # create file name and append it to a list\n",
    "    filename = f'{i}.png'\n",
    "    filenames.append(filename)\n",
    "    \n",
    "    # save frame\n",
    "    plt.savefig(filename)\n",
    "    plt.close()\n",
    "\n",
    "def pil_grid(images, max_horiz=np.iinfo(int).max):\n",
    "    n_images = len(images)\n",
    "    n_horiz = min(n_images, max_horiz)\n",
    "    h_sizes, v_sizes = [0] * n_horiz, [0] * (n_images // n_horiz)\n",
    "    for i, im in enumerate(images):\n",
    "        h, v = i % n_horiz, i // n_horiz\n",
    "        h_sizes[h] = max(h_sizes[h], im.size[0])\n",
    "        v_sizes[v] = max(v_sizes[v], im.size[1])\n",
    "    h_sizes, v_sizes = np.cumsum([0] + h_sizes), np.cumsum([0] + v_sizes)\n",
    "    im_grid = Image.new('RGB', (h_sizes[-1], v_sizes[-1]), color='white')\n",
    "    for i, im in enumerate(images):\n",
    "        im_grid.paste(im, (h_sizes[i % n_horiz], v_sizes[i // n_horiz]))\n",
    "    return im_grid\n",
    "\n",
    "# build gif\n",
    "with imageio.get_writer('scatterwimg.gif', mode='I') as writer:\n",
    "    for i in range(len(filenames)):\n",
    "        if i >= len(blue_composite)*0.05 and i <= len(blue_composite)*0.95:\n",
    "            imgs = [Image.fromarray(x) for x in images[i]] + [Image.open(filenames[i])]\n",
    "            pil_grid(imgs).save('temp.png')\n",
    "            image = imageio.imread('temp.png')\n",
    "            os.remove('temp.png')\n",
    "            writer.append_data(image)\n",
    "        \n",
    "# Remove files\n",
    "for filename in set(filenames):\n",
    "    os.remove(filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Blue\n",
    "average_blue = []\n",
    "for idx, (bl, cmp) in enumerate(blue_composite):\n",
    "    if idx < len(blue_composite)*0.05 or idx > len(blue_composite)*0.95:\n",
    "        point = get_intensity(bl, cmp, outer_bool=True)\n",
    "    else:\n",
    "        point = get_intensity(bl, cmp)\n",
    "    average_blue.append(point)\n",
    "\n",
    "#Red\n",
    "average_red = []\n",
    "for idx, (rd, cmp) in enumerate(red_composite):\n",
    "    if idx < len(red_composite)*0.05 or idx > len(red_composite)*0.95:\n",
    "        point = get_intensity(rd, cmp, outer_bool=True)\n",
    "    else:\n",
    "        point = get_intensity(rd, cmp)\n",
    "    average_red.append(point)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Blue\n",
    "points_blue = []\n",
    "for idx, (bl, cmp) in enumerate(blue_composite):\n",
    "    if idx < len(blue_composite)*0.05 or idx > len(blue_composite)*0.95:\n",
    "        point = get_intensity_points(bl, cmp, outer_bool=True)\n",
    "    else:\n",
    "        point = get_intensity_points(bl, cmp)\n",
    "    points_blue.append(point)\n",
    "\n",
    "#Red\n",
    "points_red = []\n",
    "for idx, (rd, cmp) in enumerate(red_composite):\n",
    "    if idx < len(red_composite)*0.05 or idx > len(red_composite)*0.95:\n",
    "        point = get_intensity_points(rd, cmp, outer_bool=True)\n",
    "    else:\n",
    "        point = get_intensity_points(rd, cmp)\n",
    "    points_red.append(point)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# bluevsred.gif\n",
    "import imageio\n",
    "import os\n",
    "from matplotlib import cm\n",
    "from matplotlib.colors import Normalize \n",
    "from scipy.interpolate import interpn\n",
    "\n",
    "def density_scatter( x , y, ax = None, sort = True, bins = 20, **kwargs )   :\n",
    "    \"\"\"\n",
    "    Scatter plot colored by 2d histogram\n",
    "    \"\"\"\n",
    "    if ax is None :\n",
    "        fig , ax = plt.subplots()\n",
    "    data , x_e, y_e = np.histogram2d( x, y, bins = bins, density = True )\n",
    "    z = interpn( ( 0.5*(x_e[1:] + x_e[:-1]) , 0.5*(y_e[1:]+y_e[:-1]) ) , data , np.vstack([x,y]).T , method = \"splinef2d\", bounds_error = False)\n",
    "\n",
    "    #To be sure to plot all data\n",
    "    z[np.where(np.isnan(z))] = 0.0\n",
    "\n",
    "    # Sort the points by density, so that the densest points are plotted last\n",
    "    if sort :\n",
    "        idx = z.argsort()\n",
    "        x, y, z = x[idx], y[idx], z[idx]\n",
    "\n",
    "    ax.scatter( x, y, c=z, marker='o', label='Outer',alpha=0.01 )\n",
    "\n",
    "    norm = Normalize(vmin = np.min(z), vmax = np.max(z))\n",
    "    cbar = fig.colorbar(cm.ScalarMappable(norm = norm), ax=ax)\n",
    "    cbar.ax.set_ylabel('Density')\n",
    "\n",
    "    return ax\n",
    "\n",
    "\n",
    "filenames = []\n",
    "for i in range(len(points_blue)):\n",
    "    fig , ax = plt.subplots()\n",
    "    plt.scatter(points_blue[i][0],points_red[i][0] , c='b', marker='o', label='Outer',alpha=0.01)\n",
    "    plt.scatter(points_blue[i][1], points_red[i][1], c='r', marker='o', label='Inner',alpha=0.01)\n",
    "    plt.legend(loc='upper left')\n",
    "    plt.xlabel('Blue')\n",
    "    plt.ylabel('Red')\n",
    "    ax = plt.gca()\n",
    "    #set log scale\n",
    "    ax.set_xscale('log')\n",
    "    ax.set_yscale('log')\n",
    "    plt.title(f'Blue vs Red Intensity at {i}')\n",
    "\n",
    "\n",
    "    # create file name and append it to a list\n",
    "    filename = f'{i}.png'\n",
    "    filenames.append(filename)\n",
    "    \n",
    "    # save frame\n",
    "    plt.savefig(filename)\n",
    "    plt.close()\n",
    "\n",
    "# build gif\n",
    "with imageio.get_writer('bluevsred.gif', mode='I') as writer:\n",
    "    for filename in filenames:\n",
    "        image = imageio.imread(filename)\n",
    "        writer.append_data(image)\n",
    "        \n",
    "# Remove files\n",
    "for filename in set(filenames):\n",
    "    os.remove(filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points_blue, points_red\n",
    "\n",
    "#Blue vs Red scatter plot for each pixel for each frame for both inner and outer\n",
    "for i in range(len(points_blue)):\n",
    "    plt.scatter(points_blue[i][0],points_red[i][0] , c='b', marker='o', label='Outer',alpha=0.01)\n",
    "    plt.scatter(points_blue[i][1], points_red[i][1], c='r', marker='o', label='Inner',alpha=0.01)\n",
    "    plt.legend(loc='upper left')\n",
    "    plt.xlabel('Blue')\n",
    "    plt.ylabel('Red')\n",
    "    plt.title('Blue vs Red')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(points_blue)):\n",
    "    plt.scatter(list(zip(*points_blue))[0][i], list(zip(*points_red))[0][i], c='r', marker='o', label='Inner',alpha=0.01)\n",
    "\n",
    "ax = plt.gca()\n",
    "ax.set_xscale('log')\n",
    "ax.set_yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(points_blue)):\n",
    "    plt.scatter(list(zip(*points_blue))[1][i], list(zip(*points_red))[1][i], c='b', marker='o', label='Outer',alpha=0.01)\n",
    "\n",
    "ax = plt.gca()\n",
    "#set log scale\n",
    "ax.set_xscale('log')\n",
    "ax.set_yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(points_blue)):\n",
    "    plt.scatter(list(zip(*points_blue))[0][i], list(zip(*points_red))[0][i], c='b', marker='o', label='Inner',alpha=0.01)\n",
    "    plt.scatter(list(zip(*points_blue))[1][i], list(zip(*points_red))[1][i], c='r', marker='o', label='Outer',alpha=0.01)\n",
    "    \n",
    "ax = plt.gca()\n",
    "#set log scale\n",
    "ax.set_xscale('log')\n",
    "ax.set_yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#outer, inner\n",
    "#Make a histogram of the pixel intensities for each frame\n",
    "for i in range(len(points_blue)):\n",
    "    if i >= len(blue_composite)*0.05 and i <= len(blue_composite)*0.95:\n",
    "        plt.hist(points_red[i][0]/points_blue[i][0], bins=100, color='b', alpha=0.5, label='Outer',density=True)\n",
    "        plt.hist(points_red[i][1]/points_blue[i][1], bins=100, color='r', alpha=0.5, label='Inner',density=True)\n",
    "        plt.legend(loc='upper right')\n",
    "        plt.xlabel('Pixel Intensity')\n",
    "        plt.ylabel('Frequency')\n",
    "        plt.title('Frequency of Pixel Intensities')\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot outer and inner intensities over time with labels\n",
    "plt.plot([i[0] for i in average_blue], label='Outer')\n",
    "plt.plot([i[1] for i in average_blue], label='Inner')\n",
    "plt.legend()\n",
    "plt.title('Blue')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot outer and inner intensities over time with labels\n",
    "plt.plot([i[0] for i in average_red], label='Outer')\n",
    "plt.plot([i[1] for i in average_red], label='Inner')\n",
    "plt.legend()\n",
    "plt.title('Red')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "blue_composite = [(blue[i], composite[i]) for i in range(len(composite))]\n",
    "red_composite = [(red[i], composite[i]) for i in range(len(composite))]\n",
    "\n",
    "#Blue\n",
    "images = []\n",
    "for idx, (bl, cmp) in enumerate(blue_composite):\n",
    "    if idx < len(blue_composite)*0.05 or idx > len(blue_composite)*0.95:\n",
    "        point = get_intensity_pictures(bl, cmp, outer_bool=True)\n",
    "    else:\n",
    "        point = get_intensity_pictures(bl, cmp)\n",
    "    images.append(point)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tifffile\n",
    "\n",
    "FRAMES = [inner for inner, outer, comp in images] + [outer for inner, outer, comp in images] +  [comp for inner, outer, comp in images]\n",
    "OUT_NAME = \"blue_2.tiff\" # Name to save to\n",
    "with tifffile.TiffWriter(OUT_NAME) as tiff:\n",
    "    for img in FRAMES:\n",
    "        tiff.save(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inner_points_blue = [j for i in points_blue for j in i[0]]\n",
    "inner_points_red = [j for i in points_red for j in i[0]]\n",
    "outer_points_blue = [j for i in points_blue for j in i[1]]\n",
    "outer_points_red = [j for i in points_red for j in i[1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Intrinsic and extrinsic noise for inner\n",
    "#intrinsic = (1/n)*sum([0.5*(blue - red)**2 for blue, red in zip(blue_intensities, red_intensities)])/(mean(blue_intensities) * mean(red_intensities))\n",
    "#extrinsic = (1/n)*sum([blue*red - mean(blue_intesity)*mean(red_intensity) for blue, red in zip(blue_intensities, red_intensities)])/(mean(blue_intensities) * mean(red_intensities))\n",
    "#total = (1/n)*sum([0.5*(blue**2 - red**2) - mean(blue_intesity)*mean(red_intensity) for blue, red in zip(blue_intensities, red_intensities)])/(mean(blue_intensities) * mean(red_intensities))\n",
    "inner_points_blue_mean = np.mean(inner_points_blue)\n",
    "inner_points_red_mean = np.mean(inner_points_red)\n",
    "inner_intrinsic = (1/len(inner_points_blue))*sum([0.5*(blue - red)**2 for blue, red in zip(inner_points_blue, inner_points_red)])/(inner_points_blue_mean * inner_points_red_mean)\n",
    "\n",
    "outer_points_blue_mean = np.mean(outer_points_blue)\n",
    "outer_points_red_mean = np.mean(outer_points_red)\n",
    "outer_intrinsic = (1/len(outer_points_blue))*sum([0.5*(blue - red)**2 for blue, red in zip(outer_points_blue, outer_points_red)])/(outer_points_blue_mean * outer_points_red_mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inner_extrinsic = (1/len(inner_points_blue))*sum([blue*red - inner_points_red_mean*inner_points_blue_mean for blue, red in zip(inner_points_blue, inner_points_red)])/(inner_points_blue_mean * inner_points_red_mean)\n",
    "outer_extrinsic = (1/len(outer_points_blue))*sum([blue*red - outer_points_red_mean*outer_points_blue_mean for blue, red in zip(outer_points_blue, outer_points_red)])/(outer_points_blue_mean * outer_points_red_mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inner_total = (1/len(inner_points_blue))*sum([0.5*(blue**2 - red**2) -  inner_points_red_mean*inner_points_blue_mean for blue, red in zip(inner_points_blue, inner_points_red)])/(inner_points_red_mean*inner_points_blue_mean )\n",
    "outer_total = (1/len(outer_points_blue))*sum([0.5*(blue**2 - red**2) -  outer_points_red_mean*outer_points_blue_mean for blue, red in zip(outer_points_blue, outer_points_red)])/(outer_points_red_mean*outer_points_blue_mean )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inner_intrinsic,inner_extrinsic,inner_total"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outer_intrinsic,outer_extrinsic,outer_total"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.10 ('base')",
   "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.10 | packaged by conda-forge | (main, Feb  1 2022, 21:27:43) \n[Clang 11.1.0 ]"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "224d2b38dd5c4bf9bb346b4b7ffc570356b5a9eca7c29d4f05461d7abc1d8f91"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}