1
function [ im,endensv,priorsumv,variancev,intrinscurv ] = chambollePockTVPrimalDual( im0,imsize,imclean,lambda ,lambda_data,iterations)
2
%GENERALNEWTONALGOTV Summary of this function goes here
3
% Detailed explanation goes here
10
[xgrid,ygrid]=meshgrid(1:imsize,1:imsize);
11
gradf=computeGradientfilter(1,2,1)
12
xgradf=[0,gradf]%[-1 1 0];
13
ygradf=[0,gradf]'%[-1 1 0]';
15
xgradf2=[gradf]%[-1 1];
16
ygradf2=[gradf ]'%[-1 1]';
18
tau=0.1;%1.0/(4.0*lambda);
20
tauspace=0.0;%tau/5.0;
21
px=zeros(imsize);%imfilter(im,[gradf,0],'conv','replicate');
22
py=zeros(imsize);%imfilter(im,[gradf,0]','conv','replicate');
24
%px=px./(max(sqrt(px.^2+py.^2),theta));
25
%py=py./(max(sqrt(px.^2+py.^2),theta));
26
im0s=scale_image(im0,-1,1);
27
ims=scale_image(im0,-1,1);
30
% px = imfilter(im0s, [-1 1], 'replicate');
32
% py = imfilter(im0s, [-1 1]', 'replicate');
33
% reprojection = max(theta, sqrt(px.^2 + py.^2));
34
% px = px./reprojection;
35
% py = py./reprojection;
42
% px=padarray(pout(:,:,1),[5,5],'replicate');
43
% py=padarray(pout(:,:,2),[5,5],'replicate');
44
% [pxx,d]=gradient(px,0.5);
45
% [d,pyy]=gradient(py,0.5);
47
% div_p=div_p(6:(end-5),6:(end-5));
48
dx=imfilter(ims,[gradf,0],'conv','replicate');
49
dy=imfilter(ims,[gradf,0]','conv','replicate');
54
% px=pstep/lambda.*dx+canonxspace0;
56
% py=pstep/lambda.*dy+canonyspace0;
57
% reprojection=max(1.0,sqrt((px).^2+(py).^2));
58
% px=px./reprojection;
59
% py=py./reprojection;
62
px=px+tau*dx*lambda;%./max(1.0,sqrt(dx2+dy2));
63
py=py+tau*dy*lambda;%./max(1.0,sqrt(dx2+dy2));
64
reprojection=max(1.0,sqrt((px).^2+(py).^2));
67
cxx=imfilter(px,xgradf,'conv','replicate');
68
cyy=imfilter(py,ygradf,'conv','replicate');
70
cyx=imfilter(py,xgradf,'conv','replicate');
71
cxy=imfilter(px,ygradf,'conv','replicate');
75
lagim=lambda_data*(ims-im0s)-lambda*div_p0;
77
imsnew=ims-tauprimal*lagim;
82
dx=imfilter(ims,[0.5 0 0.5],'conv','replicate');
83
dy=imfilter(ims,[0.5 0 0.5]','conv','replicate');
84
reprojection=max(1.0,sqrt((dx).^2+(dy).^2));
87
fx=-dy*tauspace*lambda;
88
fy=dx*tauspace*lambda;
102
ims=interp2(ims, cfx, cfy,'cubic');
105
% div_p0 = imfilter(px0, [-1 1 0], 'corr', 0)+ ...
106
% imfilter(py0, [-1 1 0]', 'corr', 0);
108
% ims=im0s+lambda*div_p;
109
% lagpx = imfilter(im0s+lambda*div_p0, [-0.5 0 0.5], 'replicate');
111
% lagpy = imfilter(im0s+lambda*div_p0, [-0.5 0 0.5]', 'replicate');
112
% norm=max(theta,sqrt(lagpx.^2+lagpy.^2));
115
% fx=-lagpy*tauspace*lambda;
116
% fy=lagpx*tauspace*lambda;
129
% px=interp2(px, cfx, cfy,'cubic');
130
% py=interp2(py, cfx, cfy,'cubic');
132
% reprojection = max(theta, sqrt(px.^2 + py.^2));
133
% px = px./reprojection;
134
% py = py./reprojection;
137
curv=(cxx-cyy).^2+(cxy+cyx).^2;
138
intrinscurv=[intrinscurv,sum(sum(curv))];
142
im=scale_image(ims,0,255);
143
imshow([uint8(scale_image(im,0,255)),uint8(im0)]);
144
dx=imfilter(imsnew,xgradf,'conv','replicate');
145
dy=imfilter(imsnew,ygradf,'conv','replicate');
146
prior=lambda.*((dx.*dx+dy.*dy)).^(0.5);
147
dataterm=lambda_data*(im0s-imsnew).^2./2;
148
endens=prior+dataterm;
150
endensv=[endensv,sum(sum(endens))];
152
priorsum=sum(sum(prior));%mean(mean(abs(cxx+cyy)));
153
priorsumv=[priorsumv,priorsum];
155
varderiv=mean(mean((im-imclean).^2));
156
variancev=[variancev,varderiv];