2 tong2711 tong2711 于 2015.06.09 11:13 提问

请问matlab里自带的deconvblind函数运用的是哪种盲解卷积的算法,求原理过程~~~

% Parse inputs to verify valid function calling syntaxes and arguments
[J,P,NUMIT,DAMPAR,READOUT,WEIGHT,sizeI,classI,sizePSF,FunFcn,FunArg] = ...
parse_inputs(varargin{:});

% 1. Prepare parameters for iterations
%
% Create indexes for image according to the sampling rate
idx = repmat({':'},[1 length(sizeI)]);

wI = max(WEIGHT.*(READOUT + J{1}),0);% at this point - positivity constraint

fw = fftn(WEIGHT);
clear WEIGHT;
DAMPAR22 = (DAMPAR.^2)/2;

% 2. L_R Iterations
%
lambda = 2*any(J{4}(:)~=0);
for k = (lambda + 1) : (lambda + NUMIT),

% 2.a Make an image and PSF predictions for the next iteration

if k > 2,% image
lambda = (J{4}(:,1).'*J{4}(:,2))/(J{4}(:,2).'*J{4}(:,2) + eps);
lambda = max(min(lambda,1),0);
% stability enforcement lambda(0,1)之间
end
Y = max(J{2} + lambda*(J{2} - J{3}),0);% image positivity constraint

if k > 2,% PSF
lambda = (P{4}(:,1).'*P{4}(:,2))/(P{4}(:,2).'*P{4}(:,2) + eps);
lambda = max(min(lambda,1),0);% stability enforcement

end
B = max(P{2} + lambda*(P{2} - P{3}),0);% PSF positivity constraint
sumPSF = sum(B(:));
B = B/(sum(B(:)) + (sumPSF==0)*eps);% normalization is a necessary constraint,
% because given only input image, the algorithm cannot even know how much
% power is in the image vs PSF. Therefore, we force PSF to satisfy this
% type of normalization: sum to one.

% 2.b Make core for the LR estimation
CC = corelucy(Y,psf2otf(B,sizeI),DAMPAR22,wI,READOUT,1,idx,[],[]);

% 2.c Determine next iteration image & apply positivity constraint
J{3} = J{2};
H = psf2otf(P{2},sizeI);
scale = real(ifftn(conj(H).*fw)) + sqrt(eps);
J{2} = max(Y.*real(ifftn(conj(H).*CC))./scale,0);
clear scale;
J{4} = [J{2}(:)-Y(:) J{4}(:,1)];
clear Y H;

% 2.d Determine next iteration PSF & apply positivity constraint + normalization
P{3} = P{2};
H = fftn(J{3});
scale = otf2psf(conj(H).*fw,sizePSF) + sqrt(eps);
P{2} = max(B.*otf2psf(conj(H).*CC,sizePSF)./scale,0);
clear CC H;

sumPSF = sum(P{2}(:));
P{2} = P{2}/(sumPSF + (sumPSF==0)*eps);

if ~isempty(FunFcn),
FunArg{1} = P{2};
P{2} = feval(FunFcn,FunArg{:});
end;
P{4} = [P{2}(:)-B(:) P{4}(:,1)];
end
clear fw wI;

% 3. Convert the right array (for cell it is first array, for notcell it is
% second array) to the original image class & output the whole
num = 1 + strcmp(classI{1},'notcell');
if ~strcmp(classI{2},'double'),
J{num} = images.internal.changeClass(classI{2},J{num});
end

if num == 2,% the input & output is NOT a cell
P = P{2};
J = J{2};
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Function: parse_inputs
function [J,P,NUMIT,DAMPAR,READOUT,WEIGHT,sizeI,classI,sizePSF,FunFcn,FunArg] ...
= parse_inputs(varargin)
%
% Outputs:
% I=J{1} the input array (could be any numeric class, 2D, 3D)
% P=P{1} the operator that distorts the ideal image
%
% Defaults:
%
NUMIT=[];NUMIT_d = 10; % Number of iterations, usually produces good
% result by 10.
DAMPAR =[];DAMPAR_d = 0;% No damping is default
WEIGHT =[]; % All pixels are of equal quality, flat-field is one
READOUT=[];READOUT_d= 0;% Zero readout noise or any other
% back/fore/ground noise associated with CCD camera.
% Or the Image is corrected already for this noise by user.
FunFcn = '';FunFcn_d = '';
FunArg = {};FunArg_d = {};
funnum = [];funnum_d = nargin+1;

narginchk(2,inf);% no constraint on max number
% because of FUN args

% First, assign the inputs starting with the cell/not cell image & PSF
switch iscell(varargin{1}) + iscell(varargin{2}),

case 0, % no-cell array is used to do a single set of iterations
classI{1} = 'notcell';

J{1} = varargin{1};% create a cell array in order to do the iterations
P{1} = varargin{2};
case 1,
error(message('images:deconvblind:IandInitpsfMustBeOfSameType'))
case 2,% input cell is used to resume the interrupted iterations or
classI{1} = 'cell';% to interrupt the iteration to resume them later
J = varargin{1};
P = varargin{2};
if length(J)~=length(P),
error(message('images:deconvblind:IandInitpsfMustBeOfSameSize'))
end
end;

% check the Image, which is the first array of the J-cell
[sizeI, sizePSF] = padlength(size(J{1}), size(P{1}));
classI{2} = class(J{1});

validateattributes(J{1},{'uint8' 'uint16' 'double' 'int16' 'single'},...
{'real' 'nonempty' 'finite'},mfilename,'I',1);

if prod(sizeI)<2,
error(message('images:deconvblind:inputImageMustHaveAtLeast2Elements'))
elseif ~isa(J{1},'double'),
J{1} = im2double(J{1});
end

% check the PSF, which is the first array of the P-cell
validateattributes(P{1},{'uint8' 'uint16' 'double' 'int16' 'single'},...
{'real' 'nonempty' 'finite' 'nonzero'},mfilename,'INITPSF',2);

if prod(sizePSF)<2,
error(message('images:deconvblind:initPSFMustHaveAtLeast2Elements'))
elseif ~isa(P{1},'double'),
P{1} = double(P{1});
end

% now since the image&PSF are OK&double, we assign the rest of the J & P cells
len = length(J);
if len == 1,% J = {I} will be reassigned to J = {I,I,0,0}
J{2} = J{1};
J{3} = 0;
J{4}(prod(sizeI),2) = 0;
P{2} = P{1};
P{3} = 0;
P{4}(prod(sizePSF),2) = 0;
elseif len ~= 4,% J = {I,J,Jm1,gk} has to have 4 or 1 arrays
error(message('images:deconvblind:inputCellsMustBe1or4ElementNumArrays'))
else % check if J,Jm1,gk are double in the input cell
if ~all([isa(J{2},'double'),isa(J{3},'double'),isa(J{4},'double')]),
error(message('images:deconvblind:ImageCellElementsMustBeDouble'))
elseif ~all([isa(P{2},'double'),isa(P{3},'double'),isa(P{4},'double')]),
error(message('images:deconvblind:psfCellElementsMustBeDouble'))
end
end;

% Second, Find out if we have a function to put additional constraints on the PSF
%

function_classes = {'inline','function_handle','char'};
idx = [];
for n = 3:nargin,
idx = strmatch(class(varargin{n}),function_classes);
if ~isempty(idx),
[FunFcn,msgStruct] = fcnchk(varargin{n}); %only works on char, making it inline
if ~isempty(msgStruct)
error(msgStruct)
end
FunArg = [{0},varargin(n+1:nargin)];
try % how this function works, just in case.
feval(FunFcn,FunArg{:});
catch ME
Ftype = {'inline object','function_handle','expression ==>'};
Ffcnstr = {' ',' ',varargin{n}};
error(message('images:deconvblind:userSuppliedFcnFailed', Ftype{ idx }, Ffcnstr{ idx }, ME.message))
end
funnum = n;
break
end
end

if isempty(idx),
FunFcn = FunFcn_d;
FunArg = FunArg_d;
funnum = funnum_d;
end

%
% Third, Assign the inputs for general deconvolution:
%
if funnum>7
error(message('images:validate:tooManyInputs',mfilename));
end

switch funnum,
case 4,% deconvblind(I,PSF,NUMIT,fun,...)
NUMIT = varargin{3};
case 5,% deconvblind(I,PSF,NUMIT,DAMPAR,fun,...)
NUMIT = varargin{3};
DAMPAR = varargin{4};
case 6,% deconvblind(I,PSF,NUMIT,DAMPAR,WEIGHT,fun,...)
NUMIT = varargin{3};
DAMPAR = varargin{4};
WEIGHT = varargin{5};
case 7,% deconvblind(I,PSF,NUMIT,DAMPAR,WEIGHT,READOUT,fun,...)
NUMIT = varargin{3};
DAMPAR = varargin{4};
WEIGHT = varargin{5};
READOUT = varargin{6};
end

% Forth, Check validity of the gen.conv. input parameters:
%
% NUMIT check number of iterations
if isempty(NUMIT),
NUMIT = NUMIT_d;
else %verify validity
validateattributes(NUMIT,{'double'},...
{'scalar' 'positive' 'integer' 'finite'},...
mfilename,'NUMIT',3);
end

% DAMPAR check damping parameter
if isempty(DAMPAR),
DAMPAR = DAMPAR_d;
elseif (numel(DAMPAR)~=1) && ~isequal(size(DAMPAR),sizeI),
error(message('images:deconvblind:damparMustBeSizeOfInputImage'));
elseif ~isa(DAMPAR,classI{2}),
error(message('images:deconvblind:damparMustBeSameClassAsInputImage'));
elseif ~strcmp(classI{2},'double'),
DAMPAR = im2double(DAMPAR);
end

if ~isfinite(DAMPAR),
error(message('images:deconvblind:damparMustBeFinite'));
end

% WEIGHT check weighting
if isempty(WEIGHT),
WEIGHT = ones(sizeI);
else
numw = numel(WEIGHT);
validateattributes(WEIGHT,{'double'},{'finite'},mfilename,'WEIGHT',5);
if (numw ~= 1) && ~isequal(size(WEIGHT),sizeI),
error(message('images:deconvblind:weightMustBeSizeOfInputImage'));
elseif numw == 1,
WEIGHT = repmat(WEIGHT,sizeI);
end;
end

% READOUT check read-out noise
if isempty(READOUT),
READOUT = READOUT_d;
elseif (numel(READOUT)~=1) && ~isequal(size(READOUT),sizeI),
error(message('images:deconvblind:readoutMustBeSizeOfInputImage'));
elseif ~isa(READOUT,classI{2}),
error(message('images:deconvblind:readoutMustBeSameClassAsInputImage'));
elseif ~strcmp(classI{2},'double'),
READOUT = im2double(READOUT);
end
if ~isfinite(READOUT),
error(message('images:deconvblind:readoutMustBeFinite'));
end;

2个回答

fdscwm12
fdscwm12   2015.06.09 11:19

我去,这个代码有点长的啊。

tong2711
tong2711   2015.06.10 09:05

这个代码就是matlab里自带的函数程序

Csdn user default icon
上传中...
上传图片
插入图片
准确详细的回答,更有利于被提问者采纳,从而获得C币。复制、灌水、广告等回答会被删除,是时候展现真正的技术了!
其他相关推荐
MATLAB中盲去卷积deconvblind函数的使用
deconvblind:使用盲解卷积的去模糊图像。[J,PSF] = deconvblind(I,INITPSF)使用最大似然算法对图像I解卷积,返回去模糊图像J和恢复的点扩散函数PSF。 生成的PSF是与INITPSF相同大小的正数组,归一化,所以它的总和增加到1。PSF的恢复受其初始猜测大小INITPSF的影响较大,而其值较小(一个数组是一个更安全的猜测)。I 可以是N维数组。为了改善恢复,可...
盲去卷积原理及在图像复原的应用
盲去卷积的原理及在图像复原的应用
matlab自带函数-盲卷积-加噪-卷积-滤波-小结
总结自网上、matlab帮助文档等,都是图像复原中经常用到的基础函数或操作。可以模拟图像降质过程和用一些经典方法盲解卷积复原的过程。 一、卷积:conv2、convn、convmtx2 卷积的计算步骤: (1)卷积核绕自己的核心元素顺时针旋转180度 (2)移动卷积核的中心元素,使它位于输入图像待处理像素的正上方 (3)在旋转后的卷积核中,将输入图像的像素值作为权重相乘 (4)第三
盲卷积要点
Blind Deconvolution Primer From home photographers to world-class astronomers, anyone who has taken a photograph, has taken a blurry picture, including NASA. At these times, it would
盲反卷积算法复原图像matlab实现
盲反卷积算法复原图像matlab实现,有着详细的注释和解释,可以让人很好读懂。
第一步,我想我总算找到用Matlab计算"卷积"以及,进行"可逆"的"解卷积"的途径了
 首先, 这会使得Fourier transform self-deconvolution 在振动光谱中的练习成为可能; 其次, 我要尝试, 这次是否把Fourier deconvolution用于temperature modulated differential scanning calorimetry曲线的解析中是否能够得到至少是近似的可用的rapid reverse和rapid
基于MATLAB并结合IBD算法的盲迭代反卷积法进行图像复原
基于MATLAB并结合IBD算法的盲迭代反卷积法进行图像复原图像复原 盲迭代 反卷积 IBD PSF估计
matlab盲解卷积去噪
针对于模糊图像复原,psf复原,加权数组去模糊,FUN复原图像
反卷积算法
Deconvolution All microscopes are limited by the laws of physics, and these laws state that when light passes through a medium, that light will bend. This is one of the most common causes of haze a
盲卷积复原算法
基于matlab的盲卷积算法复原,其复原效果明显,能够很好的抑制振铃效应