/
perform_windowed_dct_transform.m
113 lines (90 loc) · 2.23 KB
/
perform_windowed_dct_transform.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
function y = perform_windowed_dct_transform(x,w,q,n, options)
% perform_windowed_dct_transform - compute a local DCT transform
%
% Forward transform:
% MF = perform_windowed_dct_transform(M,w,q,n, options);
% Backward transform:
% M = perform_windowed_dct_transform(MF,w,q,n, options);
%
% w is the width of the window used to perform local computation.
% q is the spacing betwen each window.
%
% MF(:,:,,i,j) contains the transform around point ((i-1)*q,(j-1)*q)+1.
%
% A typical use, for an redundancy of 2x2 could be w=2*q
%
% Copyright (c) 2007 Gabriel Peyre
options.transform_type = 'dct';
y = perform_windowed_fourier_transform(x,w,q,n, options);
return;
%%%%%%%%%%% OLD CODE %%%%%%%%%%%%
options.null = 0;
if size(x,3)>1
dir = -1;
if nargin<4
% assume power of 2 size
n = q*size(x,3);
n = 2^floor(log2(n));
end
else
dir = 1;
n = size(x,1);
end
% perform sampling
t = 1:q:n-w+1;
[Y,X] = meshgrid(t,t);
p = size(X,1);
w = ceil(w/2)*2;
w1 = w/2;
t = 0:w-1;
[dY,dX] = meshgrid(t,t);
X = reshape(X,[1 1 p p]);
Y = reshape(Y,[1 1 p p]);
X = repmat( X, [w w 1 1] );
Y = repmat( Y, [w w 1 1] );
dX = repmat( dX, [1 1 p p] );
dY = repmat( dY, [1 1 p p] );
X1 = X+dX;
Y1 = Y+dY;
if 0
I = find(X1<1); X1(I) = 1-X1(I);
I = find(X1>n); X1(I) = 2*n+1-X1(I);
I = find(Y1<1); Y1(I) = 1-Y1(I);
I = find(Y1>n); Y1(I) = 2*n+1-Y1(I);
end
% build a weight function
if isfield(options, 'window_type')
window_type = options.window_type;
else
window_type = 'constant';
end
if strcmp(window_type, 'sin') || strcmp(window_type, 'sine')
t = linspace(0,1,w);
W = sin(t(:)*pi).^2;
W = W * W';
elseif strcmp(window_type, 'constant')
W = ones(w);
else
error('Unkwnown winow.');
end
I = X1 + (Y1-1)*n;
if dir==1
y = x(I) .* repmat( W, [1 1 p p] );
y = my_dct_transform( y, +1 );
else
x = my_dct_transform( x, -1 );
weight = zeros(n); y = zeros(n);
for i=1:p
for j=1:p
y(I(:,:,i,j)) = y(I(:,:,i,j)) + x(:,:,i,j);
weight(I(:,:,i,j)) = weight(I(:,:,i,j)) + W;
end
end
y = real( y./weight );
end
function y = my_dct_transform(x,dir)
for i=1:size(x,3)
for j=1:size(x,4)
y(:,:,i,j) = perform_dct_transform(x(:,:,i,j),dir);
end
end