~gerald-mwangi/+junk/ThesisMatlab

« back to all changes in this revision

Viewing changes to Matlab/StructureTensFeaturesAnal/chambollePockTVPrimalDual.m

  • Committer: gerald.mwangi at gmx
  • Date: 2017-10-05 12:14:48 UTC
  • Revision ID: gerald.mwangi@gmx.de-20171005121448-wesl3srbemcec5eg
Added code for extension of chambolle pock

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
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
 
4
 
 
5
im=im0;
 
6
endensv=[];
 
7
priorsumv=[];
 
8
variancev=[];
 
9
intrinscurv=[];
 
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]';
 
14
 
 
15
xgradf2=[gradf]%[-1 1];
 
16
ygradf2=[gradf ]'%[-1 1]';
 
17
theta=1.0;
 
18
tau=0.1;%1.0/(4.0*lambda);
 
19
tauprimal=0.1;
 
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');
 
23
 
 
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);
 
28
%   im0s=im0;
 
29
%   ims=im0;
 
30
%         px = imfilter(im0s, [-1 1], 'replicate');
 
31
%          
 
32
%         py = imfilter(im0s, [-1 1]', 'replicate');
 
33
%                     reprojection = max(theta, sqrt(px.^2 + py.^2));
 
34
%         px = px./reprojection;
 
35
%         py = py./reprojection;
 
36
for i=1:iterations
 
37
   
 
38
 
 
39
 
 
40
 
 
41
 
 
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);
 
46
%         div_p=pxx+pyy;
 
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');
 
50
            dx2=dx.*dx;
 
51
            dy2=dy.*dy;
 
52
 
 
53
 
 
54
            % px=pstep/lambda.*dx+canonxspace0;
 
55
            % 
 
56
            % py=pstep/lambda.*dy+canonyspace0;
 
57
            % reprojection=max(1.0,sqrt((px).^2+(py).^2));
 
58
            % px=px./reprojection;
 
59
            % py=py./reprojection;
 
60
            % canonx=px;
 
61
            % canony=py;
 
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));
 
65
            px=px./reprojection;
 
66
            py=py./reprojection;
 
67
                cxx=imfilter(px,xgradf,'conv','replicate');
 
68
            cyy=imfilter(py,ygradf,'conv','replicate');
 
69
 
 
70
            cyx=imfilter(py,xgradf,'conv','replicate');
 
71
            cxy=imfilter(px,ygradf,'conv','replicate');
 
72
            
 
73
        
 
74
                    div_p0 = cxx+cyy;
 
75
            lagim=lambda_data*(ims-im0s)-lambda*div_p0;
 
76
            
 
77
                 imsnew=ims-tauprimal*lagim;
 
78
        
 
79
        ims=imsnew;           
 
80
            
 
81
            
 
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));
 
85
            dx=dx./reprojection;
 
86
            dy=dy./reprojection;
 
87
                fx=-dy*tauspace*lambda;
 
88
    fy=dx*tauspace*lambda;
 
89
    
 
90
                cfx=fx+xgrid;
 
91
        cfy=fy+ygrid;
 
92
        id=cfx<1;
 
93
        cfx(id)=1;
 
94
        id=cfx>imsize;
 
95
        cfx(id)=imsize;
 
96
        id=cfy<1;
 
97
        cfy(id)=1;
 
98
        id=cfy>imsize;
 
99
        cfy(id)=imsize;
 
100
        
 
101
 
 
102
        ims=interp2(ims, cfx, cfy,'cubic');
 
103
 
 
104
 
 
105
%             div_p0 = imfilter(px0, [-1 1 0], 'corr', 0)+ ...
 
106
%                 imfilter(py0, [-1 1 0]', 'corr', 0);  
 
107
%     
 
108
%     ims=im0s+lambda*div_p;
 
109
%             lagpx = imfilter(im0s+lambda*div_p0, [-0.5 0 0.5], 'replicate');
 
110
%          
 
111
%         lagpy = imfilter(im0s+lambda*div_p0, [-0.5 0 0.5]', 'replicate');
 
112
%         norm=max(theta,sqrt(lagpx.^2+lagpy.^2));
 
113
%         lagpx=lagpx./norm;
 
114
%         lagpy=lagpy./norm;
 
115
%     fx=-lagpy*tauspace*lambda;
 
116
%     fy=lagpx*tauspace*lambda;
 
117
%     
 
118
%                 cfx=fx+xgrid;
 
119
%         cfy=fy+ygrid;
 
120
%         id=cfx<1;
 
121
%         cfx(id)=1;
 
122
%         id=cfx>imsize;
 
123
%         cfx(id)=imsize;
 
124
%         id=cfy<1;
 
125
%         cfy(id)=1;
 
126
%         id=cfy>imsize;
 
127
%         cfy(id)=imsize;
 
128
%         
 
129
%         px=interp2(px, cfx, cfy,'cubic');
 
130
%         py=interp2(py, cfx, cfy,'cubic');
 
131
    
 
132
%         reprojection = max(theta, sqrt(px.^2 + py.^2));
 
133
%         px = px./reprojection;
 
134
%         py = py./reprojection;
 
135
    
 
136
    
 
137
     curv=(cxx-cyy).^2+(cxy+cyx).^2;
 
138
     intrinscurv=[intrinscurv,sum(sum(curv))];
 
139
 
 
140
     figure(3);
 
141
 
 
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;
 
149
 
 
150
    endensv=[endensv,sum(sum(endens))];
 
151
 
 
152
    priorsum=sum(sum(prior));%mean(mean(abs(cxx+cyy)));
 
153
    priorsumv=[priorsumv,priorsum];
 
154
 
 
155
    varderiv=mean(mean((im-imclean).^2));
 
156
    variancev=[variancev,varderiv];
 
157
    drawnow
 
158
end
 
159
 
 
160
 
 
161
end
 
162