✅作者简介：热爱科研的Matlab仿真开发者，擅长毕业设计辅导、数学建模、数据处理、程序设计科研仿真。
🍎完整代码获取 定制创新 论文复现点击：Matlab科研工作室
👇 关注我领取海量matlab电子书和数学建模资料
🍊个人信条：做科研，博学之、审问之、慎思之、明辨之、笃行之，是为：博学慎思，明辨笃行。
🔥 内容介绍
一、引言
计算成像作为现代成像技术的重要领域，旨在通过算法和计算手段从测量数据中重建高质量图像。相位恢复问题是计算成像中的关键挑战之一，而基于优化的方法为解决这一问题提供了有效途径。交替方向乘子法（ADMM）因其在处理复杂优化问题时的高效性和灵活性，在相位恢复以及相关的压缩感知成像任务中得到了广泛应用。
二、相位恢复问题概述
- 相位恢复的重要性
：在许多成像系统中，探测器通常只能测量光的强度信息，而相位信息对于准确重建物体的真实结构至关重要。例如，在 X 射线晶体学、光学显微镜和天文成像等领域，恢复相位信息能够显著提高图像的分辨率和质量，帮助科学家获取更详细的样本信息。
- 相位恢复的挑战
：从仅有的强度测量数据中恢复相位是一个不适定问题，因为多个不同的相位分布可能产生相同的强度模式。这意味着存在无数个可能的解，需要通过合理的约束和优化算法来找到最符合实际情况的相位。
- ADMM 在相位恢复中的应用
- ADMM (L1)
：L1 范数常用于促进解的稀疏性。在相位恢复中，通过将图像的某些特征（如像素值的梯度）的 L1 范数作为正则化项加入目标函数，ADMM 能够在恢复相位的同时，使重建图像具有一定的稀疏特性。例如，图像中的边缘部分在 L1 范数的约束下，能够更清晰地凸显出来，有助于去除噪声和提高图像的清晰度。
- ADMM (MWPR_L1)
：MWPR（Multi - Window Phase Retrieval）结合 L1 范数的方法，通过在多个窗口内进行相位恢复，并利用 L1 范数进行约束。这种方法能够充分利用图像在不同局部区域的特征，提高相位恢复的准确性。不同窗口的设置可以捕捉到图像不同尺度的信息，从而在复杂结构的图像重建中表现出色。
- ADMM (TV)
：总变分（TV）正则化是一种常用的图像恢复技术。ADMM 结合 TV 正则化，通过最小化图像的总变分来平滑图像，同时保持图像的边缘信息。在相位恢复中，TV 正则化可以有效地抑制噪声，使得重建图像更加平滑自然，尤其适用于处理具有块状或纹理特征的图像。
- ADMM (L1)
四、压缩感知与相关算法
- 压缩感知基础
：压缩感知理论指出，对于满足一定稀疏条件的信号，可以从远少于传统奈奎斯特采样定理所需的测量数据中精确重建。在计算成像中，这意味着可以通过较少的测量获取图像信息，大大提高了数据采集效率。
- 压缩感知在成像中的算法
- Denoise
：在压缩感知成像过程中，噪声会影响图像的重建质量。去噪算法旨在从含噪的测量数据中去除噪声，提高重建图像的质量。常见的去噪方法包括基于小波变换、稀疏表示等，通过对噪声特性的分析和信号的稀疏建模，实现噪声的有效抑制。
- Twist
：Twist 算法是一种用于压缩感知成像的特殊方法，它通过对测量数据进行特定的变换和处理，利用信号的稀疏性和结构信息，实现图像的快速准确重建。该算法在处理具有特定结构的图像数据时，能够显著提高重建效率和精度。
- Fista
：快速迭代收缩阈值算法（Fista）是一种加速的梯度下降算法，常用于求解具有稀疏约束的优化问题。在压缩感知成像中，Fista 通过快速迭代更新解，能够在保证重建精度的同时，提高算法的收敛速度。它在处理大规模数据时表现出较高的效率，是一种常用的图像重建算法。
- Fista - WeightedTV
：Fista 结合加权总变分（WeightedTV）的方法，通过为不同区域的图像元素赋予不同的权重，进一步优化了总变分正则化。这种方法能够根据图像的局部特征，更加灵活地平衡平滑和边缘保持，从而在复杂图像重建中获得更好的效果。
- Fista - TV_Laplace
：Fista 与 TV 正则化结合 Laplace 先验的方法，利用 Laplace 分布对图像的稀疏性进行建模。Laplace 先验能够更好地捕捉图像中的边缘和细节信息，与 TV 正则化和 Fista 算法相结合，在重建图像时能够在抑制噪声的同时，保留更多的图像细节，提高图像的视觉质量。
- Denoise
五、ADMM 与压缩感知算法的协同
在计算成像中，ADMM 与压缩感知相关算法常常协同工作。例如，ADMM 可以用于解决相位恢复中的优化问题，而压缩感知算法则用于从少量测量数据中重建图像。通过将两者结合，可以在低采样率下实现高质量的图像重建。ADMM 的优化能力能够更好地处理相位恢复中的约束条件，而压缩感知算法的稀疏重建特性则能够充分利用数据的稀疏性，提高重建效率和准确性。
六、总结与展望
基于 ADMM 的优化方法在计算成像的相位恢复以及相关的压缩感知任务中展现出强大的能力。通过合理选择正则化项（如 L1、TV 等），ADMM 能够有效地处理复杂的优化问题，实现高质量的图像重建。同时，与压缩感知算法的协同应用进一步拓展了其在低采样率成像等领域的应用潜力。未来，随着计算成像技术的不断发展，ADMM 和相关算法有望在提高重建速度、增强对复杂场景的适应性以及结合多模态数据等方面取得更多突破，推动计算成像技术在医学、遥感、材料科学等众多领域的广泛应用。
⛳️ 运行结果
📣 部分代码
function output = CG(o,u0,op,z1,dz,lambda,deltaA,deltaB,gamma,mu,a0,ITER,alpha)
dgamma = 0.1;
dmu = 5e-3;
da0 = 1;
dITER = 500;
beta = 0.1;
q = [1,1];
dalpha = 1;
method = 0;
av = true;
i = sqrt(-1);RMSE = []; tm = 0;
if ~isreal(op) || ~isnumeric(op) || ceil(op)-floor(op)==1
disp('Error: Wrong setting of the operation, please choose op={0,1,2,3}.'); out = []; return
end
[N,M,K] = size(o); % the size of the input set of the (noisy) intensity observations
if mod(N,2)~=0 || mod(M,2)~=0
disp('Error: The size of the intensity observations must be even integers'); out = []; return
end
% if K<2
% disp('Error: It must be at least 2 intensity observations, otherwise phase retrieval doesnt make sense'); out = []; return
% end
if (nargin<9 || isempty(gamma)), gamma=dgamma; end
if (nargin<14 || isempty(alpha)), alpha=dalpha; end
if (nargin<11 || isempty(a0)), a0=da0; end
if any(~isreal(gamma)) || any(~isnumeric(gamma)) || any(gamma<=0)
disp('Error: "gamma" must be a vector of K positive real values'); out = []; return
end
if (nargin<10 || isempty(mu)), mu=dmu; end
if ~isreal(mu) || ~isnumeric(mu) || mu<0
disp('Error: "mu" must be a non-negative real value'); out = []; return
end
if (nargin<12 || isempty(ITER)), ITER=dITER; end
if q(1) < 1 || q(2) < 1, disp('Error: q must be not smaller than 1'); out = []; return ; end
Nya = ceil(q(1)*N/2)*2; Nxa = ceil(q(2)*M/2)*2; % the support of the computed Fourier images
cRMSE = false; % control parameter which control whether RMSE should be calculated
% it the true u0 is given - RMSE is computed and the output variable 'out' is the calculated object
% otherwise - RMSE=[] is empty variable and 'out' is the calculated volumetric sensor plane wave field
if ~isempty(u0) % the input object is used to estimate the reconstruction accuracy in RMSE
[Ny0,Nx0,~] = size(u0); % the size of the input object
if Ny0>Nya || Nx0>Nxa
disp('Error: The size of the object "u0" should not be larger than the size of the intensity observations!');out = []; return;
end
if mod(Ny0,2)~=0 || mod(Nx0,2)~=0
disp('Error: The size of the object "u0" must be even integers'); out = []; return
end
cRMSE = true; RMSE = zeros(2,ITER);
end
% IF mod(op,2)=0, i.e. op IS EVEN (op={2,4,6,8...}) THEN WE USE F-DDT, OTHERWISE - ASD
if method < 2
Nya = Nya+Ny0;Nxa = Nxa+Nx0;
end
coordy = Nya/2-N/2+1:Nya/2+N/2; coordx = Nxa/2-M/2+1:Nxa/2+M/2; % center part of the sensor plane wave field
bourdery = Nya/2-Ny0/2+1:Nya/2+Ny0/2; bourderx = Nxa/2-Nx0/2+1:Nxa/2+Nx0/2; % center part of the object
gamma0=0.9;
gamma1=0.05;
gamma2=0.05;
gamma3=0.01;
gammar = gamma0/K*ones(K,1);
gamma_scaler = 2;
[V,B,strs,tm]=TransformMatrix(0,z1,dz,lambda,deltaA,deltaB,Nya,Nxa,K,gammar,mu,av);
u0_iter = zeros(Nya, Nxa); a_=1/2;
u0_iter=a_+u0_iter;
for jndex = 1:ITER
PP = zeros(Nya, Nxa);
for index = 1:K
PP1 = ifft2(fft2(u0_iter).*V(:,:,index));
PP1(bourdery,bourderx) = PP1(bourdery,bourderx) - o(:,:,index);
PP = PP + ifft2(fft2(PP1).*B(:,:,index));
end
PP1 = gamma0*PP + abs(gamma1*TV(TV(u0_iter,1),1)+gamma2*TV(TV(u0_iter,2),2));
u0_iter = u0_iter-gamma3*PP1;
out = u0_iter(bourdery,bourderx);
if cRMSE
figure(5)
v3 = abs(out) - abs(u0); RMSE(1,jndex)=sqrt(mean(mean( ( v3 ).^2 )));
x1 = angle(out) - angle(u0); txmp = mean(x1(:)); x1 = x1-txmp; RMSE(2,jndex)=sqrt(mean(mean( ( x1 ).^2 )));
subplot(2,3,1), imshow(abs(out),[]), title('amplitude, ADMM'),xlabel(['RMSE=' num2str(RMSE(1,jndex))])
subplot(2,3,2), plot(RMSE(1,1:jndex)), xlabel(['\gamma=' num2str(mean(gamma)) ',\alpha_r=' num2str(alpha) ',\mu=' num2str(mu)])
drawnow
end
gamma1 = gamma1;
gamma2 = gamma2;
gamma3 = gamma3;
end
output=u0_iter;
