# Simpson College MATH 290 Spring 2021 - PowerPollen Research # Authors: Georgia Porter, Eli Marriott, Alyssa Whitham, Ethan Humble # PowerPollen Mentor Contacts: Cheng Pan, Alex Austin # Advisor: Ross Sweet # Deepnote setup !apt update !apt install ffmpeg libsm6 libxext6 -y !pip install xlrd !pip install openpyxl # Imports import math import os import re import time import scipy.stats import cv2 import imutils import numpy as np import matplotlib.pyplot as plt import seaborn as sns import pandas as pd # Function to compute the number of grains in an image. Input file name and directory. # Output is an integer count. def grain_counter(fname, direc): f = os.path.join(direc, fname) resizeHeight = 1000 im = cv2.imread(f) im = imutils.resize(im, height=resizeHeight) im = cv2.cvtColor(im, cv2.COLOR_BGR2GRAY) im = cv2.blur(im, (9, 9)) ret, binary1 = cv2.threshold(im, 80, 255, cv2.THRESH_BINARY_INV) kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (5, 1)) closed = cv2.morphologyEx(binary1, cv2.MORPH_CLOSE, kernel) contours1, _ = cv2.findContours(closed, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE) grains = 0 minimum_area = 900 connected_grain_area = 2800 temp = 0 for i in contours1: area = cv2.contourArea(i) if area > minimum_area: if area > connected_grain_area: temp = area/connected_grain_area temp = math.ceil(temp) # Assumption of rounding up grains += temp else: grains += 1 return(grains) # Vectorize the grain_counter function. vgrain_counter = np.vectorize(grain_counter) # Function to plot the original image and the image with contours of grains added. # Input is a file name and directory, output is a pair of captioned images. def grain_imager(fname, direc): f = os.path.join(direc, fname) resizeHeight = 1000 im = cv2.imread(f) im = imutils.resize(im, height=resizeHeight) im2 = im.copy() gray1 = cv2.cvtColor(im, cv2.COLOR_BGR2GRAY) blurred1 = cv2.blur(gray1, (9, 9)) ret, binary1 = cv2.threshold(blurred1, 80, 255, cv2.THRESH_BINARY_INV) kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (5, 1)) closed = cv2.morphologyEx(binary1, cv2.MORPH_CLOSE, kernel) contours1, _ = cv2.findContours(closed, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE) minimum_area = 900 connected_grain_area = 2800 temp = 0 grains = 0 for i in contours1: area = cv2.contourArea(i) if area > minimum_area: cv2.drawContours(im, [i], -1, (0, 255, 0), 5) if area > connected_grain_area: temp = area/connected_grain_area temp = math.ceil(temp) # Assumption of rounding up grains += temp else: grains += 1 text = "Pollen Grains:{grains}.".format(grains=grains) dy = 30 y0 = resizeHeight - dy for i, line in enumerate(text.split('\n')): y = y0 + i*dy counted = cv2.putText(im, line, (10, y), cv2.FONT_HERSHEY_DUPLEX, 1, (255, 0, 0), 3) titles = ["Original Image", "Counted"] images = [im2, counted] plt.figure(figsize=(20,12)) for i in range(2): plt.subplot(2,3,i+1),plt.imshow(images[i]) plt.title(titles[i]) plt.xticks([]),plt.yticks([]) return(plt.draw()) def tube_counter(fname, direc, match_im): resizeHeight = 1000 tubecount = 0 min_max_list = list() # Read comparison image and test image f = os.path.join(direc, fname) img1 = cv2.imread(match_im) img1 = imutils.resize(img1, height=resizeHeight) img2 = cv2.imread(f) img2 = imutils.resize(img2, height=resizeHeight) # Convert to grayscale gray1 = cv2.cvtColor(img1, cv2.COLOR_BGR2GRAY) gray2 = cv2.cvtColor(img2, cv2.COLOR_BGR2GRAY) # Sharpen Image gaussian_3 = cv2.GaussianBlur(gray2, (0,0), 2.0) gray2 = cv2.addWeighted(gray2, 1.5, gaussian_3, -0.25, -30) points = sns.kdeplot(gray2.ravel(), cut = 0).get_lines()[0].get_data() # pull point data from density plot combined = pd.DataFrame(np.column_stack(points),columns=['x','y']) # format point data as data frame resizeWidth = img2.shape[1] numPixels = resizeHeight * resizeWidth combined['freq'] = combined['y'] * numPixels freqCutoff = 950 # Set the cutoff value for pixel value frequency lower = (combined['freq'] < freqCutoff) & (combined.shift()['freq'] > freqCutoff) upper = (combined['freq'] < freqCutoff) & (combined.shift(-1)['freq'] > freqCutoff) min_max = combined[lower | upper] # Filter for index differences larger than one if min_max.shape[0] == 3: minPixel = min_max.iloc[0]['x'] maxPixel = min_max.iloc[1]['x'] elif min_max.shape[0] == 4: minPixel = min_max.iloc[1]['x'] maxPixel = min_max.iloc[2]['x'] else: minPixel = 50 maxPixel = 150 # Threshold Images using MIN & MAX blurred1 = cv2.blur(gray1, (9, 9)) ret, binary1 = cv2.threshold(blurred1, 80, 255, cv2.THRESH_BINARY_INV) binary2 = cv2.inRange(gray2, minPixel, maxPixel) # Find contours kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (5, 1)) closed1 = cv2.morphologyEx(binary1, cv2.MORPH_CLOSE, kernel) closed2 = cv2.morphologyEx(binary2, cv2.MORPH_CLOSE, kernel) contours1, _ = cv2.findContours(closed1,cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE) cnt1 = contours1[0] contours2, _ = cv2.findContours(closed2,cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE) # Run MatchShape function for i in contours2: perimeter = cv2.arcLength(i, True) area = cv2.contourArea(i) testreturn = cv2.matchShapes(cnt1, i, 1, 0.0) if testreturn > 1 and testreturn < 1000 and area > 400: cv2.drawContours(img2, [i], -1, (0, 255, 0), 5) tubecount += 1 return(tubecount) # Vectorize the tube_counter function. vtube_counter = np.vectorize(tube_counter) # Tube count visualization def tube_imager(fname, direc, match_im): # Variables resizeHeight = 1000 tubecount = 0 min_max_list = list() # Read comparison image and test image f = os.path.join(direc, fname) img1 = cv2.imread(match_im) img1 = imutils.resize(img1, height=resizeHeight) img2 = cv2.imread(f) img2 = imutils.resize(img2, height=resizeHeight) img2_orig = img2.copy() # Convert to grayscale gray1 = cv2.cvtColor(img1, cv2.COLOR_BGR2GRAY) gray2 = cv2.cvtColor(img2, cv2.COLOR_BGR2GRAY) # Sharpen Image gaussian_3 = cv2.GaussianBlur(gray2, (0,0), 2.0) gray2 = cv2.addWeighted(gray2, 1.5, gaussian_3, -0.25, -30) points = sns.kdeplot(gray2.ravel(), cut = 0).get_lines()[0].get_data() # pull point data from density plot combined = pd.DataFrame(np.column_stack(points),columns=['x','y']) # format point data as data frame resizeWidth = img2.shape[1] numPixels = resizeHeight * resizeWidth combined['freq'] = combined['y'] * numPixels freqCutoff = 950 # Set the cutoff value for pixel value frequency lower = (combined['freq'] < freqCutoff) & (combined.shift()['freq'] > freqCutoff) upper = (combined['freq'] < freqCutoff) & (combined.shift(-1)['freq'] > freqCutoff) min_max = combined[lower | upper] # Filter for index differences larger than one if min_max.shape[0] == 3: minPixel = min_max.iloc[0]['x'] maxPixel = min_max.iloc[1]['x'] elif min_max.shape[0] == 4: minPixel = min_max.iloc[1]['x'] maxPixel = min_max.iloc[2]['x'] else: minPixel = 50 maxPixel = 150 # Threshold Images using MIN & MAX blurred1 = cv2.blur(gray1, (9, 9)) ret, binary1 = cv2.threshold(blurred1, 80, 255, cv2.THRESH_BINARY_INV) binary2 = cv2.inRange(gray2, minPixel, maxPixel) # Find contours kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (5, 1)) closed1 = cv2.morphologyEx(binary1, cv2.MORPH_CLOSE, kernel) closed2 = cv2.morphologyEx(binary2, cv2.MORPH_CLOSE, kernel) contours1, _ = cv2.findContours(closed1,cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE) cnt1 = contours1[0] contours2, _ = cv2.findContours(closed2,cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE) # Run MatchShape function for i in contours2: perimeter = cv2.arcLength(i, True) area = cv2.contourArea(i) testreturn = cv2.matchShapes(cnt1, i, 1, 0.0) if testreturn > 1 and testreturn < 1000 and area > 400: cv2.drawContours(img2, [i], -1, (0, 255, 0), 5) tubecount += 1 text = "Tube count:{testreturn}.".format(testreturn=testreturn) dy = 30 y0 = resizeHeight - dy for i, line in enumerate(text.split('\n')): y = y0 + i*dy counted = cv2.putText(img2, line, (10, y),cv2.FONT_HERSHEY_DUPLEX, 1, (0, 0, 25), 3) # Plot step by step images in a grid titles = ["Original Image", "Counted"] images = [img2_orig, counted] plt.figure(figsize=(20,12)) for i in range(2): plt.subplot(2,3,i+1),plt.imshow(images[i]) plt.title(titles[i]) plt.xticks([]),plt.yticks([]) return(plt.draw()) #Run Grain Counter directory = '/work/Math290Spring2021/Pure' pure_data = pd.read_excel('/work/Math290Spring2021/Image_count_result.xlsx', sheet_name='Pure', engine = 'openpyxl') match = '/work/Pure (1)-Cropped.tif' df = pd.DataFrame(os.listdir(directory), columns = ['FileName']) # Create data frame of file names. df['Sample'] = df['FileName'].str.extract('\(([^)]+)') # Make column with just sample number. df = df.dropna() # Drop missing values corresponding to non-sample files, e.g., Readme df['Sample'] = df['Sample'].astype(str).astype(int) # Convert sample column to integer. df['est_grain'] = vgrain_counter(df['FileName'], directory) # Compute estimated number of grains for each sample. df = pd.merge(df, pure_data[['Sample', 'Total Pollen', 'with Tube']], on = 'Sample', how = 'left') # Merge with actual counts from Excel file. df['fracdiff_grain'] = (df['est_grain'] - df['Total Pollen']) / df['Total Pollen'] # Compute the fractional difference in estimated and true counts. df['absfracdiff_grain'] = df['fracdiff_grain'].abs() print(df['absfracdiff_grain'].describe()) # Summary statistics of the fractional differences. plt.figure() df['absfracdiff_grain'].hist(bins = 20) # Plot histogram of fractional differences. plt.figure() df.boxplot(column='absfracdiff_grain') # Print boxplot of fractional differences. grain_imager("Pure (23).tif", directory) # View sample image and counted grains. #Run Tube Counter df['est_tube'] = vtube_counter(df['FileName'], directory, match) df['fracdiff_tube'] = (df['est_tube'] - df['with Tube']) / df['with Tube'] df_wtubes = df df['absfdiff_tube'] = df['fracdiff_tube'].abs() print(df['absfdiff_tube'].describe()) # Summary statistics of the fractional differences. plt.figure() df['absfdiff_tube'].hist(bins = 20) # Plot histogram of fractional differences. plt.figure() df.boxplot(column='absfdiff_tube') # Print boxplot of fractional differences. print(df.to_string()) tube_imager('Pure (72).tif', directory, match) print(df_wtubes['fracdiff_tube'].abs().median()) print(df_wtubes['fracdiff_grain'].abs().median())