# -*- coding: utf-8 -*-
"""
Created on Mon Mar 10 09:30:05 2025

@author: Simulator
"""

import cv2
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.fftpack import fft, fftfreq
from scipy.signal import find_peaks
from scipy.ndimage import gaussian_filter

def extract_diagonal(video_path):
    cap = cv2.VideoCapture(video_path)
    frame_count = int(cap.get(cv2.CAP_PROP_FRAME_COUNT))
    ret, frame = cap.read()
    
    
    if not ret:
        print("Failed to read the video.")
        return None
    
    height, width, _ = frame.shape
    diag_length = min(height, width)  # Ensuring we stay within bounds
    
    diagonals = np.zeros((frame_count, diag_length))
    
    frame_idx = 0
    while ret:
        gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)  # Convert to grayscale
        filtered_image = gaussian_filter(gray, sigma=(0, 10))
        diag_pixels = [filtered_image[height-1-i, i] for i in range(diag_length)]  # Extract diagonal
        diagonals[frame_idx, :] = diag_pixels
        frame_idx += 1
        ret, frame = cap.read()
    
    cap.release()
    #imgor = cv2.imread(r'C:\Users\Simulator\Desktop\Quantification\counter_TileScan 9_t105_s10_RAW_ch00.jpg')
    
    
    #height, width, _ = imgor.shape
    #diag_length = min(height, width)  # Ensuring we stay within bounds
    
    #diagonals = np.zeros((1, diag_length))
    
    
    #gray = cv2.cvtColor(imgor, cv2.COLOR_BGR2GRAY)  # Convert to grayscale
    #filtered_image = gaussian_filter(gray, sigma=(0, 10))
    #diag_pixels = [filtered_image[height-1-i, i] for i in range(diag_length)]  # Extract diagonal
    #diagonals = diag_pixels[0:1250]
    
    return diagonals

def save_diagonal_image(diagonals, output_path="diagonal_output.png"):
    plt.imsave(output_path, diagonals.T, cmap='gray', origin='lower')
    print(f"Saved output image as {output_path}")

def analyze_wavelength(diagonals):
    
    # Compute FFT across diagonal to find the most periodic frame
    # frame_fft = np.abs(fft(diagonals, axis=1))
    # best_frame_idx = np.argmax(np.max(frame_fft[1:], axis=1))  # Ignore DC component 
    # Extract intensity along diagonal for the best frame
    best_frame_series = diagonals

    
    M = len(best_frame_series)
    x = np.array(range(0, M))/9.33
    y = (np.array(best_frame_series) - np.min(best_frame_series)) / (np.max(best_frame_series) - np.min(best_frame_series))
    best_frame_series = np.ravel(best_frame_series)  # flatten to 1D if needed
    yf = fft(best_frame_series)
    xf = fftfreq(M, d=1)
    
    df = pd.DataFrame({"Frequency (frames⁻¹)": xf[:M//2], "FFT Magnitude": yf[:M//2]})
    df.to_excel("FFT_Spatial_Analysis.xlsx", index=False)
    
    
    plt.figure(figsize=(8, 5))
    plt.plot(xf[:M//2], np.abs(yf[:M//2]))
    plt.xlabel("Frequency (1/px)")
    plt.ylabel("Power")
    plt.title("Fourier Spectrum of Ripples")
    plt.savefig("fft_spectrumwave.png")
    plt.show()
    
    peak_indices, _ = find_peaks(np.abs(yf[:M//2]), height=10)
    if len(peak_indices) > 0:
        dominant_wave = xf[peak_indices[0]]
        print (dominant_wave)
    
    
    
    # Plot intensity over diagonal for the most periodic frame
    plt.figure(figsize=(10, 4))
    plt.plot(x, y )
    plt.xlabel("Diistance ($\mu$m)")
    plt.ylabel("Intensity")
    plt.title("Intensity Over Diagonal in a Frame")
    plt.legend()
    plt.show()

def main():
    video_path = r'H:\Supplementary Video 14.mov'  # Change this to your video file
    diagonals = extract_diagonal(video_path)
    if diagonals is not None:
        save_diagonal_image(diagonals)
        analyze_wavelength(diagonals)

if __name__ == "__main__":
    main()