# -*- coding: utf-8 -*-
"""
Created on Wed Aug 26 15:44:24 2026

@author: chuffard
"""


import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

import os , sys

import cv2

import glob

from datetime import datetime

import os.path
import warnings

from pathlib import Path
import math





def create_bg_max_img(files, f,num_of_bg_imgs ):

    photo_max=[]

    #Create a max pixel image of the previous 20 images

    #But if within the first 20 images taken, create a max pixel image of the next 20 images that were taken afterward

    #This means we can't start calculating data until atleast 41 images have been collected.

    count=num_of_bg_imgs

    while count>0:

        if f>=num_of_bg_imgs:

            photo=cv2.imread(files[f-count])

        else:

            photo=cv2.imread(files[f+count])

        photo_value=photo[:,:,2]

        photo_value=cv2.normalize(photo_value,None,0, 255,cv2.NORM_MINMAX)

        if len(photo_max)==0:

            photo_max=photo_value.copy()

        else:

            max_values=photo_max.astype(int)-photo_value.astype(int)

            photo_max[max_values<0]=photo_value[max_values<0]

        count=count-1

    return(photo_max)



def find_circle_mask(photo_max):

    circle_markers=cv2.Canny(cv2.convertScaleAbs(photo_max),25,25)

    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(photo_max)

    cv2.circle(mask,center,radius,(255),-1)

    return(mask)



def measure_atn(mask, photo_max,files, f):

    photo=cv2.imread(files[f])

    photo_red=photo[:,:,2] ##cv2 reads images as BGR, so red is the third channel

    photo_red_mask=photo_red[mask==255]

    photo_red_mask[photo_red_mask==0]=1

    photo_max_mask=photo_max[mask==255]

    photo_max_mask[photo_max_mask==0]=1

    attenuation=np.mean(-np.log(photo_red_mask/photo_max_mask))

    file_name=os.path.basename(files[f])

    data=pd.DataFrame([[file_name,attenuation]],columns=['file_name','atn'])

    return(data)



out_path='W:/Analyses/2024 MARS deployment'


#import the files names

im_path='C:/SES_in'





fileslist = glob.glob("C:/SES_in/*.jpg")



#for some reason it's making the second slash a backslash instead of forward. Replace it

fileslist = [path.replace('\\', '/') for path in fileslist]



filesdf = pd.DataFrame(fileslist, columns=['filenames'])

if os.path.isdir(out_path)==False:

    os.mkdir(out_path)



###Crissy I added/modified this code to illiminate warnings###

filenames = filesdf[filesdf['filenames'].str.contains('Ext')]
filenames = filesdf[filesdf['filenames'].str.contains('5min')]

dates= filenames['filenames'].str.split('_').str[-1]

collect_time=filenames['filenames'].str.split('_').str[-3]



files=pd.DataFrame()

files['filenames']=filenames

files['date']=dates

files['collect_time']=collect_time

files=files.sort_values(by=['date'])

sorted_filelist = files['filenames'].tolist()



#Identify the files that needs to be processed

#Open existing data file and find the name of the last image file that was analyzed and analyze every image after that one

        

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

    data_link=open(os.path.join(out_path,os.path.basename(out_path)+'_atn.csv'),'r')

    last_line=data_link.readlines()[-1]

    data_link.close()

    for x in np.flip(np.arange(0,len(sorted_filelist))):

        if sorted_filelist[x].split('/')[-1] in last_line:

            files_to_be_analyzed=np.arange(x+1,len(sorted_filelist))

            break

#If no data file exists, then just start with the first image

else:

    files_to_be_analyzed=np.arange(0,len(sorted_filelist))



#Create an empty dataframe for new data to be recorded

num_of_bg_imgs=20

if len(sorted_filelist)>(num_of_bg_imgs*2):

    for f in files_to_be_analyzed:

        all_atn_data=pd.DataFrame(columns=['file_name','atn'])

        photo_max=create_bg_max_img(sorted_filelist,f,num_of_bg_imgs)

        mask=find_circle_mask(photo_max)

        atn_data=measure_atn(mask,photo_max,sorted_filelist,f)

        all_atn_data=pd.concat([all_atn_data,atn_data])

        #Save the data to the existing file, or create a new file

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

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

        else:

            all_atn_data.to_csv(os.path.join(out_path,os.path.basename(out_path)+'_atn.csv'))

 #       cv2.imwrite(os.path.join(out_path,'Outline_'+str(os.path.basename(sorted_filelist[f]))),example_img)
 
 
 #############################
 #%% Now convert to POC flux
 


MARS2024atn = r"W:/Analyses/2024 MARS deployment/all_atn_data.csv"
MARS2024atn = pd.read_csv(MARS2024atn)

# Extract the pattern into a new column
MARS2024atn['file_name'] = MARS2024atn['file_name'].astype(str).str.strip()

MARS2024atn[['date', 'time']] = MARS2024atn['file_name'].str.extract(r'(\d{6})-(\d{6})')

# Convert to actual datetime objects for analysis
MARS2024atn['timestamp'] = pd.to_datetime(MARS2024atn['date'] + MARS2024atn['time'], format='%y%m%d%H%M%S')

MARS2024atn.columns



MARS2024atn["DateTime"] = pd.to_datetime(MARS2024atn['timestamp'])
#convert atn to POC flux

#from Estapa et al. 2024 optical sediment trap calibration. also accounting for image and collect area since those aren't in our og atn calculation

SEStraparea = 0.5;
#imgarea = length(find(isfinite(bgr(:))))./(0.07194.^2)./(1000000.^2);  % pix / (pix/um^2) / (um^2/m^2)#this was in Meg's example code but it's not for SES

#This is adapted to pyton from matlab code Meg gave us
SESimgarea = math.pi*(25/1000/2)**2  # diameter of lower end of funnel is 25 mm, compute area in m^2
SEScolltime = 5/(60*24) #in days
MARS2024atn['F_atn']= MARS2024atn['atn']/SEScolltime*SESimgarea/SEStraparea


#These are from SES at Station M in 2024 paper
MARS2024atn['POCmodel_mg']= 10**3.4 * MARS2024atn['F_atn']**0.94
MARS2024atn['POCmodel_g']= MARS2024atn['POCmodel_mg']/1000

MARS2024atn['POCmodel_g'].mean()#(0.2753247346815628)
MARS2024atn['POCmodel_g'].plot()
plt.show()

MARS2024atn.to_csv('//Thalassa/ProjectLibrary/902305_Event_Detection_with_SES/Analyses/2024 MARS deployment/MARS2024atn_POC_red.csv')