Lucy-Richardson加速算法以及PSF计算MATLAB代码

前面我们介绍了Lucy-Richardson的加速算法，这里给出其implement和PSF的显微成像计算公式。

https://blog.csdn.net/weixin_41923961/article/details/81157082

function result=Accelaration_DeconvLucy(data,psf,iteration,method)
%-----------------------------------------------
%Inputs:
% data      single diffraction limit image
% psf        PSF
%iteration Maximum deconvolution iterations {example:20}
%method accelarate method{'NoAcce'(normal),'Acce1','Acce2','Acce3'}
%------------------------------------------------
%Output:
% result     Lucy-Richardson deconvolution result
%------------------------------------------------
% reference:
% [1].D. S. C. Biggs and M. Andrews, "Acceleration of iterative image restoration
% algorithms," Appl. Opt. 36(8), 1766â€?1775 (1997).%   Copyright  2018 Weisong Zhao et al
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
psf= psf./sum(sum(psf));
psf0=psf;
B=floor(max(size(data,1),size(data,2))/4);
data=padarray(data, [B,B] , 'replicate');
% move the PSF so that the maximum is at pixel 1 in the corner
psfm=psf;
psf = zeros(size(data,1),size(data,2));
psf(1 : size(psfm,1), 1 : size(psfm,2)) = psfm;
[~,idx] = max(psf(:));
[x,y,~] = ind2sub(size(psf), idx);
psf = circshift(psf, -[x,y,1] + 1);
otf = fftn(psf);
fprintf('Start deconvolution\n')
% Richardson Lucy update rule
rliter = @(estimate, data, otf)fftn(data./ max(ifftn(otf .* fftn(estimate)), 1e-6));
switch methodcase'NoAcce'estimate=data;estimate=max(estimate,0.001);for iter = 1:iterationif mod(iter,10)==0fprintf('iteration %d \n',iter)endestimate= estimate.*ifftn(conj(otf).*rliter(estimate,data,otf));estimate=max(estimate,1e-5);%             estimate=estimate./max(max(estimate));endresult=estimate(B:end-B, B:end-B)./max(max(estimate(B:end-B, B:end-B)));case 'Acce1'estimate=data;estimate=max(estimate,0.001);k=2;for iter = 1:iterationif mod(iter,10)==0fprintf('iteration %d \n',iter)endestimate= estimate.*(ifftn(conj(otf).*rliter(estimate,data,otf))).^k;%         estimate=estimate./max(max(estimate));estimate=max(estimate,1e-5);endresult=estimate(B:end-B, B:end-B)./max(max(estimate(B:end-B, B:end-B)));case 'Acce2'estimate=data;estimate=max(estimate,0.001);delta=2;for iter = 1:iterationif mod(iter,10)==0fprintf('iteration %d \n',iter)endestimatek=estimate;estimate= estimate.*ifftn(conj(otf).*rliter(estimate,data,otf));between_k_and_k1=estimate-estimatek;estimate=estimate+delta*between_k_and_k1;estimate=max(estimate,0.00001);estimate=estimate./max(max(estimate));endresult=estimate(B:end-B, B:end-B)./max(max(estimate(B:end-B, B:end-B)));case 'Acce3'yk=data;xk=zeros(size(data));vk=zeros(size(data));for iter = 1:iterationxk_update=xk;xk= max(yk.*real(ifftn(conj(otf).*rliter(yk,data,otf)))./real(ifftn(fftn(ones(size(data))).*otf)),0.00001);vk_update=vk;vk=max(xk-yk,0.00001);if iter==1alpha=0;elsealpha=sum(sum(vk_update.*vk))/(sum(sum(vk_update.*vk_update))+eps);alpha=max(min(alpha,1),0);endyk=max(xk+alpha*(xk-xk_update),0.00001);endresult=yk(B+1:end-B, B+1:end-B)./max(max(yk(B+1:end-B, B+1:end-B)));
end

function Ipsf=Generate_PSF(pixel,lamda,NA,n,z)
%-----------------------------------------------
%Source code for generating PSF
%pixel     pixel size {example:65*10^-9}
%lamda   wavelength {example:532*10^-9}
%NA        NA {example:1.49}
%n          number of psf {example:64}
%z          defocus length {example:1*10^-6}
%------------------------------------------------
%Output:
% Ipsf    PSF%-------------------------------------------------------------------------------------
%                    Copyright  2018 Weisong Zhao
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%-------------------------------------------------------------------------------------
%%
if nargin < 4 || isempty(n)n=32;
end
if nargin < 5 || isempty(z)z=0;
end
sin2=((1-(1-NA^2))/2);
u=8*pi*z*sin2/lamda;
h=@(r,p) 2*exp((1i*u*(p.^2))/2).*besselj(0,2*pi*r*NA/lamda.*p);
x=-n*pixel:pixel:n*pixel;
[X,Y]=meshgrid(x,x);
[~,s1]=cart2pol(X,Y);
idx=s1<=1;
IP=zeros(size(X));
k=1;
for f=1:1:size(s1)for j=1:1:size(s1)if idx(f,j)==0IP(f,j)=0;elseo=s1(idx);r=o(k);k=k+1;II=@(p)h(r,p);IP(f,j)=integral(II,0,1);endend
end
Ipsf=abs(IP.^2);
Ipsf=Ipsf./sum(sum( Ipsf));

Lucy-Richardson加速算法以及PSF计算MATLAB代码相关推荐

灰狼（GWO）算法（附完整Matlab代码，可直接复制）
尊重他人劳动成果,请勿转载! 有问题可留言或私信,看到了都会回复解答! 其他算法请参考: 1.粒子群(PSO)优化算法(附完整Matlab代码,可直接复制)https://blog.csdn.net/ ...
粒子群(PSO)算法（附完整Matlab代码，可直接复制）
在粒子群优化算法中,每个解可用一只鸟(粒子)表示,目标函数就是鸟群所需要寻找的食物源.寻找最优解的过程中,粒子包含两种行为:个体行为和群体行为. 个体行为:粒子根据自身在寻优过程中的最优解更新自己的位 ...
25-混合A星算法Hybrid_Astar路径规划MATLAB代码
资源: Hybrid-Astar(混合A星算法)路径规划MATLAB代码-电子商务文档类资源-CSDN文库主要内容: 以车辆的运动学模型为节点,以当前点到终点的Astar距离和RS距离两者最大的距离 ...
matlab圆周率计算,matlab代码求圆周率的简单算法
说起圆周率的算法很多人都会想起一大堆的无穷级数等各种表达式,但是这样的算法需要比较高的数学推理水品,而且对于很多初学者而言很难理解.程序员不能仅仅是机械的写程序,必须要真正的理解程序中每个代码的意义, ...
图像处理：边缘提取算法(边缘提取算子总结)——Matlab代码实现
边缘提取算子一阶: Roberts算子.Sobel算子.Prewitt算子.Kirsch算子.Robinson算子二阶: Laplacian算子.Canny算子.Marr-Hildreth(Lo ...
a*算法matlab代码_NSGAII多目标优化算法讲解（附MATLAB代码）
小编今天为大家讲解NSGA-II多目标优化算法,提到多目标优化,大家可能第一个就想到NSGA-II算法,今天小编就带领大家解开NSGA-II的神秘面纱. NSGA-II全称是快速非支配排序遗传算法,这 ...
Dijkstra算法和Floyd算法详解（MATLAB代码）
一.Dijkstra算法 1.算法简介 Dijkstra算法是由E.W.Dijkstra于1959年提出,又叫迪杰斯特拉算法,它应用了贪心算法模式,是目前公认的最好的求解最短路径的方法.算法解决的是有 ...
【路径规划-TSP问题】基于粒子群结合蚁群算法求解旅行商问题附matlab代码
1 内容介绍一种基于粒子群优化的蚁群算法求解TSP问题的方法.该方法在求解TSP问题时,利用粒子群优化的思想,对蚁群算法的参数取值进行优化并选择.在粒子群算法中,将蚁群算法的5个参数(q,α,β,ρ ...
多尺度结构元素形态学边缘检测算法的研究-含Matlab代码
目录一.引言二.数学形态学理论概述三.实验验证四.参考文献五.Matlab代码获取一.引言使用数字图像处理技术来解决计算机视觉.人工智能.生物遥感器视觉等领域所涉及到的图像问题时,最重要 ...

Lucy-Richardson加速算法以及PSF计算MATLAB代码

Lucy-Richardson加速算法以及PSF计算MATLAB代码相关推荐

最新文章

热门文章