# -*- 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:/SESanalysistemp")#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


# #import the files names
# im_path='C:/SESanalysistemp'
# 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:/SESanalysistemp'
files= glob.glob(os.path.join(im_path,'*.jpg'))

out_directory = 'C:/SES_MARS_out'
#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)
        particle_minor_axis.append(min_axis)
        particle_area.append(um_area)
        particle_ESD.append(ESD)
        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:/SES_MARS_out/All_particle_measurements.csv')