79 lines
1.4 KiB
Mathematica
79 lines
1.4 KiB
Mathematica
|
close all;
|
||
|
|
clear; clc;
|
||
|
|
addpath /Users/maggicl/Git/matlab2tikz/src/;
|
||
|
|
|
||
|
|
%% Load Default Img Data
|
||
|
|
load('blur_data/B.mat');
|
||
|
|
B=double(B);
|
||
|
|
load('blur_data/A.mat');
|
||
|
|
A=double(A);
|
||
|
|
|
||
|
|
% Show Image
|
||
|
|
figure
|
||
|
|
im_l=min(min(B));
|
||
|
|
im_u=max(max(B));
|
||
|
|
imshow(B,[im_l,im_u])
|
||
|
|
title('Blured Image')
|
||
|
|
|
||
|
|
% Vectorize the image (row by row)
|
||
|
|
b=B';
|
||
|
|
b=b(:);
|
||
|
|
|
||
|
|
AT = A' * A;
|
||
|
|
bt = A' * b;
|
||
|
|
|
||
|
|
IL = ichol(AT, struct('type', 'nofill', 'diagcomp', 0.01));
|
||
|
|
|
||
|
|
[x, rvec] = myCG(AT, bt, zeros(size(b)), 200, 1e-6);
|
||
|
|
[x2, ~, ~, ~, rvec2] = pcg(AT, bt, 1e-6, 200, IL, IL');
|
||
|
|
|
||
|
|
figure;
|
||
|
|
semilogy(rvec / norm(bt));
|
||
|
|
hold on;
|
||
|
|
semilogy(rvec2 / norm(bt));
|
||
|
|
hold off;
|
||
|
|
title('Residual norms over iteration (y is log)')
|
||
|
|
matlab2tikz('showInfo', false, '../res_log.tex');
|
||
|
|
|
||
|
|
|
||
|
|
X = zeros(250, 250);
|
||
|
|
for i = 0:249
|
||
|
|
for j = 1:250
|
||
|
|
X(i + 1, j) = x(i * 250 + j);
|
||
|
|
end
|
||
|
|
end
|
||
|
|
|
||
|
|
figure;
|
||
|
|
im_l=min(min(X));
|
||
|
|
im_u=max(max(X));
|
||
|
|
imshow(X,[im_l,im_u])
|
||
|
|
title('Sharp Image (myCG)')
|
||
|
|
matlab2tikz('showInfo', false, '../img_my.tex');
|
||
|
|
|
||
|
|
|
||
|
|
X2 = zeros(250, 250);
|
||
|
|
for i = 0:249
|
||
|
|
for j = 1:250
|
||
|
|
X2(i + 1, j) = x2(i * 250 + j);
|
||
|
|
end
|
||
|
|
end
|
||
|
|
|
||
|
|
figure;
|
||
|
|
im_l=min(min(X2));
|
||
|
|
im_u=max(max(X2));
|
||
|
|
imshow(X2,[im_l,im_u])
|
||
|
|
title('Sharp Image (rcg)')
|
||
|
|
matlab2tikz('showInfo', false, '../img_rcg.tex');
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
%% Validate Test values
|
||
|
|
load('test_data/A_test.mat');
|
||
|
|
load('test_data/x_test_exact.mat');
|
||
|
|
load('test_data/b_test.mat');
|
||
|
|
|
||
|
|
%res=||x^*-A^{-1}b||
|
||
|
|
res=x_test_exact-inv(A_test)*b_test;
|
||
|
|
norm(res);
|
||
|
|
%(Now do it with your CG and Matlab's PCG routine!!!)
|