%给图像加扰动，扰动方式以zernike多项式的形式加入，Image为源图像，Co为要加的zernike系数，为行向量，
%输出分别为扰动后的图像，扰动PV值，扰动RMS值
function [OUT,dis_pv,dis_rms]=Image_disturb(Image,Co)
Image=imread('1.bmp');
M=size(Image,1);N=size(Image,2);
%--------所加扰动，用zernike多项式表示--------------------
Co=[1,1,1,0.02,0.02,0.02,0.02,0.02,0.02];
N_Zer=size(Co,2);%所用zernike阶数
disturb=zeros(M,N);
Zer=zernike_poly(M,N);%1-36阶
for q=1:N_Zer 
    disturb=disturb+Co(1,q)*Zer{q};
end
dis_pv=max(max(disturb/2/pi))-min(min(disturb/2/pi));
n_fi=nnz(disturb);%nnz(A)返回A中非0元素的个数 
dis_rms=sqrt(sum(sum((disturb/2/pi).^2))/n_fi);
%--------------计算点扩散函数---------------------
pupil_r=20;%光瞳半径，单位mm
W=Hn(disturb,M,N,pupil_r);
h=fftshift(fft2(W));
p=(abs(h)).^2;
surf(abs(W));
figure
surf(p)
%-------------计算扰动后的图像--------------------
Image=double(Image)/double(max(max(Image)));
I=fftshift(fft2(Image));
P=fftshift(fft2(p));
OUT=ifftshift(ifft2(P.*I));
OUT=abs(OUT);
OUT=OUT/max(max(OUT));

figure
subplot(1,2,1)
imshow(Image);
title('原始图像')
subplot(1,2,2)
imshow(OUT);
title('加扰动后的图像')




%计算广义光瞳函数
%fi是波前畸变值
%N是参与计算矩阵大小
function Hn_pupil=Hn(fi,M,N,pupil_r)
H=exp(i*fi);
P=pupil(pupil_r,M,N);
Hn_pupil=H.*P;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%画出圆光瞳函数，输入半径，可以得到圆瞳大小
function P=pupil(R,M,N)
ix=-R:2*R/(M-1):R;
iy=R:-2*R/(N-1):-R;
[x,y]=meshgrid(ix,iy);
P=zeros(M,N);
i=find(sqrt(x.^2+y.^2)<=R);
P(i)=1;