# -*- coding: utf-8 -*-
"""
Created on Thu Apr 30 12:08:50 2026

@author: chuffard
"""



out_path='C:/SES_out'



#import the files names

im_path='C:/SES_in'

#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


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

#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 'Ext' 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= 0.0
    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_path,'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'])


    if os.path.isfile(os.path.join(out_path,os.path.basename(out_path)+'_roi.csv')):

        data.to_csv(os.path.join(out_path,os.path.basename(out_path)+'_roi.csv'), mode='a', header=False)

    else:

        data.to_csv(os.path.join(out_path,os.path.basename(out_path)+'_roi.csv'))

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_out/All_particle_measurements.csv')