close all;clear all;clc; %
iterative=300; %設迭代次數為300次吧
imagename='cameraman.tif'; %你想要提取相位的圖像名稱
phaseimage='phase.png'; %要保存的相位圖像名稱
figure(1),imshow(imagename);
imagename='cameraman.tif'; %你想要提取相位的圖像名稱
phaseimage='phase.png'; %要保存的相位圖像名稱
figure(1),imshow(imagename);
%空域輸入圖像的幅度（是已知的，也就是清晰的圖像，它的灰度就是幅值）和相位圖像（待恢復）
known_abs_spatial=imread(imagename); %作為輸入圖像的幅度，是已知的
%known_abs_spatial =rgb2gray(known_abs_spatial);%注意要用單通道圖像做實驗，如果你讀取的是彩色圖像，那就吧這行取消注釋變成灰度圖像吧
known_abs_spatial=im2double(known_abs_spatial); %將圖像灰度映射到0～1
known_abs_spatial=imread(imagename); %作為輸入圖像的幅度，是已知的
%known_abs_spatial =rgb2gray(known_abs_spatial);%注意要用單通道圖像做實驗，如果你讀取的是彩色圖像，那就吧這行取消注釋變成灰度圖像吧
known_abs_spatial=im2double(known_abs_spatial); %將圖像灰度映射到0～1
unknown_phase=known_abs_spatial; %Peppers圖像作為輸入圖像的相位，也即為待恢復的數據，
%要求它和known_abs_spatial大小一致，所以這里直接賦值就好了
unknown_phase=im2double(unknown_phase); %將圖像灰度映射到0～1
unknown_phase2=unknown_phase*2*pi; %相位范圍映射到0-2*pi
unknown_phase2(unknown_phase2>pi)=unknown_phase2(unknown_phase2>pi)-2*pi;%進一步映射至[-pi,+pi]
%要求它和known_abs_spatial大小一致，所以這里直接賦值就好了
unknown_phase=im2double(unknown_phase); %將圖像灰度映射到0～1
unknown_phase2=unknown_phase*2*pi; %相位范圍映射到0-2*pi
unknown_phase2(unknown_phase2>pi)=unknown_phase2(unknown_phase2>pi)-2*pi;%進一步映射至[-pi,+pi]
[width,length]=size(known_abs_spatial); %獲取圖像的大小
input=known_abs_spatial.*exp(i*unknown_phase2); %最終輸入圖像:幅度*e^(i*相位角度)，它是復數圖像
known_abs_fourier=abs(fft2(input)); %先將input圖像進行傅立葉變換，然后取模，就是傅氏變換后的幅度
%以下開始迭代求相位
phase_estimate=pi*rand(width,length); %這是生成了一副大小為(width*length)的圖像
%它的像素值是在[0,pi]范圍內隨機生成的。
figure(2),imshow(phase_estimate)
%以下開始迭代
for p=1:iterative
signal_estimate_spatial=known_abs_spatial.*exp(i*phase_estimate); %Step 1 構造estimated signal：還是幅度*e^(i*相位角度)變成復數形式
temp1=fft2(signal_estimate_spatial); %傅立葉變換到頻域
temp_ang=angle(temp1); %求相位弧度，它的范圍是[-pi,pi]
signal_estimate_fourier=known_abs_fourier.*exp(i*temp_ang); %Step 2 替換傅氏變換后的幅度，產生estimate Fourier transform
temp2=ifft2(signal_estimate_fourier); %Step 3 對Step 2產生的estimate Fourier transform進行傅立葉反變換，又變換到空域了
phase_estimate=angle(temp2); %Step 4:estimated phase
% IS=abs(abs(temp2)-abs(temp1)).^2;
% MSE=sum(IS(:))/256^2%計算均方誤差
end
%以上循環就是通過隨便預設一個相位圖像，在循環中不斷調整逼近真實的相位，直到滿足條件（也就是我們求的相位和真實相位非常接近的時候）
%不過這里我們只需要設定一個比較大的循環就可以了，基本上都可以滿足條件了，這個激光原理就講過了。
input=known_abs_spatial.*exp(i*unknown_phase2); %最終輸入圖像:幅度*e^(i*相位角度)，它是復數圖像
known_abs_fourier=abs(fft2(input)); %先將input圖像進行傅立葉變換，然后取模，就是傅氏變換后的幅度
%以下開始迭代求相位
phase_estimate=pi*rand(width,length); %這是生成了一副大小為(width*length)的圖像
%它的像素值是在[0,pi]范圍內隨機生成的。
figure(2),imshow(phase_estimate)
%以下開始迭代
for p=1:iterative
signal_estimate_spatial=known_abs_spatial.*exp(i*phase_estimate); %Step 1 構造estimated signal：還是幅度*e^(i*相位角度)變成復數形式
temp1=fft2(signal_estimate_spatial); %傅立葉變換到頻域
temp_ang=angle(temp1); %求相位弧度，它的范圍是[-pi,pi]
signal_estimate_fourier=known_abs_fourier.*exp(i*temp_ang); %Step 2 替換傅氏變換后的幅度，產生estimate Fourier transform
temp2=ifft2(signal_estimate_fourier); %Step 3 對Step 2產生的estimate Fourier transform進行傅立葉反變換，又變換到空域了
phase_estimate=angle(temp2); %Step 4:estimated phase
% IS=abs(abs(temp2)-abs(temp1)).^2;
% MSE=sum(IS(:))/256^2%計算均方誤差
end
%以上循環就是通過隨便預設一個相位圖像，在循環中不斷調整逼近真實的相位，直到滿足條件（也就是我們求的相位和真實相位非常接近的時候）
%不過這里我們只需要設定一個比較大的循環就可以了，基本上都可以滿足條件了，這個激光原理就講過了。
phase_estimate(phase_estimate<0)=phase_estimate(phase_estimate<0)+2*pi; %把estimate_phase從[-pi,+pi]，映射到[0,2pi]
retrieved=phase_estimate/(2*pi);%再映射到[0,1]
retrieved=phase_estimate/(2*pi);%再映射到[0,1]
% IS=abs(abs(temp2)-abs(temp1)).^2;
% MSE1=sum(IS(:))/256^2%計算均方誤差
% figure(2);plot(log10(MSE),'LineWidth',1.5);
% xlabel('Iterative number');%迭代的次數
% ylabel('Logarithm of Mean Square Error');%均方誤差的對數
figure (3)
imshow(retrieved);title('相位圖像')%顯示我們提取到的相位圖像
imwrite(retrieved,phaseimage)
%<span style="color:rgb(79,79,79);font-size:16px;font-family:'PingFang SC', 'Microsoft YaHei', SimHei, Arial, SimSun;">
%</span>