pycromanager_optoTracking.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9107532d-8682-4607-bef0-f12b13112603",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pycromanager import Bridge, Acquisition, multi_d_acquisition_events, Dataset\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import os\n",
    "import shutil\n",
    "import math\n",
    "import tifffile\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import cm\n",
    "from skimage import io, draw, morphology, filters, measure, transform, segmentation\n",
    "from scipy import spatial\n",
    "from scipy.ndimage import distance_transform_edt, center_of_mass\n",
    "from dmdUtils import *\n",
    "from imgUtils import *\n",
    "\n",
    "from IPython.display import clear_output\n",
    "\n",
    "try:\n",
    "    s1.close()\n",
    "    s1 = initialize_led_connection()\n",
    "except:\n",
    "    s1 = initialize_led_connection()\n",
    "\n",
    "bridge = Bridge()\n",
    "core = bridge.get_core()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a3c7478-767d-4818-af9e-79bea17312a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# manually focus the dmd at the TIRF plane\n",
    "\n",
    "core.set_config(\"-TIRF\", \"REFLECT-488LP\")\n",
    "core.set_exposure(200)\n",
    "v_650(s1, 40000)\n",
    "mickey_mouse()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f82da4a4-65da-4dbe-9ce1-2a6bbdf59137",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Snap mickey mouse image for record keeping\n",
    "\n",
    "core.snap_image()\n",
    "tagged_image = core.get_tagged_image()\n",
    "h, w = tagged_image.tags['Height'], tagged_image.tags['Width']\n",
    "pixels = np.reshape(tagged_image.pix, newshape=[h, w])\n",
    "\n",
    "metadata_dict = tagged_image.tags\n",
    "metadata_dict['timestamp'] = str(time.time())\n",
    "metadata_dict['stage_x_pos'] = str(core.get_x_position())\n",
    "metadata_dict['stage_y_pos'] = str(core.get_y_position())\n",
    "core.set_config(\"Focus\", \"TIPFSOffset\")\n",
    "core.get_position()\n",
    "metadata_dict['ti_zpos'] = core.get_position()\n",
    "core.set_config(\"Focus\", \"TIZDrive\")\n",
    "core.get_position()\n",
    "metadata_dict['pfs_value'] = core.get_position()\n",
    "\n",
    "txt_timestamp = time.strftime('%Y%m%d_%H%M%S', time.localtime(time.time()))\n",
    "\n",
    "image_dest = f\"C:\\\\Data\\\\Jason\\\\mickey_records\\\\mickey_{txt_timestamp}.tif\"\n",
    "tifffile.imwrite(image_dest, pixels, metadata = metadata_dict)\n",
    "    \n",
    "pattern = np.zeros((600,800), dtype = np.uint8)\n",
    "np_to_dmd(pattern)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7388b05-9a77-483e-ac1d-1d82eb3c23d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# figure out coordinate mapping between dmd pixels and image pixels\n",
    "# this can be used to estimate the affine transform\n",
    "\n",
    "print(core.get_position())\n",
    "core.set_config(\"-TIRF\", \"REFLECT-488LP\")\n",
    "core.set_exposure(200)\n",
    "\n",
    "v_650(s1, 60000)\n",
    "\n",
    "pattern = np.zeros((600,800), dtype = np.uint8)\n",
    "np_to_dmd(pattern)\n",
    "time.sleep(1)\n",
    "core.snap_image()\n",
    "tagged_image = core.get_tagged_image()\n",
    "h, w = tagged_image.tags['Height'], tagged_image.tags['Width']\n",
    "bg_pixels = np.reshape(tagged_image.pix, newshape=[h, w])\n",
    "\n",
    "recorded_peaks = []\n",
    "for dmd_coord_pair in [[150, 150], [450, 450], [450, 150], [150, 450]]:\n",
    "\n",
    "    pattern = np.zeros((600,800), dtype = np.uint8)\n",
    "    x, y = dmd_coord_pair\n",
    "    pattern[y-10:y+11,x-10:x+11] = 1\n",
    "    np_to_dmd(pattern)\n",
    "    time.sleep(1)\n",
    "    \n",
    "    core.snap_image()\n",
    "    tagged_image = core.get_tagged_image()\n",
    "    h, w = tagged_image.tags['Height'], tagged_image.tags['Width']\n",
    "    pixels = np.reshape(tagged_image.pix, newshape=[h, w])\n",
    "    \n",
    "    \n",
    "    smth_pix = filters.gaussian(pixels.astype(float) - bg_pixels.astype(float), 5, preserve_range = True)\n",
    "    ypeak, xpeak = np.unravel_index(np.argmax(smth_pix), shape = (h,w))\n",
    "    recorded_peaks.append(xpeak)\n",
    "    recorded_peaks.append(ypeak)\n",
    "    \n",
    "    plt.figure(figsize = (8,8))\n",
    "    plt.imshow(pixels)\n",
    "    plt.scatter(xpeak, ypeak, s = 500, facecolor = 'None', ec = 'k')\n",
    "    print(xpeak, ypeak)\n",
    "    plt.show()\n",
    "    \n",
    "recorded_values = open('C:\\\\Data\\\\Jason\\\\SCRIPTS\\\\calibration_coordinates.txt', 'w')\n",
    "for val in recorded_peaks:\n",
    "    recorded_values.write(str(val) + '\\n')\n",
    "recorded_values.close()\n",
    "\n",
    "pattern = np.zeros((600,800), dtype = np.uint8)\n",
    "np_to_dmd(pattern)\n",
    "v_650(s1, 0)\n",
    "\n",
    "\n",
    "\n",
    "x1, y1, x2, y2, x3, y3, x4, y4 = 150, 150, 450, 450, 450, 150, 150, 450\n",
    "src = np.array([x1, y1, x2, y2, x3, y3, x4, y4]).reshape((4,2))\n",
    "\n",
    "recorded_values = open('C:\\\\Data\\\\Jason\\\\SCRIPTS\\\\calibration_coordinates.txt', 'r')\n",
    "x1, y1, x2, y2, x3, y3, x4, y4 = [int(i.strip('\\n')) for i in recorded_values.readlines()]\n",
    "dst = np.array([x1, y1, x2, y2, x3, y3, x4, y4]).reshape((4,2))\n",
    "\n",
    "dmd_tform = tf.estimate_transform('affine', src, dst)\n",
    "dmd_tform_mat = dmd_tform.params\n",
    "\n",
    "dmd_affine_as_string = ''\n",
    "dmd_nums = [str(i) for i in dmd_tform_mat.ravel()[:6]]\n",
    "for i in dmd_nums[:6]:\n",
    "    dmd_affine_as_string += i + ';'\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea2ea3a5-0be9-45a0-a1d2-a63bf2f86c41",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.average(recorded_peaks[::2]))\n",
    "print(np.average(recorded_peaks[1::2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99bc3f75-aa9a-45c3-a125-31cc2734c528",
   "metadata": {},
   "outputs": [],
   "source": [
    "def binarize_img_thresh(input_img, thresh):\n",
    "    #corrected_image = (input_img.astype(float) - dark_avg.astype(float))/flatfield_ch1.astype(float)\n",
    "    corrected_image = input_img.astype(float)\n",
    "    smth_corrected_image = filters.gaussian(corrected_image, 2, preserve_range = True)\n",
    "    threshold = thresh\n",
    "    binary_image = smth_corrected_image > threshold\n",
    "    binary_image = segmentation.clear_border(binary_image)\n",
    "    binary_image = morphology.remove_small_holes(binary_image, 200)\n",
    "    binary_image = morphology.remove_small_objects(binary_image, 500)\n",
    "    binary_image = morphology.binary_closing(binary_image, selem = morphology.disk(5))\n",
    "    return(binary_image)\n",
    "\n",
    "def run_detection_protocol(mode, channel_list, channel_id):\n",
    "    core.set_config(mode, channel_list[channel_id])\n",
    "    core.set_exposure(200)\n",
    "    core.snap_image()\n",
    "    tagged_image = core.get_tagged_image()\n",
    "    h, w = tagged_image.tags['Height'], tagged_image.tags['Width']\n",
    "    img = np.reshape(tagged_image.pix, newshape=[h, w])\n",
    "    return(img)\n",
    "\n",
    "def test_cells_present(input_img, thresh):\n",
    "   \n",
    "    binary_image = binarize_img_thresh(input_img, thresh)\n",
    "    \n",
    "    # ignore the 100 pixels arouund the border because they are not reliable\n",
    "    D = 100\n",
    "    binary_image[:D,:] = False\n",
    "    binary_image[:,:D] = False\n",
    "    binary_image[-D:,:] = False\n",
    "    binary_image[:,-D:] = False\n",
    "    \n",
    "    lblim = measure.label(binary_image)\n",
    "    reg = measure.regionprops(lblim, intensity_image = input_img)\n",
    "   \n",
    "    area_data = []\n",
    "    polarization_data = []\n",
    "    aspect_ratio_data = []\n",
    "   \n",
    "    filtered_bin_img = np.zeros_like(binary_image)\n",
    "    \n",
    "   \n",
    "    if len(reg) > 0:\n",
    "        \n",
    "        for r in reg:\n",
    "            yg, xg = r.centroid\n",
    "            #print(xg, yg, r.area)\n",
    "            yw, xw = r.weighted_centroid\n",
    "            pol_magnitude = np.hypot(yg - yw, xg - xw)\n",
    "            polarization_data.append(pol_magnitude)\n",
    "            area_data.append(r.area)\n",
    "            if r.minor_axis_length > 0:\n",
    "                aspect_ratio_data.append(r.major_axis_length/r.minor_axis_length)\n",
    "            else:\n",
    "                aspect_ratio_data.append(np.nan)\n",
    "        area_data = np.array(area_data)\n",
    "        aspect_ratio_data = np.array(aspect_ratio_data)\n",
    "        pass_test_conditions = (area_data > 1000) & (aspect_ratio_data < 5)\n",
    "\n",
    "        if np.sum(pass_test_conditions) > 0:\n",
    "            for i in np.where(pass_test_conditions)[0] + 1:\n",
    "                filtered_bin_img[lblim == i] = 1\n",
    "           \n",
    "            return(True, filtered_bin_img)\n",
    "             \n",
    "        else:\n",
    "            return(False, filtered_bin_img)\n",
    "    else:\n",
    "        return(False, filtered_bin_img)\n",
    "   \n",
    "   \n",
    "def find_promising_cell(ch0_stack, ch1_stack, thresh):\n",
    "    # find a consistently migrating cell\n",
    "   \n",
    "    candidate_found = False\n",
    "    target_x = None\n",
    "    target_y = None\n",
    "   \n",
    "    # parameters for identifying candidates:\n",
    "    min_breathing_room = 0 # minimum distance (in pixels) between nearest cell (or edge of image)\n",
    "    min_cell_persistence = 0.8 # minimum required \"chemotactic index\"\n",
    "    min_cell_speed = 3 #minimum average frame-to-frame displacement of centroid (in pixels/frame)\n",
    "    min_cell_area = 1500  # minimum cell size needed in pixels\n",
    "    max_cell_area = 15000 # maximum cell size in pixels\n",
    "    min_cell_ch0_sig = 600 # minimum average signal (prestimulus) in channel 0\n",
    "    min_cell_ch1_sig = 0 # minimum average signal (prestimulus) in channel 1\n",
    "\n",
    "   \n",
    "    # set up a distance map to figure out how close a given coordinate is to the edge of the image\n",
    "    edge = np.zeros((1024, 1024))\n",
    "    edge[1:-1,1:-1] = 1\n",
    "    edge_dist = distance_transform_edt(edge)\n",
    "    \n",
    "    bin_stack = []\n",
    "    for t in range(len(ch1_stack)):\n",
    "        input_img = ch1_stack[t]\n",
    "        binary_img = binarize_img_thresh(input_img, thresh)\n",
    "        bin_stack.append(binary_img)\n",
    "        \n",
    "    bin_stack = np.array(bin_stack)\n",
    "    lbl_stack, num_potential_cells = measure.label(bin_stack, return_num = True)\n",
    "    \n",
    "    plt.figure(figsize = (8,8))\n",
    "    plt.subplot(aspect = 'equal')\n",
    "    cmap = cm.get_cmap('jet')\n",
    "    cvec = [cmap(i) for i in np.linspace(0,1,len(bin_stack))]\n",
    "    for i in range(len(bin_stack)):\n",
    "        plt.contour(bin_stack[i], levels = [False], colors = [cvec[i]])\n",
    "    plt.show()\n",
    "    \n",
    "    # if at least one cell is found, move on with analysis\n",
    "    if np.sum(bin_stack[-1]) > 0:\n",
    "\n",
    "        cell_position_x = []\n",
    "        cell_position_y = []\n",
    "        ch0_mean_signal = []\n",
    "        ch1_mean_signal = []\n",
    "        cell_areas = []\n",
    "        cell_speed = []\n",
    "        cell_persistence = []\n",
    "        breathing_room = []\n",
    "\n",
    "        for n in range(1, num_potential_cells + 1):\n",
    "\n",
    "            if n in lbl_stack[-1]:\n",
    "                isolated_cell = lbl_stack == n\n",
    "                X_positions = []\n",
    "                Y_positions = []\n",
    "                sizes = []\n",
    "                for t in range(len(isolated_cell)):\n",
    "                    y, x = center_of_mass(isolated_cell[t])\n",
    "                    X_positions.append(x)\n",
    "                    Y_positions.append(y)\n",
    "                    sizes.append(np.sum(isolated_cell[t]))\n",
    "\n",
    "                cell_position_x.append(X_positions[-1])\n",
    "                cell_position_y.append(Y_positions[-1])\n",
    "                cell_areas.append(sizes[-1])\n",
    "                ch0_mean_signal.append(np.average(ch0_stack[-1][isolated_cell[-1]]))\n",
    "                ch1_mean_signal.append(np.average(ch1_stack[-1][isolated_cell[-1]]))\n",
    "\n",
    "                path_length = np.sum(np.hypot(np.diff(X_positions), np.diff(Y_positions)))\n",
    "                avg_speed = np.average(np.hypot(np.diff(X_positions), np.diff(Y_positions)))\n",
    "                end_to_end_displacement = np.hypot(X_positions[-1] - X_positions[0], Y_positions[-1] - Y_positions[0])\n",
    "                persistence_index = end_to_end_displacement/path_length\n",
    "\n",
    "                cell_speed.append(avg_speed)\n",
    "                cell_persistence.append(persistence_index)\n",
    "\n",
    "        # if at least 2 cells are present:\n",
    "        if len(cell_position_x) > 1:\n",
    "            Dmat = spatial.distance.squareform(spatial.distance.pdist(np.vstack([cell_position_x, cell_position_y]).T))\n",
    "            for i in range(len(Dmat)):\n",
    "                closest_cell = np.min(np.delete(Dmat[i], i))\n",
    "                closest_edge = edge_dist[int(cell_position_y[i]), int(cell_position_x[i])]\n",
    "\n",
    "                breathing_room.append(np.min([closest_cell, closest_edge]))\n",
    "               \n",
    "        # if only 1 cell is present\n",
    "        else:\n",
    "            breathing_room.append(edge_dist[int(cell_position_y[0]), int(cell_position_x[0])])\n",
    "\n",
    "        # convert to numpy arrays for easy element-wise calculations\n",
    "        cell_position_x = np.array(cell_position_x)\n",
    "        cell_position_y = np.array(cell_position_y)\n",
    "        ch0_mean_signal = np.array(ch0_mean_signal)\n",
    "        ch1_mean_signal = np.array(ch1_mean_signal)\n",
    "        cell_areas = np.array(cell_areas)\n",
    "        cell_speed = np.array(cell_speed)\n",
    "        cell_persistence = np.array(cell_persistence)\n",
    "        breathing_room = np.array(breathing_room)\n",
    "\n",
    "        celldata = {'x': cell_position_x, 'y': cell_position_y, 'ch0_F': ch0_mean_signal, 'ch1_F': ch1_mean_signal, 'Area': cell_areas, 'Speed': cell_speed, 'Persistence': cell_persistence, 'BreathingRoom': breathing_room}\n",
    "        df = pd.DataFrame.from_dict(celldata)\n",
    "        print(df)\n",
    "\n",
    "        # test all identified cells to see if they meet minimum requirements as outlined above\n",
    "        # if more than one is found, identify the one with the most \"breathing room\"\n",
    "        cond1 = (breathing_room > min_breathing_room)\n",
    "        cond2 = (cell_persistence > min_cell_persistence)\n",
    "        cond3 = (cell_speed > min_cell_speed)\n",
    "        cond4 = (cell_areas > min_cell_area)\n",
    "        cond5 = (cell_areas < max_cell_area)\n",
    "        cond6 = (ch0_mean_signal > min_cell_ch0_sig)\n",
    "        cond7 = (ch1_mean_signal > min_cell_ch1_sig)\n",
    "        min_requirements = (cond1 & cond2 & cond3 & cond4 & cond5 & cond6 & cond7)\n",
    "\n",
    "       \n",
    "        if np.sum(min_requirements) > 0:\n",
    "            target_x = cell_position_x[min_requirements][np.argmax(breathing_room[min_requirements])]\n",
    "            target_y = cell_position_y[min_requirements][np.argmax(breathing_room[min_requirements])]\n",
    "            candidate_found = True\n",
    "        else:\n",
    "            candidate_found = False\n",
    "\n",
    "    else:\n",
    "        candidate_found = False\n",
    "       \n",
    "    return(candidate_found, target_x, target_y)\n",
    "\n",
    "\n",
    "\n",
    "def run_migration_test(mode, channel_list, num_imgs, wait_interval):\n",
    "\n",
    "    ch0_stack = []\n",
    "    ch1_stack = []\n",
    "\n",
    "    for repeat in range(num_imgs):\n",
    "        start_time = time.time()\n",
    "\n",
    "        core.set_config(mode, channel_list[0])\n",
    "        core.set_exposure(200)\n",
    "        core.snap_image()\n",
    "        tagged_image = core.get_tagged_image()\n",
    "        h, w = tagged_image.tags['Height'], tagged_image.tags['Width']\n",
    "        pixels = np.reshape(tagged_image.pix, newshape=[h, w])\n",
    "        ch0_stack.append(pixels)\n",
    "\n",
    "        core.set_config(mode, channel_list[1])\n",
    "        core.set_exposure(200)\n",
    "        core.snap_image()\n",
    "        tagged_image = core.get_tagged_image()\n",
    "        h, w = tagged_image.tags['Height'], tagged_image.tags['Width']\n",
    "        pixels = np.reshape(tagged_image.pix, newshape=[h, w])\n",
    "        ch1_stack.append(pixels)\n",
    "\n",
    "        end_time = time.time()\n",
    "        \n",
    "        remaining_time = wait_interval - (end_time - start_time)\n",
    "        if remaining_time > 0:\n",
    "            time.sleep(remaining_time)\n",
    "       \n",
    "    return(np.array(ch0_stack), np.array(ch1_stack))\n",
    "\n",
    "\n",
    "def make_data_directory(folder_base_name):\n",
    "    # create a storage folder & metadata file\n",
    "    folder_number = len([i for i in os.listdir(\"C:\\\\Data\\\\Jason\\\\python_controlled_acq\\\\\") if i[:-4] == folder_base_name])\n",
    "    folder_path = \"C:\\\\Data\\\\Jason\\\\python_controlled_acq\\\\\" + folder_base_name + '_' + str(folder_number).zfill(3)\n",
    "    while not os.path.exists(folder_path):\n",
    "        try:\n",
    "            os.mkdir(folder_path)\n",
    "        except:\n",
    "            print('Directory conflict, trying another folder name')\n",
    "            folder_number += 1\n",
    "            folder_path = \"C:\\\\Data\\\\Jason\\\\python_controlled_acq\\\\\" + folder_base_name + '_' + str(folder_number).zfill(3)\n",
    "    return(folder_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "864ce7e1-2672-48fe-832c-b2ee961a3e6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('x', core.get_x_position())\n",
    "print('y', core.get_y_position())\n",
    "print('z', core.get_position())\n",
    "print('exposure_ms', core.get_exposure())\n",
    "affine_matrix = np.vstack([np.array([float(i) for i in core.get_pixel_size_affine_as_string().split(';')]).reshape(2, 3), [0,0,1]])\n",
    "\n",
    "\n",
    "aft_im_to_stage = transform.AffineTransform(affine_matrix)\n",
    "aft_dmd_to_im = transform.AffineTransform(dmd_tform_mat)\n",
    "\n",
    "print(affine_matrix)\n",
    "print(dmd_tform_mat)\n",
    "\n",
    "dmd_scalar_factor = np.average(aft_dmd_to_im.scale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b88280c4-490f-4667-92cc-e26ea8614fbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# collect current coords\n",
    "# MANUALLY NAVIGATE TO CENTER OF WELL\n",
    "center_y = core.get_y_position()\n",
    "center_x = core.get_x_position()\n",
    "\n",
    "print(center_x, center_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ea8adf2-36ac-4f25-b777-de783e469539",
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate a circular-shaped set of imaging stage coordinates \n",
    "R = 7\n",
    "\n",
    "rr, cc = draw.disk((0, 0), R)\n",
    "unique_rs = np.sort(np.unique(rr))\n",
    "\n",
    "rr_snake = []\n",
    "cc_snake = []\n",
    "\n",
    "flip = -1\n",
    "for i in unique_rs:\n",
    "    for j in cc[rr == i]:\n",
    "        rr_snake.append(i)\n",
    "        cc_snake.append(flip * j)\n",
    "    flip *= -1\n",
    "\n",
    "relative_y = 200 * np.array(cc_snake)\n",
    "relative_x = 200 * np.array(rr_snake)\n",
    "absolute_x = center_x + np.array(relative_x)\n",
    "absolute_y = center_y + np.array(relative_y)\n",
    "position_num_ids = range(len(absolute_x))\n",
    "print(f'Set up to image {len(position_num_ids)} points\\n')\n",
    "\n",
    "\n",
    "plt.figure(figsize = (8,8))\n",
    "plt.subplot(aspect = 'equal')\n",
    "plt.plot(absolute_x, absolute_y, color = 'k', ls = '--')\n",
    "plt.scatter(absolute_x, absolute_y)\n",
    "\n",
    "theta = np.linspace(0, 2 * math.pi, 100)\n",
    "plt.plot(3200 * np.cos(theta) + center_x, 3200 * np.sin(theta) + center_y)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df1bb738-f265-45d5-950c-8c002d452e5f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "DMD_CENTER_X = 740\n",
    "DMD_CENTER_Y = 600\n",
    "\n",
    "\n",
    "pattern = np.zeros((600,800), dtype = np.uint8)\n",
    "np_to_dmd(pattern)\n",
    "v_650(s1, 0)\n",
    "\n",
    "NUM_REPEATS = 3\n",
    "\n",
    "num_images_to_acquire = 60\n",
    "wait_time_interval = 6\n",
    "folder_base_name = '220314_KWC_35-52-97_simpleStimScreen'\n",
    "mode = '-TIRF'\n",
    "channel_list = [\"561-488LP\", \"640-488LP\", \"REFLECT-488LP\"]\n",
    "\n",
    "zlevel = 0\n",
    "track_cell = True\n",
    "segmentation_channel = 1\n",
    "intensity_thresh = 650\n",
    "\n",
    "# set up ROI with side length 2*D\n",
    "D = 200\n",
    "ROI_xmin = DMD_CENTER_X - D\n",
    "ROI_xmax = DMD_CENTER_X + D\n",
    "ROI_ymin = DMD_CENTER_Y - D\n",
    "ROI_ymax = DMD_CENTER_Y + D\n",
    "\n",
    "\n",
    "\n",
    "for repeat in range(NUM_REPEATS):\n",
    "    \n",
    "    for i in position_num_ids:\n",
    "        x_R = np.nan\n",
    "        y_R = np.nan\n",
    "        x_L = np.nan\n",
    "        y_L = np.nan\n",
    "\n",
    "        # move to position\n",
    "        print('Position' + str(i + 1).zfill(3) + ' of ' + str(len(position_num_ids)).zfill(3))\n",
    "        core.set_xy_position(np.double(absolute_x[i]), np.double(absolute_y[i]))\n",
    "        core.wait_for_system()\n",
    "        time.sleep(1)\n",
    "\n",
    "        # run cell detection protocol\n",
    "        detection_image = run_detection_protocol(mode, channel_list, segmentation_channel)\n",
    "        cells_present_condition, segmented_proposed_cells = test_cells_present(detection_image, intensity_thresh)\n",
    "\n",
    "        # display info for user\n",
    "        plt.figure(figsize = (16,8))\n",
    "\n",
    "        # current_position\n",
    "        plt.subplot(121, aspect = 'equal')\n",
    "        plt.plot(absolute_x, absolute_y, color = 'k', ls = '--', zorder = 0)\n",
    "        plt.scatter(absolute_x, absolute_y, zorder = 1)\n",
    "        plt.scatter(absolute_x[i], absolute_y[i], s = 200, lw = 2, ec = 'k', facecolor = 'None')\n",
    "        \n",
    "        # image and segmentation\n",
    "        plt.subplot(122)\n",
    "        plt.imshow(detection_image)\n",
    "        plt.colorbar()\n",
    "        plt.contour(segmented_proposed_cells, levels = [False], colors = 'w', linestyles = '--')\n",
    "        plt.title('Cell Presence Test')\n",
    "        plt.show()\n",
    "\n",
    "        # if preliminary test was passed, run secondary test\n",
    "        if cells_present_condition:\n",
    "\n",
    "            print('Cell likely present, starting migration test')\n",
    "            \n",
    "            \n",
    "            ch0_data, ch1_data = run_migration_test(mode, channel_list, num_imgs = 5, wait_interval = 5) #DEFAULT IS 5, 5\n",
    "            continue_imaging, x_in_img, y_in_img = find_promising_cell(ch0_data, ch1_data, intensity_thresh)\n",
    "\n",
    "            plt.figure(figsize = (8,8))\n",
    "            plt.imshow(np.max(ch1_data, axis = 0))\n",
    "            plt.scatter(x_in_img, y_in_img, facecolor = 'None', ec = 'w', s = 500)\n",
    "            plt.title('Max Projection over Time')\n",
    "            plt.show()\n",
    "\n",
    "            stage_x_sequence = []\n",
    "            stage_y_sequence = []\n",
    "            \n",
    "            if continue_imaging:\n",
    "                print('Migration test passed! Moving on to opto experiment')\n",
    "                \n",
    "                # move promising cell roughly to center of dmd field\n",
    "                move_from = aft_im_to_stage([x_in_img, y_in_img])[0]\n",
    "                move_to = aft_im_to_stage([DMD_CENTER_X, DMD_CENTER_Y])[0]\n",
    "                DX, DY = move_from - move_to\n",
    "                new_x = core.get_x_position() + DX\n",
    "                new_y = core.get_y_position() + DY\n",
    "                core.set_xy_position(new_x, new_y)\n",
    "                \n",
    "                # determine assay type\n",
    "                brightness = 50000 #np.random.choice([0, 15000, 30000, 45000])\n",
    "                variable_options = ['none', 'global', 'front', 'side', 'back', 'lateral90', 'lateral45']\n",
    "                variable_choice = len(os.listdir(\"C:\\\\Data\\\\Jason\\\\python_controlled_acq\\\\\"))%len(variable_options)\n",
    "                assay_choice = variable_options[variable_choice]\n",
    "                print(f'Running {assay_choice}-type assay')\n",
    "                folder_path = make_data_directory(folder_base_name + '_' + assay_choice + '_' + str(brightness) + '_Intensity')\n",
    "                \n",
    "                ch0_stack = []\n",
    "                bin_stack = []\n",
    "                lght_bin_stack = []\n",
    "\n",
    "                experiment_start_time = time.time()\n",
    "                \n",
    "                lost_cell_count = 0\n",
    "                \n",
    "                continuous_imaging = True\n",
    "                for t in range(num_images_to_acquire):\n",
    "                    \n",
    "                    if lost_cell_count <= 3:\n",
    "\n",
    "                        start_time = time.time()\n",
    "                        description_dict = {'status': 'observing'}\n",
    "\n",
    "                        # collect images in sequence (separated from other processing for the sake of speed)\n",
    "                        tiff_collection = []\n",
    "                        md_collection = []\n",
    "                        for channel in range(len(channel_list)):\n",
    "                            core.set_config(mode, channel_list[channel])\n",
    "                            core.set_exposure(200)\n",
    "                            core.snap_image()\n",
    "                            tagged_image = core.get_tagged_image()\n",
    "                            metadata_dict = tagged_image.tags\n",
    "                            metadata_dict['elapsed_time_s'] = str(time.time() - experiment_start_time)\n",
    "                            metadata_dict['stage_x_pos'] = str(core.get_x_position())\n",
    "                            metadata_dict['stage_y_pos'] = str(core.get_y_position())\n",
    "                            metadata_dict['dmd_affine_transform'] = dmd_affine_as_string\n",
    "                            metadata_dict['roi_xmin_xmax_ymin_ymax'] = str(ROI_xmin) + ';' + str(ROI_xmax) + ';' + str(ROI_ymin) + ';' + str(ROI_ymax)\n",
    "                            tiff_collection.append(tagged_image)\n",
    "                            md_collection.append(metadata_dict)\n",
    "\n",
    "                        # now save those images and run analysis\n",
    "                        for channel in range(len(channel_list)):\n",
    "                            img = tiff_collection[channel]\n",
    "                            metadata_dict = img.tags\n",
    "                            h, w = metadata_dict['Height'], metadata_dict['Width']\n",
    "                            pixels = np.reshape(img.pix, newshape=[h, w])\n",
    "\n",
    "                            image_dest = folder_path + f\"\\\\img_channel{str(channel).zfill(3)}_position{str(i).zfill(3)}_time{str(t).zfill(9)}_z{str(zlevel).zfill(3)}.tif\"\n",
    "                            tifffile.imwrite(image_dest, pixels[ROI_ymin:ROI_ymax, ROI_xmin:ROI_xmax], metadata = md_collection[channel])\n",
    "\n",
    "                            if channel == segmentation_channel:\n",
    "\n",
    "                                ch0_stack.append(pixels)\n",
    "\n",
    "                                binary_image = binarize_img_thresh(pixels, intensity_thresh)\n",
    "\n",
    "                                isolated_cell = get_closest_cell(binary_image, DMD_CENTER_X, DMD_CENTER_Y)\n",
    "\n",
    "                                if np.sum(isolated_cell) > 0:\n",
    "                                    possible_y_center, possible_x_center = center_of_mass(isolated_cell)\n",
    "                                    dist_from_target = np.hypot(possible_y_center - DMD_CENTER_Y, possible_x_center - DMD_CENTER_X)\n",
    "                                    if dist_from_target > 100:\n",
    "                                        isolated_cell = np.zeros_like(isolated_cell)\n",
    "\n",
    "                                bin_stack.append(isolated_cell)\n",
    "\n",
    "                                # update optogenetic stimulation status\n",
    "                                if (t >= 10):\n",
    "\n",
    "                                    # calculate angle for opto stim\n",
    "                                    if track_cell:\n",
    "                                        dx_track = np.average(np.diff(stage_x_sequence[:10]))\n",
    "                                        dy_track = -np.average(np.diff(stage_y_sequence[:10])) # y sign is flipped between image and stage coords\n",
    "                                        theta = np.arctan2(dy_track, dx_track)\n",
    "\n",
    "                                    else:\n",
    "                                        theta = get_move_polar_angle(bin_stack[7], bin_stack[10])\n",
    "\n",
    "                                    # set brightness\n",
    "                                    brightness = 50000\n",
    "                                    pattern = np.zeros((600,800))\n",
    "\n",
    "                                    # check if cell is being tracked and update points if possible\n",
    "                                    if np.sum(bin_stack[-1]) > 0:\n",
    "                                        \n",
    "                                        if assay_choice == 'none':\n",
    "                                            pass\n",
    "                                        \n",
    "                                        elif assay_choice == 'global':\n",
    "                                            pattern = np.ones((600,800))\n",
    "                                            \n",
    "                                        elif assay_choice == 'front':\n",
    "                                            x_R, y_R, x_L, y_L = get_lateral_points(bin_stack[-1], theta - math.pi/2) # -math.pi/2 for front, +math.pi/2 for back\n",
    "                                            x_in_dmd, y_in_dmd = dmd_from_img(x_L, y_L)\n",
    "                                            rr_disk, cc_disk = draw.disk((y_in_dmd, x_in_dmd), 20, shape = (600, 800))\n",
    "                                            pattern[rr_disk, cc_disk] = 1\n",
    "                                            \n",
    "                                        elif assay_choice == 'side':\n",
    "                                            x_R, y_R, x_L, y_L = get_lateral_points(bin_stack[-1], theta) # -math.pi/2 for front, +math.pi/2 for back\n",
    "                                            x_in_dmd, y_in_dmd = dmd_from_img(x_L, y_L)\n",
    "                                            rr_disk, cc_disk = draw.disk((y_in_dmd, x_in_dmd), 20, shape = (600, 800))\n",
    "                                            pattern[rr_disk, cc_disk] = 1\n",
    "                                            \n",
    "                                        elif assay_choice == 'back':\n",
    "                                            x_R, y_R, x_L, y_L = get_lateral_points(bin_stack[-1], theta + math.pi/2) # -math.pi/2 for front, +math.pi/2 for back\n",
    "                                            x_in_dmd, y_in_dmd = dmd_from_img(x_L, y_L)\n",
    "                                            rr_disk, cc_disk = draw.disk((y_in_dmd, x_in_dmd), 20, shape = (600, 800))\n",
    "                                            pattern[rr_disk, cc_disk] = 1\n",
    "                                            \n",
    "                                        elif assay_choice == 'lateral90':\n",
    "                                            x_R, y_R, x_L, y_L = get_lateral_points(bin_stack[-1], theta) # -math.pi/2 for front, +math.pi/2 for back\n",
    "                                            x_in_dmd, y_in_dmd = dmd_from_img(x_L, y_L)\n",
    "                                            rr_disk, cc_disk = draw.disk((y_in_dmd, x_in_dmd), 20, shape = (600, 800))\n",
    "                                            pattern[rr_disk, cc_disk] = 1\n",
    "                                            x_in_dmd, y_in_dmd = dmd_from_img(x_R, y_R)\n",
    "                                            rr_disk, cc_disk = draw.disk((y_in_dmd, x_in_dmd), 20, shape = (600, 800))\n",
    "                                            pattern[rr_disk, cc_disk] = 1\n",
    "                                            \n",
    "                                        elif assay_choice == 'lateral45':\n",
    "                                            x_R, y_R, x_L, y_L = get_lateral_points(bin_stack[-1], theta - math.pi/2 + math.pi/4) # -math.pi/2 for front, +math.pi/2 for back\n",
    "                                            x_in_dmd, y_in_dmd = dmd_from_img(x_L, y_L)\n",
    "                                            rr_disk, cc_disk = draw.disk((y_in_dmd, x_in_dmd), 20, shape = (600, 800))\n",
    "                                            pattern[rr_disk, cc_disk] = 1\n",
    "                                            x_R, y_R, x_L, y_L = get_lateral_points(bin_stack[-1], theta - math.pi/2 - math.pi/4)\n",
    "                                            x_in_dmd, y_in_dmd = dmd_from_img(x_L, y_L)\n",
    "                                            rr_disk, cc_disk = draw.disk((y_in_dmd, x_in_dmd), 20, shape = (600, 800))\n",
    "                                            pattern[rr_disk, cc_disk] = 1\n",
    "\n",
    "                                        np_to_dmd(pattern)\n",
    "                                        v_650(s1, brightness)\n",
    "                                        \n",
    "                                        status = 'tracking cell'\n",
    "                                        lost_cell_count = 0\n",
    "                                        \n",
    "                                    else:\n",
    "                                        status = 'lost cell'\n",
    "                                        lost_cell_count += 1\n",
    "                                        \n",
    "\n",
    "                                    R = np.round(dmd_scalar_factor * 20, 1)\n",
    "                                    et = np.round(time.time() - experiment_start_time, 3)\n",
    "                                    description_dict = {'status': status, 'assay': assay_choice, 'theta': theta, 'Intensity': str(brightness),'image_x': np.round(x_R, 2),'image_y': np.round(y_R, 2), 'image_r': R, 'elapsed_time': et}\n",
    "                                    print(description_dict)\n",
    "\n",
    "\n",
    "\n",
    "                        # save binary image and predicted dmd binary image\n",
    "                        bin_im_path = folder_path + f\"\\\\img_channel{str(len(channel_list)).zfill(3)}_position{str(i).zfill(3)}_time{str(t).zfill(9)}_z{str(zlevel).zfill(3)}.tif\"\n",
    "                        tifffile.imwrite(bin_im_path, isolated_cell.astype(np.int16)[ROI_ymin:ROI_ymax, ROI_xmin:ROI_xmax])\n",
    "\n",
    "                        stage_x_sequence.append(core.get_x_position())\n",
    "                        stage_y_sequence.append(core.get_y_position())\n",
    "\n",
    "                        dmd_im_path = folder_path + f\"\\\\img_channel{str(len(channel_list) + 1).zfill(3)}_position{str(i).zfill(3)}_time{str(t).zfill(9)}_z{str(zlevel).zfill(3)}.tif\"\n",
    "                        dmd_im = predict_DMD_image()\n",
    "                        tifffile.imwrite(dmd_im_path, dmd_im.astype(np.int16)[ROI_ymin:ROI_ymax, ROI_xmin:ROI_xmax], metadata = description_dict)\n",
    "                        lght_bin_stack.append(dmd_im)\n",
    "\n",
    "                        end_time = time.time()\n",
    "                        remaining_time = wait_time_interval - (end_time - start_time)\n",
    "                        if remaining_time > 0:\n",
    "                            time.sleep(remaining_time)\n",
    "\n",
    "                        # if tracking, move cell just before the loop starts over\n",
    "                        if track_cell:    \n",
    "\n",
    "                            current_cell_pos_y, current_cell_pos_x = center_of_mass(isolated_cell)\n",
    "                            move_from = aft_im_to_stage([current_cell_pos_x, current_cell_pos_y])[0]\n",
    "                            move_to = aft_im_to_stage([DMD_CENTER_X, DMD_CENTER_Y])[0]\n",
    "                            DX, DY = move_from - move_to\n",
    "                            if np.hypot(DX, DY) < 100:\n",
    "                                new_x = core.get_x_position() + DX\n",
    "                                new_y = core.get_y_position() + DY\n",
    "                                core.set_xy_position(new_x, new_y)\n",
    "                    else:\n",
    "                        print('Lost cell, removing data')\n",
    "                        try:\n",
    "                            shutil.rmtree(folder_path)\n",
    "                            continuous_imaging = False\n",
    "                        except OSError as e:\n",
    "                            print(\"Error: %s - %s.\" % (e.filename, e.strerror))\n",
    "                \n",
    "                if continuous_imaging:\n",
    "                    # verify presence of opto tools in blue and green channels\n",
    "                    core.set_config('-TIRF', '488')\n",
    "                    core.set_exposure(200)\n",
    "                    core.snap_image()\n",
    "                    tagged_image = core.get_tagged_image()\n",
    "                    metadata_dict = tagged_image.tags\n",
    "                    metadata_dict['elapsed_time_s'] = str(time.time() - experiment_start_time)\n",
    "                    metadata_dict['stage_x_pos'] = str(core.get_x_position())\n",
    "                    metadata_dict['stage_y_pos'] = str(core.get_y_position())\n",
    "                    metadata_dict['dmd_affine_transform'] = dmd_affine_as_string\n",
    "                    metadata_dict['roi_xmin_xmax_ymin_ymax'] = str(ROI_xmin) + ';' + str(ROI_xmax) + ';' + str(ROI_ymin) + ';' + str(ROI_ymax)\n",
    "                    h, w = metadata_dict['Height'], metadata_dict['Width']\n",
    "                    pixels = np.reshape(tagged_image.pix, newshape=[h, w])\n",
    "                    image_dest = folder_path + f\"\\\\img_channel{str(len(channel_list) + 2).zfill(3)}_position{str(i).zfill(3)}_time{str(t).zfill(9)}_z{str(zlevel).zfill(3)}.tif\"\n",
    "                    tifffile.imwrite(image_dest, pixels[ROI_ymin:ROI_ymax, ROI_xmin:ROI_xmax], metadata = md_collection[channel])\n",
    "\n",
    "                    core.set_config('-TIRF', '405')\n",
    "                    core.set_exposure(200)\n",
    "                    core.snap_image()\n",
    "                    tagged_image = core.get_tagged_image()\n",
    "                    metadata_dict = tagged_image.tags\n",
    "                    metadata_dict['elapsed_time_s'] = str(time.time() - experiment_start_time)\n",
    "                    metadata_dict['stage_x_pos'] = str(core.get_x_position())\n",
    "                    metadata_dict['stage_y_pos'] = str(core.get_y_position())\n",
    "                    metadata_dict['dmd_affine_transform'] = dmd_affine_as_string\n",
    "                    metadata_dict['roi_xmin_xmax_ymin_ymax'] = str(ROI_xmin) + ';' + str(ROI_xmax) + ';' + str(ROI_ymin) + ';' + str(ROI_ymax)\n",
    "                    h, w = metadata_dict['Height'], metadata_dict['Width']\n",
    "                    pixels = np.reshape(tagged_image.pix, newshape=[h, w])\n",
    "                    image_dest = folder_path + f\"\\\\img_channel{str(len(channel_list) + 3).zfill(3)}_position{str(i).zfill(3)}_time{str(t).zfill(9)}_z{str(zlevel).zfill(3)}.tif\"\n",
    "                    tifffile.imwrite(image_dest, pixels[ROI_ymin:ROI_ymax, ROI_xmin:ROI_xmax], metadata = md_collection[channel])\n",
    "\n",
    "\n",
    "                \n",
    "                # turn off lights and pattern for fresh restart\n",
    "                v_650(s1, 0)\n",
    "                pattern = np.zeros((600,800), dtype = np.uint8)\n",
    "                np_to_dmd(pattern)\n",
    "\n",
    "                # display for user\n",
    "                if track_cell:\n",
    "                    stage_x_sequence = np.array(stage_x_sequence)\n",
    "                    stage_y_sequence = np.array(stage_y_sequence)\n",
    "\n",
    "                    plt.figure(figsize = (6,6))\n",
    "                    plt.subplot(aspect = 'equal')\n",
    "                    plt.plot(stage_x_sequence - stage_x_sequence[0], stage_y_sequence - stage_y_sequence[0])\n",
    "                    plt.title(f'Position: {position_num_ids[i]}')\n",
    "                    plt.show()\n",
    "\n",
    "\n",
    "                    plt.figure(figsize = (8,8))\n",
    "                    plt.subplot(aspect = 'equal')\n",
    "                    cmap = cm.get_cmap('jet')\n",
    "                    cvec = [cmap(timestep) for timestep in np.linspace(0,1,len(bin_stack))]\n",
    "\n",
    "                    for i in range(len(bin_stack)):\n",
    "                        for ctrs in measure.find_contours(bin_stack[i], level = False):\n",
    "                            y_segment_imcoords, x_segment_imcoords = ctrs.T\n",
    "                            y_segment_stgcoords, x_segment_stgcoords = aft_im_to_stage(np.array([y_segment_imcoords, x_segment_imcoords]).T).T\n",
    "\n",
    "                            x_translated = -x_segment_stgcoords + stage_x_sequence[i] - stage_x_sequence[0]\n",
    "                            y_translated = -y_segment_stgcoords + stage_y_sequence[i] - stage_y_sequence[0]\n",
    "\n",
    "                            plt.plot(x_translated, y_translated, color = cvec[i])\n",
    "\n",
    "                        for ctrs in measure.find_contours(lght_bin_stack[i], level = False):\n",
    "                            y_segment_imcoords, x_segment_imcoords = ctrs.T\n",
    "                            y_segment_stgcoords, x_segment_stgcoords = aft_im_to_stage(np.array([y_segment_imcoords, x_segment_imcoords]).T).T\n",
    "\n",
    "                            x_translated = -x_segment_stgcoords + stage_x_sequence[i] - stage_x_sequence[0]\n",
    "                            y_translated = -y_segment_stgcoords + stage_y_sequence[i] - stage_y_sequence[0]\n",
    "\n",
    "                            plt.plot(x_translated, y_translated, color = cvec[i], ls = '--')\n",
    "\n",
    "                    plt.show()\n",
    "\n",
    "            else:\n",
    "                print('Likely no migrating cells here')\n",
    "        else:\n",
    "            print('Likely no cells here')\n",
    "\n",
    "        print()\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "734ac293-ea80-4cea-9c1e-9c20cd93fce5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}