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