1% RES = reconSFpyr(PYR, INDICES, LEVS, BANDS, TWIDTH)
2%
3% Reconstruct image from its steerable pyramid representation, in the Fourier
4% domain, as created by buildSFpyr.
5%
6% PYR is a vector containing the N pyramid subbands, ordered from fine
7% to coarse.  INDICES is an Nx2 matrix containing the sizes of
8% each subband.  This is compatible with the MatLab Wavelet toolbox.
9%
10% LEVS (optional) should be a list of levels to include, or the string
11% 'all' (default).  0 corresonds to the residual highpass subband.
12% 1 corresponds to the finest oriented scale.  The lowpass band
13% corresponds to number spyrHt(INDICES)+1.
14%
15% BANDS (optional) should be a list of bands to include, or the string
16% 'all' (default).  1 = vertical, rest proceeding anti-clockwise.
17%
18% TWIDTH is the width of the transition region of the radial lowpass
19% function, in octaves (default = 1, which gives a raised cosine for
20% the bandpass filters).
21
22%%% MODIFIED VERSION, 7/04, uses different lookup table for radial frequency!
23
24% Eero Simoncelli, 5/97.
25
26function res = reconSFpyr(pyr, pind, levs, bands, twidth)
27
28%%------------------------------------------------------------
29%% DEFAULTS:
30
31if (exist('levs') ~= 1)
32  levs = 'all';
33end
34
35if (exist('bands') ~= 1)
36  bands = 'all';
37end
38
39if (exist('twidth') ~= 1)
40  twidth = 1;
41elseif (twidth <= 0)
42  fprintf(1,'Warning: TWIDTH must be positive.  Setting to 1.\n');
43  twidth = 1;
44end
45
46%%------------------------------------------------------------
47
48nbands = spyrNumBands(pind);
49
50maxLev =  1+spyrHt(pind);
51if strcmp(levs,'all')
52  levs = [0:maxLev]';
53else
54  if (any(levs > maxLev) | any(levs < 0))
55    error(sprintf('Level numbers must be in the range [0, %d].', maxLev));
56  end
57  levs = levs(:);
58end
59
60if strcmp(bands,'all')
61  bands = [1:nbands]';
62else
63  if (any(bands < 1) | any(bands > nbands))
64    error(sprintf('Band numbers must be in the range [1,3].', nbands));
65  end
66  bands = bands(:);
67end
68
69%----------------------------------------------------------------------
70
71dims = pind(1,:);
72ctr = ceil((dims+0.5)/2);
73
74[xramp,yramp] = meshgrid( ([1:dims(2)]-ctr(2))./(dims(2)/2), ...
75    ([1:dims(1)]-ctr(1))./(dims(1)/2) );
76angle = atan2(yramp,xramp);
77log_rad = sqrt(xramp.^2 + yramp.^2);
78log_rad(ctr(1),ctr(2)) =  log_rad(ctr(1),ctr(2)-1);
79log_rad  = log2(log_rad);
80
81%% Radial transition function (a raised cosine in log-frequency):
82[Xrcos,Yrcos] = rcosFn(twidth,(-twidth/2),[0 1]);
83Yrcos = sqrt(Yrcos);
84YIrcos = sqrt(abs(1.0 - Yrcos.^2));
85
86if (size(pind,1) == 2)
87  if (any(levs==1))
88    resdft = fftshift(fft2(pyrBand(pyr,pind,2)));
89  else
90    resdft = zeros(pind(2,:));
91  end
92else
93  resdft = reconSFpyrLevs(pyr(1+prod(pind(1,:)):size(pyr,1)), ...
94      pind(2:size(pind,1),:), ...
95      log_rad, Xrcos, Yrcos, angle, nbands, levs, bands);
96end
97
98lo0mask = pointOp(log_rad, YIrcos, Xrcos(1), Xrcos(2)-Xrcos(1), 0);
99resdft = resdft .* lo0mask;
100
101%% residual highpass subband
102if any(levs == 0)
103  hi0mask = pointOp(log_rad, Yrcos, Xrcos(1), Xrcos(2)-Xrcos(1), 0);
104  hidft = fftshift(fft2(subMtx(pyr, pind(1,:))));
105  resdft = resdft + hidft .* hi0mask;
106end
107
108res = real(ifft2(ifftshift(resdft)));
109