重建代码-叠层-傅里叶得层重建-傅里叶叠层-

%%initializationcloseall;clear;clc;%%parameterslamuda=0.632;%波长umD_led=4*1000;%LED间距h_err=86*1000;%LED与样品间距（校正前）ledMM=15;ledNN=15;%15X15的LED阵列ledM=ledMM;ledN=ledNN;k_lamuda=2*pi/lamuda;%波数pixel_size=6.5;%相机像素尺寸mag=4;%显微镜放大倍率NA=0.1;%物镜数值孔径M=128;N=128;%拍摄到的图像大小D_pixel=pixel_size/mag;%像面像素尺寸kmax=NA*k_lamuda;%物镜数值孔径对应的最大波数半径%重构后图像相比原始图像放大倍率1+2*D_pixel*sin(theta)/lamdaceil向上取整单个维度!!MAGimg=ceil(1+2*D_pixel*7*D_led/sqrt((7*D_led)^2+h_err^2)/lamuda);%重构后图像相比原始图像放大倍率MM=M*MAGimg;NN=N*MAGimg;%重构后图像大小Niter0=10;%初始迭代次数Niter1=50;%225幅图像迭代次数Cuv=8;%模拟退火法的初始搜索步长%%misalignmentobjdx=835.7472;%样品中所选的小区域的位置objdy=835.7472;theta=5;%*(2*rand(1)-1);%LED板的旋转误差leddx=1000;%*(2*rand(1)-1);%LED板的横向位置误差leddy=1000;%*(2*rand(1)-1);h=h_err+1000;%*(2*rand(1)-1);%LED板和样品的真实间距disp('iterthetaled_dxled_dyh分别为');disp(vpa([0,theta,leddx,leddy,h],4));%%objectI0=im2double(imread('I0.bmp'));P0=im2double(imread('P0.bmp'));I0=imresize(I0,[MM,NN]);%物体的真实强度P0=imresize(P0,[MM,NN]);%物体的真实相位o1=sqrt(I0).*exp(sqrt(-1).*P0);%物体的复振幅O1=fftshift(fft2(o1));%fftshift的作用正是让正半轴部分和负半轴部分的图像分别关于各自的中心对称。因为直接用fft得出的数据与频率不是对应的，fftshift可以纠正过来。fft是离散傅立叶变换；fft(X);fft(X,N);fft(X,[],DIM)，fft2是2维离散傅立叶变换；fft2(X);fft2(X,MROWS,NCOLS)%%频域坐标[Fx1,Fy1]=meshgrid(-(N/2):(N/2-1),-(M/2):(M/2-1));Fx1=Fx1./(N*D_pixel).*(2*pi);%原始图像的频域坐标Fy1=Fy1./(M*D_pixel).*(2*pi);%原始图像的频域坐标Fx2=Fx1.*Fx1;Fy2=Fy1.*Fy1;Fxy2=Fx2+Fy2;Pupil0=zeros(M,N);Pupil0(Fxy2<=(kmax^2))=1;%物镜在频域的孔径[Fxx1,Fyy1]=meshgrid(-(NN/2):(NN/2-1),-(MM/2):(MM/2-1));Fxx1=Fxx1(1,:)./(N*D_pixel).*(2*pi);%重构图像的频域坐标Fyy1=Fyy1(:,1)./(M*D_pixel).*(2*pi);%重构图像的频域坐标mesh(Pupil0)%%每个LED在频域对应的像素坐标ledpos_N=2;ledpos_err=zeros(ledMM,ledNN,2);ledpos_true=zeros(ledMM,ledNN,2);ledpos_sum=zeros(2,ledMM*ledNN,ledpos_N);k=zeros(ledMM,ledNN);fori=1:ledMMforj=1:ledNNk(i,j)=j+(i-1)*ledNN;m=i-(ledMM+1)/2;n=j-(ledNN+1)/2;x=objdx-(m*D_led);y=objdy-(n*D_led);u=k_lamuda*x/sqrt(x^2+y^2+h_err^2);v=k_lamuda*y/sqrt(x^2+y^2+h_err^2);Fx1_temp=abs(Fxx1-u.*ones(1,NN));ledpos_err(i,j,1)=find(Fx1_temp==min(Fx1_temp));Fy1_temp=abs(Fyy1-v.*ones(MM,1));ledpos_err(i,j,2)=find(Fy1_temp==min(Fy1_temp));%不知道LED板有误差时的每个LED在频域对应的像素坐标m=[(i-(ledMM+1)/2),(j-(ledNN+1)/2)]*[cos(theta*pi/180);sin(theta*pi/180)];%%%%n=[(i-(ledMM+1)/2),(j-(ledNN+1)/2)]*[-sin(theta*pi/180);cos(theta*pi/180)];%%%%%%x=objdx-(m*D_led+leddx);%%%y=objdy-(n*D_led+leddy);%%%u=k_lamuda*x/sqrt(x^2+y^2+h^2);v=k_lamuda*y/sqrt(x^2+y^2+h^2);Fx1_temp=abs(Fxx1-u.*ones(1,NN));ledpos_true(i,j,1)=find(Fx1_temp==min(Fx1_temp));Fy1_temp=abs(Fyy1-v.*ones(MM,1));ledpos_true(i,j,2)=find(Fy1_temp==min(Fy1_temp));%已知LED板的误差时的每个LED在频域对应的像素坐标ledpos_sum(1,k(i,j),1)=ledpos_err(i,j,1);ledpos_sum(2,k(i,j),1)=ledpos_err(i,j,2);ledpos_sum(1,k(i,j),2)=ledpos_true(i,j,1);ledpos_sum(2,k(i,j),2)=ledpos_true(i,j,2);endend%%imagerecordingnoise=0.4;%%百分制噪声zhendang=1;Isum=zeros(M,N,ledMM*ledNN);fori=1:ledMMforj=1:ledNNuo=ledpos_true(i,j,1);vo=ledpos_true(i,j,2);O1P0=O1((vo-M/2):(vo-1+M/2),(uo-N/2):(uo-1+N/2))./(MAGimg^2).*Pupil0;o10=ifft2(fftshift(O1P0));oI10=abs(o10).^2;st=zhendang.*rand+(1-zhendang/2);ifst<0st=0.0001;endoI10=st*oI10;ifi<=10&&i>=6&&j<=10&&j>=6elses=mean2(oI10)*noise/0.8;oI10=oI10+randn(size(oI10)).*s;oI10(oI10<(s))=0;%噪声淹没endIsum(:,:,k(i,j))=oI10;%225幅拍摄到的图像imshow(Isum(:,:,k(i,j)),[])pause(0.000001)endendI_mid=Isum(:,:,(ledMM*ledNN+1)/2);%中央LED照明时对应的图像%%生成由中心绕圈向外的迭代顺序（更新顺序）ord_ijsum=zeros(ledMM,ledNN);ord_isum=zeros(1,ledM*ledN);ord_jsum=zeros(1,ledM*ledN);ord_ii=(ledMM+1)/2;ord_jj=(ledNN+1)/2;ord_isum(1,1)=ord_ii;ord_jsum(1,1)=ord_jj;ord_ijsum(ord_ii,ord_jj)=1;led_num=1;direction=0;while(min(min(ord_ijsum))==0)led_num=led_num+1;direction2=direction+1;ord_ii2=round(ord_ii+sin(pi/2*direction2));ord_jj2=round(ord_jj+cos(pi/2*direction2));if(ord_ijsum(ord_ii2,ord_jj2)==0)direction=direction2;endord_ii=round(ord_ii+sin(pi/2*direction));ord_jj=round(ord_jj+cos(pi/2*direction));ord_isum(1,led_num)=ord_ii;ord_jsum(1,led_num)=ord_jj;ord_ijsum(ord_ii,ord_jj)=1;end%%不知道LED板的误差时的错误重构结果%oI=imresize(I_mid,MAGimg,'bilinear');%figure;imshow(oI,[]);%oP=ones(MM,NN).*0.5;%o=sqrt(oI).*exp(sqrt(-1).*oP);%O=fftshift(fft2(o));%Pupil=Pupil0;%PupilSUM=O.*0;%%c=zeros(1,ledM*ledN);%Irmse=zeros(1,Niter1);%Prmse=Irmse;%Ico=Irmse;Pco=Ico;%error_now=0;%error_bef=inf;%%alpha=1;beta=1;%foriter=1:Niter1%forled_num=1:1:(ledM*ledN)%%i=ord_isum(led_num);%j=ord_jsum(led_num);%%uo=ledpos_err(i,j,1);%vo=ledpos_err(i,j,2);%%OP_bef=O((vo-M/2):(vo-1+M/2),(uo-N/2):(uo-1+N/2))./(MAGimg^2).*Pupil;%o_bef=ifft2(fftshift(OP_bef));%oI_bef=abs(o_bef).^2;%%oI_cap=Isum(:,:,k(i,j));%o_aft=sqrt(oI_cap)./sqrt(oI_bef).*o_bef;%OP_aft=fftshift(fft2(o_aft));%%OP_diff=(OP_aft-OP_bef).*(MAGimg^2);%guodu=O((vo-M/2):(vo-1+M/2),(uo-N/2):(uo-1+N/2));%%O((vo-M/2):(vo-1+M/2),(uo-N/2):(uo-1+N/2))=guodu+...%alpha*abs(Pupil).*conj(Pupil).*OP_diff./max(max(abs(Pupil)))./(abs(Pupil).^2+1);%Pupil=Pupil+...%beta*(abs(OP_bef).*conj(OP_bef)).*OP_diff./max(max(abs(guodu)))./(abs(OP_bef).^2+1000).*Pupil0;%%Pupil=Pupil+...%%beta*(abs(guodu).*conj(guodu)).*OP_diff./max(max(abs(guodu)))./(abs(guodu).^2+eps).*Pupil0;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(iter==1)%PupilSUM((vo-M/2):(vo-1+M/2),(uo-N/2):(uo-1+N/2))=Pupil0+...%PupilSUM((vo-M/2):(vo-1+M/2),(uo-N/2):(uo-1+N/2)).*(ones(M,N)-Pupil0);%else%O=O.*PupilSUM;%end%%c(k(i,j))=sum(sum(oI_bef))/sum(sum(oI_cap));%%Isum(:,:,k(i,j))=Isum(:,:,k(i,j)).*c(k(i,j));%error_now=error_now+sum(sum((abs(o_bef)-sqrt(oI_cap)).^2));%%end%%if((error_bef-error_now)/error_bef<0.01)%%%Reducethestepsizewhennosufficientprogressismade%alpha=alpha/2;%beta=beta/2;%%StoptheiterationwhenAlphaislessthan0.001(convergenced)%if(alpha<0.001)%alpha=0;%end%%end%%ifiter>Niter0%error_bef=error_now;%end%%%o=ifft2(fftshift(O));%oI=abs(o).^2;%oP=angle(o);%Ico(iter)=corr2(I0,oI);%Pco(iter)=corr2(P0,oP);%disp([iterIco(iter)Pco(iter)]);%oP=oP-mean2(oP)+mean2(P0);%oI=oI-mean2(oI)+mean2(I0);%Idif=I0-oI;%Pdif=P0-oP;%Irmse(iter)=sqrt(mean2(Idif.^2));%Prmse(iter)=sqrt(mean2(Pdif.^2));%disp([iterIrmse(iter)Prmse(iter)]);%%if(alpha==0)%break;end%%end%%pinpu=log(O.*conj(O)+1);%figure;imshow(pinpu,[],'border','tight','initialmagnification','fit');%set(gcf,'Position',[0,0,384,384]);%axisnormal;%figure;imshow(oI,[],'border','tight','initialmagnification','fit');%set(gcf,'Position',[0,0,384,384]);%axisnormal;%figure;imshow(oP,[],'border','tight','initialmagnification','fit');colormap(hot)%set(gcf,'Position',[0,0,384,384]);%axisnormal;%%已知LED板的误差时的正确重构结果oI=imresize(I_mid,MAGimg,'bilinear');figure;imshow(oI,[])title('oI')oP=ones(MM,NN).*0.5;%%%相位初始化o=sqrt(oI).*exp(sqrt(-1).*oP);O=fftshift(fft2(o));Pupil=Pupil0;PupilSUM=O.*0;alpha=1;beta=1;foriter=1:6forled_num=1:1:(ledM*ledN)i=ord_isum(led_num);j=ord_jsum(led_num);uo=ledpos_true(i,j,1);vo=ledpos_true(i,j,2);OP_bef=O((vo-M/2):(vo-1+M/2),(uo-N/2):(uo-1+N/2))./(MAGimg^2).*Pupil;o_bef=ifft2(fftshift(OP_bef));oI_bef=abs(o_bef).^2;oI_cap=Isum(:,:,k(i,j));o_aft=sqrt(oI_cap)./sqrt(oI_bef).*o_bef;OP_aft=fftshift(fft2(o_aft));OP_diff=(OP_aft-OP_bef).*(MAGimg^2);guodu=O((vo-M/2):(vo-1+M/2),(uo-N/2):(uo-1+N/2));O((vo-M/2):(vo-1+M/2),(uo-N/2):(uo-1+N/2))=guodu+...alpha*abs(Pupil).*conj(Pupil).*OP_diff./max(max(abs(Pupil)))./(abs(Pupil).^2+1);Pupil=Pupil+...beta*(abs(OP_bef).*conj(OP_bef)).*OP_diff./max(max(abs(guodu)))./(abs(OP_bef).^2+1000).*Pupil0;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(iter==1)PupilSUM((vo-M/2):(vo-1+M/2),(uo-N/2):(uo-1+N/2))=Pupil0+...PupilSUM((vo-M/2):(vo-1+M/2),(uo-N/2):(uo-1+N/2)).*(ones(M,N)-Pupil0);elseO=O.*PupilSUM;endendendpinpu=log(O.*conj(O)+1);figure;imshow(pinpu,[],'border','tight','initialmagnification','fit');set(gcf,'Position',[0,0,384,384]);axisnormal;o=ifft2(fftshift(O));oI=abs(o).^2;oP=angle(o);figure;imshow(oI,[],'border','tight','initialmagnification','fit');set(gcf,'Position',[0,0,384,384]);axisnormal;figure;imshow(oP,[],'border','tight','initialmagnification','fit');colormap(hot)set(gcf,'Position',[0,0,384,384]);axisnormal;Ico=corr2(oI,I0);%%%求相关系数Pco=corr2(oP,P0);disp([IcoPco]);oP=oP-mean2(oP)+mean2(P0);oI=oI-mean2(oI)+mean2(I0);Idif=I0-oI;Pdif=P0-oP;Irmse=sqrt(mean2(Idif.^2));Prmse=sqrt(mean2(Pdif.^2));disp([IrmsePrmse]);%%SC-FPM重构过程theta_cor=0;%LED校正参数初始化leddx_cor=0;leddy_cor=0;h_cor=h_err;ledpos_cor=ledpos_err;%LED频域坐标初始化Irmse=zeros(1,Niter0+Niter1);Prmse=Irmse;Ico=Prmse;Pco=Ico;c=zeros(1,ledM*ledN);alpha=1;beta=1;Isum0=Isum;%error=zeros(1,Niter0+Niter1);tic;foriter=1:(Niter0+Niter1)error_now=0;if(iter<=(Niter0+1))%每次初始迭代25幅图后的O和P要初始化Pupil=Pupil0;oI=imresize(I_mid,MAGimg,'bilinear');oP=ones(MM,NN).*0.5;o=sqrt(oI).*exp(sqrt(-1).*oP);O=fftshift(fft2(o));%O=O.*PupilSUM;endif(iter<=Niter0)led_seq=1:1:49;%初始迭代25幅图（初始迭代49幅图更稳定一些）elseled_seq=1:1:(ledM*ledN);%之后迭代225幅图endif((rem(iter,2)==1)&&(iter>1)&&(Cuv>2))||((iter>(Niter0+2))&&(Cuv>1))%((iter>(Niter0/2+1))&&(Cuv>2))||((iter>(Niter0+1))&&(Cuv>1))%随着迭代逐步减小模拟退火法的搜索步长Cuv=Cuv/2;endforled_num=led_seqi=ord_isum(led_num);j=ord_jsum(led_num);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%模拟退火法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%uo=ledpos_cor(i,j,1);vo=ledpos_cor(i,j,2);OP_bef=O((vo-M/2):(vo-1+M/2),(uo-N/2):(uo-1+N/2))./(MAGimg^2).*Pupil;o_bef=ifft2(fftshift(OP_bef));oI_bef=abs(o_bef).^2;oI_cap=Isum(:,:,k(i,j));IdiffMIN=mean2((oI_bef-oI_cap).*(oI_bef-oI_cap));%LED原始频域坐标对应的图像之差uoMIN0=uo;voMIN0=vo;uoMIN=uo;voMIN=vo;uoshift=round((rand(1,8)+ones(1,8)).*[-1,-1,-1,0,1,1,1,0].*Cuv);%向周围8个方向随机搜索voshift=round((rand(1,8)+ones(1,8)).*[-1,0,1,1,1,0,-1,-1].*Cuv);%向周围8个方向随机搜索forNshift=1:size(uoshift,2)uo=uoMIN0+uoshift(Nshift);vo=voMIN0+voshift(Nshift);OP_bef=O((vo-M/2):(vo-1+M/2),(uo-N/2):(uo-1+N/2))./(MAGimg^2).*Pupil;o_bef=ifft2(fftshift(OP_bef));oI_bef=abs(o_bef).^2;Idiff=mean2((oI_bef-oI_cap).*(oI_bef-oI_cap));if(Idiff<IdiffMIN)%找到图像之差最小时对应的LED频域坐标平移量uoMIN=uo;voMIN=vo;IdiffMIN=Idiff;endendledpos_cor(i,j,1)=uoMIN;%更新LED频域坐标（25个或225个）ledpos_cor(i,j,2)=voMIN;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%模拟退火法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%更新O和P%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%uo=ledpos_cor(i,j,1);vo=ledpos_cor(i,j,2);OP_bef=O((vo-M/2):(vo-1+M/2),(uo-N/2):(uo-1+N/2))./(MAGimg^2).*Pupil;o_bef=ifft2(fftshift(OP_bef));oI_bef=abs(o_bef).^2;oI_cap=Isum(:,:,k(i,j));if((mean2(oI_cap)>0.1)&&(mean2(oI_bef)>0.1))||((mean2(oI_cap)<0.1)&&(mean2(oI_bef)<0.1))o_aft=sqrt(oI_cap)./sqrt(oI_bef).*o_bef;elseo_aft=o_bef;endOP_aft=fftshift(fft2(o_aft));OP_diff=(OP_aft-OP_bef).*(MAGimg^2);guodu=O((vo-M/2):(vo-1+M/2),(uo-N/2):(uo-1+N/2));O((vo-M/2):(vo-1+M/2),(uo-N/2):(uo-1+N/2))=guodu+...alpha*abs(Pupil).*conj(Pupil).*OP_diff./max(max(abs(Pupil)))./(abs(Pupil).^2+1);if(iter>Niter0)Pupil=Pupil+...beta*(abs(OP_bef).*conj(OP_bef)).*OP_diff./max(max(abs(guodu)))./(abs(OP_bef).^2+1000).*Pupil0;endifiter>=Niter0+2error_now=error_now+sum(sum((abs(o_bef)-sqrt(oI_cap)).^2));endifrem(iter,2)==1&&iter<=Niter0+1ifi==8&&j==8c(k(i,j))=sum(sum(oI_bef))/sum(sum(oI_cap));elsec(k(i,j))=sum(sum(oI_bef))/sum(sum(oI_cap))/c(k(8,8));Isum(:,:,k(i,j))=Isum(:,:,k(i,j)).*c(k(i,j));endelseifiter<=Niter0Isum=Isum0;endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%更新O和P%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%非线性回归%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Rm=zeros(1,(ledM*ledN));Rn=zeros(1,(ledM*ledN));Ru=zeros(1,(ledM*ledN));Rv=zeros(1,(ledM*ledN));forled_num=led_seqi=ord_isum(led_num);j=ord_jsum(led_num);Rm(1,led_num)=i-(ledMM+1)/2;Rn(1,led_num)=j-(ledNN+1)/2;Ru(1,led_num)=Fyy1(ledpos_cor(i,j,1));Rv(1,led_num)=Fxx1(ledpos_cor(i,j,2));endR=[Rm',Rn',Ru',Rv'];Q0=[theta_cor,leddx_cor,leddy_cor,h_cor];F=@(Q,R)((k_lamuda.*(objdx-(D_led*cos(pi/180*Q(1)).*R(:,1)+D_led*sin(pi/180*Q(1)).*R(:,2)+Q(2)))./...sqrt((objdx-(D_led*cos(pi/180*Q(1)).*R(:,1)+D_led*sin(pi/180*Q(1)).*R(:,2)+Q(2))).^2+...(objdy-(-D_led*sin(pi/180*Q(1)).*R(:,1)+D_led*cos(pi/180*Q(1)).*R(:,2)+Q(3))).^2+Q(4).^2)-R(:,3)).^2+...(k_lamuda.*(objdy-(-D_led*sin(pi/180*Q(1)).*R(:,1)+D_led*cos(pi/180*Q(1)).*R(:,2)+Q(3)))./...sqrt((objdx-(D_led*cos(pi/180*Q(1)).*R(:,1)+D_led*sin(pi/180*Q(1)).*R(:,2)+Q(2))).^2+...(objdy-(-D_led*sin(pi/180*Q(1)).*R(:,1)+D_led*cos(pi/180*Q(1)).*R(:,2)+Q(3))).^2+Q(4).^2)-R(:,4)).^2);Q=nlinfit(R,zeros(size(R,1),1),F,Q0);%非线性回归分析disp(vpa([iter,Q],4));theta_cor=Q(1);leddx_cor=Q(2);leddy_cor=Q(3);h_cor=Q(4);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%非线性回归%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%led在频域的表示ledpos_N=ledpos_N+1;ledpos_sum(1,:,ledpos_N)=ledpos_sum(1,:,ledpos_N-1);ledpos_sum(2,:,ledpos_N)=ledpos_sum(2,:,ledpos_N-1);fori=1:ledMMforj=1:ledNNm=[(i-(ledMM+1)/2),(j-(ledNN+1)/2)]*[cos(theta_cor*pi/180);sin(theta_cor*pi/180)];n=[(i-(ledMM+1)/2),(j-(ledNN+1)/2)]*[-sin(theta_cor*pi/180);cos(theta_cor*pi/180)];x=objdx-(m*D_led+leddx_cor);y=objdy-(n*D_led+leddy_cor);u=k_lamuda*x/sqrt(x^2+y^2+h_cor^2);v=k_lamuda*y/sqrt(x^2+y^2+h_cor^2);Fx1_temp=abs(Fxx1-u.*ones(1,NN));ledpos_cor(i,j,1)=find(Fx1_temp==min(Fx1_temp));Fy1_temp=abs(Fyy1-v.*ones(MM,1));ledpos_cor(i,j,2)=find(Fy1_temp==min(Fy1_temp));%更新LED频域坐标（225个）ledpos_sum(1,k(i,j),ledpos_N-1)=ledpos_cor(i,j,1);ledpos_sum(2,k(i,j),ledpos_N-1)=ledpos_cor(i,j,2);endend%%%%%自适应步进ifiter>=Niter0+2&&((error_bef-error_now)/error_bef<0.01)%Reducethestepsizewhennosufficientprogressismadealpha=alpha/2;beta=beta/2;%StoptheiterationwhenAlphaislessthan0.001(convergenced)if(alpha<0.001)alpha=0;endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%