find_bolt_and_ellipse.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%matplotlib notebook\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import cv2\n",
    "import numpy as np\n",
    "from skimage import color, img_as_ubyte\n",
    "import os, sys\n",
    "from scipy.spatial.distance import pdist\n",
    "import csv\n",
    "from skimage.transform import hough_ellipse\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#Find centroids of bolts near edge of circles in find_pmts() using findContours()\n",
    "\n",
    "#Output is a list of bolt locations for a given circle \n",
    "def ret_centres(image, params) :\n",
    "    out = []\n",
    "    x = params[0]\n",
    "    y = params[1]\n",
    "    r = params[2]\n",
    "    L = int(r*r_mult)\n",
    "\n",
    "    gray_image = image[y-L:y+L, x-L:x+L]\n",
    "    ret,thresh = cv2.threshold(gray_image,127,255,0)\n",
    "    gray_image = color.gray2rgb(img_as_ubyte(gray_image))\n",
    "    \n",
    "\n",
    "    # find contours in the binary image\n",
    "    contours, hierarchy = cv2.findContours(thresh,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)\n",
    "\n",
    "        #Loop over bolts and find centroid\n",
    "    for c in contours:\n",
    "    #moments\n",
    "        M = cv2.moments(c)\n",
    "        area = cv2.contourArea(c)\n",
    "        \n",
    "        if (area <=2) or (area>= 500) :  # skip ellipses smaller then 10x10\n",
    "            continue\n",
    "            \n",
    "           # calculate x,y coordinate of center, breaks if m00 = 0 output location in large image (offset)\n",
    "        elif M[\"m00\"] != 0:\n",
    "            cX = int(M[\"m10\"] / M[\"m00\"])\n",
    "            cY = int(M[\"m01\"] / M[\"m00\"])\n",
    "            out.append([cX + (x-L), cY + (y-L)])\n",
    "             #   print(area)\n",
    " #           cv2.circle(gray_image,(cX, cY),10,(0,255,255),2)\n",
    "        else:\n",
    "            continue\n",
    "    \n",
    "  #  cv2.imshow('img', gray_image)\n",
    "  #  cv2.waitKey(0)\n",
    "    return out\n",
    "\n",
    "\n",
    "#Orders bolts correctly\n",
    "def number_bolts(circle, bolt_ring) :\n",
    "    out = []\n",
    "    circ_x = circle[0]\n",
    "    circ_y = circle[1]\n",
    "    circ_r = circle[2]\n",
    "    L = int(circ_r*r_mult)  \n",
    "      \n",
    "    nodes = np.asarray(bolt_ring)\n",
    "    for i in np.linspace(0, 2*3.14159, num=60) : \n",
    "        search_x = int(circ_x + L*np.sin(i))\n",
    "        search_y = int(circ_y - L*np.cos(i))\n",
    "        dist = np.sum((nodes - [search_x,search_y])**2, axis=1)\n",
    "        closest = np.argmin(dist)\n",
    "            \n",
    "        if bolt_ring[closest] not in out :\n",
    "            out.append(bolt_ring[closest])\n",
    "    return out\n",
    "\n",
    "\n",
    "#Diameters of bolts\n",
    "def bolt_avg_dia(circle, bolt_ring) :\n",
    "    dias = []\n",
    "    num = []\n",
    "    lengths = []\n",
    "    num.append(len(bolt_ring))\n",
    "    \n",
    "    #debug\n",
    "\n",
    "    \n",
    "    \n",
    "    circ_x = circle[0]\n",
    "    circ_y = circle[1]\n",
    "    circ_r = circle[2]\n",
    "    L = int(circ_r*r_mult)          \n",
    "\n",
    "#    test_image = image[y-L:y+L, x-L:x+L]\n",
    "#    ret,thresh = cv2.threshold(test_image,127,255,0)\n",
    "#    test_image = color.gray2rgb(img_as_ubyte(test_image))\n",
    "    \n",
    "\n",
    "    centre = np.asarray([circ_x, circ_y]) \n",
    "    nodes = np.asarray(bolt_ring)\n",
    "    \n",
    "   # print(f'centre:{centre}')\n",
    "    #print(f'nodes:{nodes}')\n",
    "    #print(f'bolt_ring:{bolt_ring}')\n",
    "    for i in np.linspace(0, 3.14159, num=60) : \n",
    "        x_1 = int(circ_x + L*np.sin(i))\n",
    "        y_1 = int(circ_y - L*np.cos(i))\n",
    "\n",
    "        dist_1 = np.sum((nodes - [x_1,y_1])**2, axis=1)\n",
    "        \n",
    "            \n",
    "        dia_1 = bolt_ring[np.argmin(dist_1)]\n",
    "        dist_2 = np.sum((nodes - dia_1)**2, axis=1)\n",
    "        \n",
    "        dia_2 = bolt_ring[np.argmax(dist_2)]\n",
    "        \n",
    "        rad_1 = centre - dia_1\n",
    "        rad_2 = centre - dia_2\n",
    "        #print(i, np.argmin(dist_1), np.argmax(dist_2), dia_1)\n",
    "        angle = np.arccos(np.dot(rad_1, rad_2) / (np.linalg.norm(rad_1) * np.linalg.norm(rad_2)))\n",
    "        if (angle > 3) and ([dia_1, dia_2] not in dias) :\n",
    "            #print(dia_1[0])\n",
    "            dias.append([dia_1, dia_2])\n",
    "            lengths.append(np.sqrt((dia_1[0]-dia_2[0])**2 + (dia_1[1]-dia_2[1])**2))\n",
    "            #print(angle)        \n",
    "    return([dias, num, lengths])\n",
    "    \n",
    "#plot histo of results\n",
    "def eval_reco(odir, info, img) :\n",
    "    #dias = info[0]\n",
    "    num_bolts = info[0]\n",
    "    lengths = info[1]\n",
    "    \n",
    "  #  print(num_bolts[0])\n",
    "   ## print(lengths[0])\n",
    "    #ax1 = plt.subplot(212)\n",
    "    #ax1.margins(1)           # Default margin is 0.05, value 0 means fit\n",
    "    #ax1.imshow(img)\n",
    "\n",
    "    ax2 = plt.subplot(221)\n",
    "    ax2.hist(num_bolts, np.arange(10,25,1), edgecolor='black', linewidth=0.75)\n",
    "    ax2.set_xlabel('Number of bolts')\n",
    "    ax2.set_ylabel('PMTs')\n",
    "    ax2.margins(0.1)\n",
    "    \n",
    "    ax3 = plt.subplot(222)\n",
    "    ax3.hist(lengths, edgecolor='black', linewidth=0.75)\n",
    "    ax3.set_xlabel('Percentage uncertainty in diameter')\n",
    "    ax3.set_ylabel('Bolt Pairs')\n",
    "    ax3.margins(0.1)\n",
    "    #plt.show()\n",
    "    \n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Approx pmt locations\n",
    "\n",
    "def find_pmts(pmt_img, bolt_img, pmt_locs) :\n",
    "\n",
    "    out = []\n",
    "\n",
    "    #r = 210 \n",
    "    circles = []\n",
    "    ellipses = []\n",
    "   # print(pmt_locs)\n",
    "    r = int(0.5*np.average(np.diff( np.unique(np.array(pmt_locs['pmt_pos'])[:,1]).astype('float64') ))*r_mult) # MULTIPLIER SCALES R IT IS SLIGHTLY TOO BIG AS TAKEN FROM GAPS\n",
    "    for i in range(len(pmt_locs['pmt_pos'])) :\n",
    "            posn = pmt_locs['pmt_pos'][i]\n",
    "            label = pmt_locs['pmt_label'][i]\n",
    "            \n",
    "            if posn[0] > 300 and posn[1] > 300 :\n",
    "                \n",
    "\n",
    "                #print(f'PMT: {posn} {label}')\n",
    "                \n",
    "                one_pmt_img = pmt_img[int(posn[0])-r:int(posn[0])+r, int(posn[1])-r:int(posn[1])+r]\n",
    "             #   one_pmt_bolt_img_orig = bolt_img[int(posn[0])-r:int(posn[0])+r, int(posn[1])-r:int(posn[1])+r]\n",
    "                \n",
    "\n",
    "                #ellipse fit :\n",
    "                '''\n",
    "                edges = cv2.Canny(gray_pmt_img,0,100)\n",
    "                print([i, j])\n",
    "                result = hough_ellipse(edges, accuracy=50, threshold=100, min_size=r, max_size=2*r)\n",
    "                result.sort(order='accumulator')\n",
    "                if len(result) == 0 :\n",
    "                    continue\n",
    "                best = list(result[-1])\n",
    "                yc, xc, a, b = [int(round(x)) for x in best[1:5]]\n",
    "                circles.append([int(j)+200+yc, int(i)+100+xc, int(a), int(b)])\n",
    "                '''\n",
    "                #blob\n",
    "                '''\n",
    "               # detector = cv2.SimpleBlobDetector()\n",
    "                # Setup SimpleBlobDetector parameters.\n",
    "                params = cv2.SimpleBlobDetector_Params()\n",
    "                #params.minThreshold = 250\n",
    "                #params.filterByColor = True\n",
    "                #params.blobColor = 0\n",
    "                # Filter by Area.\n",
    "                params.filterByArea = True\n",
    "                params.minArea = 500\n",
    "                # Filter by Circularity\n",
    "                params.filterByCircularity = True\n",
    "                params.minCircularity = 0.1\n",
    "\n",
    "                # Create a detector with the parameters\n",
    "                detector = cv2.SimpleBlobDetector_create(params)\n",
    "                # Detect blobs.\n",
    "                keypoints = detector.detect(inv_pmt_img)                \n",
    "                im_with_keypoints = cv2.drawKeypoints(inv_pmt_img, keypoints, np.array([]), (0,0,255), cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)\n",
    "                print(keypoints)\n",
    "                # Show blobs\n",
    "                cv2.imshow(\"Keypoints\", im_with_keypoints)\n",
    "                cv2.waitKey(0)\n",
    "                '''\n",
    "                contours, hierarchy = cv2.findContours(one_pmt_img,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)\n",
    "                area = 0\n",
    "                coords = []\n",
    "                \n",
    "                #debug:draw contours\n",
    "                #cont_pmt_img = cv2.drawContours(one_pmt_img, contours, -1, (255,255,0), 3)\n",
    "                #cv2.imshow('img', cont_pmt_img)\n",
    "                #cv2.waitKey(0)\n",
    "                \n",
    "                \n",
    "                for c in contours:\n",
    "                #moments\n",
    "                    new_M = cv2.moments(c)\n",
    "                    new_area = cv2.contourArea(c)\n",
    "\n",
    "                    if (new_area <=300) or (new_area>= 100000) :  # skip ellipses smaller then 10x10\n",
    "                        continue\n",
    "\n",
    "                       # calculate x,y coordinate of center, breaks if m00 = 0 output location in large image (offset)\n",
    "\n",
    "                    elif new_M[\"m00\"] != 0:\n",
    "                        cX = int(new_M[\"m10\"] / new_M[\"m00\"])\n",
    "                        cY = int(new_M[\"m01\"] / new_M[\"m00\"])\n",
    "                                 #   print(area)\n",
    "\n",
    "                    else:\n",
    "                        continue\n",
    "\n",
    "                    if new_area >= area :    #TODO:BETTER CHECK?\n",
    "                        area = new_area\n",
    "                        coords = [cX, cY, r]\n",
    "                       # print(area)\n",
    "                    \n",
    "                #Fit an ellipse  \n",
    "                #cut image to center of centroid of pmt\n",
    "                #print(coords)\n",
    "                #TODO: posn coords are wrtong way round\n",
    "                \n",
    "                #TODO: Check by drawing contours\n",
    "                if len(coords) > 0 :\n",
    "                    x_cor_circ = int(posn[1]) - r\n",
    "                    y_cor_circ = int(posn[0]) - r\n",
    "\n",
    "                    x_cor_ell = int(coords[1]) - r\n",
    "                    y_cor_ell = int(coords[0]) - r\n",
    "\n",
    "                    #print([x_cor_circ, x_cor_ell])\n",
    "                    #corrected to center image on centre of blob\n",
    "                    one_pmt_bolt_img = bolt_img[int(coords[1])+y_cor_circ-r:int(coords[1])+y_cor_circ+r, int(coords[0])+x_cor_circ-r:int(coords[0])+x_cor_circ+r]\n",
    "\n",
    "                   # cv2.imshow('bolt', one_pmt_bolt_img)\n",
    "                   # cv2.imshow('pmt', one_pmt_img)\n",
    "                   #' cv2.imshow('bold', one_pmt_bolt_img_orig)\n",
    "                    #cv2.waitKey(0)\n",
    "                \n",
    "\n",
    "                    #print([i, j, coords, r, int(i)+100-r, int(j)+200-r])\n",
    "                    circles.append([int(coords[0])  + x_cor_circ, int(coords[1]) + y_cor_circ, int(r), label])\n",
    "                    ellipse_params =  fit_to_bolts(one_pmt_bolt_img, coords) #[1,1,1,1,1] #   \n",
    "                    ellipses.append([ellipse_params[0]+ x_cor_circ, ellipse_params[1]+y_cor_circ, ellipse_params[2], ellipse_params[3], ellipse_params[4], label, ellipse_params[5]])\n",
    "\n",
    "            #    circles.append([int(j)+200, int(i)+100, int(r)])\n",
    "    return circles, ellipses\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Main Bolt Finder\n",
    "#TODO : MARGIN AT BOTTOM OF IMAGE\n",
    "def find_centroids(idir, odir, img_name, guess) :\n",
    "    #pmt_guess = open(f'/home/dm3315/Documents/SK/PMT_learning/images/ring_reco/four_gaps/guess_meth/pmt_loc.txt', \"r\")\n",
    "\n",
    "    feature_info = open(f'{odir}coords/{img_name}.txt', \"w\")\n",
    "    \n",
    "    img = cv2.imread(f'{idir}{img_name}.png')\n",
    "\n",
    "    #Set up separate masks for bolts and pmts and an output image\n",
    "    train_labels = img\n",
    "    lower_pmts = np.array([254, 0, 0], dtype = \"uint16\")\n",
    "    upper_pmts = np.array([255,0,0], dtype = \"uint16\")\n",
    "    lower_bolts = np.array([0, 0, 254], dtype = \"uint16\")\n",
    "    upper_bolts = np.array([0,0,255], dtype = \"uint16\")\n",
    "\n",
    "    pmts_gray = cv2.inRange(train_labels, lower_pmts, upper_pmts)\n",
    "    bolts_gray = cv2.inRange(train_labels, lower_bolts, upper_bolts)\n",
    "    img = cv2.cvtColor(bolts_gray,cv2.COLOR_GRAY2BGR) #used for making output image\n",
    "\n",
    "    #print('images made')\n",
    "\n",
    "    #Calls HoughCircles() to find multiple PMTs in an image\n",
    "\n",
    "    #Outpts list of circle parameters\n",
    "\n",
    "\n",
    "\n",
    "    \n",
    "    #find circles in image\n",
    "    \n",
    "    #TODO: Do circles then filter out bads, then do ellipse to save time\n",
    "    circles, ellipses = find_pmts(pmts_gray, bolts_gray, guess) ######### TODO:DO ELLIPSE FIT IN HERE RETURNED CIRCLES ARE INSRTEAD ELLIPSES\n",
    "    bolt_centres = []\n",
    "\n",
    "    print(f'number of circles:{len(circles)}, ellipses:{len(ellipses)}')\n",
    "    #loop over circles and find bolts, end up with list of circles, each a  list of bolts (no circle parameters are saved )\n",
    "    bolt_no = 0\n",
    "    pmt_reco = 0\n",
    "    \n",
    "    #info contains list of circles, each a list of bolt centres, and number of bolts\n",
    "    info = []\n",
    "    info.append([])\n",
    "    info.append([])\n",
    "    dia_err = []\n",
    "    \n",
    "    \n",
    "    \n",
    "    #output ellipse fits\n",
    "    for i in range(len(circles)) :\n",
    "        \n",
    "        circ = circles[i]\n",
    "        ell = ellipses[i]\n",
    "        \n",
    "       # print([[circ[0], circ[1]], [ell[0], ell[1]]])\n",
    "        cv2.circle(img,(circ[0],circ[1]),circ[2],(0,100,0),2)\n",
    "\n",
    "        bolt_ring = ret_centres(bolts_gray, circ) # list of bolt centres\n",
    "        #temp_bolts_gray = bolts_gray.copy()\n",
    "\n",
    "       # cv2.imshow('img2', img)\n",
    "     #   cv2.waitKey(0)\n",
    "        if len(bolt_ring) < 3 :\n",
    "           # cv2.circle(img,(circ[0],circ[1]),circ[2],(0,255,0),2)\n",
    "           # feature_info.write(f'00    {circ[1]}    {circ[0]}\\n')\n",
    "            continue\n",
    "            \n",
    "        bolts_ord = number_bolts(circ, bolt_ring) # ordered list of bolt centres\n",
    "        dias = bolt_avg_dia(circ, bolts_ord)      #lis of bolt centres, diameters and number of bolts in circle\n",
    "\n",
    "\n",
    "        avg_dia = 0\n",
    "        avg_centre = 0\n",
    "        num = len(dias[2])\n",
    "        \n",
    "      #  if num<3:\n",
    "      #          cv2.circle(img,(circ[0],circ[1]),circ[2],(0,100,0),2)\n",
    "      #          continue\n",
    "        pmt_reco += 1\n",
    "        \n",
    "        #Draw diameters here\n",
    "        #for j in dias[0] : \n",
    "        #    print(j, j[0])\n",
    "            #avg_dia += np.linalg.norm(np.array(np.array(j[1]) - np.array(j[0]))) / num\n",
    "            #avg_centre += \n",
    "           # cv2.line(img, (j[0][0], j[0][1]), (j[1][0], j[1][1]), (0, 255, 0), thickness=2, lineType=8)\n",
    "          #  cv2.imshow('img', img)\n",
    "           # cv2.waitKey(0)\n",
    "        #Draw hough circles on img\n",
    "        #if num !=0 :\n",
    "            #avg_dia = sum_dia/num\n",
    "\n",
    "            \n",
    "        cv2.ellipse(img, (int(ell[0]), int(ell[1])), (int(ell[2]), int(ell[3])), int(ell[4]), 0, 360, (0, 255, 255), 2)\n",
    "        cv2.putText(img, f'{circ[3]}', (int(ell[0]), int(ell[1])), cv2.FONT_HERSHEY_SIMPLEX, 1, (255,255,255), 2)\n",
    "\n",
    "        feature_info.write(f'00:{circ[0]}:{circ[1]}:{circ[2]}:{ell[0]}:{ell[1]}:{ell[2]}:{ell[3]}:{ell[4]}:{circ[3]}:{ell[6]}\\n')\n",
    "\n",
    "        avg_dia = np.average(dias[2])\n",
    "        dia_err = np.asarray(dias[2])/avg_dia\n",
    "        info[0].extend(dias[1])\n",
    "        info[1].extend((dia_err - 1)*100)\n",
    "        \n",
    "\n",
    "        bolt_no = 0\n",
    "    #draw bolt locations\n",
    "        for k in bolts_ord :\n",
    "            bolt_no += 1\n",
    "#        cv2.circle(img, (j[1], j[0]), 5, (0,0,255), -1)\n",
    "            cv2.putText(img, f'{bolt_no}', (k[0], k[1]), cv2.FONT_HERSHEY_SIMPLEX, 1, (255,255,255), 1)\n",
    "            feature_info.write(f'{bolt_no}:{k[0]}:{k[1]}:{circ[3]}\\n')\n",
    "\n",
    "    #print(info)\n",
    "    #make a plot of bolt numbers and diameters in whole image\n",
    "    print(f'number of pmts = {len(circles)}')\n",
    "    print(f'number of pmts reconstructed = {pmt_reco}')\n",
    "    print(f'percentage PMTs = {pmt_reco}')\n",
    "\n",
    "    print(f'number of bolt pairs = {np.sum(info[0])/2}')\n",
    "    print(f'number of diameters = {len(info[1])}')\n",
    "    eval_reco(odir, info, img)\n",
    "    cv2.imwrite(f'{odir}plots/{img_name}_labelled.png', img )\n",
    "    feature_info.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Ellipse fitter - run on small bolt image\n",
    "\n",
    "def fit_to_bolts(img, init_coords) :\n",
    "\n",
    "\n",
    "    #print(circle)\n",
    "    circ_x = init_coords[0]\n",
    "    circ_y = init_coords[1]\n",
    "    circ_r = init_coords[2]\n",
    "    L = int(circ_r)          \n",
    "    \n",
    "    x_cor_ell = circ_x - circ_r\n",
    "    y_cor_ell = circ_y-circ_r\n",
    "    \n",
    "    best_param = [x_cor_ell, y_cor_ell, 10, 10, 0]\n",
    "    draw_temp_img = img.copy()\n",
    "    centre = np.asarray([circ_x, circ_y])\n",
    "\n",
    "    # gray_image = bolts_gray[circ_y-L:circ_y+L, circ_x-L:circ_x+L]\n",
    "\n",
    "   # cv2.imshow('img', img)\n",
    "   # cv2.waitKey(0)\n",
    "    \n",
    "    #First loop fits a circle and its centre \n",
    "    score = 0\n",
    "    for i in range(int(circ_r*0.4)) : #r1\n",
    "        for j in range(40) : # xcen\n",
    "            for k in range(40) :# ycen\n",
    "                new_r1 = circ_r*0.8 + i\n",
    "                new_r2 = circ_r*0.8 + i\n",
    "\n",
    "                mult = 2\n",
    "\n",
    "                new_cen = (L - 20*mult + j*mult, L - 20*mult + k*mult)\n",
    "                new_r = (int(new_r1), int(new_r2))\n",
    "                new_ang = 0\n",
    "                circ_mask = np.zeros(img.shape,dtype=np.uint8)\n",
    "\n",
    "               # print(new_cen, new_r, new_ang, 0, 360, 1, 1)\n",
    "                cv2.ellipse(circ_mask, new_cen, new_r, new_ang, 0, 360, color=1 , thickness=1)\n",
    "                fit = img*circ_mask\n",
    "                new_score = np.sum(fit)\n",
    "                \n",
    "                \n",
    "                if new_score >= score :\n",
    "                    score = new_score\n",
    "                    best_param = [new_cen[0]+x_cor_ell, new_cen[1]+y_cor_ell, new_r1, new_r2, new_ang, score ]\n",
    "\n",
    "#############################################################################\n",
    "#second loop does the ellipse \n",
    "    for j in range(int(circ_r*0.8)) : #r2\n",
    "        for m in [0, 45] : # angle\n",
    "                #print(i)\n",
    "                new_r1 = best_param[2]\n",
    "                new_r2 = circ_r*0.6+ j\n",
    "\n",
    "                new_r = (int(new_r1), int(new_r2))\n",
    "                new_ang = m\n",
    "                new_cen = (best_param[0]-x_cor_ell, best_param[1]-y_cor_ell)\n",
    "                circ_mask = np.zeros(img.shape,dtype=np.uint8)\n",
    "                cv2.ellipse(circ_mask, new_cen, new_r, new_ang, 0, 360, 1, 1)\n",
    "                fit = img*circ_mask\n",
    "\n",
    "                #Debug\n",
    "                #cv2.imshow('img', fit)\n",
    "                #cv2.waitKey(0)   \n",
    "                \n",
    "                init_rad = circ_r\n",
    "                new_score = np.sum(fit)\n",
    "                if new_score >= score :\n",
    "                    score = new_score\n",
    "                    best_param = [new_cen[0]+x_cor_ell, new_cen[1]+y_cor_ell, new_r1, new_r2, new_ang, score ] #TODO:DIFFERENT\n",
    "    #debug\n",
    "    #print(best_param)\n",
    "    #print(int(best_param[0]-x_cor_ell))\n",
    "    #cv2.ellipse(draw_temp_img, (int(best_param[0]-x_cor_ell), int(best_param[1]-y_cor_ell)), (int(best_param[2]), int(best_param[3])), int(best_param[4]), 0, 360, (255, 255, 255), 2)\n",
    "  #  cv2.imwrite(f'{odir}plots/temp/{circ_x}.png', draw_temp_img )\n",
    "\n",
    "    return best_param\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run below here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['794_pred_substacks_far_unet.19']\n",
      "processing image: 794\n",
      "number of circles:54, ellipses:54\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-2-a54089b8c529>:105: RuntimeWarning: invalid value encountered in arccos\n",
      "  angle = np.arccos(np.dot(rad_1, rad_2) / (np.linalg.norm(rad_1) * np.linalg.norm(rad_2)))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of pmts = 54\n",
      "number of pmts reconstructed = 52\n",
      "percentage PMTs = 52\n",
      "number of bolt pairs = 375.0\n",
      "number of diameters = 202\n",
      "time taken = 1310.2389328479767 seconds\n"
     ]
    }
   ],
   "source": [
    "import time \n",
    "start_time = time.time()\n",
    "\n",
    "dataset = 'substacks_far'\n",
    "seg_files = 'ring_avg_all'\n",
    "method = 'manual_intervention'\n",
    "supervision = 'supervised'\n",
    "info_filename = 'fixed_pmt_loc_label'\n",
    "\n",
    "#img_name = 'SING0046'\n",
    "#img_path = f'/home/dm3315/Documents/SK/PMT_learning/images/datasets/{dataset}/output/output/{img_name}.png'\n",
    "idir = f'/home/dm3315/Documents/SK/PMT_learning/images/datasets/{dataset}/output/{seg_files}/'\n",
    "odir = f'{idir}feature_info/{supervision}/'\n",
    "#img_path = '/home/dm3315/Documents/SK/PMT_learning/images/datasets/bolts_pmt/output/output/large_test_out/4_pred_noverlay.png'\n",
    "\n",
    "f = []\n",
    "for (dirpath, dirnames, filenames) in os.walk(f'{idir}'):\n",
    "    f.extend(filenames)\n",
    "    break\n",
    "\n",
    "#f = [x[:-4] for x in f]\n",
    "#f = f[:]\n",
    "f = ['794_pred_substacks_far_unet.19']\n",
    "print(f)\n",
    "#f = ['B7730726_pred_substacks_far_unet.19']\n",
    "\n",
    "guess_info = {}\n",
    "with open (f'/home/dm3315/Documents/SK/PMT_learning/images/datasets/{dataset}/output/{seg_files}/smodule/{method}/{info_filename}.txt', 'r') as info:\n",
    "    for row in csv.reader(info,delimiter=':') :\n",
    "        if str(row[0]) not in guess_info:\n",
    "            guess_info[str(row[0])] = {'pmt_pos':[], 'pmt_label':[]}\n",
    "        guess_info[str(row[0])]['pmt_pos'].append([float(row[1]), float(row[2])])\n",
    "        guess_info[str(row[0])]['pmt_label'].append([float(row[3]), float(row[4]), float(row[5]), float(row[6])])\n",
    "        \n",
    "#print(guess_info)\n",
    "\n",
    "global r_mult #for cutting windows around PMTs\n",
    "r_mult = float(1.0)\n",
    "\n",
    "for i in f :\n",
    "    if i[:3] in guess_info :\n",
    "        print(f'processing image: {i[:3]}')\n",
    "        #print(f'{guess_info[str(i[:4])]}')\n",
    "        find_centroids(idir, odir, i, guess_info[str(i[:3])])\n",
    "        plt.close()\n",
    "    else:\n",
    "        print(f'PMT information not provided for {i[:3]}')\n",
    "        \n",
    "print(f'time taken = {time.time() - start_time} seconds')\n",
    "\n",
    "\n",
    "cv2.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}