Advertisement
Guest User

Untitled

a guest
Dec 3rd, 2019
108
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 3.23 KB | None | 0 0
  1. %%lab 2
  2. %% 1
  3. m=100;
  4. n=100;
  5. sztitan = [215,126];
  6. b=sqrt(m*n);
  7. sz = [m,n];
  8. q1 = @(k)[0 -1 0;-1 4+k^2 -1;0 -1 0];
  9. q2 = @(k)[0 -0.1 -10;-0.1 20.4+k^2 -0.1;-10 -0.1 0];
  10. q3 = @(k)[0 -1 -2;-1 8+k^2 -1;-2 -1 0];
  11. k = 0.01;
  12. %imagesc(gmrfprec(sz,q2(k)));
  13. Q = gmrfprec(sz,q3(k)); %q3 var bra med 2 eller 3 neighbour niv?er f?r data i labben
  14. %% 1.a)
  15. s = Q \ sparse(m*n/2+m/2,1,1,size(Q,1),1); %covariance for colon, when inverting Q
  16. imagesc(reshape(s,b,[]));
  17. colorbar
  18. %% 1.b)
  19. mu = 5;
  20. R = chol(Q);
  21. x = mu + R \ randn(size(R,1),1);
  22. imagesc(reshape(x,b,[]))
  23. colorbar
  24. %% 1.c) higher order neigbourhoods
  25. Q=Q^2;
  26. s = Q \ sparse(m*n/2+m/2,1,1,size(Q,1),1);
  27. imagesc(reshape(s,b,[]));
  28. %% 2 Interpolation
  29. load lab3.mat
  30. %% Image to restore
  31. imagesc(reshape(xmiss,100,[]))
  32. %% Construct the observation matrix based on the location of the known pixels:
  33. log_risk = log(Insurance(:,2)./Insurance(:,1));
  34. Y=logrisk(known);
  35.  
  36. A = sparse(1:length(Y), find(known), 1,length(Y), numel(xmiss));
  37. spy(A)
  38.  
  39. %% Given the observation matrix, a vector of covariates for all grid points,
  40. %and precision matrices for x and ?; the joint matrices can be constructed
  41. %as
  42.  
  43. Aall = [A ones(length(Y),1)];
  44. Qbeta = 1e-6 * speye(1); %use SPEYE to create a SPARSE-matrix
  45. Qall = blkdiag(Q, Qbeta);
  46.  
  47. %% c) Compute the posterior precision matrix
  48. %assume a very small observation-uncertainty
  49. Qeps = 1e5 * speye(length(Y));
  50. %posterior precisionx
  51. Qxy = Qall + Aall'*Qeps*Aall;
  52. spy(Qxy);
  53. %% reordering
  54. p = amd(Qxy);
  55. Qxy = Qxy(p,p);
  56. Aall = Aall(:,p);
  57. spy(Qxy);
  58. %% Compute conditional expectation ... and do inverse reordering
  59.  
  60. % E_xy = inv(Qxy)\Aall'*Qeps*Y;
  61. R=chol(Qxy);
  62. E_xy = R\(R'\(Aall'*Qeps*Y));
  63.  
  64. E_xy(p) = E_xy;
  65. E_zy = [speye(size(Q,1)) ones(size(Q,1),1)]*E_xy;
  66. subplot(2,2,1);
  67. imagesc(reshape(xmiss,100,[]));
  68. colorbar
  69. subplot(2,2,2);
  70. imagesc(reshape(E_zy,100,[]));
  71. colorbar
  72. subplot(2,2,3);
  73. imagesc(reshape(xtrue,100,[]));
  74. colorbar
  75. subplot(2,2,4);
  76. x_rec = reshape(E_zy,100,[]);
  77. imagesc(x_rec-xtrue);
  78. colorbar
  79. %%
  80. titan = imread('titan.jpg');
  81. imshow(titan);
  82. %% reconstruction of image from titan
  83. k = 0.001;
  84. Q = gmrfprec(sztitan,q1(k));
  85. Q=Q*Q;
  86. miss=0.9;
  87. titan = imread('titan.jpg');
  88. titan=double(titan)/255;
  89. xtrue = titan;
  90. known=(rand(size(titan))>miss);
  91. dknown=double(known);
  92. xmiss=xtrue.*dknown;
  93. Y=xmiss(known);
  94. A = sparse(1:length(Y), find(known), 1,length(Y), numel(xmiss));
  95. Aall = [A ones(length(Y),1)];
  96. Qbeta = 1e-6 * speye(1);
  97. Qall = blkdiag(Q, Qbeta);
  98. %assume a very small observation-uncertainty
  99. Qeps = 1e5 * speye(length(Y));
  100. %posterior precision
  101. Qxy = Qall + Aall'*Qeps*Aall;
  102. p = amd(Qxy);
  103. Qxy = Qxy(p,p);
  104. Aall = Aall(:,p);
  105. R=chol(Qxy);
  106. E_xy = R\(R'\(Aall'*Qeps*Y));
  107. E_xy(p) = E_xy;
  108. E_zy = [speye(size(Q,1)) ones(size(Q,1),1)]*E_xy;
  109. subplot(2,3,1);
  110. % imagesc(reshape(xmiss,215,[]));
  111. imshow(reshape(xmiss,215,[]));
  112. colorbar
  113. subplot(2,3,2);
  114. % imagesc(reshape(E_zy,215,[]));
  115. imshow(reshape(E_zy,215,[]));
  116. colorbar
  117. subplot(2,3,3);
  118. % imagesc(reshape(xtrue,215,[]));
  119. imshow(reshape(xtrue,215,[]));
  120. colorbar
  121. subplot(2,3,4);
  122. x_rec = reshape(E_zy,215,[]);
  123. plot(x_rec-xtrue,'xk');
  124. subplot(2,3,5)
  125. s = Q \ sparse(m*n/2+m/2,1,1,size(Q,1),1); %covariance for colon, when inverting Q
  126. imagesc(reshape(s,215,[]));
  127. colorbar
  128.  
  129.  
  130.  
  131. %%
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement