Color_metrics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Standard libraries\n",
    "import os\n",
    "from scipy.optimize import curve_fit, root\n",
    "from sklearn.metrics import r2_score, auc\n",
    "from scipy import stats\n",
    "from scipy.signal import deconvolve\n",
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "import pandas as pd\n",
    "\n",
    "# Plotting and image libraries\n",
    "from matplotlib import pyplot as plt, style\n",
    "import matplotlib as mpl\n",
    "import matplotlib.image as mpimg\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "# Custom libraries\n",
    "import frgtools.misc as frgm\n",
    "from frgtools.curveprocessing import gaussian\n",
    "\n",
    "# Other utilities\n",
    "from natsort import natsorted\n",
    "from tqdm import tqdm\n",
    "import ruptures as rpt\n",
    "\n",
    "# OpenCV\n",
    "import cv2\n",
    "\n",
    "import cv2\n",
    "import numpy as np\n",
    "from tqdm import tqdm\n",
    "import plotly.express as px\n",
    "\n",
    "# Jupyter notebook settings\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "from scipy.ndimage import gaussian_filter\n",
    "\n",
    "# Style settings\n",
    "style.use('ggplot')  # or plt.style.use('ggplot')\n",
    "mpl.rcParams.update(mpl.rcParamsDefault)\n",
    "\n",
    "# Color choice\n",
    "color_choice = ['#F8DDA4', '#C8D5B9', '#68B0AB', '#4A7C59'][::-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def modified_z_score(intensity):\n",
    "    median_int = np.median(intensity)\n",
    "    mad_int = np.median([np.abs(intensity -median_int)])\n",
    "    modified_z_scores = 0.6745 * (intensity- median_int) / mad_int\n",
    "    return modified_z_scores\n",
    "\n",
    "def fixer(y,m=3, th=7):\n",
    "    try:\n",
    "        th = 7 # binarization threshold. \n",
    "        spikes = abs(np.array(modified_z_score(np.diff(y)))) > th\n",
    "        y_out = y.copy() # So we don’t overwrite y\n",
    "        for i in np.arange(len(spikes)):\n",
    "            if spikes[i] != 0: # If we have an spike in position i\n",
    "                w = np.arange(i-m,i+1+m) # we select 2 m + 1 points around our spike\n",
    "                w2 = w[spikes[w] == 0] # From such interval, we choose the ones which are not spikes\n",
    "                y_out[i] = np.mean(y[w2]) # and we average their values\n",
    "    except:\n",
    "        m = 1\n",
    "        th = 7 # binarization threshold. \n",
    "        spikes = abs(np.array(modified_z_score(np.diff(y)))) > th\n",
    "        y_out = y.copy() # So we don’t overwrite y\n",
    "        for i in np.arange(len(spikes)):\n",
    "            if spikes[i] != 0: # If we have an spike in position i\n",
    "                w = np.arange(i-m,i+1+m) # we select 2 m + 1 points around our spike\n",
    "                w2 = w[spikes[w] == 0] # From such interval, we choose the ones which are not spikes\n",
    "                y_out[i] = np.mean(y[w2]) # and we average their values        \n",
    " \n",
    "    return y_out\n",
    "\n",
    "RGB_SCALE = 255\n",
    "CMYK_SCALE = 100\n",
    "\n",
    "\n",
    "def rgb_to_cmyk(r, g, b):\n",
    "    if (r, g, b) == (0, 0, 0):\n",
    "        # black\n",
    "        return 0, 0, 0, CMYK_SCALE\n",
    "\n",
    "    # rgb [0,255] -> cmy [0,1]\n",
    "    c = 1 - r / RGB_SCALE\n",
    "    m = 1 - g / RGB_SCALE\n",
    "    y = 1 - b / RGB_SCALE\n",
    "\n",
    "    # extract out k [0, 1]\n",
    "    min_cmy = min(c, m, y)\n",
    "    c = (c - min_cmy) / (1 - min_cmy)\n",
    "    m = (m - min_cmy) / (1 - min_cmy)\n",
    "    y = (y - min_cmy) / (1 - min_cmy)\n",
    "    k = min_cmy\n",
    "\n",
    "    # rescale to the range [0,CMYK_SCALE]\n",
    "    return c * CMYK_SCALE, m * CMYK_SCALE, y * CMYK_SCALE, k * CMYK_SCALE\n",
    "# test = df['img'][6][300,:,:,:]/255\n",
    "\n",
    "def rgb_to_cmyk(rgb):\n",
    "\tr, g, b = rgb[...,0], rgb[...,1], rgb[...,2]\n",
    "\tk = 1 - np.max(rgb, axis=-1)\n",
    "\tc = (1-r-k)/(1-k)\n",
    "\tm = (1-g-k)/(1-k)\n",
    "\ty = (1-b-k)/(1-k)\n",
    "\treturn np.dstack([c, m, y, k])\n",
    "\n",
    "def rgb_to_cmyk_2(rgb):\n",
    "\twith np.errstate(invalid='ignore', divide='ignore'):\n",
    "\t\tK = 1 - np.max(test, axis=2)\n",
    "\t\tC = (1-test[...,0] - K)/(1-K)\n",
    "\t\tM = (1-test[...,1] - K)/(1-K)\n",
    "\t\tY = (1-test[...,2] - K)/(1-K)\n",
    "\treturn np.dstack([C, M, Y, K])\n",
    "\n",
    "def sigmoid(x, L ,x0, k, b):\n",
    "    y = L / (1 + np.exp(-k*(x-x0))) + b\n",
    "    return (y)\n",
    "\n",
    "def remove_inf(x, y):\n",
    "    mask = ~np.isinf(x) & ~np.isinf(y)\n",
    "    x = x[mask]\n",
    "    y = y[mask]\n",
    "    return x, y\n",
    "\n",
    "def linear(t, b, m):\n",
    "    y = m*t + b\n",
    "    return (y)\n",
    "\n",
    "def avrami(t, k, n):\n",
    "    y = 1 - np.exp(-k*(t)**n)\n",
    "    return (y)\n",
    "\n",
    "def JMAK(t, e, k, n):\n",
    "    R = 8.314 #J/mol/K\n",
    "    T = 273.15+230 #K\n",
    "    t0=0\n",
    "    y = 1 - np.exp(-(k*np.exp(-e/(R*T))*(t-t0))**n)\n",
    "    return (y)\n",
    "\n",
    "def nan_helper(y):\n",
    "    return np.isnan(y), lambda z: z.nonzero()[0]\n",
    "\n",
    "\n",
    "# Define the RGB to CMYK conversion function\n",
    "def rgb_to_cmyk(rgb):\n",
    "    r, g, b = rgb[...,0], rgb[...,1], rgb[...,2]\n",
    "    k = 1 - np.max(rgb, axis=-1)\n",
    "    c = (1-r-k)/(1-k)\n",
    "    m = (1-g-k)/(1-k)\n",
    "    y = (1-b-k)/(1-k)\n",
    "    return np.dstack([c, m, y, k])\n",
    "\n",
    "# Function to calculate the K (black) mean for a video matrix\n",
    "# def calculate_K_mean(vm):\n",
    "#     # Convert each frame to CMYK and calculate the mean of the K channel\n",
    "#     return [1-np.mean(rgb_to_cmyk(vm[frame, :, :, :]/255)[:,:,3]) for frame in range(vm.shape[0])]\n",
    "\n",
    "# Apply the function to the 'VM' column of the DataFrame\n",
    "# all_samples_df['K_mean'] = all_samples_df['VM'].apply(calculate_K_mean)\n",
    "\n",
    "\n",
    "\n",
    "def calculate_K_mean(vm):\n",
    "    K_means = []\n",
    "    for frame in range(vm.shape[0]):\n",
    "        # Convert the frame to CMYK\n",
    "        cmyk_image = rgb_to_cmyk(vm[frame, :, :, :]/255)\n",
    "        \n",
    "        # Extract the K channel\n",
    "        K_channel = cmyk_image[:,:,3]\n",
    "        \n",
    "        # Calculate the 2.5th and 97.5th percentiles (95% CI)\n",
    "        ci_low, ci_high = np.percentile(K_channel, [2.5, 97.5])\n",
    "        \n",
    "        # Mask the values outside the CI\n",
    "        masked_K_channel = np.ma.masked_outside(K_channel, ci_low, ci_high)\n",
    "        \n",
    "        # Calculate the mean of the values within the CI\n",
    "        K_mean = np.ma.mean(masked_K_channel)\n",
    "        \n",
    "        K_means.append(K_mean)\n",
    "        \n",
    "    return K_means\n",
    "\n",
    "\n",
    "def calculate_cv(values):\n",
    "    mean_value = np.nanmean(values)\n",
    "    std_dev = np.nanstd(values)\n",
    "    cv = (std_dev / mean_value) * 100 if mean_value != 0 else 0\n",
    "    return cv\n",
    "\n",
    "\n",
    "\n",
    "def calculate_auc(values):\n",
    "    normalized_values = values / np.nanmax(values)\n",
    "    return auc(time_axis, normalized_values)\n",
    "\n",
    "\n",
    "def calculate_k_vm(vm):\n",
    "    return np.array([rgb_to_cmyk(frame) for frame in vm])\n",
    "\n",
    "def calculate_K_delta(values):\n",
    "    return values[0] - values[-1]\n",
    "\n",
    "\n",
    "\n",
    "def process_video(fids, output_path, export_video=True, return_matrix=False, apply_stabilization=True, fps=60, interval=10, crop_x0=500, crop_x1=3850, crop_y0=500, crop_y1=5000):\n",
    "    Full_frame = slice(crop_x0, crop_x1), slice(crop_y0, crop_y1) #TODO change this w.r.t. your video size, this crops the video, basically crop this to your hot plate or whatever\n",
    "\n",
    "\n",
    "    size = (Full_frame[1].stop-Full_frame[1].start, Full_frame[0].stop-Full_frame[0].start)\n",
    "    fourcc = cv2.VideoWriter_fourcc(*'XVID')\n",
    "    \n",
    "    out = None\n",
    "    if export_video:\n",
    "        out = cv2.VideoWriter(output_path, fourcc, fps, frameSize=size)\n",
    "        if not out.isOpened():\n",
    "            print(\"Error: VideoWriter didn't open\")\n",
    "            return None\n",
    "    \n",
    "    ref_img = cv2.imread(fids[0])\n",
    "    if ref_img is None:\n",
    "        print(\"Error reading reference image\")\n",
    "        return None\n",
    "    \n",
    "    ref_img_cropped = ref_img[Full_frame]\n",
    "    ref_img_cropped_gray = cv2.cvtColor(ref_img_cropped, cv2.COLOR_BGR2GRAY)\n",
    "    \n",
    "    mask = np.ones_like(ref_img_cropped_gray, dtype=np.uint8)\n",
    "    \n",
    "    num_frames = (len(fids) - 1) // interval\n",
    "    video_matrix = np.zeros((num_frames, size[1], size[0], 3), dtype=np.uint8) if return_matrix else None\n",
    "    \n",
    "    counter = 0\n",
    "    \n",
    "    for n in tqdm(range(0, num_frames * interval, interval)):\n",
    "        temp = cv2.imread(fids[n])\n",
    "        if temp is None:\n",
    "            continue\n",
    "        \n",
    "        cropped_temp = temp[Full_frame]\n",
    "        cropped_temp_gray = cv2.cvtColor(cropped_temp, cv2.COLOR_BGR2GRAY)\n",
    "        \n",
    "        if apply_stabilization:\n",
    "            warp_matrix = np.eye(2, 3, dtype=np.float32)\n",
    "            criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 5000, 1e-6)\n",
    "            _, warp_matrix = cv2.findTransformECC(ref_img_cropped_gray, cropped_temp_gray, warp_matrix, cv2.MOTION_EUCLIDEAN, criteria, inputMask=mask, gaussFiltSize=5)\n",
    "            aligned_temp = cv2.warpAffine(cropped_temp, warp_matrix, (size[0], size[1]), flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP)\n",
    "            \n",
    "        # if apply_stabilization:\n",
    "        #     warp_matrix = np.eye(2, 3, dtype=np.float32)  # Initialize the warp matrix as an identity matrix\n",
    "        #     criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 5000, 1e-6)\n",
    "            \n",
    "        #     # Use cv2.MOTION_AFFINE instead of cv2.MOTION_EUCLIDEAN\n",
    "        #     _, warp_matrix = cv2.findTransformECC(ref_img_cropped_gray, cropped_temp_gray, warp_matrix, cv2.MOTION_AFFINE, criteria, inputMask=mask, gaussFiltSize=5)\n",
    "            \n",
    "        #     aligned_temp = cv2.warpAffine(cropped_temp, warp_matrix, (size[0], size[1]), flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP)\n",
    "\n",
    "        else:\n",
    "            aligned_temp = cropped_temp\n",
    "        \n",
    "        if export_video:\n",
    "            out.write(aligned_temp)\n",
    "        \n",
    "        if return_matrix:\n",
    "            video_matrix[counter] = aligned_temp\n",
    "            counter += 1\n",
    "    \n",
    "    if export_video:\n",
    "        out.release()\n",
    "    \n",
    "    cv2.destroyAllWindows()\n",
    "    \n",
    "    if return_matrix:\n",
    "        return video_matrix\n",
    "    \n",
    "\n",
    "\n",
    "def apply_perspective_transform(frame, M, cols, rows):\n",
    "    return cv2.warpPerspective(frame, M, (cols, rows))\n",
    "\n",
    "\n",
    "def fit_data_and_find_x_intercept(K_mean, frame_axis, start=2000, stop=3000, sigma=15):\n",
    "    if K_mean is None or np.isnan(K_mean).all():\n",
    "        return None\n",
    "\n",
    "    K_mean_array = np.array(K_mean)\n",
    "    K_mean_array = K_mean_array[~np.isnan(K_mean_array)]\n",
    "    # sigma = 15\n",
    "\n",
    "    smoothed_norm_k_mean_values = gaussian_filter(K_mean_array, sigma) / np.nanmax(K_mean_array)\n",
    "    \n",
    "    frame_axis_array = np.array(frame_axis)[:len(K_mean_array)]\n",
    "    mask = (frame_axis_array >= start) & (frame_axis_array <= stop)\n",
    "    fit_x = frame_axis_array[mask]\n",
    "    fit_y = 1 - smoothed_norm_k_mean_values[mask]\n",
    "\n",
    "    coefficients = np.polyfit(fit_x, fit_y, 1)\n",
    "    slope, intercept = coefficients\n",
    "\n",
    "    x_intercept = -intercept / slope if slope != 0 else None\n",
    "\n",
    "    return x_intercept"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get your data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_folder_0 = '/Users/deniz/Library/CloudStorage/SynologyDrive-MyDrive/Characterization_RAWDATA/Go_Pro/20230925_Inorganic_T33/combined'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fids_0 = []\n",
    "for f in frgm.listdir(image_folder_0, display = False):\n",
    "    if 'control' in f and '_rgb' not in f:\n",
    "        continue\n",
    "    fids_0.append(f)\n",
    "fids_0 = natsorted(fids_0)\n",
    "print(len(fids_0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first figure out the cropping to the full spread of samples range\n",
    "temp = cv2.imread(fids_0[1000])\n",
    "plt.imshow(temp[500:3850, 500:5000])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_path = '/Users/deniz/Library/CloudStorage/SynologyDrive-MyDrive/Inorganic_PSK/GoPro_analysis/video_export/T33/T33.avi'\n",
    "video_matrix = process_video(fids_0, output_path, fps=30, interval=100, export_video=False, return_matrix=True, apply_stabilization=True, crop_x0=500, crop_x1=3850, crop_y0=500, crop_y1=5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "video_matrix.size*video_matrix.itemsize/1e9 # GB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import h5py\n",
    "\n",
    "# # Save the video_matrix as h5 if over 4GB\n",
    "# with h5py.File('T33.h5', 'w') as f:\n",
    "#     f.create_dataset('video_matrix', data=video_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import h5py\n",
    "\n",
    "# # Load the video_matrix\n",
    "# with h5py.File('b29-b31_VM.h5', 'r') as f:\n",
    "#     video_matrix = f['video_matrix'][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import pickle as pkl\n",
    "\n",
    "# with open('b29-b31_VM.pkl', 'wb') as f:\n",
    "#     pkl.dump(video_matrix, f)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import pickle as pkl\n",
    "\n",
    "# with open('b29-b31_VM.pkl', 'rb') as f:\n",
    "#     video_matrix = pkl.load(f)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import cv2\n",
    "\n",
    "# def write_video_from_matrix(video_matrix, output_path, fps=60):\n",
    "#     # Get the shape of the video_matrix\n",
    "#     num_frames, height, width, _ = video_matrix.shape\n",
    "    \n",
    "#     # Define the codec and create VideoWriter object\n",
    "#     fourcc = cv2.VideoWriter_fourcc(*'XVID')\n",
    "#     out = cv2.VideoWriter(output_path, fourcc, fps, (width, height))\n",
    "    \n",
    "#     if not out.isOpened():\n",
    "#         print(\"Error: VideoWriter didn't open\")\n",
    "#         return\n",
    "    \n",
    "#     # Write each frame to the video\n",
    "#     for i in range(num_frames):\n",
    "#         out.write(video_matrix[i])\n",
    "    \n",
    "#     # Release the VideoWriter\n",
    "#     out.release()\n",
    "#     print(f\"Video written to {output_path}\")\n",
    "\n",
    "# # Example usage\n",
    "# output_path = '/Users/deniz/Library/CloudStorage/SynologyDrive-MyDrive/Tandem_project/GoPro_Analysis/20230814_B29/output_video/B29_r2_stabilized_transformed.avi'\n",
    "# write_video_from_matrix(video_matrix, output_path)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "# Define the points in the original image (corners of the original image)\n",
    "rows, cols = video_matrix.shape[1:3]\n",
    "pts1 = np.float32([[769, 79], [3796, 35], [387, 3295], [4328, 3274]])  # top left, top right, bottom left, bottom right\n",
    "\n",
    "# Define the magnitude of the transformation (e.g., 0.05 for a 5% stretch at the top)\n",
    "magnitude = 0.01\n",
    "\n",
    "# Define where those points will be in the transformed image (stretching the top)\n",
    "pts2 = np.float32([[0, -rows * magnitude], [cols - 1, -rows * magnitude], [0, rows - 1], [cols - 1, rows - 1]])\n",
    "\n",
    "# Compute the perspective transformation matrix\n",
    "M = cv2.getPerspectiveTransform(pts1, pts2)\n",
    "\n",
    "# Get the first original image from video_matrix\n",
    "original_image = video_matrix[0, :, :, 1]\n",
    "\n",
    "# Apply the transformation to the first frame to see what it looks like\n",
    "transformed_image = apply_perspective_transform(original_image, M, cols, rows)\n",
    "\n",
    "# Display the original image using Plotly Express\n",
    "fig_original = px.imshow(original_image)\n",
    "fig_original.show()\n",
    "\n",
    "# Display the transformed image using Plotly Express\n",
    "fig_transformed = px.imshow(transformed_image)\n",
    "fig_transformed.show()\n",
    "\n",
    "# # If the transformed image looks good, uncomment the following lines to apply the transformation to the whole video_matrix\n",
    "for i in tqdm(range(video_matrix.shape[0])):\n",
    "    video_matrix[i, :, :, :] = apply_perspective_transform(video_matrix[i, :, :, :], M, cols, rows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "# Create a blank matrix with placeholders for 64 samples in a format similar to center_points_0\n",
    "\n",
    "\n",
    "# Create an 8x8 blank matrix with placeholders (None) for each sample\n",
    "center_points_0 = np.array([\n",
    "    [(1612, 389), (2747, 389)],\n",
    "    [(1612, 1221), (2747, 1221)],\n",
    "    [(1612, 2053), (2747, 2053)],\n",
    "    [(1612, 2896), (2747, 2896)],\n",
    "\n",
    "])\n",
    "\n",
    "# sample_names_0 = np.arange(63, -1, -1).reshape(8, 8)\n",
    "\n",
    "sample_names_0 = np.array([ [10, 110],\n",
    "                            [5, 15],\n",
    "                            [1, 11],\n",
    "                            [0, 10],])\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display the original or transformed image\n",
    "ax, fig = plt.subplots(figsize=(5, 5))\n",
    "plt.imshow(transformed_image, cmap='gray')\n",
    "\n",
    "# Loop through sample_names_0 and center_points_0 to annotate the image\n",
    "for i in range(center_points_0.shape[0]):\n",
    "    for j in range(center_points_0.shape[1]):\n",
    "        x, y = center_points_0[i, j]\n",
    "        sample_name = sample_names_0[i, j]\n",
    "        \n",
    "        # Annotate the image\n",
    "        plt.text(x, y, str(sample_name), color='white', fontsize=10, ha='center', va='center')\n",
    "\n",
    "# Show the image with annotations\n",
    "plt.axis('off')\n",
    "# make it tight\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# drop noisy section of data if not reasonably repairable\n",
    "\n",
    "# video_matrix = np.concatenate([video_matrix[:2700], video_matrix[2901:]], axis=0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the window size (70 pixels in x and y)\n",
    "window_size = 250\n",
    "num_frames = video_matrix.shape[0]\n",
    "\n",
    "# Create a slice matrix to hold the slices\n",
    "slice_shape = (2 * window_size, 2 * window_size, video_matrix.shape[3])\n",
    "slice_matrix = np.empty((*sample_names_0.shape, *slice_shape))\n",
    "\n",
    "# Extract the slices from the video_matrix using the center points\n",
    "for i in range(sample_names_0.shape[0]):\n",
    "    for j in range(sample_names_0.shape[1]):\n",
    "        x, y = center_points_0[i, j]\n",
    "\n",
    "        # Boundary checks\n",
    "        x_min = max(0, x - window_size)\n",
    "        x_max = min(video_matrix.shape[2], x + window_size)\n",
    "        y_min = max(0, y - window_size)\n",
    "        y_max = min(video_matrix.shape[1], y + window_size)\n",
    "\n",
    "        temp_slice = video_matrix[0, y_min:y_max, x_min:x_max, :]\n",
    "\n",
    "        # Padding in case the slice is smaller than the window\n",
    "        pad_x_min = window_size - (x - x_min)\n",
    "        pad_x_max = window_size + (x_max - x)\n",
    "        pad_y_min = window_size - (y - y_min)\n",
    "        pad_y_max = window_size + (y_max - y)\n",
    "\n",
    "        slice_matrix[i, j, pad_y_min:pad_y_max, pad_x_min:pad_x_max, :] = temp_slice\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "# Create a list to hold the data for each sample\n",
    "samples_data = []\n",
    "\n",
    "# Iterate through the samples in the slice_matrix\n",
    "for i in tqdm(range(slice_matrix.shape[0])):\n",
    "    for j in range(slice_matrix.shape[1]):\n",
    "        # Extract the sample name\n",
    "        sample_name = sample_names_0[i, j]\n",
    "\n",
    "        # Create a new video matrix for the current sample\n",
    "        sample_video_matrix = np.empty((num_frames, *slice_shape))\n",
    "\n",
    "        # Extract the slices for the current sample from the video_matrix\n",
    "        for frame in range(num_frames):\n",
    "            x, y = center_points_0[i, j]\n",
    "            if np.isnan(x) or np.isnan(y):\n",
    "                continue # Skip if x or y is NaN\n",
    "            x, y = int(x), int(y)\n",
    "\n",
    "            # Check that the coordinates are within valid bounds\n",
    "            if x - window_size < 0 or x + window_size >= video_matrix.shape[2] or y - window_size < 0 or y + window_size >= video_matrix.shape[1]:\n",
    "                continue # Skip if coordinates are too close to the edges\n",
    "\n",
    "            sample_video_matrix[frame] = video_matrix[frame, y-window_size:y+window_size, x-window_size:x+window_size, :]\n",
    "\n",
    "        # Append the current sample's data to the list\n",
    "        samples_data.append((sample_name, sample_video_matrix))\n",
    "\n",
    "\n",
    "\n",
    "# Create a DataFrame from the list of samples data\n",
    "all_samples_df = pd.DataFrame(samples_data, columns=['sample', 'VM'])\n",
    "\n",
    "# You now have a DataFrame (all_samples_df) that contains the sample names and video matrices for all samples\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_samples_df['K_VM'] = all_samples_df['VM'].apply(calculate_k_vm)\n",
    "\n",
    "all_samples_df['RGB_mean'] = all_samples_df['VM'].apply(\n",
    "    lambda vm: [np.mean(vm[frame, :, :, :], axis=(0, 1)) for frame in range(vm.shape[0])]\n",
    ")\n",
    "time_axis = [frame * 120*100 / 60**2 for frame in range(len(all_samples_df.loc[0, 'RGB_mean']))] # in hours: 120 is the camera picture interval in seconds * 10 for the interval of 10\n",
    "frame_axis = [frame for frame in range(len(all_samples_df.loc[0, 'RGB_mean']))]\n",
    "\n",
    "# Apply the function to the 'VM' column of the DataFrame\n",
    "all_samples_df['K_mean'] = all_samples_df['VM'].apply(calculate_K_mean)\n",
    "# all_samples_df['cv'] = all_samples_df['RGB_mean'].apply(calculate_cv)\n",
    "all_samples_df['auc'] = all_samples_df['K_mean'].apply(calculate_auc)\n",
    "\n",
    "all_samples_df['K_delta'] = all_samples_df['K_mean'].apply(calculate_K_delta)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_samples_df['x_intercept'] = all_samples_df['K_mean'].apply(\n",
    "    fit_data_and_find_x_intercept, \n",
    "    args=(frame_axis, 5, 18, 2)\n",
    ")\n",
    "# fitting needs work "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Specify the sample values you're interested in\n",
    "# samples_to_nan = [47, 46, 45, 36, 17]\n",
    "\n",
    "# # Set all other column values to NaN for these samples\n",
    "# all_samples_df.loc[all_samples_df['sample'].isin(samples_to_nan), 'VM':'K_delta'] = np.nan\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.ticker import MultipleLocator\n",
    "from scipy.ndimage import gaussian_filter\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "filtered_df = all_samples_df.dropna(subset=['K_mean'])\n",
    "\n",
    "sigma = 1\n",
    "\n",
    "# Your existing code for plotting\n",
    "for index, row in filtered_df.iterrows():\n",
    "    k_mean_values = row['K_mean']\n",
    "    smoothed_norm_k_mean_values = gaussian_filter(k_mean_values, sigma) / np.nanmax(gaussian_filter(k_mean_values, sigma)) - (np.nanmin(gaussian_filter(k_mean_values, sigma)) / np.nanmax(gaussian_filter(k_mean_values, sigma)))\n",
    "    smoothed_norm_k_mean_values = smoothed_norm_k_mean_values/np.nanmax(smoothed_norm_k_mean_values)\n",
    "\n",
    "    # smoothed_norm_k_mean_values = gaussian_filter(k_mean_values, sigma) / np.nanmax(k_mean_values)\n",
    "    plt.plot(frame_axis, 1-smoothed_norm_k_mean_values, alpha=1)\n",
    "\n",
    "\n",
    "# add a vertical line\n",
    "# plt.axvline(x=2700, color='k', linestyle='--', linewidth=1)\n",
    "# plt.axvline(x=2900, color='k', linestyle='--', linewidth=1)\n",
    "# Set major ticks at intervals of 24\n",
    "ax = plt.gca()\n",
    "# ax.xaxis.set_major_locator(MultipleLocator(24))\n",
    "\n",
    "# Set minor ticks at intervals of 6\n",
    "# ax.xaxis.set_minor_locator(MultipleLocator(6))\n",
    "\n",
    "plt.xlim(0, 20)\n",
    "# plt.ylim(0, .6)\n",
    "plt.xlabel('Frame')\n",
    "plt.ylabel('1-cmyk_AverageK_95ci')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "#the next resonable step would be to fit the linear portion of these 'S' shape curves to get an onset value to further dimensionality reduce.\n",
    "## Note samples that did not degrade by end of test look odd, e.g. max value will not be at end of measurement, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extra Scratch work"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from matplotlib.ticker import MultipleLocator\n",
    "# from scipy.ndimage import gaussian_filter\n",
    "# import matplotlib.pyplot as plt\n",
    "# import numpy as np\n",
    "\n",
    "# filtered_df = all_samples_df.dropna(subset=['K_mean'])\n",
    "\n",
    "# sigma = 15\n",
    "\n",
    "# # Your existing code for plotting\n",
    "# for index, row in filtered_df.iterrows():\n",
    "#     k_mean_values = row['K_mean']\n",
    "#     smoothed_norm_k_mean_values = gaussian_filter(k_mean_values, sigma) / np.nanmax(k_mean_values)\n",
    "    \n",
    "#     # Gather data for fitting between frames 2000 and 3000\n",
    "#     frame_axis_np = np.array(frame_axis)\n",
    "#     mask = (frame_axis_np >= 2000) & (frame_axis_np <= 3000)\n",
    "#     fit_x = frame_axis_np[mask]\n",
    "#     fit_y = (1 - smoothed_norm_k_mean_values)[mask]\n",
    "\n",
    "#     # Fit the data with a line\n",
    "#     coefficients = np.polyfit(fit_x, fit_y, 1)\n",
    "#     slope, intercept = coefficients\n",
    "\n",
    "#     # Calculate the x-intercept (where y=0)\n",
    "#     x_intercept = -intercept / slope if slope != 0 else None\n",
    "    \n",
    "#     # Plot the fitted line\n",
    "#     fit_line_y = slope * fit_x + intercept\n",
    "#     plt.plot(fit_x, fit_line_y, '-', label=f'Fit for curve {index}')\n",
    "\n",
    "# # Customize plot\n",
    "# plt.xlabel('Frame')\n",
    "# plt.ylabel('1-cmyk_AverageK_95ci')\n",
    "# # plt.xlim(0,4000)\n",
    "# # plt.legend()\n",
    "# plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from matplotlib.ticker import MultipleLocator\n",
    "# from scipy.ndimage import gaussian_filter\n",
    "# import matplotlib.pyplot as plt\n",
    "# import numpy as np\n",
    "\n",
    "# filtered_df = all_samples_df.dropna(subset=['K_mean'])\n",
    "# samples_to_keep = sample_names_0[sample_names_0 >= 48]\n",
    "# samples_to_keep = samples_to_keep.flatten()\n",
    "# filtered_df = all_samples_df[all_samples_df['sample'].isin(samples_to_keep)]\n",
    "\n",
    "# sigma = 15\n",
    "\n",
    "# # Your existing code for plotting\n",
    "# for index, row in filtered_df.iterrows():\n",
    "#     k_mean_values = row['K_mean']\n",
    "#     smoothed_norm_k_mean_values = gaussian_filter(k_mean_values, sigma) / np.nanmax(k_mean_values)\n",
    "    \n",
    "#     # Color lines differently based on sample name\n",
    "#     if row['sample'] in range(48, 57):\n",
    "#         plt.plot(time_axis, smoothed_norm_k_mean_values, alpha=1, color='blue')\n",
    "#     else:\n",
    "#         plt.plot(time_axis, smoothed_norm_k_mean_values, alpha=1, color='red')\n",
    "\n",
    "# # Add a vertical line\n",
    "# # plt.axvline(x=2700, color='k', linestyle='--', linewidth=1)\n",
    "# # plt.axvline(x=2900, color='k', linestyle='--', linewidth=1)\n",
    "# # Set major ticks at intervals of 24\n",
    "# ax = plt.gca()\n",
    "# ax.xaxis.set_major_locator(MultipleLocator(24))\n",
    "\n",
    "# # Set minor ticks at intervals of 6\n",
    "# ax.xaxis.set_minor_locator(MultipleLocator(6))\n",
    "\n",
    "# plt.xlim(0, 140)\n",
    "# plt.xlabel('Time (Hrs)')\n",
    "# plt.ylabel('cmyk_AverageK_95ci')\n",
    "# plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import matplotlib.pyplot as plt\n",
    "\n",
    "# # Create a time axis (in hours) assuming each frame is 20 minutes apart\n",
    "# time_axis = [frame * 20 / 60 for frame in range(len(all_samples_df.loc[0, 'RGB_mean']))] \n",
    "\n",
    "# # Iterate through the DataFrame and plot the RGB values for each sample\n",
    "# for index, row in all_samples_df.iterrows():\n",
    "#     # Extract the RGB_mean values for the current sample\n",
    "#     rgb_mean_values = row['RGB_mean']\n",
    "\n",
    "#     # Extract the R, G, and B values\n",
    "#     r_values = [rgb[0] for rgb in rgb_mean_values]\n",
    "#     g_values = [rgb[1] for rgb in rgb_mean_values]\n",
    "#     b_values = [rgb[2] for rgb in rgb_mean_values]\n",
    "\n",
    "#     # Plot the R, G, and B values over time\n",
    "#     plt.plot(time_axis, r_values, label='R - ', color='red', alpha=0.5)\n",
    "#     plt.plot(time_axis, g_values, label='G - ', color='green', alpha=0.5)\n",
    "#     plt.plot(time_axis, b_values, label='B - ', color='blue', alpha=0.5)\n",
    "\n",
    "# plt.xlabel('Time (hours)')\n",
    "# plt.ylabel('RGB Value')\n",
    "# # plt.title('RGB Values Over Time for All Samples')\n",
    "# # plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') # Position the legend outside the plot\n",
    "# plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extra: To make a color bar out of images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# all_samples_df['deg_image_vector'] = ''\n",
    "# all_samples_df['time_vector'] = ''\n",
    "\n",
    "# # test_frame_0 = df['img'][0][20]\n",
    "# # test_frame_1 = df['img'][0][21]\n",
    "\n",
    "# # test_combined_frame = np.concatenate((test_frame_0, test_frame_1), axis=1)fi\n",
    "# index = 0\n",
    "# exp = 1\n",
    "# for index in range(7, 8,1):\n",
    "\n",
    "#     fig, ax = plt.subplots(figsize=(200,120), constrained_layout=False) #(200,120)\n",
    "#     start = 5\n",
    "\n",
    "#     for frame in range(start, 200, 1):\n",
    "#         test_frame = np.stack(np.round(all_samples_df['K_mean'][index][frame],3)*np.ones((500,500,1)))\n",
    "#         test_frame_array = np.round(all_samples_df['K_mean'][index][frame],3)\n",
    "#         sampled_time_array = time_axis[frame]\n",
    "\n",
    "#         if frame == start:\n",
    "#             test_combined_frame = test_frame\n",
    "#             test_combined_frame_array = test_frame_array\n",
    "#             sampled_time_combined_array = sampled_time_array\n",
    "\n",
    "\n",
    "#         if frame != start:\n",
    "#             test_combined_frame = np.concatenate((test_combined_frame, test_frame), axis=1)\n",
    "#             test_combined_frame_array = np.append(test_combined_frame_array, test_frame_array)\n",
    "#             sampled_time_combined_array = np.append(sampled_time_combined_array, sampled_time_array)\n",
    "\n",
    "\n",
    "\n",
    "#     ax.axis('off')\n",
    "#     # plt.tight_layout()\n",
    "\n",
    "#     print(all_samples_df['sample'][index])\n",
    "#     ax.imshow(test_combined_frame, cmap='Greys', vmin=0, vmax=1)\n",
    "#     # plt.imsave('test_combined_frame.png', test_combined_frame)\n",
    "#     print(index)\n",
    "#     # plt.title(df['name'][index])\n",
    "#     plt.show()\n",
    "\n",
    "#     all_samples_df['deg_image_vector'][index] = test_combined_frame_array\n",
    "#     all_samples_df['time_vector'][index] = sampled_time_combined_array\n",
    "# del test_combined_frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.7.9 ('ds')",
   "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.7.9"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "c86dfbc17468201bf110f641e56a1df81e3c9871b3a7a52a78db9d9bb9a33d6a"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}