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
|
% tvdantzig_example.m
%
% Test out tvdantzig code (TV minimization with bounded residual correlation).
%
% Written by: Justin Romberg, Caltech
% Email: jrom@acm.caltech.edu
% Created: October 2005
%
% use implicit, matrix-free algorithms ?
largescale = 0;
path(path, './Optimization');
path(path, './Measurements');
path(path, './Data');
% test image = 32x32 piece of cameraman's arm
load camera
I = camera(81:112,37:68);
n = 32;
N = n*n;
I = I/norm(I(:));
I = I - mean(I(:));
x = reshape(I,N,1);
% num obs
K = 300;
% permutation P and observation set OMEGA
P = randperm(N)';
q = randperm(N/2-1)+1;
OMEGA = q(1:K/2)';
% measurement matrix
if (largescale)
A = @(z) A_f(z, OMEGA, P);
At = @(z) At_f(z, N, OMEGA, P);
% obsevations
b = A(x);
% initial point
x0 = At(b);
else
FT = 1/sqrt(N)*fft(eye(N));
A = sqrt(2)*[real(FT(OMEGA,:)); imag(FT(OMEGA,:))];
A = [1/sqrt(N)*ones(1,N); A];
At = [];
% observations
b = A*x;
% initial point
x0 = A'*b;
end
epsilon = 5e-3;
tvI = sum(sum(sqrt([diff(I,1,2) zeros(n,1)].^2 + [diff(I,1,1); zeros(1,n)].^2 )));
disp(sprintf('Original TV = %.3f', tvI));
time0 = clock;
xp = tvdantzig_logbarrier(x0, A, At, b, epsilon, 1e-3, 5, 1e-8, 1500);
Ip = reshape(xp, n, n);
disp(sprintf('Total elapsed time = %f secs\n', etime(clock,time0)));
|
