# -*- coding: utf-8 -*-
"""
Created on Tue Mar  7 12:35:03 2023

@author: chuffard
"""

#!/Users/cdurkin/opt/anaconda3/bin/python3 #why this folder? Do I need to have all the images in this folder?

# import os library
import os
 
# change the current directory
# to specified directory
os.chdir(r"C:/Users/chuffard/Documents/MATLAB")#Folder where the images are

os.getcwd()

import numpy as np #why add the "as np" part?
import pandas as pd
import matplotlib.pyplot as plt
import os , sys
import scipy.ndimage as ndi
import cv2 #why does this not have an "as ... part?
from skimage import measure
import glob

pip install opencv-python 

#import the files names
im_path='C:/Users/chuffard/Documents/MATLAB'
files= glob.glob(os.path.join(im_path,'*.jpg'))


#calculate the mean of the previous ten images
files=np.sort(files)
photo_mean=None #why is this none?
photo_max=None #why is this none?
count=0

#for this notebook, just need to choose the test image.  Here it is image number 51 in the list
file_count=51
while count < 10:
    file_count=file_count-1
    photo=cv2.imread(files[file_count])
    photo_hsv=cv2.cvtColor(photo,cv2.COLOR_BGR2HSV)
    photo_value=photo_hsv[:,:,2]
    photo_value=cv2.normalize(photo_value,None,0, 255,cv2.NORM_MINMAX)
    if np.array(photo_mean).any() == None:
        photo_mean=photo_value/10
        photo_max=photo_value.copy()
    else:
        photo_mean=photo_mean + photo_value/10
        max_values=photo_max.astype(int)-photo_value.astype(int)
        photo_max[max_values<0]=photo_value[max_values<0]
        
    count=count+1
photo_mean=photo_mean.astype(int)
plt.figure()
plt.imshow(photo_mean,cmap=plt.cm.gray)
plt.title('Mean of last 10 images')
plt.figure()
plt.imshow(photo_max,cmap=plt.cm.gray)


photo1=cv2.imread(files[51])
photo1_hsv=cv2.cvtColor(photo1,cv2.COLOR_BGR2HSV)
#HSV = Hue, Saturation, Value.  We are using the value which is basically greyscale
photo1_value=photo1_hsv[:,:,2]
plt.figure()
plt.imshow(photo1_value,cmap=plt.cm.gray)
plt.title('Example image to process, gray value')

photo1_norm=cv2.normalize(photo1_value,None,0, 255,cv2.NORM_MINMAX)
photo_denoise1=photo1_norm-photo_max.astype(int)
if np.min(photo_denoise1)<0:
    photo_denoise=photo_denoise1+abs(np.min(photo_denoise1))
else:
    photo_denoise=photo_denoise1.copy()
photo_denoise=cv2.convertScaleAbs(photo_denoise)
photo_denoise=cv2.normalize(photo_denoise,None,0, 255,cv2.NORM_MINMAX)
plt.figure()
plt.imshow(photo_denoise,cmap=plt.cm.gray)
plt.title('Image with background noise and scratches removed')

photo1_norm=cv2.normalize(photo1_value,None,0, 255,cv2.NORM_MINMAX)
photo_denoise1=photo1_norm-photo_mean
if np.min(photo_denoise1)<0:
    photo_denoise=photo_denoise1+abs(np.min(photo_denoise1))
else:
    photo_denoise=photo_denoise1.copy()
#photo_denoise=cv2.convertScaleAbs(photo_denoise)
photo_denoise=cv2.normalize(photo_denoise,None,0, 255,cv2.NORM_MINMAX)
plt.figure()
plt.imshow(photo_denoise,cmap=plt.cm.gray)
plt.title('Image with background noise and scratches removed')

#This is a median blur of the mean image, and inverted
photo_blur=cv2.medianBlur(cv2.convertScaleAbs(photo_mean),251)
photo_blur=abs(photo_blur.astype(int)-255)
plt.figure()
plt.imshow(photo_blur,cmap=plt.cm.gray)
plt.title("Inverted, median blur of mean image")

#Create an image with no background
photo_nobg=photo_denoise-photo_blur
if np.min(photo_nobg)<0:
    photo_nobg=photo_nobg+abs(np.min(photo_nobg))


plt.figure()
plt.imshow(photo1_value.astype(int),cmap=plt.cm.gray)
plt.title('Origianl greyscale image of interest')
photo_nobg=cv2.normalize(photo_nobg,None,0, 255,cv2.NORM_MINMAX)
plt.figure()
plt.imshow(photo_nobg.astype(int),cmap=plt.cm.gray)
plt.title('Background removed image of interest')


circle_markers=np.zeros_like(photo_mean)

#Find the edges in the blurred image
circle_markers=cv2.Canny(cv2.convertScaleAbs(photo_blur),6,6)
##Also use the edges in the mean image.
circle_markers3=cv2.Canny(cv2.convertScaleAbs(photo_mean),10,10)
circle_markers[circle_markers3==255]=255
#To help connect the circle, mask out the darkest pixels in the blured image too
circle_markers[photo_blur<np.mean(photo_blur[500:1500,500:1500])]=255

#Close the circle
kernel = np.ones((5, 5), np.uint8)
circle_marker_dilation=cv2.dilate(circle_markers,kernel,1)  
contours, hierarchy = cv2.findContours(circle_marker_dilation, cv2.RETR_EXTERNAL , cv2.CHAIN_APPROX_NONE)
circle_markers2=np.zeros_like(circle_markers)
for cnt in contours:
    cv2.drawContours(circle_markers2,[cnt],-1,255,-1)
circle_markers2=cv2.erode(circle_markers2,kernel,1)  

plt.figure()
plt.imshow(circle_markers2,cmap=plt.cm.gray)
plt.title('filled edges')

#Now have it find the circle for the largest object
contours,_= cv2.findContours(cv2.convertScaleAbs(circle_markers2), cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)
areas = [cv2.contourArea(c) for c in contours]
sorted_areas = np.sort(areas)
cnt=contours[areas.index(sorted_areas[-1])] 
(x,y),radius = cv2.minEnclosingCircle(cnt)
center = (int(x),int(y))
radius = int(radius)-50
plt.figure()
mask=np.zeros_like(photo_nobg)
cv2.circle(mask,center,radius,(255),-1)
plt.imshow(mask,cmap=plt.cm.gray)
plt.title('automatic circle fitting to largest object in above image')

#Use adaptive thesholding to find potential particles
markers=cv2.adaptiveThreshold(cv2.convertScaleAbs(photo_nobg), 255,cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY_INV, 801, 3)

#Also find the very darkest pixels to help close some of the gaps in the adaptive threshold. Label those dark pixels=50
dark_markers=np.zeros_like(markers)
dark_markers[photo_nobg<0.4*np.median(photo_mean[mask==255])]=50

#Add those dark pixels into the adaptive thesholded pixels and label them all value=255
markers[dark_markers==50]=255
markers[mask==0]=0

#Get rid of smallest particle noise
kernel = np.ones((5, 5), np.uint8)   
marker_erosion=cv2.erode(markers,kernel,1)  
marker_erosion=cv2.erode(marker_erosion,kernel,1)  
marker_dilation=cv2.dilate(marker_erosion,kernel,1)  
marker_dilation=cv2.dilate(marker_dilation,kernel,1)  

#Fill in the shapes and label them with value=100
contours, hierarchy = cv2.findContours(marker_dilation, cv2.RETR_EXTERNAL , cv2.CHAIN_APPROX_NONE)
markers2=np.zeros_like(markers)
for cnt in contours:
    cv2.drawContours(markers2,[cnt],-1,100,-1)
#Identify edges and label them with value = 50
edges=cv2.Canny(cv2.convertScaleAbs(photo_nobg),50,80)
edges[edges==255]=100

#Sum together the markers, dark pixels, and edges so that each particle has values ranging between 50 and 200, depending on how it was detected
combined=markers2+edges+dark_markers
plt.figure()
plt.imshow(markers2,cmap=plt.cm.gray)
plt.figure()
plt.imshow(combined,cmap=plt.cm.gray)

num_labels, label_image = cv2.connectedComponents(markers2)
sharp_particles=[]
sharp_img=np.zeros_like(label_image)
for n in np.arange(1,num_labels):
    pixels=combined[label_image==n]
    if len(pixels[pixels>=150])/len(pixels)>0:
        sharp_particles.append(n)
        sharp_img[label_image==n]=255
            
detected_img=sharp_img.astype('ubyte')
label_image_detect = cv2.connectedComponentsWithStats(detected_img)
plt.figure()
plt.imshow(detected_img,cmap=plt.cm.gray)

particle_edges=cv2.Canny(detected_img,100,200,5)
kernel = np.ones((5, 5), np.uint8)
edge_dilation=cv2.dilate(particle_edges,kernel,1)  

example_img=photo1.copy()
example_img[edge_dilation==255]=[127,0,255]
plt.figure()
plt.imshow(example_img)

properties = measure.regionprops(label_image_detect[1])

count=0
particle_area = []
particle_ESD = []
particle_perimeter = []
particle_major_axis= []
particle_minor_axis= []
particle_numberID = []
file_name = []
bounding_box = []
scale= (0.07194)
scale_area = scale**2 
file='test'
for x in properties:
        count = count +1
        px_area=x.area
        um_area = px_area / scale_area
        min_axis = x.minor_axis_length / scale
        maj_axis = x.major_axis_length / scale
        perim = x.perimeter / scale
        ESD = 2*(np.sqrt(um_area/np.pi))
        particle_perimeter.append(perim)
        particle_major_axis.append(maj_axis)
        particle_minor_axis.append(min_axis)
        particle_area.append(um_area)
        particle_ESD.append(ESD)
        particle_numberID.append(count)
        file_name.append(file)
        bounding_box.append(str([x.bbox[0],x.bbox[2],x.bbox[1],x.bbox[3]]))

#Input all data into a dataframe.#
data = pd.DataFrame(np.stack((particle_numberID, particle_area, particle_ESD, particle_minor_axis, particle_major_axis, particle_perimeter, file_name, bounding_box),-1),columns=['Number','Area','ESD','minor_length','major_length','perimeter','file_name', 'bounding_box'])


#Start here

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os , sys
import scipy.ndimage as ndi
import cv2
from skimage import measure
import glob

im_path='C:/Users/chuffard/Documents/SES TEST/subset test'
files= glob.glob(os.path.join(im_path,'*.jpg'))

out_directory = 'C:/Users/chuffard/Documents/SES_TEST'
#try out_directory = 'C:/Users/chuffard/Documents/SES_TEST' ok taht worked

count=0
particle_area = []
particle_ESD = []
particle_perimeter = []
particle_major_axis= []
particle_minor_axis= []
particle_numberID = []
file_name = []
bounding_box = []

files= (glob.glob(os.path.join(im_path,'*.jpg')))
Ext_files=[]

for x in files:
    if '_SESExt_' in x:
        Ext_files.append(x)
files=np.sort(np.array(Ext_files))

file_count=int(len(files))-1

for f in files:
    photo_mean=None
    count=0
    file_count2=file_count
    #Create a mean image of the previous 3 images
    while count<10:
        file_count2=file_count2-1
        photo=cv2.imread(files[file_count2])
        photo_hsv=cv2.cvtColor(photo,cv2.COLOR_BGR2HSV)
        photo_value=photo_hsv[:,:,2]
        photo_value=cv2.normalize(photo_value,None,0, 255,cv2.NORM_MINMAX)
        if np.array(photo_mean).any() == None:
            photo_mean=photo_value/10
        else:
            photo_mean=photo_mean + photo_value/10
        count=count+1
    
    photo_mean=photo_mean.astype(int)
   
   #Read the image and subtract the mean to get the image without background
    photo=cv2.imread(files[file_count])
    photo_hsv=cv2.cvtColor(photo,cv2.COLOR_BGR2HSV)
    photo_value=photo_hsv[:,:,2]
    photo_norm=cv2.normalize(photo_value,None,0,255,cv2.NORM_MINMAX)
    photo_denoise1=photo_norm-photo_mean
    if np.min(photo_denoise1)<0:
        photo_denoise=photo_denoise1+abs(np.min(photo_denoise1))
    else:
        photo_denoise=photo_denoise1.copy()
    photo_denoise=cv2.normalize(photo_denoise,None,0,255,cv2.NORM_MINMAX)
    photo_blur=cv2.medianBlur(cv2.convertScaleAbs(photo_mean),251)
    photo_blur=abs(photo_blur.astype(int)-255)
    photo_nobg=photo_denoise-photo_blur
    if np.min(photo_nobg)<0:
        photo_nobg=photo_nobg+abs(np.min(photo_nobg))
    photo_nobg=cv2.normalize(photo_nobg,None,0, 255,cv2.NORM_MINMAX)
        
    circle_markers=np.zeros_like(photo_mean)
    circle_markers=cv2.Canny(cv2.convertScaleAbs(photo_blur),6,6)
    circle_markers3=cv2.Canny(cv2.convertScaleAbs(photo_mean),10,10)
    circle_markers[circle_markers3==255]=255
    circle_markers[photo_blur<np.mean(photo_blur[500:1500,500:1500])]=255
    kernel = np.ones((5, 5), np.uint8)
    circle_marker_dilation=cv2.dilate(circle_markers,kernel,1)  
    contours, hierarchy = cv2.findContours(circle_marker_dilation, cv2.RETR_EXTERNAL , cv2.CHAIN_APPROX_NONE)
    circle_markers2=np.zeros_like(circle_markers)
    for cnt in contours:
        cv2.drawContours(circle_markers2,[cnt],-1,255,-1)
    circle_markers2=cv2.erode(circle_markers2,kernel,1)  
    
    contours,_= cv2.findContours(cv2.convertScaleAbs(circle_markers2), cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)
    areas = [cv2.contourArea(c) for c in contours]
    sorted_areas = np.sort(areas)
    cnt=contours[areas.index(sorted_areas[-1])] 
    (x,y),radius = cv2.minEnclosingCircle(cnt)
    center = (int(x),int(y))
    radius = int(radius)-50
    mask=np.zeros_like(cv2.convertScaleAbs(photo_nobg))
    cv2.circle(mask,center,radius,(255),-1)
    
    markers=cv2.adaptiveThreshold(cv2.convertScaleAbs(photo_nobg), 255,cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY_INV,801, 3)
    dark_markers=np.zeros_like(markers)
    dark_markers[photo_nobg<0.4*np.mean(photo_mean[mask==255])]=50
    markers[dark_markers==50]=255   
    markers[mask==0]=0
    kernel = np.ones((5, 5), np.uint8)
    marker_erosion=cv2.erode(markers,kernel,1) 
    marker_erosion=cv2.erode(marker_erosion,kernel,1) 
    marker_dilation=cv2.dilate(marker_erosion,kernel,1)  
    marker_dilation=cv2.dilate(marker_dilation,kernel,1)  
    contours, hierarchy = cv2.findContours(marker_dilation, cv2.RETR_EXTERNAL , cv2.CHAIN_APPROX_NONE)
    markers2=np.zeros_like(markers)
    for cnt in contours:
        cv2.drawContours(markers2,[cnt],-1,100,-1)
    edges=cv2.Canny(cv2.convertScaleAbs(photo_nobg),50,80)
    edges[edges==255]=100
    combined=markers2+edges+dark_markers
    
    num_labels, label_image = cv2.connectedComponents(markers2)
    sharp_particles=[]
    sharp_img=np.zeros_like(label_image)
    for n in np.arange(1,num_labels):
        pixels=combined[label_image==n]
        if len(pixels[pixels>=150])/len(pixels)>0: #also saying that atleast 1% of the pixels have to be infocus edges
            sharp_particles.append(n)
            sharp_img[label_image==n]=255
            
    detected_img=sharp_img.astype('ubyte')
    label_image_detect = cv2.connectedComponentsWithStats(detected_img)

    #Create image with detected particles outlined for quality checking
    particle_edges=cv2.Canny(detected_img,100,200,5)
    kernel = np.ones((5, 5), np.uint8)
    edge_dilation=cv2.dilate(particle_edges,kernel,1)  
    example_img=photo.copy()
    example_img[edge_dilation==255]=[127,0,255]
    cv2.imwrite(os.path.join(out_directory,'Outline_'+str(os.path.basename(files[file_count]))),example_img)

    #Record measurements of detected particles and save in csv file
    scale= 0.07194
    scale_area = scale**2   #square pixels per square micron
    properties = measure.regionprops(label_image_detect[1])
    for x in properties:
        count = count +1
        px_area=x.area
        um_area = px_area / scale_area
        min_axis = x.minor_axis_length / scale
        maj_axis = x.major_axis_length / scale
        perim = x.perimeter / scale
        ESD = 2*(np.sqrt(um_area/np.pi))
        particle_perimeter.append(perim)
        particle_major_axis.append(maj_axis)diam
        particle_minor_axis.append(min_axis)talkative
        particle_area.append(um_area)
        particle_ESD.append(ESD)nyari apa
        particle_numberID.append(count)
        file_name.append(str(os.path.basename(files[file_count])))
        bounding_box.append(str([x.bbox[0],x.bbox[2],x.bbox[1],x.bbox[3]]))
    file_count=file_count-1

data = pd.DataFrame(np.stack((particle_numberID, particle_area, particle_ESD, particle_minor_axis, particle_major_axis, particle_perimeter, file_name, bounding_box),-1),columns=['Number','Area','ESD','minor_length','major_length','perimeter','file_name', 'bounding_box'])
data.to_csv('C:/Users/chuffard/Documents/MATLAB/All_particle_measurements.csv')