对光学图像还原与分析的修正与补充（以及对MWORKS的补充）

主要有以下几个可能的原因导致拟合结果全为1：

Airy函数的定义问题：在径向距离r=0处，Airy函数的计算会出现除以零的情况，可能导致数值不稳定。
初始参数猜测不合理：不合理的初始参数可能导致拟合算法无法正确收敛。
数据预处理问题：如果CCD数据未正确归一化或存在噪声，可能影响拟合结果。
拟合目标函数定义错误：拟合强度时，可能需要直接拟合强度而不是幅度。

下面，我将提供一个经过修正的完整Matlab代码，解决上述问题，并添加三维光强图的可视化。同时，我会详细解释

% 清空环境变量，关闭所有图形窗口
clear;
close all;
clc;

%% 1. 读取和预处理CCD数据
% 读取CCD图像数据
I_raw = imread('ccd_data.png'); % 请将'ccd_data.png'替换为您的实际文件名

% 将数据转换为双精度类型以便计算
I = double(I_raw);

% 如果图像是彩色的，转换为灰度图
if size(I, 3) == 3
    I = rgb2gray(I);
end

% 显示原始光强度图像
figure;
imshow(I, []);
title('原始CCD光强度数据');

%% 2. 可视化光强度的三维图
figure;
surf(I, 'EdgeColor', 'none');
colormap jet;
colorbar;
view(2); % 从上方观察
title('原始CCD光强度三维图');

%% 3. 计算光场幅值和傅里叶变换
% 计算光场幅值（假设相位信息可忽略）
field_amplitude = sqrt(I);

% 执行二维傅里叶变换并移频
F = fftshift(fft2(field_amplitude));

% 计算傅里叶变换的强度
F_intensity = abs(F).^2;

% 显示空间频谱的三维图
figure;
surf(log(F_intensity + 1), 'EdgeColor', 'none');
colormap jet;
colorbar;
view(2); % 从上方观察
title('空间频谱（光场的傅里叶变换）三维图');

%% 4. 定义径向坐标和展开数据用于拟合
% 获取图像尺寸
[rows, cols] = size(I);

% 生成网格坐标
[X, Y] = meshgrid(1:cols, 1:rows);

% 计算中心坐标
centerX = cols / 2;
centerY = rows / 2;

% 计算每个像素到中心的径向距离
R = sqrt((X - centerX).^2 + (Y - centerY).^2);

% 将矩阵展开为列向量
R_flat = R(:);
I_flat = I(:);

%% 5. 定义Airy斑的理论模型函数（直接拟合光强度）
% Airy斑强度分布函数
airy_function = @(params, r) params(1) * (2 * besselj(1, params(2) * r) ./ (params(2) * r)).^2;

% 处理r=0的情况，使用极限值
airy_function = @(params, r) params(1) * ( (2 * besselj(1, params(2) * r) ./ (params(2) * r + (r==0)) ).^2 .* (r ~=0) + 1 * (r ==0));

% 参数初始猜测值 [I0, k]
initial_params = [max(I_flat), pi/100]; % 适当调整初始k值

%% 6. 应用最小二乘法进行拟合
% 为了提高拟合效率，选择中心区域的数据进行拟合
% 定义拟合半径阈值（根据数据选择合适值）
fit_radius = min(centerX, centerY) / 2; % 只使用中心一半的数据

% 筛选拟合数据
fit_indices = R_flat <= fit_radius;
R_fit = R_flat(fit_indices);
I_fit = I_flat(fit_indices);

% 定义最小二乘法的目标函数
objective_function = @(params) airy_function(params, R_fit) - I_fit;

% 设置最小二乘法选项
options = optimoptions('lsqnonlin', 'Display', 'iter', 'TolFun',1e-6,'TolX',1e-6, 'MaxIterations', 1000, 'MaxFunctionEvaluations', 2000);

% 执行非线性最小二乘拟合
[params_fit, resnorm, residuals] = lsqnonlin(objective_function, initial_params, [], [], options);

% 显示拟合结果
disp('拟合参数：');
disp(['I0 = ', num2str(params_fit(1))]);
disp(['k = ', num2str(params_fit(2))]);

%% 7. 计算拟合后的强度分布并显示
% 计算拟合后的强度分布
I_model = airy_function(params_fit, R);

% 显示拟合的强度图像
figure;
surf(I_model, 'EdgeColor', 'none');
colormap jet;
colorbar;
view(2); % 从上方观察
title('拟合后的强度分布三维图');

% 对比原始和拟合的强度分布
figure;
plot(R_fit, I_fit, 'b.', 'MarkerSize', 10);
hold on;
plot(R_fit, airy_function(params_fit, R_fit), 'r-', 'LineWidth', 2);
xlabel('径向距离');
ylabel('强度');
legend('观测强度', '拟合强度');
title('强度分布对比');
grid on;
hold off;

%% 8. 还原原始物体信息
% 提取拟合参数k
k_fit = params_fit(2);

% 已知参数（根据实际情况替换）
lambda = 632.8e-9; % 波长（米），例如He-Ne激光器
focal_length = 0.5; % 焦距（米），请根据实际情况替换

% 计算孔径半径
aperture_radius = (1.22 * lambda * focal_length) / (k_fit);

% 显示估计的孔径半径
fprintf('估计的孔径半径：%.6e 米\n', aperture_radius);

%% 9. 数据滤波和噪声处理（可选）
% 创建高斯滤波器
H = fspecial('gaussian', size(I), 10); % 高斯滤波器，标准差为10

% 对原始图像进行滤波
I_filtered = imfilter(I, H, 'replicate');

% 重新计算光场幅值和傅里叶变换
field_amplitude_filtered = sqrt(I_filtered);
F_filtered = fftshift(fft2(field_amplitude_filtered));
F_intensity_filtered = abs(F_filtered).^2;

% 显示滤波后的空间频谱三维图
figure;
surf(log(F_intensity_filtered + 1), 'EdgeColor', 'none');
colormap jet;
colorbar;
view(2); % 从上方观察
title('滤波后的空间频谱三维图');

%% 10. 结束语
disp('数据分析完成。');

运行步骤

准备数据：
- 确保您的CCD数据文件（例如ccd_data.png）位于Matlab的当前工作目录下，或者提供完整的文件路径。
设置参数：
- 根据您的实验条件，调整lambda（波长）和focal_length（焦距）的值。
运行脚本：
- 将上述代码保存为analysis_corrected.m，然后在Matlab命令窗口中运行：

检查结果：

脚本将依次显示原始光强图、空间频谱三维图、拟合后的强度分布三维图以及拟合对比图。
在命令窗口中，您将看到估计的孔径半径。

MWORK代码

Julia代码，用于读取CCD采集的光强度数据，进行傅里叶变换，应用最小二乘法拟合，并还原衍射物体的信息。请确保您的Julia环境中安装了所需的包，并根据实际情况调整文件名和参数。

确保安装了以下Julia包：

using Pkg
Pkg.add(["Images", "ImageIO", "Plots", "FFT", "Optim", "SpecialFunctions", "FileIO"])

# 清空环境变量
workspace()

# 导入必要的包
using Images
using ImageIO
using Plots
using FFT
using Optim
using SpecialFunctions
using FileIO

# 设置绘图后端（可选，根据需求选择，例如 :gr()）
# plotly() # 使用 Plotly 后端
gr() # 使用 GR 后端

# 关闭所有图形窗口（在Julia中无需显式关闭）

## 1. 读取和预处理CCD数据
# 读取CCD图像数据
I_raw = load("hag.png") # 请将'hag.png'替换为您的实际文件名

# 如果图像是彩色的，转换为灰度图
if ndims(I_raw) == 3
    I = Gray.(I_raw)
else
    I = I_raw
end

# 将数据转换为浮点类型以便计算
I = Float64.(I)

# 显示原始光强度图像
plot()
heatmap(I, color=:gray, title="原始CCD光强度数据", axis=false)
gui() # 显示图像窗口

## 2. 计算光场幅值和傅里叶变换
# 计算光场幅值（假设相位信息可忽略）
field_amplitude = sqrt.(I)

# 执行二维傅里叶变换
F = fftshift(fft.(fft.(field_amplitude, dims=1), dims=2))

# 计算傅里叶变换的强度
F_intensity = abs.(F).^2

# 显示空间频谱
plot()
heatmap(log.(F_intensity .+ 1), color=:gray, title="空间频谱（光场的傅里叶变换）", axis=false)
gui()

## 3. 定义径向坐标和展开数据用于拟合
# 获取图像尺寸
rows, cols = size(I)

# 生成网格坐标
X = repeat(1:cols, inner=rows)
Y = repeat(1:rows, outer=cols)
X = reshape(X, rows, cols)
Y = reshape(Y, rows, cols)

# 计算中心坐标
centerX = cols / 2
centerY = rows / 2

# 计算每个像素到中心的径向距离
R = sqrt.((X .- centerX).^2 .+ (Y .- centerY).^2)

# 将矩阵展开为向量
R_flat = vec(R)
I_flat = vec(I)

## 4. 定义Airy斑的理论模型函数
# Airy斑强度分布函数
function airy_function(params, r)
    I0, k = params
    # 处理r=0的情况，避免除以零
    intensity = zeros(length(r))
    for i in 1:length(r)
        if r[i] == 0
            intensity[i] = I0
        else
            intensity[i] = I0 * (2 * besselj(1, k * r[i]) / (k * r[i]))^2
        end
    end
    return intensity
end

## 5. 应用最小二乘法进行拟合
# 参数初始猜测值 [I0, k]
I0_initial = maximum(I_flat)
k_initial = π / 100.0
initial_params = [I0_initial, k_initial]

# 定义最小二乘法的目标函数
function objective(params, r, I_obs)
    return airy_function(params, r) - I_obs
end

# 为了提高拟合效率，选择中心区域的数据进行拟合
# 定义拟合半径阈值（根据数据选择合适值）
fit_radius = min(centerX, centerY) / 2.0 # 只使用中心一半的数据

# 筛选拟合数据
fit_indices = findall(r -> r <= fit_radius, R_flat)
R_fit = R_flat[fit_indices]
I_fit = I_flat[fit_indices]

# 定义优化问题
optimize_objective(params) = objective(params, R_fit, I_fit)

# 执行非线性最小二乘拟合
result = optimize(optimize_objective, initial_params, NelderMead())

# 提取拟合参数
params_fit = Optim.minimizer(result)
I0_fit, k_fit = params_fit

# 显示拟合结果
println("拟合参数：")
println("I0 = ", I0_fit)
println("k = ", k_fit)

## 6. 计算拟合后的强度分布并显示
# 计算拟合后的强度分布
I_model = airy_function(params_fit, R)

# 显示拟合的强度图像
plot()
heatmap(I_model, color=:jet, title="拟合后的强度分布", axis=false)
gui()

# 对比原始和拟合的强度分布
using Statistics

# 只对拟合区域进行绘图
plot(R_fit, I_fit, seriestype=:scatter, label="观测强度", markersize=2)
plot!(R_fit, airy_function(params_fit, R_fit), label="拟合强度", linewidth=2)
xlabel!("径向距离")
ylabel!("强度")
title!("强度分布对比")
legend()
gui()

## 7. 还原原始物体信息
# 已知参数（根据实际情况替换）
lambda = 632.8e-9 # 波长（米），例如He-Ne激光器
focal_length = 0.5 # 焦距（米），请根据实际情况替换

# 计算孔径半径
aperture_radius = (1.22 * lambda * focal_length) / k_fit

# 显示估计的孔径半径
println("估计的孔径半径：", aperture_radius, " 米")

## 8. 数据滤波和噪声处理（可选）
using ImageFiltering

# 创建高斯滤波器
σ = 10.0
I_filtered = imfilter(I, Kernel.gaussian(Float64, σ))

# 重新计算光场幅值和傅里叶变换
field_amplitude_filtered = sqrt.(I_filtered)
F_filtered = fftshift(fft.(fft.(field_amplitude_filtered, dims=1), dims=2))
F_intensity_filtered = abs.(F_filtered).^2

# 显示滤波后的空间频谱
plot()
heatmap(log.(F_intensity_filtered .+ 1), color=:gray, title="滤波后的空间频谱", axis=false)
gui()

## 9. 结束语
println("数据分析完成。")

代码说明与关键点

包导入与安装：
- Images 和 ImageIO：用于读取和处理图像。
- Plots：用于绘图和可视化。
- FFT：进行傅里叶变换。
- Optim：进行优化和拟合。
- SpecialFunctions：提供Bessel函数等特殊函数。
- ImageFiltering：进行图像滤波。
读取和预处理CCD数据：
- 使用 load 函数读取图像。
- 检查图像是否为彩色，如果是则转换为灰度图。
- 将图像数据转换为浮点数类型，便于后续计算。
傅里叶变换：
- 计算光场幅值的平方根，假设相位信息可忽略。
- 使用 fft 进行二维傅里叶变换，并使用 fftshift 将零频移至中心。
- 计算傅里叶变换的强度，并使用对数尺度进行可视化。
径向坐标和数据展开：
- 生成网格坐标，计算每个像素到中心的径向距离。
- 将二维矩阵展开为一维向量，便于拟合。
定义Airy斑模型：
- 定义了一个处理 r = 0 情况的 airy_function 函数，以避免除以零错误。
最小二乘法拟合：
- 选择中心区域的数据进行拟合，以减少边缘噪声的影响。
- 使用 Optim.jl 的 optimize 函数进行非线性最小二乘拟合。这里使用了 NelderMead 优化方法，适用于非线性优化问题。
- 提取拟合参数 I0_fit 和 k_fit。
拟合结果可视化：
- 显示拟合后的强度分布。
- 绘制观测强度与拟合强度的对比图，检查拟合效果。
还原原始物体信息：
- 使用拟合得到的 k_fit 计算孔径半径，基于衍射理论公式。
数据滤波和噪声处理（可选）：
- 使用高斯滤波器对原始图像进行平滑处理，减少噪声影响。
- 重新计算傅里叶变换，并显示滤波后的空间频谱。
结束语：
- 输出分析完成的提示。

关键注意事项

Airy函数的定义：
- 在 airy_function 中，特别处理了 r = 0 的情况，直接赋予峰值强度 I0，避免除以零错误。
初始参数的选择：
- 初始猜测值 I0_initial 和 k_initial 需要根据实际数据进行调整。过小或过大的初始值可能导致优化失败或陷入局部最优。
优化方法的选择：
- 使用了 NelderMead 方法，适用于非线性优化。如果拟合效果不理想，可以尝试其他优化方法，如 BFGS 或 TrustRegion。
数据筛选：
- 为了提高拟合的准确性，仅使用中心区域的数据进行拟合，避免边缘区域的噪声影响。
可视化：
- 使用 heatmap 和 surf 的替代方法进行二维图像的可视化。
- 使用 Plots 包的 gui() 函数显示图像窗口。
滤波器的选择：
- 高斯滤波器的标准差 σ 可以根据实际噪声水平进行调整。较大的 σ 会导致图像更平滑，但也可能损失细节。

常见问题排查

拟合结果全为1：
- 原因：可能是 airy_function 在 r = 0 处没有正确处理，导致所有点被赋予相同的值。
- 解决方案：确保 airy_function 正确处理 r = 0 的情况，如代码中所示。
优化未收敛：
- 原因：初始参数猜测不合理，或数据噪声过大。
- 解决方案：调整初始参数，或对数据进行更多的滤波和预处理。
图像显示问题：
- 原因：缺少适当的绘图后端或图像数据格式不正确。
- 解决方案：确保安装并正确使用 Plots 包的后端，如 GR 或 Plotly。