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
|
% tvqc_quantization_example.m
%
% Takes random measurements of the 'boats' image, quantizes, and
% reconstructs using quadractically constrained TV.
%
% Written by: Justin Romberg, Caltech
% Email: jrom@acm.caltech.edu
% Created: October 2005
%
path(path, './Optimization');
path(path, './Measurements');
path(path, './Data');
% boats
n = 256;
N = n*n;
load boats
img = boats;
I = img/norm(img(:));
I = I - mean(I(:));
x = reshape(I,N,1);
% num obs
K = 25000;
% permutation P and observation set OMEGA
P = randperm(N)';
q = randperm(N/2-1)+1;
OMEGA = q(1:K/2)';
% measurement implicit matrices
A = @(z) A_f(z, OMEGA, P);
At = @(z) At_f(z, N, OMEGA, P);
% take measurements
Ax = A(x);
% quantization codebook
mxax = max(abs(Ax));
bins = 10;
C = 2*mxax*(linspace(1/(2*bins),1-1/(2*bins),bins)-1/2);
% quantize
b = zeros(K,1);
symbols = zeros(K,1);
for kk = 1:K
[mn,mni] = min(abs(Ax(kk)-C));
b(kk) = C(mni);
symbols(kk) = mni;
end
% error (just for future reference)
e = b - Ax;
% the error should be roughly this size
epsilon = (2*mxax/bins)/sqrt(12)*sqrt(K);
% initial point = min l2
x0 = At(b);
tvI = sum(sum(sqrt([diff(I,1,2) zeros(n,1)].^2 + [diff(I,1,1); zeros(1,n)].^2 )));
disp(sprintf('original TV = %8.3f', tvI));
time0 = clock;
xp = tvqc_logbarrier(x0, A, At, b, epsilon, 1e-1, 5);
Ip = reshape(xp, n, n);
disp(sprintf('Total elapsed time = %f secs\n', etime(clock,time0)));
figure(1); imagesc(I); colormap(gray); colorbar;
title('Original image');
figure(2); imagesc(Ip); colormap(gray); colorbar;
title('Reconstructed image');
|
