{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Analyse images to determine asymmetry of comets ##\n", "\n", "Look at images, find good comets using three thresholds and determine their asymmetry." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Some useful routines\n", "\n", "Convert RGBA images to RGB and RGB to greyscale. Ceiling division. Analysis routines...\n", "\n", "Note that must set negative = 1 in call to rgb2grey in processImage for MC data, negative = 0 for real images." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Date and time 2021-06-03 14:06:43.249654\n", " \n", "Date and time 2021-06-03 14:06:43.641679\n", "Time since last check is 0:00:00.392025\n" ] } ], "source": [ "import datetime\n", "now = datetime.datetime.now()\n", "print(\"Date and time \",str(now))\n", "#\n", "import numpy as np\n", "#\n", "def rgba2rgb(rgba, background = (255, 255, 255)):\n", " '''\n", " Function to convert RGBA images into RGB format. Input RGBA image (and background); output RGB image.\n", " '''\n", " rows, cols, chans = rgba.shape\n", " #\n", " debug = False\n", " #\n", " if debug:\n", " print(\"------------------------------------------------------------------------------------\")\n", " print(\"Running rgba2rgb\")\n", " if chans == 4:\n", " print(\"RGBA image\")\n", " elif chans == 3:\n", " print(\"RGB image\")\n", " return rgba\n", " else:\n", " print(\"Channel number is\",chans)\n", " sys.exit()\n", " else:\n", " assert chans == 4, 'RGBA image must have 4 channels.'\n", " #\n", " rgb = np.zeros((rows, cols, 3), dtype = 'float32')\n", " r, g, b, a = rgba[:,:,0], rgba[:,:,1], rgba[:,:,2], rgba[:,:,3]\n", " #\n", " a = np.asarray(a, dtype='float32')/255.0\n", " #\n", " R, G, B = background\n", " #\n", " rgb[:, :, 0] = r*a + (1.0 - a)*R\n", " rgb[:, :, 1] = g*a + (1.0 - a)*G\n", " rgb[:, :, 2] = b*a + (1.0 - a)*B\n", " #\n", " return np.asarray(rgb, dtype = np.uint8)\n", "#\n", "#\n", "#\n", "def rgb2grey(rgb, negative = 0, withHists = False):\n", " '''\n", " Convert RGB image to greyscale. Input RGB (and flag indicating negative required), output greyscale image.\n", " '''\n", " rows, cols, chans = rgb.shape\n", " #\n", " debugHere = False\n", " #\n", " if debugHere:\n", " print(\"------------------------------------------------------------------------------------\")\n", " print(\"Running rgb2grey\")\n", " if chans == 4:\n", " print(\"RGBA image\")\n", " elif chans == 3:\n", " print(\"RGB image\")\n", " elif chans == 1:\n", " print(\"Greyscale image\")\n", " return rgb\n", " else:\n", " print(\"Channel number is\",chans)\n", " sys.exit()\n", " else:\n", " assert chans == 3, 'RGB image must have 3 channels.'\n", " #\n", " grey = np.zeros((rows, cols), dtype = 'float32')\n", " #\n", " r, g, b = rgb[:,:,0], rgb[:,:,1], rgb[:,:,2]\n", " #\n", " grey[:, :] = (0.2125*(r*negative + (negative - 1.0)*r) + \n", " 0.7154*(g*negative + (negative - 1.0)*g) + \n", " 0.0721*(b*negative + (negative - 1.0)*b))\n", " if withHists:\n", " print(\" \")\n", " print(\"Max intensities red\",np.amax(r),\"blue\",np.amax(b),\"green\",np.amax(g))\n", " print(\" \")\n", " nBins = int(256/8)\n", " nBins = 256\n", " plt.figure(figsize = (10, 9))\n", " plt.subplot(4, 1, 1)\n", " plt.hist(np.ravel(r), bins = nBins, color = 'r')\n", " plt.yscale(\"log\")\n", " plt.xlim(-1, 256)\n", " plt.subplot(4, 1, 2)\n", " plt.hist(np.ravel(g), bins = nBins, color = 'g')\n", " plt.yscale(\"log\")\n", " plt.xlim(-1, 256)\n", " plt.subplot(4, 1, 3)\n", " plt.hist(np.ravel(b), bins = nBins, color = 'b')\n", " plt.yscale(\"log\")\n", " plt.xlim(-1, 256)\n", " plt.subplot(4, 1, 4)\n", " plt.hist(np.ravel(grey), bins = nBins, color = 'k')\n", " plt.yscale(\"log\")\n", " plt.xlim(-1, 256)\n", " plt.tight_layout()\n", " plt.show()\n", " #\n", " return np.asarray(grey, dtype = np.uint8)\n", "#\n", "#\n", "#\n", "def ceilDiv(a, b):\n", " '''\n", " Return a//b rounded up.\n", " '''\n", " ceiling = -(-a//b)\n", " return ceiling\n", "#\n", "#\n", "#\n", "def processImage(imgRaw, hotCut, withHists):\n", " #\n", " import sys\n", " import numpy as np\n", " import scipy.ndimage as ndimage\n", " import matplotlib.pyplot as plt\n", " %matplotlib inline\n", " # \n", " # Number of rows and columns. Depth is 3 for RGB, 4 for RGBA image. A is opacity (alpha)\n", " nRows = imgRaw.shape[0] \n", " nCols = imgRaw.shape[1]\n", " nDepth = imgRaw.shape[2]\n", " #\n", " nThresh = len(thresh)\n", " img = np.zeros((nRows, nCols))\n", " imgThr = np.zeros((nRows, nCols, nThresh))\n", " #\n", " # Determine image format and process accordingly\n", " if nDepth == 4:\n", " imgRGB = rgba2rgb(imgRaw)\n", " imgGrey = rgb2grey(imgRGB, 0, withHists)\n", " elif nDepth == 3:\n", " imgRGB = imgRaw\n", " imgGrey = rgb2grey(imgRGB, 1, withHists)\n", " elif nDepth == 1:\n", " imgRGB = imgRaw\n", " imgGrey = imgRaw\n", " else:\n", " print(\" \")\n", " print(\"Unexpected image depth\",nDepth)\n", " sys.stop()\n", " #\n", " # Remove \"hot\" pixels\n", " hotPixels = imgGrey > hotCut\n", " imgGrey[hotPixels] = 0\n", " #\n", " # Rescale \n", " imgGrey = np.asarray(255/hotCut*imgGrey, dtype = np.uint8)\n", " #\n", " if withHists:\n", " print(\" \")\n", " print(\"Maximum pixel value before removing hot pixels\",np.amax(imgGrey))\n", " print(\"Minimum pixel value \",np.amin(imgGrey))\n", " print(\"Type of raw image file is\",imgRaw.dtype)\n", " print(\"Type of greyscale image file is\",imgGrey.dtype)\n", " print(\"Number of rows\",nRows,\"of columns\",nCols,\"of pixels\",nRows*nCols,\"and depth\",nDepth)\n", " print(\"Hot pixel cut\",hotCut)\n", " print(\"Number of hot pixels\",np.sum(hotPixels))\n", " print(\"Maximum pixel value after removing hot pixels\",np.amax(imgGrey))\n", " print(\"Minimum pixel value \",np.amin(imgGrey))\n", " #\n", " return imgGrey, nRows, nCols\n", "#\n", "#\n", "#\n", "def findWheels(imgGrey, threshold):\n", " #\n", " if not hasattr(findWheels, \"firstCall\"):\n", " findWheels.firstCall = True\n", " print(\"------------------------------------------------------------------------------------\")\n", " print(\"Running findWheels\")\n", " if debug:\n", " findWheels.firstCall = True\n", " print(\"------------------------------------------------------------------------------------\")\n", " print(\"Running findWheels\")\n", " #\n", " # Create image with cluster threshold and run watershed algorithm\n", " clusImg = imgGrey > threshold\n", " if debug:\n", " plt.figure(figsize = (pltX, pltY))\n", " plt.imshow(clusImg);\n", " colDotsClus = ndimage.watershed_ift(clusImg.astype(np.uint8), markers)\n", " #\n", " # Remove \"isolated\" markers (i.e. markers in regions where no cluster found)\n", " colDotsClus[rMark, cMark] = colDotsClus[rMark + 1, cMark]\n", " #\n", " # Find the value of the marker of the last identified regions\n", " mMaxClus = np.amax(colDotsClus)\n", " #\n", " # Positions of marker values\n", " boolClus = colDotsClus >= mStart\n", " nFoundClus = len(np.unique(colDotsClus[boolClus]))\n", " rMarkerClus = np.zeros(nFoundClus)\n", " rMarkerClus = np.unique(colDotsClus[boolClus])\n", " #\n", " # Select the clusters, first pass (determine accepted number). Must exclude tiny clusters (noise) and any \n", " # really large clusters (background regions in image).\n", " nClus = 0\n", " cMarkerClus = np.zeros(nFoundClus)\n", " maxInDotClus = 0\n", " for nR in range(0, nFoundClus):\n", " boolClus = colDotsClus == rMarkerClus[nR]\n", " nHereCl = np.sum(boolClus)\n", " if nHereCl < minClusPixels or nHereCl > maxClusPixels:\n", " continue\n", " maxInDotClus = max(maxInDotClus, nHereCl)\n", " cMarkerClus[nClus] = rMarkerClus[nR]\n", " nClus += 1\n", " #\n", " # Information on clusters\n", " nInCluster = np.zeros(nClus).astype(int)\n", " iClusSum = np.zeros(nClus)\n", " #\n", " # Information on pixels in clusters\n", " lClus = np.sum(clusImg)\n", " #\n", " # Safe size of arrays would be number of pixels in image. Try to reduce size by using number of pixels in clusters\n", " # The factor lFact can be used to expand array sizes\n", " lFact = 4\n", " indexCl = np.zeros(lFact*lClus)\n", " cPixelsCl = np.zeros(lFact*lClus).astype(int)\n", " rPixelsCl = np.zeros(lFact*lClus).astype(int)\n", " iPixelsCl = np.zeros(lFact*lClus)\n", " #\n", " # Temporary information \n", " cPixelsHere = np.zeros(lClus).astype(int)\n", " rPixelsHere = np.zeros(lClus).astype(int)\n", " iPixelsHere = np.zeros(lClus)\n", " #\n", " # Find pixels in clusters, determine positions and plot \n", " if findWheels.firstCall:\n", " fig = plt.figure(figsize = (pltX, pltY))\n", " ax = fig.add_subplot(1, 1, 1)\n", " plt.title(\"Clusters with threshold \" + str(threshold), fontsize = 12)\n", " plt.xlabel('x pixel', fontsize = 12)\n", " plt.ylabel('y pixel', fontsize = 12)\n", " #\n", " # Figure control\n", " xOffMax = 2\n", " yOffMax = 2\n", " mSize = 0.001\n", " nCol = 0\n", " #\n", " nLastCl = 0\n", " #\n", " for nC in range(0, nClus):\n", " #\n", " # Clusters\n", " boolClus = colDotsClus == cMarkerClus[nC]\n", " nHereCl = np.sum(boolClus).astype(int)\n", " #\n", " nInCluster[nC] = nHereCl\n", " #\n", " # Indices of x and y pixels\n", " rPixelsHere, cPixelsHere = np.where(boolClus)\n", " #\n", " # Intensities in pixels\n", " iPixelsHere = imgGrey[rPixelsHere, cPixelsHere]\n", " iClusSum[nC] = np.sum(iPixelsHere)\n", " #\n", " indexCl[nLastCl:nLastCl + nHereCl] = nC*np.ones(nHereCl)\n", " cPixelsCl[nLastCl:nLastCl + nHereCl] = cPixelsHere[:]\n", " rPixelsCl[nLastCl:nLastCl + nHereCl] = rPixelsHere[:]\n", " iPixelsCl[nLastCl:nLastCl + nHereCl] = iPixelsHere[:]\n", " #\n", " if findWheels.firstCall:\n", " plt.scatter(cPixelsHere, rPixelsHere, s = mSize, c = colorTab[nCol], marker = 'o')\n", " rLab = np.amax(rPixelsHere[0:nHereCl]) + yOffMax\n", " cLab = np.amax(cPixelsHere[0:nHereCl]) + xOffMax\n", " plt.text(cLab, rLab, str(nC), color = colorTab[nCol])\n", " nCol = nCol + 1\n", " if nCol > nColTab - 1:\n", " nCol = 0\n", " #\n", " nLastCl = nLastCl + nHereCl\n", " #\n", " #\n", " if findWheels.firstCall:\n", " plt.xlim(-0.05*nCols, 1.05*nCols)\n", " plt.ylim(-0.05*nRows, 1.05*nRows)\n", " plt.grid(color = 'g')\n", " print(\" \")\n", " plt.show()\n", " #\n", " findWheels.firstCall = False\n", " #\n", " return nInCluster, rPixelsCl, cPixelsCl, iPixelsCl\n", "#\n", "#\n", "#\n", "def edgeFinderIn(picture, edgeWidth, useDiag = True):\n", " '''\n", " Return array containing pixels in edges (of width edgeWidth) of input (thresholded) image.\n", " The edges are inside the original image. The flag useDiag ensures \"corner\" pixels are selected.\n", " '''\n", " #\n", " shiftR = edgeWidth\n", " shiftC = edgeWidth\n", " nRows, nCols = picture.shape\n", " edges = np.full((nRows, nCols), False)\n", " imgShift = np.full((nRows, nCols), False)\n", " #\n", " # Right edge\n", " imgShift[0:nRows, 0:nCols - shiftC] = picture[0:nRows, shiftC:nCols]\n", " imgShift[0:nRows, nCols - shiftC:nCols] = False\n", " edges = np.logical_and(picture, np.logical_not(imgShift))\n", " #\n", " # Left edge \n", " imgShift[0:nRows, shiftC:nCols] = picture[0:nRows, 0:nCols - shiftC]\n", " imgShift[0:nRows, 0:shiftC] = False\n", " edges = np.logical_or(edges, np.logical_and(picture, np.logical_not(imgShift)))\n", " #\n", " # Lower edge \n", " imgShift[shiftR:nRows, 0:nCols] = picture[0:nRows - shiftR, 0:nCols]\n", " imgShift[0:shiftR, 0:nCols] = False \n", " edges = np.logical_or(edges, np.logical_and(picture, np.logical_not(imgShift)))\n", " #\n", " # Upper edge \n", " imgShift[0:nRows - shiftR, 0:nCols] = picture[shiftR:nRows, 0:nCols] \n", " imgShift[nRows - shiftR:nRows, 0:nCols] = False\n", " edges = np.logical_or(edges, np.logical_and(picture, np.logical_not(imgShift)))\n", " #\n", " if useDiag:\n", " #\n", " # Left upper edge\n", " imgShift[0:nRows, 0:nCols - shiftC] = picture[0:nRows, shiftC:nCols]\n", " imgShift[0:nRows, nCols - shiftC:nCols] = False\n", " imgShift[shiftR:nRows, 0:nCols] = imgShift[0:nRows - shiftR, 0:nCols]\n", " imgShift[0:shiftR, 0:nCols] = False \n", " edges = np.logical_or(edges, np.logical_and(picture, np.logical_not(imgShift)))\n", " #\n", " # Left lower edge\n", " imgShift[0:nRows, 0:nCols - shiftC] = picture[0:nRows, shiftC:nCols]\n", " imgShift[0:nRows, nCols - shiftC:nCols] = False\n", " imgShift[0:nRows - shiftR, 0:nCols] = imgShift[shiftR:nRows, 0:nCols] \n", " imgShift[nRows - shiftR:nRows, 0:nCols] = False\n", " edges = np.logical_or(edges, np.logical_and(picture, np.logical_not(imgShift)))\n", " #\n", " # Right upper edge \n", " imgShift[0:nRows, shiftC:nCols] = picture[0:nRows, 0:nCols - shiftC]\n", " imgShift[0:nRows, 0:shiftC] = False\n", " imgShift[shiftR:nRows, 0:nCols] = imgShift[0:nRows - shiftR, 0:nCols]\n", " imgShift[0:shiftR, 0:nCols] = False \n", " edges = np.logical_or(edges, np.logical_and(picture, np.logical_not(imgShift)))\n", " #\n", " # Right lower edge \n", " imgShift[0:nRows, shiftC:nCols] = picture[0:nRows, 0:nCols - shiftC]\n", " imgShift[0:nRows, 0:shiftC] = False\n", " imgShift[0:nRows - shiftR, 0:nCols] = imgShift[shiftR:nRows, 0:nCols] \n", " imgShift[nRows - shiftR:nRows, 0:nCols] = False\n", " edges = np.logical_or(edges, np.logical_and(picture, np.logical_not(imgShift)))\n", " #\n", " return edges\n", "#\n", "#\n", "#\n", "def expander(picture, edgeWidth, useDiag = True):\n", " '''\n", " Return array containing thresholded regions expanded by band of width edgeWidth. The flag useDiag ensures\n", " \"corner\" pixels are included correctly.\n", " '''\n", " #\n", " shiftR = 1\n", " shiftC = 1\n", " nRows, nCols = picture.shape\n", " edges = np.full((nRows, nCols), False)\n", " edgeSum = np.full((nRows, nCols), False)\n", " imgShift = np.full((nRows, nCols), False)\n", " #\n", " for nE in range(0, edgeWidth):\n", " #\n", " # Left edge\n", " imgShift[0:nRows, 0:nCols - shiftC] = picture[0:nRows, shiftC:nCols]\n", " imgShift[0:nRows, nCols - shiftC:nCols] = False\n", " edges = np.logical_and(np.logical_not(picture), imgShift)\n", " #\n", " # Right edge \n", " imgShift[0:nRows, shiftC:nCols] = picture[0:nRows, 0:nCols - shiftC]\n", " imgShift[0:nRows, 0:shiftC] = False\n", " edges = np.logical_or(edges, np.logical_and(np.logical_not(picture), imgShift))\n", " #\n", " # Upper edge \n", " imgShift[shiftR:nRows, 0:nCols] = picture[0:nRows - shiftR, 0:nCols]\n", " imgShift[0:shiftR, 0:nCols] = False \n", " edges = np.logical_or(edges, np.logical_and(np.logical_not(picture), imgShift))\n", " #\n", " # Lower edge \n", " imgShift[0:nRows - shiftR, 0:nCols] = picture[shiftR:nRows, 0:nCols] \n", " imgShift[nRows - shiftR:nRows, 0:nCols] = False\n", " edges = np.logical_or(edges, np.logical_and(np.logical_not(picture), imgShift))\n", " #\n", " if useDiag:\n", " #\n", " # Left upper edge\n", " imgShift[0:nRows, 0:nCols - shiftC] = picture[0:nRows, shiftC:nCols]\n", " imgShift[0:nRows, nCols - shiftC:nCols] = False\n", " imgShift[shiftR:nRows, 0:nCols] = imgShift[0:nRows - shiftR, 0:nCols]\n", " imgShift[0:shiftR, 0:nCols] = False \n", " edges = np.logical_or(edges, np.logical_and(np.logical_not(picture), imgShift))\n", " #\n", " # Left lower edge\n", " imgShift[0:nRows, 0:nCols - shiftC] = picture[0:nRows, shiftC:nCols]\n", " imgShift[0:nRows, nCols - shiftC:nCols] = False\n", " imgShift[0:nRows - shiftR, 0:nCols] = imgShift[shiftR:nRows, 0:nCols] \n", " imgShift[nRows - shiftR:nRows, 0:nCols] = False\n", " edges = np.logical_or(edges, np.logical_and(np.logical_not(picture), imgShift))\n", " #\n", " # Right upper edge \n", " imgShift[0:nRows, shiftC:nCols] = picture[0:nRows, 0:nCols - shiftC]\n", " imgShift[0:nRows, 0:shiftC] = False\n", " imgShift[shiftR:nRows, 0:nCols] = imgShift[0:nRows - shiftR, 0:nCols]\n", " imgShift[0:shiftR, 0:nCols] = False \n", " edges = np.logical_or(edges, np.logical_and(np.logical_not(picture), imgShift))\n", " #\n", " # Right lower edge \n", " imgShift[0:nRows, shiftC:nCols] = picture[0:nRows, 0:nCols - shiftC]\n", " imgShift[0:nRows, 0:shiftC] = False\n", " imgShift[0:nRows - shiftR, 0:nCols] = imgShift[shiftR:nRows, 0:nCols] \n", " imgShift[nRows - shiftR:nRows, 0:nCols] = False\n", " edges = np.logical_or(edges, np.logical_and(np.logical_not(picture), imgShift))\n", " #\n", " picture = np.logical_or(picture, edges)\n", " edgeSum = np.logical_or(edges, edgeSum)\n", " #\n", " return picture, edgeSum\n", "#\n", "#\n", "#\n", "def clusExpand(nInCluster, rPixelsCl, cPixelsCl, widthEx):\n", " #\n", " debugHere = False\n", " #\n", " if not hasattr(clusExpand, \"firstCall\"):\n", " clusExpand.firstCall = True\n", " print(\"------------------------------------------------------------------------------------\")\n", " print(\"Running clusExpand\")\n", " if debug:\n", " clusExpand.firstCall = True\n", " print(\"------------------------------------------------------------------------------------\")\n", " print(\"Running clusExpand\")\n", " #\n", " nClus = len(nInCluster)\n", " #\n", " # Temporary information \n", " lClus = np.amax(nInCluster)\n", " lFact = 1\n", " cPixelsHere = np.zeros(lFact*lClus).astype(int)\n", " rPixelsHere = np.zeros(lFact*lClus).astype(int)\n", " thisPic = np.zeros((nRows, nCols))\n", " thisPicEx = np.zeros((nRows, nCols))\n", " thisEdge = np.zeros((nRows, nCols))\n", " #\n", " # Information on expanded cluster\n", " nInClusEx = np.zeros(nClus).astype(int)\n", " iClusExSum = np.zeros(nClus)\n", " #\n", " # Information on pixels in expanded cluster\n", " lFact = 30\n", " indexClEx = np.zeros(lFact*lClus)\n", " cPixelsClEx = np.zeros(lFact*lClus).astype(int)\n", " rPixelsClEx = np.zeros(lFact*lClus).astype(int)\n", " iPixelsClEx = np.zeros(lFact*lClus)\n", " #\n", " # Information on edge\n", " nInEdge = np.zeros(nClus).astype(int)\n", " iEdgeSum = np.zeros(nClus)\n", " #\n", " # Information on pixels in expanded edge\n", " lFact = 30\n", " indexEdge = np.zeros(lFact*lClus)\n", " cPixelsEdge = np.zeros(lFact*lClus).astype(int)\n", " rPixelsEdge = np.zeros(lFact*lClus).astype(int)\n", " iPixelsEdge = np.zeros(lFact*lClus)\n", " #\n", " # Pictures before and after expanding and edges (all clusters)\n", " clusPic = np.zeros((nRows, nCols))\n", " clusPicEx = np.zeros((nRows, nCols))\n", " edgePic = np.zeros((nRows, nCols))\n", " #\n", " if clusExpand.firstCall:\n", " fig = plt.figure(figsize = (pltX, pltY))\n", " ax = fig.add_subplot(1, 1, 1)\n", " plt.title(\"Expanded clusters\", fontsize = 12)\n", " plt.xlabel('x pixel', fontsize = 12)\n", " plt.ylabel('y pixel', fontsize = 12)\n", " #\n", " # Figure control\n", " xOffMax = 2\n", " yOffMax = 2\n", " xOffMin = 30\n", " yOffMin = 30\n", " mSize = 0.001\n", " nCol = 0\n", " #\n", " nLastCl = 0\n", " nLastClEx = 0\n", " nLastEdge = 0\n", " for nC in range(0, nClus):\n", " #\n", " # Clusters\n", " nHereCl = nInCluster[nC]\n", " #\n", " # Indices of x and y pixels\n", " thisPic.fill(0) \n", " thisPic[rPixelsCl[nLastCl:nLastCl + nHereCl], cPixelsCl[nLastCl:nLastCl + nHereCl]] = 1\n", " clusPic += thisPic\n", " #\n", " nLastCl = nLastCl + nHereCl\n", " #\n", " # Expand cluster\n", " thisExp, thisEdge = expander(thisPic, widthEx, useDiag)\n", " nHereClEx = np.sum(thisExp).astype(int)\n", " nInClusEx[nC] = nHereClEx\n", " #\n", " # Indices of x and y pixels expanded cluster\n", " rPixelsHereEx, cPixelsHereEx = np.where(thisExp)\n", " thisPicEx.fill(0)\n", " thisPicEx[rPixelsHereEx, cPixelsHereEx] = 1\n", " clusPicEx += thisPicEx\n", " #\n", " # Intensities in pixels expanded cluster\n", " iPixelsHereEx = imgGrey[rPixelsHereEx, cPixelsHereEx]\n", " iClusExSum[nC] = np.sum(iPixelsHereEx)\n", " #\n", " indexClEx[nLastClEx:nLastClEx + nHereClEx] = nC*np.ones(nHereClEx)\n", " cPixelsClEx[nLastClEx:nLastClEx + nHereClEx] = cPixelsHereEx[:]\n", " rPixelsClEx[nLastClEx:nLastClEx + nHereClEx] = rPixelsHereEx[:]\n", " iPixelsClEx[nLastClEx:nLastClEx + nHereClEx] = iPixelsHereEx[:]\n", " #\n", " if clusExpand.firstCall:\n", " plt.scatter(cPixelsClEx[nLastClEx:nLastClEx + nHereClEx], rPixelsClEx[nLastClEx:nLastClEx + nHereClEx], \n", " s = mSize, c = colorTab[nCol], marker = 'o')\n", " rLab = np.amax(rPixelsClEx[nLastClEx:nLastClEx + nHereClEx]) + yOffMax\n", " cLab = np.amax(cPixelsClEx[nLastClEx:nLastClEx + nHereClEx]) + xOffMax\n", " plt.text(cLab, rLab, str(nC), color = colorTab[nCol])\n", " nCol = nCol + 1\n", " if nCol > nColTab - 1:\n", " nCol = 0 \n", " #\n", " nLastClEx = nLastClEx + nHereClEx\n", " #\n", " # Expanded edges\n", " nHereEdge = np.sum(thisEdge).astype(int)\n", " nInEdge[nC] = nHereEdge\n", " #\n", " # Indices of x and y pixels edge\n", " rPixelsHereEdge, cPixelsHereEdge = np.where(thisEdge)\n", " thisEdge.fill(0)\n", " thisEdge[rPixelsHereEdge, cPixelsHereEdge] = 1\n", " edgePic += thisEdge\n", " #\n", " # Intensities in pixels edge\n", " iPixelsHereEdge = imgGrey[rPixelsHereEdge, cPixelsHereEdge]\n", " iEdgeSum[nC] = np.sum(iPixelsHereEdge)\n", " #\n", " indexEdge[nLastEdge:nLastEdge + nHereEdge] = nC*np.ones(nHereEdge)\n", " cPixelsEdge[nLastEdge:nLastEdge + nHereEdge] = cPixelsHereEdge[:]\n", " rPixelsEdge[nLastEdge:nLastEdge + nHereEdge] = rPixelsHereEdge[:]\n", " iPixelsEdge[nLastEdge:nLastEdge + nHereEdge] = iPixelsHereEdge[:]\n", " nLastEdge = nLastEdge + nHereEdge\n", " #\n", " #\n", " if clusExpand.firstCall:\n", " plt.xlim(-0.05*nCols, 1.05*nCols)\n", " plt.ylim(-0.05*nRows, 1.05*nRows)\n", " plt.grid(color = 'g')\n", " print(\" \")\n", " plt.show()\n", " #\n", " clusExpand.firstCall = False\n", " #\n", " if debugHere:\n", " print(\" \")\n", " plt.figure(figsize = (pltX, pltY))\n", " plt.imshow(clusPic);\n", " plt.show()\n", " print(\" \")\n", " plt.figure(figsize = (pltX, pltY))\n", " plt.imshow(clusPicEx);\n", " plt.show()\n", " print(\" \")\n", " plt.figure(figsize = (pltX, pltY))\n", " plt.imshow(edgePic);\n", " plt.show()\n", " #\n", " return nInClusEx, rPixelsClEx, cPixelsClEx, iPixelsClEx, nInEdge, rPixelsEdge, cPixelsEdge, iPixelsEdge\n", "#\n", "#\n", "#\n", "def findRims(nInClus, rPixelsCl, cPixelsCl):\n", " #\n", " if not hasattr(findRims, \"firstCall\"):\n", " findRims.firstCall = True\n", " print(\"------------------------------------------------------------------------------------\")\n", " print(\"Running findRims\")\n", " if debug:\n", " findRims.firstCall = True\n", " print(\"------------------------------------------------------------------------------------\")\n", " print(\"Running findRims\")\n", " #\n", " useDiag = True\n", " #\n", " nClus = len(nInClus)\n", " lClus = np.amax(nInClus).astype(int)\n", " #\n", " # Information on cluster edge\n", " nInClusEd = np.zeros(nClus).astype(int)\n", " iClusEdSum = np.zeros(nClus)\n", " #\n", " # Information on pixels in cluster edge\n", " lFact = 3\n", " indexClEd = np.zeros(lFact*lClus)\n", " cPixelsClEd = np.zeros(lFact*lClus)\n", " rPixelsClEd = np.zeros(lFact*lClus)\n", " iPixelsClEd = np.zeros(lFact*lClus)\n", " #\n", " # Pictures for display/edge finding\n", " thisPic = np.zeros((nRows, nCols))\n", " thisEdge = np.zeros((nRows, nCols))\n", " thisEdgePic = np.zeros((nRows, nCols))\n", " #\n", " # Edge width required\n", " width = 1\n", " #\n", " if findRims.firstCall:\n", " fig = plt.figure(figsize = (pltX, pltY))\n", " ax = fig.add_subplot(1, 1, 1)\n", " plt.title(\"Clusters and their edges\", fontsize = 12)\n", " plt.xlabel('x pixel', fontsize = 12)\n", " plt.ylabel('y pixel', fontsize = 12)\n", " #\n", " # Figure control\n", " xOffMax = 2\n", " yOffMax = 2\n", " xOffMin = 30\n", " yOffMin = 30\n", " nCol = 0\n", " mSize = 0.001\n", " #\n", " nLastCl = 0\n", " nLastClEd = 0\n", " for nC in range(0, nClus):\n", " #\n", " # Clusters\n", " nHereCl = nInClus[nC]\n", " #\n", " # Indices of x and y pixels\n", " thisPic.fill(0)\n", " thisPic[rPixelsCl[nLastCl:nLastCl + nHereCl], cPixelsCl[nLastCl:nLastCl + nHereCl]] = 1\n", " #\n", " if findRims.firstCall:\n", " plt.scatter(cPixelsCl[nLastCl:nLastCl + nHereCl], rPixelsCl[nLastCl:nLastCl + nHereCl],\n", " s = mSize, c = colorTab[nCol], marker = 'o')\n", " rLab = np.amax(rPixelsCl[nLastCl:nLastCl + nHereCl]) + yOffMax\n", " cLab = np.amax(cPixelsCl[nLastCl:nLastCl + nHereCl]) + xOffMax\n", " plt.text(cLab, rLab, str(nC), color = colorTab[nCol])\n", " #\n", " nCol = nCol + 1\n", " if nCol > nColTab - 1:\n", " nCol = 0\n", " #\n", " nLastCl = nLastCl + nHereCl\n", " #\n", " # Edges\n", " thisEdge = edgeFinderIn(thisPic, width, useDiag)\n", " nHereClEd = np.sum(thisEdge).astype(int)\n", " #\n", " nInClusEd[nC] = nHereClEd\n", " #\n", " # Indices of x and y pixels\n", " rPixelsHereEd, cPixelsHereEd = np.where(thisEdge)\n", " thisEdgePic[rPixelsHereEd, cPixelsHereEd] = 1\n", " #\n", " # Intensities in pixels\n", " iPixelsHereEd = imgGrey[rPixelsHereEd, cPixelsHereEd]\n", " iClusEdSum[nC] = np.sum(iPixelsHereEd)\n", " #\n", " indexClEd[nLastClEd:nLastClEd + nHereClEd] = nC*np.ones(nInClusEd[nC])\n", " cPixelsClEd[nLastClEd:nLastClEd + nHereClEd] = cPixelsHereEd[0:nHereClEd]\n", " rPixelsClEd[nLastClEd:nLastClEd + nHereClEd] = rPixelsHereEd[0:nHereClEd]\n", " iPixelsClEd[nLastClEd:nLastClEd + nHereClEd] = iPixelsHereEd[0:nHereClEd]\n", " #\n", " if findRims.firstCall:\n", " plt.scatter(cPixelsHereEd, rPixelsHereEd, s = mSize, c = 'k', marker = 'o')\n", " rLab = np.amin(rPixelsHereEd[0:nHereClEd]) - yOffMin\n", " cLab = np.amin(cPixelsHereEd[0:nHereClEd]) - xOffMin\n", " plt.text(cLab, rLab, str(nC), color = 'k')\n", " #\n", " nLastClEd = nLastClEd + nHereClEd\n", " #\n", " if findRims.firstCall:\n", " plt.xlim(-0.05*nCols, 1.05*nCols)\n", " plt.ylim(-0.05*nRows, 1.05*nRows)\n", " plt.grid(color = 'g')\n", " print(\" \")\n", " plt.show()\n", " #\n", " findRims.firstCall = False\n", " #\n", " return nInClusEd, rPixelsClEd, cPixelsClEd, iPixelsClEd\n", "#\n", "#\n", "#\n", "def selCometsThree(nInClus, rPixCl, cPixCl, nInClusEd, rPixClEd, cPixClEd,\n", " nInHeadL, rPixHdL, cPixHdL, nInHdEdL, rPixHdEdL, cPixHdEdL,\n", " nInHeadH, rPixHdH, cPixHdH, nInHdEdH, rPixHdEdH, cPixHdEdH):\n", " #\n", " UseLowHead = False\n", " #\n", " if not hasattr(selCometsThree, \"firstCall\"):\n", " selCometsThree.firstCall = True\n", " print(\"------------------------------------------------------------------------------------\")\n", " print(\"Running selCometsThree\")\n", " print(\"UseLowHead\",UseLowHead)\n", " #\n", " if debug:\n", " selCometsThree.firstCall = True\n", " print(\"------------------------------------------------------------------------------------\")\n", " print(\"Running selCometsThree\")\n", " print(\"UseLowHead\",UseLowHead)\n", " #\n", " nClus = len(nInClus)\n", " lClus = np.amax(nInClus)\n", " lClusEd = np.amax(nInClusEd)\n", " nHeadL = len(nInHeadL)\n", " lHeadL = np.amax(nInHeadL)\n", " nHeadH = len(nInHeadH)\n", " lHeadH = np.amax(nInHeadH)\n", " #\n", " minRimPnts0 = 100\n", " minWheelPnts0 = 500\n", " rMin0 = 5\n", " rMax0 = nRows - rMin0\n", " cMin0 = 5\n", " cMax0 = nCols - cMin0\n", " rWidMin0 = 20\n", " rWidMax0 = 200\n", " cWidMin0 = 50\n", " cWidMax0 = 500\n", " maxRtoW = 0.1\n", " #\n", " if selCometsThree.firstCall:\n", " #\n", " # Figure text positions\n", " rTextRow = 3\n", " rTextCol = 3\n", " wTextRow = 28\n", " wTextCol = 28\n", " #\n", " nCol = 0\n", " #\n", " fig = plt.figure(figsize=(7, 5))\n", " plt.title(\"Selected clusters\", fontsize = 12)\n", " plt.xlabel('Column', fontsize = 12)\n", " plt.ylabel('Row', fontsize = 12)\n", " #\n", " pnts_num = np.zeros(nClus).astype(int)\n", " pnts_row = np.zeros((nClus, lClus)).astype(int)\n", " pnts_col = np.zeros((nClus, lClus)).astype(int)\n", " mean_row = np.zeros(nClus)\n", " mean_col = np.zeros(nClus)\n", " #\n", " edgs_num = np.zeros(nClus).astype(int)\n", " edgs_row = np.zeros((nClus, lClusEd)).astype(int)\n", " edgs_col = np.zeros((nClus, lClusEd)).astype(int)\n", " #\n", " head_num = np.zeros(nHeadL).astype(int)\n", " head_row = np.zeros((nHeadL, lHeadL)).astype(int)\n", " head_col = np.zeros((nHeadL, lHeadL)).astype(int)\n", " #\n", " sumWheel = np.zeros(nClus)\n", " sumHead = np.zeros(nHeadL)\n", " #\n", " mean_row = np.zeros(nHeadL)\n", " mean_col = np.zeros(nHeadL)\n", " #\n", " iSelRim = 0\n", " nLastCl = 0\n", " nLastClEd = 0\n", " for nC in range(0, nClus):\n", " #\n", " # Clusters\n", " nHereCl = nInClus[nC]\n", " #\n", " # Cluster edges\n", " nHereClEd = nInClusEd[nC] \n", " rMinClEd = np.amin(rPixClEd[nLastClEd:nLastClEd + nHereClEd])\n", " cMinClEd = np.amin(cPixClEd[nLastClEd:nLastClEd + nHereClEd])\n", " rMaxClEd = np.amax(rPixClEd[nLastClEd:nLastClEd + nHereClEd])\n", " cMaxClEd = np.amax(cPixClEd[nLastClEd:nLastClEd + nHereClEd])\n", " rWidClEd = rMaxClEd - rMinClEd\n", " cWidClEd = cMaxClEd - cMinClEd \n", " #\n", " rimToWheel = nInClusEd[nC]/nInClus[nC]\n", " #\n", " if (nHereClEd > minRimPnts0 and nHereCl > minWheelPnts0 and\n", " rMinClEd > rMin0 and rMaxClEd < rMax0 and cMinClEd > cMin0 and cMaxClEd < cMax0 and \n", " rWidClEd > rWidMin0 and rWidClEd < rWidMax0 and cWidClEd > cWidMin0 and cWidClEd < cWidMax0 and\n", " rimToWheel < maxRtoW):\n", " #\n", " # Look at low threshold heads\n", " nHeadsInClusL = 0\n", " nLastHd = 0\n", " nLastHdEd = 0\n", " for nHL in range(0, nHeadL):\n", " #\n", " nHereHd = nInHeadL[nHL]\n", " #\n", " # Low head edges\n", " nHereHdEd = nInHdEdL[nHL]\n", " rMinHdEd = np.amin(rPixHdEdL[nLastHdEd:nLastHdEd + nHereHdEd])\n", " cMinHdEd = np.amin(cPixHdEdL[nLastHdEd:nLastHdEd + nHereHdEd])\n", " rMaxHdEd = np.amax(rPixHdEdL[nLastHdEd:nLastHdEd + nHereHdEd])\n", " cMaxHdEd = np.amax(cPixHdEdL[nLastHdEd:nLastHdEd + nHereHdEd])\n", " rWidHdEd = rMaxHdEd - rMinHdEd\n", " cWidHdEd = cMaxHdEd - cMinHdEd\n", " #\n", " # Find any low threshold head rims completely within the cluster rim\n", " if (rMinHdEd > rMinClEd and rMaxHdEd < rMaxClEd and \n", " cMinHdEd > cMinClEd and cMaxHdEd < cMaxClEd):\n", " #\n", " nHeadsInClusL += 1\n", " nLastHdSelL = nLastHd\n", " nLastHdEdSelL = nLastHdEd\n", " nHereHdSelL = nHereHd\n", " nHereHdEdSelL = nHereHdEd\n", " #\n", " nLastHd = nLastHd + nHereHd\n", " nLastHdEd = nLastHdEd + nHereHdEd\n", " #\n", " # End of loop over high heads \n", " #\n", " # Look at high threshold heads\n", " nHeadsInClusH = 0\n", " nLastHd = 0\n", " nLastHdEd = 0\n", " for nHH in range(0, nHeadH):\n", " #\n", " # High heads\n", " nHereHd = nInHeadH[nHH]\n", " #\n", " # High head edges\n", " nHereHdEd = nInHdEdH[nHH]\n", " rMinHdEd = np.amin(rPixHdEdH[nLastHdEd:nLastHdEd + nHereHdEd])\n", " cMinHdEd = np.amin(cPixHdEdH[nLastHdEd:nLastHdEd + nHereHdEd])\n", " rMaxHdEd = np.amax(rPixHdEdH[nLastHdEd:nLastHdEd + nHereHdEd])\n", " cMaxHdEd = np.amax(cPixHdEdH[nLastHdEd:nLastHdEd + nHereHdEd])\n", " rWidHdEd = rMaxHdEd - rMinHdEd\n", " cWidHdEd = cMaxHdEd - cMinHdEd\n", " #\n", " # Find any high threshold head rims completely within the cluster rim\n", " if (rMinHdEd > rMinClEd and rMaxHdEd < rMaxClEd and \n", " cMinHdEd > cMinClEd and cMaxHdEd < cMaxClEd):\n", " #\n", " nHeadsInClusH += 1\n", " nLastHdSelH = nLastHd\n", " nLastHdEdSelH = nLastHdEd\n", " nHereHdSelH = nHereHd\n", " nHereHdEdSelH = nHereHdEd\n", " #\n", " nLastHd = nLastHd + nHereHd\n", " nLastHdEd = nLastHdEd + nHereHdEd\n", " #\n", " # End of loop over high heads \n", " #\n", " if nHeadsInClusL == 1 and nHeadsInClusH == 1:\n", " # \n", " # This rim accepted\n", " pnts_num[iSelRim] = nHereCl\n", " pnts_row[iSelRim, 0:nHereCl] = rPixCl[nLastCl:nLastCl + nHereCl].astype(int)\n", " pnts_col[iSelRim, 0:nHereCl] = cPixCl[nLastCl:nLastCl + nHereCl].astype(int)\n", " #\n", " edgs_num[iSelRim] = nHereClEd\n", " edgs_row[iSelRim, 0:nHereClEd] = rPixClEd[nLastClEd:nLastClEd + nHereClEd].astype(int)\n", " edgs_col[iSelRim, 0:nHereClEd] = cPixClEd[nLastClEd:nLastClEd + nHereClEd].astype(int)\n", " #\n", " if UseLowHead:\n", " head_num[iSelRim] = nHereHdSelL\n", " head_row[iSelRim, 0:nHereHdSelL] = rPixHdL[nLastHdSelL:nLastHdSelL + nHereHdSelL].astype(int)\n", " head_col[iSelRim, 0:nHereHdSelL] = cPixHdL[nLastHdSelL:nLastHdSelL + nHereHdSelL].astype(int)\n", " else:\n", " head_num[iSelRim] = nHereHdSelH\n", " head_row[iSelRim, 0:nHereHdSelH] = rPixHdH[nLastHdSelH:nLastHdSelH + nHereHdSelH].astype(int)\n", " head_col[iSelRim, 0:nHereHdSelH] = cPixHdH[nLastHdSelH:nLastHdSelH + nHereHdSelH].astype(int)\n", " #\n", " # Sum cluster intensities\n", " sumWheel[iSelRim] = np.sum(imgGrey[pnts_row[iSelRim, 0:pnts_num[iSelRim]], pnts_col[iSelRim, 0:pnts_num[iSelRim]]])\n", " #\n", " # Head intensities and mean head position\n", " sumHead[iSelRim] = np.sum(imgGrey[head_row[iSelRim, 0:head_num[iSelRim]], head_col[iSelRim, 0:head_num[iSelRim]]])\n", " mean_row[iSelRim] = (np.sum(imgGrey[head_row[iSelRim, 0:head_num[iSelRim]], \n", " head_col[iSelRim, 0:head_num[iSelRim]]]*\n", " head_row[iSelRim, 0:head_num[iSelRim]])/sumHead[iSelRim])\n", " mean_col[iSelRim] = (np.sum(imgGrey[head_row[iSelRim, 0:head_num[iSelRim]], \n", " head_col[iSelRim, 0:head_num[iSelRim]]]*\n", " head_col[iSelRim, 0:head_num[iSelRim]])/sumHead[iSelRim])\n", " #\n", " if selCometsThree.firstCall:\n", " plt.plot(pnts_col[iSelRim, 0:pnts_num[iSelRim]], pnts_row[iSelRim, 0:pnts_num[iSelRim]],\n", " color = colorTab[nCol])\n", " plt.plot(mean_col[iSelRim], mean_row[iSelRim], marker = '+', color = colorTab[nCol])\n", " row_label = np.amax(pnts_row[iSelRim, 0:pnts_num[iSelRim]]) + rTextRow\n", " col_label = np.amax(pnts_col[iSelRim, 0:pnts_num[iSelRim]]) + rTextCol\n", " plt.text(col_label, row_label, str(nC), color = colorTab[nCol])\n", " nCol = nCol + 1\n", " if nCol > nColTab - 1:\n", " nCol = 0\n", " #\n", " iSelRim += 1\n", " #\n", " # End of if statement, rim accepted\n", " #\n", " # End of if statement, cluster in fiducial region \n", " nLastCl = nLastCl + nHereCl\n", " nLastClEd = nLastClEd + nHereClEd\n", " #\n", " # End of loop over clusters\n", " nRimOut = iSelRim\n", " #\n", " if selCometsThree.firstCall:\n", " plt.xlim(-0.05*nCols, 1.05*nCols)\n", " plt.ylim(-0.05*nRows, 1.05*nRows)\n", " plt.grid(color = 'g')\n", " print(\" \")\n", " plt.show()\n", " print(\" \")\n", " print(\"Number of selected clusters (nRimOut)\",nRimOut)\n", " #\n", " # Remove overlaps\n", " #\n", " overlap = np.zeros(nRimOut).astype(bool)\n", " #\n", " for rimA in range(0, nRimOut - 1):\n", " rowA = edgs_row[rimA, 0:edgs_num[rimA]]\n", " colA = edgs_col[rimA, 0:edgs_num[rimA]]\n", " for rimB in range(rimA + 1, nRimOut):\n", " rowB = edgs_row[rimB, 0:edgs_num[rimB]]\n", " colB = edgs_col[rimB, 0:edgs_num[rimB]]\n", " for nP in range(0, edgs_num[rimA]):\n", " dist2 = (rowA[nP] - rowB)**2 + (colA[nP] - colB)**2\n", " if np.amin(dist2) < 1:\n", " overlap[rimA] = True\n", " overlap[rimB] = True\n", " break\n", " #\n", " #\n", " #\n", " if selCometsThree.firstCall:\n", " print(\" \")\n", " print(\"Remove the overlaps:\",np.linspace(0, nRimOut - 1, nRimOut)[overlap])\n", " #\n", " nCol = 0\n", " #\n", " fig = plt.figure(figsize = (pltX, pltY))\n", " plt.title(\"Rims no overlaps\", fontsize = 12)\n", " plt.xlabel('Column', fontsize = 12)\n", " plt.ylabel('Row', fontsize = 12)\n", " #\n", " l_json = nRimOut - int(np.sum(overlap))\n", " #\n", " j_num = np.zeros(l_json).astype(int)\n", " j_row = np.zeros((l_json, lHeadL))\n", " j_col = np.zeros((l_json, lHeadL))\n", " #\n", " w_num = np.zeros(l_json).astype(int)\n", " w_row = np.zeros((l_json, lClus))\n", " w_col = np.zeros((l_json, lClus))\n", " wm_row = np.zeros(l_json)\n", " wm_col = np.zeros(l_json)\n", " #\n", " h_num = np.zeros(l_json).astype(int)\n", " h_row = np.zeros((l_json, lHeadL))\n", " h_col = np.zeros((l_json, lHeadL))\n", " #\n", " j_rim = 0\n", " for iRim in range(0, nRimOut):\n", " if overlap[iRim]:\n", " continue\n", " j_num[j_rim] = edgs_num[iRim]\n", " j_row[j_rim, 0:j_num[j_rim]] = edgs_row[iRim, 0:edgs_num[iRim]]\n", " j_col[j_rim, 0:j_num[j_rim]] = edgs_col[iRim, 0:edgs_num[iRim]]\n", " #\n", " w_num[j_rim] = pnts_num[iRim]\n", " w_row[j_rim, 0:w_num[j_rim]] = pnts_row[iRim, 0:pnts_num[iRim]]\n", " w_col[j_rim, 0:w_num[j_rim]] = pnts_col[iRim, 0:pnts_num[iRim]]\n", " wm_row[j_rim] = mean_row[iRim]\n", " wm_col[j_rim] = mean_col[iRim]\n", " #\n", " h_num[j_rim] = head_num[iRim]\n", " h_row[j_rim, 0:h_num[j_rim]] = head_row[iRim, 0:head_num[iRim]]\n", " h_col[j_rim, 0:h_num[j_rim]] = head_col[iRim, 0:head_num[iRim]]\n", " #\n", " if selCometsThree.firstCall:\n", " plt.plot(pnts_col[iRim, 0:pnts_num[iRim]], pnts_row[iRim, 0:pnts_num[iRim]],\n", " color = colorTab[nCol])\n", " plt.plot(mean_col[iRim], mean_row[iRim], marker = '+', color = colorTab[nCol])\n", " row_label = np.amax(pnts_row[iRim, 0:pnts_num[iRim]]) + rTextRow\n", " col_label = np.amax(pnts_col[iRim, 0:pnts_num[iRim]]) + rTextCol\n", " plt.text(col_label, row_label, str(iRim), color = colorTab[nCol])\n", " nCol += 1\n", " if nCol > nColTab - 1:\n", " nCol = 0\n", " j_rim += 1\n", " #\n", " n_rim = j_rim\n", " if selCometsThree.firstCall:\n", " plt.xlim(-0.05*nCols, 1.05*nCols)\n", " plt.ylim(-0.05*nRows, 1.05*nRows)\n", " plt.grid(color = 'g')\n", " print(\" \")\n", " plt.show()\n", " #\n", " fig = plt.figure(figsize = (pltX, pltY))\n", " plt.title(\"Final rims\", fontsize = 12)\n", " plt.xlabel('Column', fontsize = 12)\n", " plt.ylabel('Row', fontsize = 12)\n", " #\n", " for j_rim in range(0, n_rim):\n", " plt.plot(j_col[j_rim, 0:j_num[j_rim]], j_row[j_rim, 0:j_num[j_rim]],\n", " color = colorTab[nCol])\n", " plt.plot(wm_col[j_rim], wm_row[j_rim], marker = '+', color = colorTab[nCol])\n", " row_label = np.amax(j_row[j_rim, 0:j_num[j_rim]]) + rTextRow\n", " col_label = np.amax(j_col[j_rim, 0:j_num[j_rim]]) + rTextCol\n", " plt.text(col_label, row_label, str(j_rim), color = colorTab[nCol])\n", " nCol += 1\n", " if nCol > nColTab - 1:\n", " nCol = 0\n", " #\n", " plt.xlim(-0.05*nCols, 1.05*nCols)\n", " plt.ylim(-0.05*nRows, 1.05*nRows)\n", " plt.grid(color = 'g')\n", " print(\" \")\n", " plt.show()\n", " #\n", " selCometsThree.firstCall = False\n", " #\n", " return n_rim, w_num, w_row, w_col, wm_row, wm_col, head_num, head_row, head_col\n", "#\n", "#\n", "#\n", "def calcQuads(nWheels, pnts_num, pnts_row, pnts_col, mean_row, mean_col, head_num, head_row, head_col):\n", " #\n", " if not hasattr(calcQuads, \"firstCall\"):\n", " calcQuads.firstCall = True\n", " print(\"------------------------------------------------------------------------------------\")\n", " print(\"Running calcQuads\")\n", " if debug:\n", " calcQuads.firstCall = True\n", " print(\"------------------------------------------------------------------------------------\")\n", " print(\"Running calcQuads\")\n", " #\n", " lWheels = 10*np.amax(pnts_num)\n", " #\n", " pnts_num_UL = np.zeros(nWheels).astype(int)\n", " pnts_row_UL = np.zeros((nWheels, lWheels)).astype(int)\n", " pnts_col_UL = np.zeros((nWheels, lWheels)).astype(int)\n", " #\n", " pnts_num_UR = np.zeros(nWheels).astype(int)\n", " pnts_row_UR = np.zeros((nWheels, lWheels)).astype(int)\n", " pnts_col_UR = np.zeros((nWheels, lWheels)).astype(int)\n", " #\n", " pnts_num_LL = np.zeros(nWheels).astype(int)\n", " pnts_row_LL = np.zeros((nWheels, lWheels)).astype(int)\n", " pnts_col_LL = np.zeros((nWheels, lWheels)).astype(int)\n", " #\n", " pnts_num_LR = np.zeros(nWheels).astype(int)\n", " pnts_row_LR = np.zeros((nWheels, lWheels)).astype(int)\n", " pnts_col_LR = np.zeros((nWheels, lWheels)).astype(int)\n", " #\n", " head_num_UL = np.zeros(nWheels).astype(int)\n", " head_row_UL = np.zeros((nWheels, lWheels)).astype(int)\n", " head_col_UL = np.zeros((nWheels, lWheels)).astype(int)\n", " #\n", " head_num_UR = np.zeros(nWheels).astype(int)\n", " head_row_UR = np.zeros((nWheels, lWheels)).astype(int)\n", " head_col_UR = np.zeros((nWheels, lWheels)).astype(int)\n", " #\n", " head_num_LL = np.zeros(nWheels).astype(int)\n", " head_row_LL = np.zeros((nWheels, lWheels)).astype(int)\n", " head_col_LL = np.zeros((nWheels, lWheels)).astype(int)\n", " #\n", " head_num_LR = np.zeros(nWheels).astype(int)\n", " head_row_LR = np.zeros((nWheels, lWheels)).astype(int)\n", " head_col_LR = np.zeros((nWheels, lWheels)).astype(int)\n", " #\n", " sumULwheel = np.zeros(nWheels)\n", " sumURwheel = np.zeros(nWheels)\n", " sumLLwheel = np.zeros(nWheels)\n", " sumLRwheel = np.zeros(nWheels)\n", " #\n", " sumULhead = np.zeros(nWheels)\n", " sumURhead = np.zeros(nWheels)\n", " sumLLhead = np.zeros(nWheels)\n", " sumLRhead = np.zeros(nWheels)\n", " #\n", " if calcQuads.firstCall:\n", " #\n", " # Text positions\n", " wTextRow = 2\n", " wTextCol = 2\n", " #\n", " fig = plt.figure(figsize=(7, 5))\n", " plt.title(\"Wheel quadrants\", fontsize = 12)\n", " plt.xlabel('Column', fontsize = 12)\n", " plt.ylabel('Row', fontsize = 12)\n", " #\n", " nCol = 0\n", " #\n", " for nW in range(0, nWheels):\n", " #\n", " # Wheels\n", " colHere = pnts_col[nW, 0:pnts_num[nW]]\n", " rowHere = pnts_row[nW, 0:pnts_num[nW]]\n", " #\n", " boolUL = np.logical_and(colHere < mean_col[nW], rowHere > mean_row[nW]) \n", " pnts_num_UL[nW] = np.sum(boolUL)\n", " pnts_row_UL[nW, 0:pnts_num_UL[nW]] = rowHere[boolUL]\n", " pnts_col_UL[nW, 0:pnts_num_UL[nW]] = colHere[boolUL]\n", " #\n", " boolUR = np.logical_and(colHere > mean_col[nW], rowHere > mean_row[nW]) \n", " pnts_num_UR[nW] = np.sum(boolUR)\n", " pnts_row_UR[nW, 0:pnts_num_UR[nW]] = rowHere[boolUR]\n", " pnts_col_UR[nW, 0:pnts_num_UR[nW]] = colHere[boolUR]\n", " #\n", " boolLL = np.logical_and(colHere < mean_col[nW], rowHere < mean_row[nW])\n", " pnts_num_LL[nW] = np.sum(boolLL)\n", " pnts_row_LL[nW, 0:pnts_num_LL[nW]] = rowHere[boolLL]\n", " pnts_col_LL[nW, 0:pnts_num_LL[nW]] = colHere[boolLL]\n", " #\n", " boolLR = np.logical_and(colHere > mean_col[nW], rowHere < mean_row[nW]) \n", " pnts_num_LR[nW] = np.sum(boolLR)\n", " pnts_row_LR[nW, 0:pnts_num_LR[nW]] = rowHere[boolLR]\n", " pnts_col_LR[nW, 0:pnts_num_LR[nW]] = colHere[boolLR]\n", " #\n", " # Heads\n", " colHere = head_col[nW, 0:head_num[nW]]\n", " rowHere = head_row[nW, 0:head_num[nW]]\n", " #\n", " boolUL = np.logical_and(colHere < mean_col[nW], rowHere > mean_row[nW])\n", " head_num_UL[nW] = np.sum(boolUL)\n", " head_row_UL[nW, 0:head_num_UL[nW]] = rowHere[boolUL]\n", " head_col_UL[nW, 0:head_num_UL[nW]] = colHere[boolUL]\n", " #\n", " boolUR = np.logical_and(colHere > mean_col[nW], rowHere > mean_row[nW])\n", " head_num_UR[nW] = np.sum(boolUR)\n", " head_row_UR[nW, 0:head_num_UR[nW]] = rowHere[boolUR]\n", " head_col_UR[nW, 0:head_num_UR[nW]] = colHere[boolUR]\n", " #\n", " boolLL = np.logical_and(colHere < mean_col[nW], rowHere < mean_row[nW]) \n", " head_num_LL[nW] = np.sum(boolLL)\n", " head_row_LL[nW, 0:head_num_LL[nW]] = rowHere[boolLL]\n", " head_col_LL[nW, 0:head_num_LL[nW]] = colHere[boolLL]\n", " #\n", " boolLR = np.logical_and(colHere > mean_col[nW], rowHere < mean_row[nW]) \n", " head_num_LR[nW] = np.sum(boolLR)\n", " head_row_LR[nW, 0:head_num_LR[nW]] = rowHere[boolLR]\n", " head_col_LR[nW, 0:head_num_LR[nW]] = colHere[boolLR]\n", " #\n", " sumULwheel[nW] = np.sum(imgGrey[pnts_row_UL[nW, 0:pnts_num_UL[nW]], pnts_col_UL[nW, 0:pnts_num_UL[nW]]])\n", " sumURwheel[nW] = np.sum(imgGrey[pnts_row_UR[nW, 0:pnts_num_UR[nW]], pnts_col_UR[nW, 0:pnts_num_UR[nW]]])\n", " sumLLwheel[nW] = np.sum(imgGrey[pnts_row_LL[nW, 0:pnts_num_LL[nW]], pnts_col_LL[nW, 0:pnts_num_LL[nW]]])\n", " sumLRwheel[nW] = np.sum(imgGrey[pnts_row_LR[nW, 0:pnts_num_LR[nW]], pnts_col_LR[nW, 0:pnts_num_LR[nW]]])\n", " #\n", " sumULhead[nW] = np.sum(imgGrey[head_row_UL[nW, 0:head_num_UL[nW]], head_col_UL[nW, 0:head_num_UL[nW]]])\n", " sumURhead[nW] = np.sum(imgGrey[head_row_UR[nW, 0:head_num_UR[nW]], head_col_UR[nW, 0:head_num_UR[nW]]])\n", " sumLLhead[nW] = np.sum(imgGrey[head_row_LL[nW, 0:head_num_LL[nW]], head_col_LL[nW, 0:head_num_LL[nW]]])\n", " sumLRhead[nW] = np.sum(imgGrey[head_row_LR[nW, 0:head_num_LR[nW]], head_col_LR[nW, 0:head_num_LR[nW]]])\n", " #\n", " # Wheels\n", " if calcQuads.firstCall:\n", " size = 0.01\n", " plt.scatter(pnts_col_UL[nW, 0:pnts_num_UL[nW]], pnts_row_UL[nW, 0:pnts_num_UL[nW]], \n", " s = size, color = 'r')\n", " plt.scatter(pnts_col_UR[nW, 0:pnts_num_UR[nW]], pnts_row_UR[nW, 0:pnts_num_UR[nW]], \n", " s = size, color = 'b')\n", " plt.scatter(pnts_col_LL[nW, 0:pnts_num_LL[nW]], pnts_row_LL[nW, 0:pnts_num_LL[nW]], \n", " s = size, color = 'orange')\n", " plt.scatter(pnts_col_LR[nW, 0:pnts_num_LR[nW]], pnts_row_LR[nW, 0:pnts_num_LR[nW]], \n", " s = size, color = 'm')\n", " #\n", " # Heads\n", " plt.scatter(head_col_UL[nW, 0:head_num_UL[nW]], head_row_UL[nW, 0:head_num_UL[nW]], \n", " s = size, color = 'w')\n", " plt.scatter(head_col_UR[nW, 0:head_num_UR[nW]], head_row_UR[nW, 0:head_num_UR[nW]], \n", " s = size, color = 'c')\n", " plt.scatter(head_col_LL[nW, 0:head_num_LL[nW]], head_row_LL[nW, 0:head_num_LL[nW]], \n", " s = size, color = 'yellow')\n", " plt.scatter(head_col_LR[nW, 0:head_num_LR[nW]], head_row_LR[nW, 0:head_num_LR[nW]], \n", " s = size, color = 'g')\n", " #\n", " row_label = np.amax(pnts_row_UR[nW, 0:pnts_num[nW]]) + wTextRow\n", " col_label = np.amax(pnts_col_UR[nW, 0:pnts_num[nW]]) + wTextCol\n", " plt.text(col_label, row_label, str(nW), color = 'k')\n", " #\n", " if calcQuads.firstCall:\n", " plt.plot(mean_col, mean_row, marker = '+', linestyle = '', color = 'k')\n", " #\n", " plt.xlim(-0.05*nCols, 1.05*nCols)\n", " plt.ylim(-0.05*nRows, 1.05*nRows)\n", " plt.grid(color = 'g')\n", " print(\" \")\n", " plt.show()\n", " #\n", " calcQuads.firstCall = False\n", " #\n", " return sumULwheel, sumLLwheel, sumURwheel, sumLRwheel, sumULhead, sumLLhead, sumURhead, sumLRhead\n", "#\n", "then = now\n", "now = datetime.datetime.now()\n", "print(\" \")\n", "print(\"Date and time\",str(now))\n", "print(\"Time since last check is\",str(now - then))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Check directory contents " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Volume in drive C is OS\n", " Volume Serial Number is 8A8B-1A1F\n", "\n", " Directory of C:\\Users\\green\\OneDrive\\OneDocuments\\Liverpool\\LivDat\\ImageAnalysis\\Images\\MCdata\\pBreak090\n", "\n", "06/05/2021 15:49