-
Notifications
You must be signed in to change notification settings - Fork 5
/
angularthresh.m
115 lines (106 loc) · 3.62 KB
/
angularthresh.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
function [vwtinv,vw,vwt,recer,kilme]=...
angularthresh(v,tperc,J,iface,wav,precon,xver)
% [vwtinv,vw,vwt,recer,kilme]=...
% ANGULARTHRESH(v,tperc,J,iface,wav,precon,xver)
%
% Hard thresholding of a model in cubed-sphere coordinates
%
% INPUT:
%
% v The NxNx6 input model
% tperc Truncation level as percentile [defaulted]
% J Maximum scale in both directions
% iface Numbered face if you don't want to do them all [0]
% wav 'D2', 'D4' or 'D6' identifying the basis
% precon [1|0 1|0] toggles preconditioning in either dimension
% xver 1 Provides optional truth and fact checking
% 0 Doesn't [default]
%
% OUTPUT:
%
% vwtinv The NxNx6 output model after thresholded reconstruction
% vw The original wavelet and scaling coefficients
% vwt The thresholded wavelet and scaling coefficients
% recer An array with the following diagnostics:
% [1] The reconstruction error l2-NORM, in percent
% [2] The number of zeroed wavelet coefficients, in percent
% [3] The actual truncation value being used
% [4] The truncation level as a percentile
% [5] The l1 norm of the thresholded coefficients, in percent
% [6] The total variation of the reconstruction, in percent
% kilme The logicals identifying where the wavelets are killed
%
% Last modified by fjsimons-at-alum.mit.edu, 1/21/2011
defval('iface',0)
defval('tperc',85)
defval('precon',[1 1])
defval('xver',0);
if iface>0
v=v(:,:,iface);
end
switch wav
case 'D2'
% The 2-tap Daubechies wavelet transform
vw=angularD2WT(v,[J J],'forward',1);
case 'D4'
% The 4-tap Daubechies wavelet transform
vw=angularD4WT(v,[J J],precon,'forward',1);
case 'D6'
% The 6-tap Daubechies wavelet transform
vw=angularD6WT(v,[J J],precon,'forward',1);
end
% Test Ignace's code as far as D4 and D6 are concerned
whereitsat=fullfile(getenv('MFILES'),'ignaceloris');
if xver==1 && ~strcmp(wav,'D2')
disp(sprintf('\n%s\n',...
'Extra verification using Ignace Loris'' code'))
switch wav
case 'D4'
addpath(whereitsat)
difer([vw-angularD4MRWT(v,[J J],precon,'forward')]/prod(size(v)),8);
rmpath(whereitsat)
% Make sure to check my own inverse transform
difer([v-angularD4WT(vw,[J J],precon,'inverse',1)]/prod(size(v)),8);
case 'D6'
addpath(whereitsat)
difer([vw-angularD6MRWT(v,[J J],precon,'forward')]/prod(size(v)),8);
rmpath(whereitsat)
% Make sure to check my own inverse transform
difer([v-angularD6WT(vw,[J J],precon,'inverse',1)]/prod(size(v)),8);
end
disp(sprintf('\n%s\n',...
'Done with extra verification using Loris'' code'))
end
% Explicit and absolute zeroing of the MAGNITUDE
trunx=prctile(abs(vw(:)),tperc);
kilme=abs(vw)<trunx;
% With this truncation, how well do you approximate Africa?
vwt=vw;
vwt(kilme)=0;
% In other words, we still need the inversion
switch wav
case 'D2'
vwtinv=angularD2WT(vwt,[J J],'inverse',1);
case 'D4'
vwtinv=angularD4WT(vwt,[J J],precon,'inverse',1);
case 'D6'
vwtinv=angularD6WT(vwt,[J J],precon,'inverse',1);
end
% Now compute the Total Variation
% Take a quick look:
% plotoncube(vwtinv,'3D',[],[],[],[],[],[],0)
tv0=0; tv=0;
for i=1:size(v,3)
[Dx,Dy]=gradient(v(:,:,i));
% Take note, this is like ell_1 of the gradient!
tv0=tv0+sum(sum(sqrt([Dx.^2+Dy.^2])));
[Dx,Dy]=gradient(vwtinv(:,:,i));
tv =tv +sum(sum(sqrt([Dx.^2+Dy.^2])));
end
% The relative reconstruction error, etc, see the help
recer=[norm(vwtinv(:)-v(:))/norm(v(:))*100 ...
sum(kilme(:))/prod(size(kilme))*100 ...
trunx ...
tperc ...
norm(vwt(:),1)/norm(vw(:),1)*100 ...
tv/tv0*100];