wtnd - Discrete N-D Wavelet Transform.
WTND(X,H,G,SCALES) calculates the N-D wavelet transform of matrix X,
where N is the number of dimensions of X. The second argument H is the
lowpass filter and the third argument G the highpass filter. WTND
generalizes the algorithm for UviWave WT to N dimensions.
For a vector, the output contains the coefficients of the DWT ordered from
the low pass residue at SCALE to the coefficients at the lowest scale, as
the following example ilustrates:
Output vector (k=3):
[------|------|------------|------------------------]
| | | |
| | | `-> 1st scale coefficients
| | `-----------> 2nd scale coefficients
| `--------------------> 3rd scale coefficients
`----------------> Low pass residue at 3rd scale
A matrix can be transformed over more than one dimension. For a 2D matrix,
the output would look like this: for every scale, the lowpass residue is
placed at the top-left corner of the corresponding subimage, the
horizontal high frequency band at the top-right, the vertical high
frequency band at the bottom-left and the diagonal high frequency band at
the bottom-right. Every successive wavelet subimage substitutes the
residue of the previous scale.
Example with 2 scales:
_______
2nd scale substituting the --> |_|_| | <-- 1st scale
first scale lowpass residue |_|_|___| horiz. detail
| | |
1st scale ----> |___|___| <-- 1st scale
vertical detail diagonal detail
Notes 1: Dimensions to transform
--------------------------------
If the passed matrix X has N dimensions, then the default behaviour is to
transform all N dimensions. It is possible to select a subset of
dimensions to transform, using the P_DIMS input (see below). For example
we could choose to do a 2D transform on a 3D matrix.
WTND(X,H,G,SCALES,DEL1,DEL2) calculates the N-D wavelet
transform of matrix X, but also allows the user to change the alignment of
the outputs with respect to the input signal. This effect is achieved by
setting to DEL1 and DEL2 the delays of H and G respectively. The default
values of DEL1 and DEL2 are calculated using the function WTCENTER.
Note 2: Using P_DIMS
--------------------
WTND(X,H,G,SCALES,DEL1,DEL2,P_DIMS) adds the ability to wt over a subset
of the matrix dimensions. P_DIMS should be a vector of flag values, one
value per dimension of X; each dimension with a flag value of 1 is wt'ed,
dimensions of size 1 or with a 0 flag value are ignored. Thus, if we have
a 3D matrix, the default value for P_DIMS will be [1 1 1]; if you want to
transform just over the rows of this matrix, you could override the
default by passing a value of P_DIMS of [0 1 0]; P_DIMS can also be a
scalar, in which case it specifies that the first P_DIMS dimensions
should be processed; for example, for a 3D matrix, a P_DIMS value of 2
will expand to a P_DIMS vector of [1 1 0].
Note 3: Different dimensions, different filters
-----------------------------------------------
Each dimension can have its own wavelet filter; the input filters h and
g can be given as cell arrays, each cell containing a filter for the
matching processed dimension. For example, let use imagine you are passing
a 3D matrix, but only want to transform the first 2 dimensions. You want
a different filter for the first and second dimensions. You would use a
P_DIMS of [1 1 0] (see above), and pass a two element cell array in h and
in g, with the h filter for dimension 2 in h{2}, etc. Similarly the
delays DEL1 and DEL2 can be vectors, with one delay for each processed
dimension.
Returns
Y - transformed matrix
SCALES - number of scales transformed
(SCALES can be set in routine if not passed as input)
See also: IWT, WT2D, IWT2D, WTCENTER, ISPLIT.
Based on wt.m from UviWave 3.0, with thanks - see below
$Id: wtnd.m,v 1.1 2004/09/26 04:00:24 matthewbrett Exp $0001 function [y, scales] = wtnd(x, h, g, scales, del1, del2, p_dims) 0002 % wtnd - Discrete N-D Wavelet Transform. 0003 % 0004 % WTND(X,H,G,SCALES) calculates the N-D wavelet transform of matrix X, 0005 % where N is the number of dimensions of X. The second argument H is the 0006 % lowpass filter and the third argument G the highpass filter. WTND 0007 % generalizes the algorithm for UviWave WT to N dimensions. 0008 % 0009 % For a vector, the output contains the coefficients of the DWT ordered from 0010 % the low pass residue at SCALE to the coefficients at the lowest scale, as 0011 % the following example ilustrates: 0012 % 0013 % Output vector (k=3): 0014 % 0015 % [------|------|------------|------------------------] 0016 % | | | | 0017 % | | | `-> 1st scale coefficients 0018 % | | `-----------> 2nd scale coefficients 0019 % | `--------------------> 3rd scale coefficients 0020 % `----------------> Low pass residue at 3rd scale 0021 % 0022 % A matrix can be transformed over more than one dimension. For a 2D matrix, 0023 % the output would look like this: for every scale, the lowpass residue is 0024 % placed at the top-left corner of the corresponding subimage, the 0025 % horizontal high frequency band at the top-right, the vertical high 0026 % frequency band at the bottom-left and the diagonal high frequency band at 0027 % the bottom-right. Every successive wavelet subimage substitutes the 0028 % residue of the previous scale. 0029 % 0030 % Example with 2 scales: 0031 % _______ 0032 % 2nd scale substituting the --> |_|_| | <-- 1st scale 0033 % first scale lowpass residue |_|_|___| horiz. detail 0034 % | | | 0035 % 1st scale ----> |___|___| <-- 1st scale 0036 % vertical detail diagonal detail 0037 % 0038 % Notes 1: Dimensions to transform 0039 % -------------------------------- 0040 % 0041 % If the passed matrix X has N dimensions, then the default behaviour is to 0042 % transform all N dimensions. It is possible to select a subset of 0043 % dimensions to transform, using the P_DIMS input (see below). For example 0044 % we could choose to do a 2D transform on a 3D matrix. 0045 % 0046 % WTND(X,H,G,SCALES,DEL1,DEL2) calculates the N-D wavelet 0047 % transform of matrix X, but also allows the user to change the alignment of 0048 % the outputs with respect to the input signal. This effect is achieved by 0049 % setting to DEL1 and DEL2 the delays of H and G respectively. The default 0050 % values of DEL1 and DEL2 are calculated using the function WTCENTER. 0051 % 0052 % Note 2: Using P_DIMS 0053 % -------------------- 0054 % 0055 % WTND(X,H,G,SCALES,DEL1,DEL2,P_DIMS) adds the ability to wt over a subset 0056 % of the matrix dimensions. P_DIMS should be a vector of flag values, one 0057 % value per dimension of X; each dimension with a flag value of 1 is wt'ed, 0058 % dimensions of size 1 or with a 0 flag value are ignored. Thus, if we have 0059 % a 3D matrix, the default value for P_DIMS will be [1 1 1]; if you want to 0060 % transform just over the rows of this matrix, you could override the 0061 % default by passing a value of P_DIMS of [0 1 0]; P_DIMS can also be a 0062 % scalar, in which case it specifies that the first P_DIMS dimensions 0063 % should be processed; for example, for a 3D matrix, a P_DIMS value of 2 0064 % will expand to a P_DIMS vector of [1 1 0]. 0065 % 0066 % Note 3: Different dimensions, different filters 0067 % ----------------------------------------------- 0068 % 0069 % Each dimension can have its own wavelet filter; the input filters h and 0070 % g can be given as cell arrays, each cell containing a filter for the 0071 % matching processed dimension. For example, let use imagine you are passing 0072 % a 3D matrix, but only want to transform the first 2 dimensions. You want 0073 % a different filter for the first and second dimensions. You would use a 0074 % P_DIMS of [1 1 0] (see above), and pass a two element cell array in h and 0075 % in g, with the h filter for dimension 2 in h{2}, etc. Similarly the 0076 % delays DEL1 and DEL2 can be vectors, with one delay for each processed 0077 % dimension. 0078 % 0079 % Returns 0080 % Y - transformed matrix 0081 % SCALES - number of scales transformed 0082 % (SCALES can be set in routine if not passed as input) 0083 % 0084 % See also: IWT, WT2D, IWT2D, WTCENTER, ISPLIT. 0085 % 0086 % Based on wt.m from UviWave 3.0, with thanks - see below 0087 % 0088 % $Id: wtnd.m,v 1.1 2004/09/26 04:00:24 matthewbrett Exp $ 0089 0090 %-------------------------------------------------------- 0091 % wt is Copyright (C) 1994, 1995, 1996, by Universidad de Vigo 0092 % 0093 % Uvi_Wave is free software; you can redistribute it and/or modify it 0094 % under the terms of the GNU General Public License as published by the 0095 % Free Software Foundation; either version 2, or (at your option) any 0096 % later version. 0097 % 0098 % Uvi_Wave is distributed in the hope that it will be useful, but WITHOUT 0099 % ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 0100 % FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 0101 % for more details. 0102 % 0103 % You should have received a copy of the GNU General Public License 0104 % along with Uvi_Wave; see the file COPYING. If not, write to the Free 0105 % Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. 0106 % 0107 % Authors: Sergio J. Garcia Galan 0108 % Cristina Sanchez Cabanelas 0109 % e-mail: Uvi_Wave@tsc.uvigo.es 0110 %-------------------------------------------------------- 0111 0112 % Programmers' notes 0113 % 0114 % We do the wt by iteratively switching the matrix dimension we want to 0115 % process to be the first dimension (columns), then filtering. The algorithm 0116 % goes up to 4D, but could be extended further by adding the relevant lines 0117 % to the switch statements, or made general by using subsref commands: 0118 % [n_dims = 2; d = subsref(c, struct('type', '()', 'subs', repmat({:}, 1, 0119 % n_dims)))]; I haven't used the subsref commands here because they can be 0120 % considerably slower than direct [c(:,:)] type referencing 0121 0122 % ----------------------------------- 0123 % CHECK PARAMETERS AND OPTIONS 0124 % ----------------------------------- 0125 0126 if nargin < 3 0127 error('Need at least matrix to wt, and wt filters'); 0128 end 0129 if nargin < 4 0130 scales = []; 0131 end 0132 if nargin < 5 0133 del1 = []; 0134 end 0135 if nargin < 6 0136 del2 = []; 0137 end 0138 if xor(isempty(del1), isempty(del2)) 0139 error('You cannot specify only one of the two delays'); 0140 end 0141 if nargin < 7 0142 p_dims = []; 0143 end 0144 0145 % ------------------------ 0146 % Process input parameters 0147 % ------------------------ 0148 0149 % Input image dimensions 0150 sz = size(x); 0151 n_dims = length(sz); 0152 0153 % Process p_dims 0154 p_d_l = prod(size(p_dims)); 0155 if p_d_l == 0 % empty 0156 p_dims = ones(1, n_dims); 0157 elseif p_d_l == 1 0158 p_dims = [ones(1, p_dims) zeros(1, n_dims-p_dims)]; 0159 elseif p_d_l < n_dims 0160 p_dims = [p_dims(:)' ones(1, n_dims-p_d_l)]; 0161 elseif p_d_l > n_dims 0162 error('p_dims array has too many values for passed matrix'); 0163 end 0164 p_dims = p_dims & sz > 1; 0165 if ~any(p_dims), y = x; return, end 0166 n_p_dims = sum(p_dims); 0167 0168 % make filters into cells, and expand to N-D if necessary 0169 if ~iscell(h), h = {h}; end 0170 if prod(size(h))==1, h = repmat(h, 1, n_p_dims); end 0171 if ~iscell(g), g = {g}; end 0172 if prod(size(g))==1, g = repmat(g, 1, n_p_dims); end 0173 0174 %-------------------------- 0175 % DELAY CALCULATION 0176 %-------------------------- 0177 0178 if isempty(del1) 0179 % Calculate delays as the C.O.E. of the filters 0180 for d = 1:n_p_dims 0181 dlp(d) = wtcenter( h{d} ); 0182 dhp(d) = wtcenter( g{d} ); 0183 end 0184 % difference between them must be even 0185 dhp = dhp + (rem(dhp-dlp,2)~=0); 0186 else 0187 % Other experimental filter delays can be forced from the arguments 0188 o = ones(1, n_p_dims); 0189 dlp=del1; if prod(size(dlp))==1, dlp = o * dlp; end 0190 dhp=del2; if prod(size(dhp))==1, dhp = o * dhp; end 0191 end 0192 0193 % get output and scales dimensions 0194 [wt_sz sc_in_sz sc_out_sz scales] = wt_dims(sz, scales, p_dims); 0195 0196 % Offsets 0197 offsets = wt_sz - sz; 0198 0199 %------------------------------ 0200 % START THE ALGORITHM 0201 %------------------------------ 0202 0203 % For doing tricksy permute thing in loop below (see help wt_dim_shifts) 0204 [shifts end_shift] = wt_dim_shifts(p_dims, sz); 0205 dims = [1:n_dims]'; 0206 0207 % t is the input matrix to be wt'ed 0208 t = x; 0209 0210 % y is the output matrix 0211 y = zeros(wt_sz(:)'); 0212 0213 % For every scale (iteration)... 0214 for sc=1:scales 0215 % get new data block 0216 in_sz = sc_in_sz(sc, :); 0217 if sc > 1 0218 % following assumes UviWave ordering 0219 switch n_dims 0220 case 2 0221 t = t(1:in_sz(1), 1:in_sz(2)); 0222 case 3 0223 t = t(1:in_sz(1), 1:in_sz(2), 1:in_sz(3)); 0224 case 4 0225 t = t(1:in_sz(1), 1:in_sz(2), 1:in_sz(3), 1:in_sz(4)); 0226 otherwise 0227 error('Not implemented - but go ahead and do it if you need it'); 0228 end 0229 end 0230 0231 for d = 1:n_p_dims 0232 0233 % move next dimension to column dimension 0234 s = shifts(d); 0235 if s, t = permute(t, circshift(dims, s)); end 0236 0237 % all the hard work is in C 0238 t = do_wtx(t, h{d}, g{d}, dlp(d), dhp(d)); 0239 0240 end 0241 0242 % and back to original shape 0243 if end_shift 0244 t = permute(t, circshift(dims, end_shift)); 0245 end 0246 0247 % reset offsets 0248 out_sz = sc_out_sz(sc, :); 0249 offsets = offsets - (out_sz - in_sz); 0250 op1 = offsets + 1; 0251 0252 % set data block into output 0253 % following assumes UviWave ordering 0254 if sc == 1 & ~any(offsets) % for speed, do the simplest case 0255 y = t; 0256 else 0257 switch n_dims 0258 case 2 0259 y(op1(1):out_sz(1) + offsets(1), ... 0260 op1(2):out_sz(2) + offsets(2)) = t; 0261 case 3 0262 y(op1(1):out_sz(1) + offsets(1), ... 0263 op1(2):out_sz(2) + offsets(2), ... 0264 op1(3):out_sz(3) + offsets(3)) = t; 0265 case 4 0266 y(op1(1):out_sz(1) + offsets(1), ... 0267 op1(2):out_sz(2) + offsets(2), ... 0268 op1(3):out_sz(3) + offsets(3), ... 0269 op1(4):out_sz(4) + offsets(4)) = t; 0270 otherwise 0271 error('Not implemented'); 0272 end 0273 end 0274 end 0275 0276 %------------------------------ 0277 % END OF THE ALGORITHM 0278 %------------------------------ 0279