-
Notifications
You must be signed in to change notification settings - Fork 1
/
adom.m
173 lines (152 loc) · 4.29 KB
/
adom.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
% This is the official code of ADOM for stripe noise removal in remote sensing images (RSI).
% ADOM: ADMM-Based Optimization Model for Stripe Noise Removal in Remote Sensing Image
% IEEE Access
% 09/25/2023
% Namwon Kim ([email protected])
function [s,bag]=adom(o,opts)
Dx=defDDxt; %x-direction finite function
Dy=defDDyt; %y-direction finite function
Dt=defDDt;
[m,n]=size(o);
s=zeros(m,n);
p1=zeros(m,n);
p2=p1;
p3=p1;
groupweight=ones(m,n);
r=ones(m,n);
lambda1=opts.lambda1;
lambda2=opts.lambda2;
beta1 = opts.beta1;
beta2 = opts.beta2;
beta3 = opts.beta3;
tol=opts.tol;
maxiter=opts.maxitr;
gamma = 0.5*(1+sqrt(5));
cir_fft = beta1*(abs(psf2otf([1; -1], [m,n]))).^2 + ...
beta3*(abs(psf2otf([1, -1], [m,n]))).^2 + beta2;
Dxs=Dx(o-s); % d_x(s)
Dys=Dy(s); % d_y(s)
Dx_o=Dx(o);
a=zeros(m,n);
b = a;
c = a;
bag = zeros(1);
a1 = 1;
a2 = 1;
cal_tolerance=1;
iter=1;
sp = s;
while cal_tolerance>tol && iter<maxiter
% update weights control
weight1 = (a2+cal_tolerance+eps)/(a1-cal_tolerance+eps);
% starting point control
if iter <= 10
a1 = a2;
a2 = (1+sqrt(1+(4*((a1)^2))))/2;
% limit coefficient
weightsz=weight1;
else
a1 = a2;
a2 = (1+sqrt(1+(2*((a1)^2))))/2;
if a1 ~= a2
weightsz=(a1)/(a2+eps);
end
end
% step size control
s = s + ((a1-0.1)/(a2))*(s-sp);
p1 = weightsz*p1;
p2 = weightsz*p2;
p3 = weightsz*p3;
V1=Dys+p1/beta1;
V2=Dxs+p2/beta2;
% a-subproblem L1 norm
a=sign(V1).*max(abs(V1)-1/beta1,0);
% b-subproblem weighted L1 norm
b=sign(V2).*max(abs(V2)-(weight1.*lambda1)/beta2,0);
% c-subproblem
% weighted L2,1 norm with contextual information
% weighted group sparse
for i=1:n
% update weights control
% contextual information
mean_s = sum(abs(s(:,i)))/n;
if i==n
idx=abs(s(:,end-1)-s(:,end))<mean_s;
elseif i==1
idx=abs(s(:,i+1)-s(:,i))<mean_s;
else
idx=abs(s(:,i+1)-s(:,i))<mean_s | abs(s(:,i-1)-s(:,i))<mean_s;
end
groupweight(idx,i) = 1/sqrt((norm(s(:,i)+cal_tolerance,'fro')*2+eps));
r(:,i)=s(:,i)+p3(:,i)/beta3;
c(:,i)=r(:,i).*max(norm(r(:,i),'fro')-(groupweight(:,i).*lambda2)/beta3,0)/(norm(r(:,i),'fro')+eps);
end
% s-subproblem
sp=s;
t1 = beta1*a-p1;
t2 = beta2*Dx_o-beta2*b+p2;
t3 = beta3*c-p3;
sr = Dt(t2,t1)+t3;
s = real(ifft2(fft2(sr)./(cir_fft+eps)));
u1=o-s;
u2=o-sp;
% update Dxs, Dys
Dxs=Dx(o-s);
Dys=Dy(s);
% update p
% The golden ratio
% 1.618 ¡Ö (sqrt(5)+1)/2
p1=p1+gamma*beta1*(Dys-a);
p2=p2+gamma*beta2*(Dxs-b);
p3=p3+gamma*beta3*(s-c);
% Formal definition (relative error)
% eta = epsilon/abs(v)=||(v-v_approx)/v||_2
cal_tolerance = norm(u1-u2,'fro')/norm(u2,'fro');
bag(iter) = cal_tolerance;
if cal_tolerance< tol
break;
end
iter = iter+1;
end % loop
end % main
%%%%--------------Subfunction-------------%%%%
function Dx=defDDxt
Dx=@(U)ForwardD1(U);
end
function Dy=defDDyt
Dy=@(U)ForwardD2(U);
end
function Dt=defDDt
Dt= @(X,Y) Dive(X,Y);
end
% ^^^^
% |||||...||
% |||||...||
% .
% .
% .
% |||||...||
% Input: discrete data
function Dux=ForwardD1(U)
Dux=[diff(U,1,2),U(:,1)-U(:,end)];
end
% -----...--
% <
% -----...--
% .
% .
% .
% -----...--
% -----...--
% -----...--
% .
% .
% .
% -----...--
function Duy=ForwardD2(U)
Duy=[diff(U,1,1);U(1,:)-U(end,:)];
end
function DtXY = Dive(X,Y)
DtXY = [X(:,end) - X(:, 1), -diff(X,1,2)]; % x-directional
DtXY = DtXY + [Y(end,:) - Y(1, :); -diff(Y,1,1)]; % x-directional + y-directional
end