THESE/codes/wave/allcode/afb3D.m

   1 function [lo, hi] = afb3D(x, af1, af2, af3)
   2
   3 % 3D Analysis Filter Bank
   4 %
   5 % USAGE:
   6 %    [lo, hi] = afb3D(x, af1, af2, af3);
   7 % INPUT:
   8 %    x - N1 by N2 by N3 array matrix, where
   9 %        1) N1, N2, N3 all even
  10 %        2) N1 >= 2*length(af1)
  11 %        3) N2 >= 2*length(af2)
  12 %        4) N3 >= 2*length(af3)
  13 %    afi - analysis filters for dimension i
  14 %       afi(:, 1) - lowpass filter
  15 %       afi(:, 2) - highpass filter
  16 % OUTPUT:
  17 %    lo - lowpass subband
  18 %    hi{d}, d = 1..7 - highpass subbands
  19 % EXAMPLE:
  20 %    x = rand(32,64,16);
  21 %    [af, sf] = farras;
  22 %    [lo, hi] = afb3D(x, af, af, af);
  23 %    y = sfb3D(lo, hi, sf, sf, sf);
  24 %    err = x - y;
  25 %    max(max(max(abs(err))))
  26 %
  27 % WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
  28 % http://taco.poly.edu/WaveletSoftware/
  29
  30 if nargin < 3
  31    af2 = af1;
  32    af3 = af1;
  33 end
  34
  35 % filter along dimension 1
  36 [L, H] = afb3D_A(x, af1, 1);
  37
  38 % filter along dimension 2
  39 [LL LH] = afb3D_A(L, af2, 2);
  40 [HL HH] = afb3D_A(H, af2, 2);
  41
  42 % filter along dimension 3
  43 [LLL LLH] = afb3D_A(LL, af3, 3);
  44 [LHL LHH] = afb3D_A(LH, af3, 3);
  45 [HLL HLH] = afb3D_A(HL, af3, 3);
  46 [HHL HHH] = afb3D_A(HH, af3, 3);
  47
  48 lo    = LLL;
  49 hi{1} = LLH;
  50 hi{2} = LHL;
  51 hi{3} = LHH;
  52 hi{4} = HLL;
  53 hi{5} = HLH;
  54 hi{6} = HHL;
  55 hi{7} = HHH;
  56
  57
  58 % LOCAL FUNCTION
  59
  60 function [lo, hi] = afb3D_A(x, af, d)
  61
  62 % 3D Analysis Filter Bank
  63 % (along one dimension only)
  64 %
  65 % [lo, hi] = afb3D_A(x, af, d);
  66 % INPUT:
  67 %    x - N1xN2xN2 matrix, where min(N1,N2,N3) > 2*length(filter)
  68 %           (Ni are even)
  69 %    af - analysis filter for the columns
  70 %    af(:, 1) - lowpass filter
  71 %    af(:, 2) - highpass filter
  72 %    d - dimension of filtering (d = 1, 2 or 3)
  73 % OUTPUT:
  74 %     lo, hi - lowpass, highpass subbands
  75 %
  76 % % Example
  77 % x = rand(32,64,16);
  78 % [af, sf] = farras;
  79 % d = 2;
  80 % [lo, hi] = afb3D_A(x, af, d);
  81 % y = sfb3D_A(lo, hi, sf, d);
  82 % err = x - y;
  83 % max(max(max(abs(err))))
  84
  85 lpf = af(:, 1);     % lowpass filter
  86 hpf = af(:, 2);     % highpass filter
  87
  88 % permute dimensions of x so that dimension d is first.
  89 p = mod(d-1+[0:2], 3) + 1;
  90 x = permute(x, p);
  91
  92 % filter along dimension 1
  93 [N1, N2, N3] = size(x);
  94 L = size(af, 1)/2;
  95 x = cshift3D(x, -L, 1);
  96 lo = zeros(L+N1/2, N2, N3);
  97 hi = zeros(L+N1/2, N2, N3);
  98
  99 for k = 1:N3
 100    lo(:, :, k) = upfirdn(x(:, :, k), lpf, 1, 2);
 101 end
 102 lo(1:L, :, :) = lo(1:L, :, :) + lo([1:L]+N1/2, :, :);
 103 lo = lo(1:N1/2, :, :);
 104
 105 for k = 1:N3
 106    hi(:, :, k) = upfirdn(x(:, :, k), hpf, 1, 2);
 107 end
 108 hi(1:L, :, :) = hi(1:L, :, :) + hi([1:L]+N1/2, :, :);
 109 hi = hi(1:N1/2, :, :);
 110
 111 % permute dimensions of x (inverse permutation)
 112 lo = ipermute(lo, p);
 113 hi = ipermute(hi, p);
 114