1.2.HDMI_HOUGH_Detection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "aebe2113-cc61-428c-84f3-733e97e39bfb",
   "metadata": {},
   "source": [
    "# HOUGH Detection\n",
    "\n",
    "In this notebook, we convert the data from each tile into an image, process this image and detect the fiducial circles.\n",
    "\n",
    "The fiducial circles are detected using a HOUGH transform and we identify 2 narrow sizes to increase accuracy where possible.\n",
    "\n",
    "The centroids of the fiducials are used to correct for the global positions of the tiles across a Nova-ST chip."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d4fc2f8-4478-45ab-9beb-fc94e7bfed40",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pysam\n",
    "import collections as c\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import pickle\n",
    "import numpy as np\n",
    "import cv2 as cv\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import time\n",
    "import os\n",
    "import time\n",
    "import datetime\n",
    "from multiprocessing import Pool\n",
    "import PIL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59c13b9b-24dd-458f-9654-9053fa39c885",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from scipy.sparse import coo_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c22e5da0-d2bf-41d8-b552-a0b563d27389",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "TILE_DIR = 'HDMI_Tiles_Data' # The output data folder from 1.1.HDMI_Processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7eb025bb-09a6-4595-acce-f7dcfec2f89b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def process_tile(info, return_images=False):\n",
    "    tiles_dir, tile_id = info\n",
    "    try:\n",
    "        with open(f'{tiles_dir}/{tile_id}_barcodes.pickle', 'rb') as tile_fh:\n",
    "            tile_data = pickle.load(tile_fh)\n",
    "    except FileNotFoundError:\n",
    "        return (f'{tile_id}', None, None)\n",
    "\n",
    "    # Tile data structure: key: barcode, values: (x, y)\n",
    "    xs = [int(x[0]) for x in tile_data.values()]\n",
    "    ys = [int(x[1]) for x in tile_data.values()]\n",
    "    max_x = max(xs)\n",
    "    max_y = max(ys)\n",
    "\n",
    "    # Create a matrix of of zeros in the size of the tile\n",
    "    arr_orig = np.zeros((max_x + 1, max_y + 1) , dtype=np.uint8)\n",
    "    for x, y in zip(xs, ys):\n",
    "        arr_orig[x, y] = 255\n",
    "\n",
    "    # Shrink the data by binning for easier processing. Lower numbers = higher resolution but more computation time\n",
    "    IMG_BINSIZE = 25\n",
    "    m = arr_orig.shape[0]\n",
    "    n = arr_orig.shape[1]\n",
    "    trim_arr = arr_orig[:(m // IMG_BINSIZE) * IMG_BINSIZE, :(n // IMG_BINSIZE) * IMG_BINSIZE]\n",
    "    arr_re = trim_arr.reshape(m // IMG_BINSIZE, IMG_BINSIZE, n // IMG_BINSIZE, IMG_BINSIZE)\n",
    "    arr_binned = arr_re.sum(3).sum(1)\n",
    "    arr = arr_binned // (arr_binned.max() / 255)\n",
    "\n",
    "    arr = arr.astype(np.uint8)\n",
    "\n",
    "    # Process the image to make fiducials stand out\n",
    "    inverted = cv.threshold(arr, 0, 255, cv.THRESH_BINARY_INV)[1]\n",
    "    nonoise = cv.fastNlMeansDenoising(inverted, None, h=100)\n",
    "    thresholded2 = cv.threshold(nonoise, 128, 255, cv.THRESH_BINARY)[1]\n",
    "\n",
    "    # Run two rounds of Hough transforms to identify the nested fiducials \n",
    "    circles1 = cv.HoughCircles(thresholded2,cv.HOUGH_GRADIENT,1,20,\n",
    "                                param1=50,param2=20,minRadius=40,maxRadius=80)\n",
    "    circles2 = cv.HoughCircles(thresholded2,cv.HOUGH_GRADIENT,1,20,\n",
    "                                param1=50,param2=20,minRadius=15,maxRadius=30)\n",
    "\n",
    "    # If we can't identify all of the circes from a set, be a bit more lenient on detection\n",
    "    if circles1.shape[1] != 8 or circles2.shape[1] !=8:        \n",
    "        circles1 = cv.HoughCircles(thresholded2,cv.HOUGH_GRADIENT,1,20,\n",
    "                                    param1=50,param2=30,minRadius=30,maxRadius=70)\n",
    "        \n",
    "        \n",
    "        circles2 = cv.HoughCircles(thresholded2,cv.HOUGH_GRADIENT,1,20,\n",
    "                                    param1=50,param2=30,minRadius=10,maxRadius=30)\n",
    "        if circles1.shape[1] != 8 or circles2.shape[1] !=8:\n",
    "            return (f'{tile_id}', None, None)\n",
    "\n",
    "    centers1 = circles1[0, :, :2]\n",
    "    centroid1 = centers1.mean(0)\n",
    "    centers2 = circles2[0, :, :2]\n",
    "    centroid2 = centers2.mean(0)\n",
    "\n",
    "    if not return_images:\n",
    "        return (f'{tile_id}', None, {\n",
    "            'circles1': circles1,\n",
    "            'circles2': circles2,\n",
    "            'centroid1': centroid1,\n",
    "            'centroid2': centroid2\n",
    "        })\n",
    "        \n",
    "    # Recreate images if they are wanted for QC\n",
    "    cimg = cv.cvtColor(thresholded2, cv.COLOR_GRAY2BGR)\n",
    "    circles = np.uint16(np.around(circles1))\n",
    "    for i in circles[0,:]:\n",
    "        cv.circle(cimg,(i[0],i[1]),i[2],(0,255,0),2)\n",
    "        cv.circle(cimg,(i[0],i[1]),2,(0,0,255),3)\n",
    "    circles = np.uint16(np.around(circles2))\n",
    "    for i in circles[0,:]:\n",
    "        cv.circle(cimg,(i[0],i[1]),i[2],(255,0,0),2)\n",
    "        cv.circle(cimg,(i[0],i[1]),2,(0,255,0),3)\n",
    "    \n",
    "    if circles1[0].shape[0] == 8 and circles2[0].shape[0] == 8:\n",
    "        cv.line(cimg, (int(centroid1[0]) - 20, int(centroid1[1]) -20), (int(centroid1[0]) + 20, int(centroid1[1]) + 20), (0, 255, 0), 2)\n",
    "        cv.line(cimg, (int(centroid2[0]) - 20, int(centroid2[1]) +20), (int(centroid2[0]) + 20, int(centroid2[1]) - 20), (255, 0, 0), 2)\n",
    "    elif circles1[0].shape[0] == 8:\n",
    "        cv.circle(cimg, (int(centroid1[0]), int(centroid1[1])), 4, (0, 255, 0))\n",
    "    elif circles2[0].shape[0] == 8:\n",
    "        cv.circle(cimg, (int(centroid2[0]), int(centroid2[1])), 4, (255, 0, 0))\n",
    "\n",
    "    return (f'{tile_id}', {\n",
    "        'raw': arr,\n",
    "        'thresh1': inverted,\n",
    "        'denoise': nonoise,\n",
    "        'thresh2': thresholded2,\n",
    "        'anno': cimg\n",
    "    }, {\n",
    "        'circles1': circles1,\n",
    "        'circles2': circles2,\n",
    "        'centroid1': centroid1,\n",
    "        'centroid2': centroid2\n",
    "    })\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ac35bc8",
   "metadata": {},
   "outputs": [],
   "source": [
    "tiles = os.listdir(TILE_DIR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ed7cddc",
   "metadata": {},
   "outputs": [],
   "source": [
    "tiles = ['_'.join(x.split('_')[:2]) for x in tiles if not 'subset' in x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c2fc8a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_tile_ids = []\n",
    "for tile in tiles:\n",
    "    all_tile_ids.append((TILE_DIR, tile))\n",
    "\n",
    "# Multiprocesses tiles in parallel, adjust processes depending on the system configuration\n",
    "with Pool(processes=72) as pool:\n",
    "    results = pool.map(process_tile, all_tile_ids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3cbdd82",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Did we miss any tiles? There are tiles where fiducial detection failed\n",
    "for lane in range(1, 5):\n",
    "    for surface in range(1, 3):\n",
    "        for swath in range(1, 7):\n",
    "            for tile in range(1, 79):\n",
    "                final_tileid = f'{lane}_{surface}{swath}{tile:02}'\n",
    "                if not os.path.exists(f'{TILE_DIR}/{final_tileid}_barcodes.pickle'):\n",
    "                    print(f'No file for lane {lane} - {final_tileid}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e573f3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ci_info = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7df3bcc3-d1a6-4056-ba7e-bb421a753283",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for x in results:\n",
    "    ci_info[x[0]] = x[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59d3c685-4122-490a-bce5-47902466c213",
   "metadata": {},
   "source": [
    "Save the results to disk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d633ede0-9b42-42bd-9141-4c7d73955b33",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "with open(f'{TILE_DIR}/circle_info.pickle', 'wb') as fh:\n",
    "    pickle.dump(ci_info, fh)"
   ]
  }
 ],
 "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.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}