% 基本參數
lambda = 632.8e-6;    % 波長（mm）
k = 2*pi/lambda;      % 波數
z = 0.3;              % 記錄距離（mm）
pix_size = 8e-3;      % 像素尺寸（mm）
N = 1024;             % 全息圖分辨率
L = N*pix_size;       % 全息圖尺寸（mm）

obj = imread('lena.jpg');
obj_gray = rgb2gray(obj);
[M,N] = size(obj_gray);
obj_amp = im2double(obj_gray)/255;  % 振幅分佈
obj_phase = 2*pi*rand(M,N);         % 隨機相位分佈

% 物光場與參考光干涉
H = exp(1j*k*z/(2*z)) * exp(-1j*pi*lambda*z/(2*z)*(ones(M,N)));  % 參考光
hologram = obj_amp .* exp(1j*obj_phase) + H;  % 干涉場
hologram = hologram ./ (abs(hologram) + eps);  # 歸一化
imwrite(abs(hologram), 'hologram.tif');       # 保存全息圖

% 提取包裹相位
wrapped_phase = angle(hologram);

% 相位展開（最小二乘法）
[unwrapped_phase, ~] = lsq_unwrap(wrapped_phase, 2*pi);

% 可視化
figure;
subplot(1,2,1); imshow(wrapped_phase, []); title('包裹相位');
subplot(1,2,2); imshow(unwrapped_phase, []); title('展開相位');

function [unwrapped, reliability] = lsq_unwrap(wrapped, pixel_size)
    % 輸入：wrapped - 包裹相位（rad）
    % 輸出：unwrapped - 展開相位（rad）
    [rows, cols] = size(wrapped);
    reliability = ones(rows, cols);
    
    % 構建線性方程組 Ax = b
    A = zeros(2*rows*cols, 2*rows*cols);
    b = zeros(2*rows*cols, 1);
    
    % 填充矩陣（基於相位梯度約束）
    k = 0;
    for i = 1:rows
        for j = 1:cols
            idx = (i-1)*cols + j;
            if i < rows
                A(2*idx-1,:) = [1, -exp(1j*k*pixel_size), 0, 0];
                b(2*idx-1) = wrapped(i,j) - wrapped(i+1,j);
            end
            if j < cols
                A(2*idx,:) = [1, 0, -exp(1j*k*pixel_size), 0];
                b(2*idx) = wrapped(i,j) - wrapped(i,j+1);
            end
        end
    end
    
    % 最小二乘求解
    x = A\b;
    unwrapped = reshape(x(1:2:end), rows, cols);
    reliability = abs(x(2:2:end));
end

function recon = fresnel_reconstruct(hologram, z, lambda, pix_size)
    [N, ~] = size(hologram);
    k = 2*pi/lambda;

    % 預處理
    H = fftshift(fft2(hologram));

    % 頻域傳播因子
    fx = (-N/2:N/2-1)*(1/(N*pix_size));
    [FX, FY] = meshgrid(fx, fx);
    Hiz = exp(1j*k*z/(2*z)) .* exp(-1j*pi*lambda*z*(FX.^2 + FY.^2));

    # 逆變換
    recon = fftshift(ifft2(ifftshift(H .* Hiz)));
    recon = recon ./ (abs(recon) + eps);
end

function clean_recon = remove_zero_order(recon)
    [c, s] = contourlet(recon, 'max_level', 3);
    c(1,:) = 0; c(2,:) = 0;  % 抑制低頻分量
    clean_recon = icontourlet(c, s);
end

% 衍射重建
recon = fresnel_reconstruct(hologram, z, lambda, pix_size);

% 零級像消除
clean_recon = remove_zero_order(recon);

% 顯示結果
figure;
subplot(1,2,1); imshow(abs(recon), []); title('原始重建像');
subplot(1,2,2); imshow(abs(clean_recon), []); title('去零級像後');

% 分層物體建模
layers = load('microchip.mat');  % 包含各層相位信息
for layer = 1:size(layers, 3)
    hologram_layer = generate_hologram(layers(:,:,layer));
    recon_layer = fresnel_reconstruct(hologram_layer, z, lambda, pix_size);
    figure; imshow(recon_layer); title(['Layer ', num2str(layer)]);
end

% 模擬振動相位調製
vibration_phase = 0.1*sin(2*pi*50*t);  % 50Hz振動
hologram_vib = hologram .* exp(1j*vibration_phase);
recon_vib = fresnel_reconstruct(hologram_vib, z, lambda, pix_size);

% 計算重建質量指標
def calculate_metrics(original, recon):
    psnr = psnr(recon, original);
    ssim_val = ssim(recon, original);
    mae = mean(abs(recon(:) - original(:)));
    return psnr, ssim_val, mae

original = imread('reference.png');
metrics = calculate_metrics(original, recon_clean);
disp(['PSNR: ', num2str(metrics(1)), ' dB']);
disp(['SSIM: ', num2str(metrics(2))]);
disp(['MAE: ', num2str(metrics(3))]);

% 多波長重建
lambda_red = 650e-6;  lambda_green = 532e-6;  lambda_blue = 450e-6;
recon_red = fresnel_reconstruct(hologram_red, z, lambda_red, pix_size);
recon_green = fresnel_reconstruct(hologram_green, z, lambda_green, pix_size);
recon_blue = fresnel_reconstruct(hologram_blue, z, lambda_blue, pix_size);
recon_rgb = cat(3, recon_red, recon_green, recon_blue);

% 加載預訓練模型
net = load('phase_unwrap_unet.mat');
unwrapped_phase = predict(net, wrapped_phase);