# -*- coding: utf-8 -*-
"""
Created on Mon Mar 10 06:11:22 2025

@author: Simulator
"""

import cv2
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
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()
    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):
    frame_count, diag_length = diagonals.shape
    

    
    # Compute FFT across frames to find the most periodic pixel
    pixel_fft = np.abs(fft(diagonals[1:,0:1200], axis=0))
    best_pixel_idx = np.argmax(np.max(pixel_fft[1:], axis=0))  # Ignore DC component
    
    # # 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 time series for the best pixel
    best_pixel_series = diagonals[:, best_pixel_idx]


    
    N = len(best_pixel_series)
    x = np.array(range(0, N))/0.33
    y = (np.array(best_pixel_series) - np.min(best_pixel_series)) / (np.max(best_pixel_series) - np.min(best_pixel_series))
    

    
    
    # Plot intensity over time for the most periodic pixel
    plt.figure(figsize=(10, 4))
    plt.plot(x, y)
    #plt.plot(best_pixel_series, label=f'Pixel {best_pixel_idx}')
    plt.xlabel("Time (minutes)")
    plt.ylabel("Intensity")
    plt.title("Intensity Over Time for a Pixel")
    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()