### Digital Image Processing

- Digital Image Processing---Wavelets
- 1. Introduction
- 2. Related Work
- 3. Approach
- 4. Experimental Setup
- 4.1 Haar filters, scaling, and wavelet functions.
- 4.2 A simple FWT using Haar filters.
- 4.3 Comparing the execution times ofwavefast and wavedec2.
- 4.4 Wavelet Toolbox functions for manipulating transform decomposition vector c.
- 4.5 Manipulating c with wavecut and wavecopy.
- 4.6 Transform coefficient display using wavedisplay.
- 4.7 Comparing the execution times of waveback and waverec2.
- 4.8 Wavelet directionality and edge detection.
- 4.9 Wavelet-based image smoothing or blurring.
- 4.10 Progressive reconstruction.
- 5. Conclusion
- 6. Appendix

# Digital Image Processing—Wavelets

**Abstract**

## 1. Introduction

## 2. Related Work

#### 2.1 Preliminary knowledge

Where CF is the center frequency of the wavelet ,s is the wavelet scale and dalta T is the sampling interval.

2. Bump Wavelet

3. Analytic Morlet Wavelet

#### 2.2 Background

for

and orthonormal because

b. The set of functions that can be represented as a series expansion of φj.k at low scales or resolutions (i.e., small j) is contained within those that can be represented at higher scales.

c.The only function that can be represented at every scale is f(x)=0.

d.Any function can be represented width arbitrary precision as j →∞.

## 3. Approach

#### 3.1 The Fast Wavelet Transform

###### 3.1.1 FWTs Using the Wavelet Toolbox

[Lo_D, Hi_D, Lo_R, Hi_RJ = wfilters(wname)

[F1,F2] = wfilters(‘wname’,‘type’)

- Lo_D, the decomposition low-pass filter
- Hi_D, the decomposition high-pass filter
- Lo_R, the reconstruction low-pass filter
- Hi_R, the reconstruction high-pass filter

Lo_D and Hi_D | (Decomposition filters) | If ‘type’ = ‘d’ |
---|---|---|

Lo_R and Hi_R | (Reconstruction filters) | If ‘type’ = ‘r’ |

Lo_D and Lo_R | (Low-pass filters) | If ‘type’ = ‘l’ |

Hi_D and Hi_R | (High-pass filters) | If ‘type’ = ‘h’ |

Wavelet | wfamily | wname |
---|---|---|

Haar | ‘haar’ | ‘haar’ |

Daubechies | ‘db’ | ‘db1’ or ‘haar’, ‘db2’, … ,‘db10’, … , ‘db45’ |

Coiflets | ‘coif’ | ‘coif1’, … , ‘coif5’ |

Symlets | ‘sym’ | ‘sym2’, … , ‘sym8’, … ,‘sym45’ |

Discrete Meyer | ‘dmey’ | ‘dmey’ |

Biorthogonal | ‘bior’ | ‘bior1.1’, ‘bior1.3’,… ‘bior6.8’ |

Reverse Biorthogonal | ‘rbio’ | ‘rbio1.1’, ‘rbio1.3’, … ‘rbio6.8’ |

In order to print the description information of wfamily(see table 7-2) in the command window of MATLAB, enter it at the MALTAB prompt:

waveinfo(wfamily)

[phi, psi, xval] = wavefun(wname, iter)

[C,S] = wavedec2(X,N,Lo D,Hi D)

[C, S] = wavedec2(X, N, wname)

STATUS | Description |
---|---|

‘sym’ | The image is extended by mirror reflecting it across its borders. This is the normal default mode |

'zpd ’ | The image is extended by padding with a value of 0. |

‘spd’ , 'sp1 ’ | The image is extended by first-order derivative extrapolation-or padding with a linear extension of the outmost two border values. |

'spO ’ | The image is extended by periodic padding. |

'ppd ’ | The image is extended by mirror reflecting it across its borders. This is the normal default mode |

‘per’ | Description The image is extended by mirror reflecting it across its borders. This is the normal default mode. The image is extended by padding with a value of 0. The image is extended by first-order derivative extrapolation-or padding with a linear extension of the outmost two border values. The image is extended by extrapolating the border values-that is, by boundary value replication. The image is extended by periodic padding. The image i extended by periodic padding after it has been padded (if necessary) to an even size using ‘spO’ extension. |

###### 3.1.2 FWTs without the Wavelet Toolbox

function [varargout] = wavefilter(wname, type)

function [c, s] = wavefast(x, n, varargin)

function nc = addcoefs(c, x, pages)

function [y, keep] = symextend(x, fl, pages)

function y = symconv(x, h, type, fl, keep, pages)

#### 3.2 Working with Wavelet Decomposition Structures

Then

###### 3.2.1 Editing Wavelet Decomposition Coefficients without the Wavelet Toolbox

function [varargout] = wavework(opcode, type, c, s, n, x)

y = reshape(x,m,n)

function [nc, y] = wavecut(type, c, s, n)

function y = wavecopy(type, c, s, n)

function nc = wavepaste(type, c, s, n, x)

###### 3.2.2 Displaying Wavelet Decomposition Coefficients

function w = wavedisplay(c, s, scale, border)

#### 3.3 The Inverse Fast Wavelet Transform

g = waverec2(C, S, wname)

g = waverec2(C, S, Lo_R, Hi_R)

function [varargout] = waveback(c, s, varargin)

Function waveback uses the equivalent computation :

#### 3.4 Wavelets in Image Processing

1. Compute the two-dimensional wavelet transform of an image.

2. Alter the transform coefficients.

3. Compute the inverse transform.

## 4. Experimental Setup

#### 4.1 Haar filters, scaling, and wavelet functions.

```
[Lo_D, Hi_D, Lo_R, Hi_R] = wfilters( 'haar')
```

The code operation effect is as follows:

```
waveinfo( 'haar');
```

```
[phi, psi, xval] = wavefun( 'haar', 10);
xaxis = zeros(size(xval));
subplot(1,2,1);
plot(xval, phi, 'k', xval, xaxis, '--k');
axis([0 1 -1.5 1.5]); axis square;
title('Haar Scaling Function');
subplot(1,2,2); plot(xval, psi, 'k', xval, xaxis, '--k');
axis([0 1 -1.5 1.5]); axis square;
title('Haar Wavelet Function');
```

#### 4.2 A simple FWT using Haar filters.

```
f = magic(4);
```

```
[c1,s1]=wavedec2(f,1,'haar')
```

When the single-scale transform described above is extended to two scales, we get

```
[c2, s2] = wavedec2(f, 2, 'haar')
```

#### 4.3 Comparing the execution times ofwavefast and wavedec2.

```
function [ratio, maxdiff] = fwtcompare(f, n, wname)
% Get transform and computation time for wavedec2.
w1 = @() wavedec2(f, n, wname);
reftime = timeit(w1);
% Get transform and computation time for wavefast.
w2 = @() wavefast(f, n, wname);
t2 = timeit(w2);
% Compare the results.
ratio = t2 / reftime;
maxdiff = abs(max(w1() - w2()));
```

```
f = imread('D:\Dip\Chapter7\sun.jpg');
[ratio, maxdifference] = fwtcompare(f, 5, 'db4')
```

#### 4.4 Wavelet Toolbox functions for manipulating transform decomposition vector c.

```
f = magic(8);
[c1,s1]=wavedec2(f,3,'haar');
approx = appcoef2(c1,s1,'haar');
horizdet2 = detcoef2( 'h', c1 , s1, 2);
newc1 = wthcoef2('h', c1 , s1, 2);
newhorizdet2 = detcoef2('h', newc1 , s1 , 2)
```

#### 4.5 Manipulating c with wavecut and wavecopy.

```
f = magic(8);
[c1,s1]=wavedec2(f,3,'haar');
approx = wavecopy('a',c1,s1);
horizdet2 = wavecopy( 'h', c1 , s1, 2);
[newc1,horizdet2] = wavecut('h', c1 , s1, 2);
newhorizdet2 = wavecopy('h', newc1 , s1 , 2)
```

Note that all extracted matrices are identical to those of the previous .

#### 4.6 Transform coefficient display using wavedisplay.

```
f = imread('D:\Dip\Chapter7\Fig0704.tif');
[c,s] = wavefast(f, 2, 'db4');title('(a) Automatic scaling; ');
wavedisplay(c, s); title('(a) Automatic scaling; ');
figure; wavedisplay(c, s, 8);title('(b) additional scaling by 8; ');
figure; wavedisplay(c, s,-8) ; title('(c) absolute values scaled by 8. ');
```

#### 4.7 Comparing the execution times of waveback and waverec2.

```
function [ratio, maxdiff] = ifwtcompare(f, n, wname)
[c1, s1] = wavedec2(f, n, wname);
w1 = @() waverec2(c1, s1, wname);
reftime = timeit(w1);
% Compute the transform and get output and computation time for
% waveback.
[c2, s2] = wavefast(f, n, wname);
w2 = @() waveback(c2, s2, wname);
t2 = timeit(w2);
% Compare the results.
ratio = t2 / reftime;
diff = double(w1()) - w2();
maxdiff = abs(max(diff(:)));
```

```
f = imread('D:\Dip\Chapter7\sun.jpg');
[ratio, maxdifference] = ifwtcompare(f, 5, 'db4');
```

#### 4.8 Wavelet directionality and edge detection.

```
f = imread('D:\Dip\Chapter7\Fig0707(a).tif');
subplot(2, 2, 1);imshow(f); title('(a) A simple test image;');
[c, s] = wavefast(f, 1, 'sym4');
subplot(2, 2, 2); imshow(wavedisplay(c, s, -6)) ; title('(b) its wavelet transform;');
[ nc, y] = wavecut ( 'a' , c, s) ;
subplot(2, 2, 3); imshow(wavedisplay(nc, s, -6)) ; title('(c) the transform modified by zeroing all approximation coefficients;');
edges = abs (waveback ( nc, s, 'sym4' ) ) ;
subplot(2, 2, 4); imshow(mat2gray(edges)); title('( d) the edge image resulting from computing the absolute value of the inverse transform.');
```

#### 4.9 Wavelet-based image smoothing or blurring.

```
f = imread('D:\Dip\Chapter7\Fig0707(a).tif');
figure;imshow(f);title('(a)the original image');
[c, s] = wavefast(f, 4, 'sym4'); %Fast wavelet decomposition is performed on scale 4
figure;
wavedisplay(c, s, 20) ; %Reconstruction on 20 scales
title('(b)sym4 four wavelet transform results');
[c,g8] = wavezero(c, s,1,'sym4'); %This function will automatically draw a window and set the size detail to 0
title('(c)Set the details of scale 1 to 0 after the inverse transformation');
[c, g8] = wavezero(c, s, 2, 'sym4');
title('(d)Set the details of scale 2 to 0 after the inverse transformation');
[c, g8] = wavezero(c, s, 3, 'sym4');
title('(e)Set the details of scale 3 to 0 after the inverse transformation');
[c, g8] = wavezero(c, s, 4, 'sym4');
title('(f)Set the details of scale 4 to 0 after the inverse transformation');
```

#### 4.10 Progressive reconstruction.

```
f = rgb2gray(imread('D:\Dip\Chapter7\sun.jpg'));
%figure;imshow(f);title('(d)the original image');
figure;[c, s] = wavefast(f, 4, 'jpeg9.7'); %4 Scale wavelet decomposition
wavedisplay(c,s,8); %Reconstruction on 20 scales
title('(a)4 scaling wavelet decomposition coefficient (reconstruction)');
f=wavecopy('a',c,s);%Wavelet decomposition myopia coefficient matrix is obtained
figure,imshow(mat2gray(f));
title('(b)Wavelet image of myopia coefficient');
[c s]=waveback(c,s,'jpeg9.7',1);
f=wavecopy('a',c,s);
figure,imshow(mat2gray(f));%The reconstructed myopic wavelet coefficient image
title('(c)The image was reconstructed 1 time');
[c s]=waveback(c,s,'jpeg9.7',1);
f=wavecopy('a',c,s);
figure,imshow(mat2gray(f));%The reconstructed myopic wavelet coefficient image
title('(d)The image was reconstructed twice');
[c s]=waveback(c,s,'jpeg9.7',1);
f=wavecopy('a',c,s);
figure,imshow(mat2gray(f));%The reconstructed myopic wavelet coefficient image
title('(e)The image was reconstructed 3 times');
[c s]=waveback(c,s,'jpeg9.7',1);
f=wavecopy('a',c,s);
figure,imshow(mat2gray(f));%The reconstructed myopic wavelet coefficient image
title('(f)The image was reconstructed 4 times');
```

## 5. Conclusion

## 6. Appendix

**Example 1**:

Write the code:

```
% Set the name of the wavelet
wname = 'db5';
%wname = 'bior2.6';
%wname = 'rbio2.6';
% Calculate four filters for a given wavelet
[Lo_D,Hi_D,Lo_R,Hi_R] = wfilters(wname);
subplot(221); stem(Lo_D);
title('Decomposition low-pass filter');
subplot(222); stem(Hi_D);
title('Decomposition high-pass filter');
subplot(223); stem(Lo_R);
title('Reconstruction low-pass filter');
subplot(224); stem(Hi_R);
title('Reconstruction high-pass filter');
xlabel('The four filters for db5')
```

The code operation effect is as follows:

**Example 2**

Write the code:

```
wfamily = 'haar';
waveinfo(wfamily);
```

The code operation effect is as follows:

**Example 3**

Write the code::

```
wname = 'rbio2.6';
[phil,psil,phi2,psi2rxval] = wavefun(wname,2)
```

The code operation effect is as follows:

**Example 4**

```
function [varargout] = wavefilter(wname, type)
% Check the input and output arguments.
error(nargchk(1, 2, nargin));
if (nargin == 1 && nargout ~= 4) || (nargin == 2 && nargout ~= 2)
error('Invalid number of output arguments.');
end
if nargin == 1 && ~ischar(wname)
error('WNAME must be a string.');
end
if nargin == 2 && ~ischar(type)
error('TYPE must be a string.');
end
% Create filters for the requested wavelet.
switch lower(wname)
case {'haar', 'db1'}
ld = [1 1]/sqrt(2); hd = [-1 1]/sqrt(2);
lr = ld; hr = -hd;
case 'db4'
ld = [-1.059740178499728e-002 3.288301166698295e-002 ...
3.084138183598697e-002 -1.870348117188811e-001 ...
-2.798376941698385e-002 6.308807679295904e-001 ...
7.148465705525415e-001 2.303778133088552e-001];
t = (0:7);
hd = ld; hd(end:-1:1) = cos(pi * t) .* ld;
lr = ld; lr(end:-1:1) = ld;
hr = cos(pi * t) .* ld;
case 'sym4'
ld = [-7.576571478927333e-002 -2.963552764599851e-002 ...
4.976186676320155e-001 8.037387518059161e-001 ...
2.978577956052774e-001 -9.921954357684722e-002 ...
-1.260396726203783e-002 3.222310060404270e-002];
t = (0:7);
hd = ld; hd(end:-1:1) = cos(pi * t) .* ld;
lr = ld; lr(end:-1:1) = ld;
hr = cos(pi * t) .* ld;
case 'bior6.8'
ld = [0 1.908831736481291e-003 -1.914286129088767e-003 ...
-1.699063986760234e-002 1.193456527972926e-002 ...
4.973290349094079e-002 -7.726317316720414e-002 ...
-9.405920349573646e-002 4.207962846098268e-001 ...
8.259229974584023e-001 4.207962846098268e-001 ...
-9.405920349573646e-002 -7.726317316720414e-002 ...
4.973290349094079e-002 1.193456527972926e-002 ...
-1.699063986760234e-002 -1.914286129088767e-003 ...
1.908831736481291e-003];
hd = [0 0 0 1.442628250562444e-002 -1.446750489679015e-002 ...
-7.872200106262882e-002 4.036797903033992e-002 ...
4.178491091502746e-001 -7.589077294536542e-001 ...
4.178491091502746e-001 4.036797903033992e-002 ...
-7.872200106262882e-002 -1.446750489679015e-002 ...
1.442628250562444e-002 0 0 0 0];
t = (0:17);
lr = cos(pi * (t + 1)) .* hd;
hr = cos(pi * t) .* ld;
case 'jpeg9.7'
ld = [0 0.02674875741080976 -0.01686411844287495 ...
-0.07822326652898785 0.2668641184428723 ...
0.6029490182363579 0.2668641184428723 ...
-0.07822326652898785 -0.01686411844287495 ...
0.02674875741080976];
hd = [0 -0.09127176311424948 0.05754352622849957 ...
0.5912717631142470 -1.115087052456994 ...
0.5912717631142470 0.05754352622849957 ...
-0.09127176311424948 0 0];
t = (0:9);
lr = cos(pi * (t + 1)) .* hd;
hr = cos(pi * t) .* ld;
otherwise
error('Unrecognizable wavelet name (WNAME).');
end
% Output the requested filters.
if (nargin == 1)
varargout(1:4) = {ld, hd, lr, hr};
else
switch lower(type(1))
case 'd'
varargout = {ld, hd};
case 'r'
varargout = {lr, hr};
otherwise
error('Unrecognizable filter TYPE.');
end
end
```

**Example 5**

```
function [c, s] = wavefast(x, n, varargin)
error(nargchk(3, 4, nargin));
if nargin == 3
if ischar(varargin{1})
[lp, hp] = wavefilter(varargin{1}, 'd');
else
error('Missing wavelet name.');
end
else
lp = varargin{1}; hp = varargin{2};
end
% Get the filter length, 'lp', input array size, 'sx', and number of
% pages, 'pages', in extended 2-D array x.
fl = length(lp); sx = size(x); pages = size(x, 3);
if ((ndims(x) ~= 2) && (ndims(x) ~= 3)) || (min(sx) < 2) ...
|| ~isreal(x) || ~isnumeric(x)
error('X must be a real, numeric 2-D or 3-D matrix.');
end
if (ndims(lp) ~= 2) || ~isreal(lp) || ~isnumeric(lp) ...
|| (ndims(hp) ~= 2) || ~isreal(hp) || ~isnumeric(hp) ...
|| (fl ~= length(hp)) || rem(fl, 2) ~= 0
error(['LP and HP must be even and equal length real, ' ...
'numeric filter vectors.']);
end
if ~isreal(n) || ~isnumeric(n) || (n < 1) || (n > log2(max(sx)))
error(['N must be a real scalar between 1 and ' ...
'log2(max(size((X))).']);
end
% Init the starting output data structures and initial approximation.
c = []; s = sx(1:2);
app = cell(pages, 1);
for i = 1:pages
app{i} = double(x(:, :, i));
end
% For each decomposition ...
for i = 1:n
% Extend the approximation symmetrically.
[app, keep] = symextend(app, fl, pages);
% Convolve rows with HP and downsample. Then convolve columns
% with HP and LP to get the diagonal and vertical coefficients.
rows = symconv(app, hp, 'row', fl, keep, pages);
coefs = symconv(rows, hp, 'col', fl, keep, pages);
c = addcoefs(c, coefs, pages);
s = [size(coefs{1}); s];
coefs = symconv(rows, lp, 'col', fl, keep, pages);
c = addcoefs(c, coefs, pages);
% Convolve rows with LP and downsample. Then convolve columns
% with HP and LP to get the horizontal and next approximation
% coeffcients.
rows = symconv(app, lp, 'row', fl, keep, pages);
coefs = symconv(rows, hp, 'col', fl, keep, pages);
c = addcoefs(c, coefs, pages);
app = symconv(rows, lp, 'col', fl, keep, pages);
end
% Append the final approximation structures.
c = addcoefs(c, app, pages);
s = [size(app{1}); s];
if ndims(x) > 2
s(:, 3) = size(x, 3);
end
%-------------------------------------------------------------------%
function nc = addcoefs(c, x, pages)
% Add 'pages' array coefficients to the wavelet decomposition vector.
nc = c;
for i = pages:-1:1
nc = [x{i}(:)' nc];
end
%-------------------------------------------------------------------%
function [y, keep] = symextend(x, fl, pages)
% Compute the number of coefficients to keep after convolution and
% downsampling. Then extend the 'pages' arrays of x in both
% dimensions.
y = cell(pages, 1);
for i = 1:pages
keep = floor((fl + size(x{i}) - 1) / 2);
y{i} = padarray(x{i}, [(fl - 1) (fl - 1)], 'symmetric', 'both');
end
%-------------------------------------------------------------------%
function y = symconv(x, h, type, fl, keep, pages)
% For the 'pages' 2-D arrays in x, convolve the rows or columns with
% h, downsample, and extract the center section since symmetrically
% extended.
y = cell(pages, 1);
for i = 1:pages
if strcmp(type, 'row')
y{i} = conv2(x{i}, h);
y{i} = y{i}(:, 1:2:end);
y{i} = y{i}(:, fl / 2 + 1:fl / 2 + keep(2));
else
y{i} = conv2(x{i}, h');
y{i} = y{i}(1:2:end, :);
y{i} = y{i}(fl / 2 + 1:fl / 2 + keep(1), :);
end
end
```

Example 6

```
function [varargout] = wavework(opcode, type, c, s, n, x)
error(nargchk(4, 6, nargin));
if (ndims(c) ~= 2) || (size(c, 1) ~= 1)
error('C must be a row vector.');
end
if (ndims(s) ~= 2) || ~isreal(s) || ~isnumeric(s) || ...
((size(s, 2) ~= 2) && (size(s, 2) ~= 3))
error('S must be a real, numeric two- or three-column array.');
end
elements = prod(s, 2); % Coefficient matrix elements.
if (length(c) < elements(end)) || ...
~(elements(1) + 3 * sum(elements(2:end - 1)) >= elements(end))
error(['[C S] must form a standard wavelet decomposition ' ...
'structure.']);
end
if strcmpi(opcode(1:3), 'pas') && nargin < 6
error('Not enough input arguments.');
end
if nargin < 5
n = 1; % Default level is 1.
end
nmax = size(s, 1) - 2; % Maximum levels in [C, S].
aflag = (lower(type(1)) == 'a');
if ~aflag && (n > nmax)
error('N exceeds the decompositions in [C, S].');
end
switch lower(type(1)) % Make pointers into C.
case 'a'
nindex = 1;
start = 1; stop = elements(1); ntst = nmax;
case {'h', 'v', 'd'}
switch type
case 'h', offset = 0; % Offset to details.
case 'v', offset = 1;
case 'd', offset = 2;
end
nindex = size(s, 1) - n; % Index to detail info.
start = elements(1) + 3 * sum(elements(2:nmax - n + 1)) + ...
offset * elements(nindex) + 1;
stop = start + elements(nindex) - 1;
ntst = n;
otherwise
error('TYPE must begin with "a", "h", "v", or "d".');
end
switch lower(opcode) % Do requested action.
case {'copy', 'cut'}
y = c(start:stop); nc = c;
y = reshape(y, s(nindex, :));
if strcmpi(opcode(1:3), 'cut')
nc(start: stop) = 0; varargout = {nc, y};
else
varargout = {y};
end
case 'paste'
if numel(x) ~= elements(end - ntst)
error('X is not sized for the requested paste.');
else
nc = c; nc(start:stop) = x(:); varargout = {nc};
end
otherwise
error('Unrecognized OPCODE.');
end
```

Example 7

```
function [nc, y] = wavecut(type, c, s, n)
error(nargchk(3, 4, nargin));
if nargin == 4
[nc, y] = wavework('cut', type, c, s, n);
else
[nc, y] = wavework('cut', type, c, s);
end
function y = wavecopy(type, c, s, n)
error(nargchk(3, 4, nargin));
if nargin == 4
y = wavework('copy', type, c, s, n);
else
y = wavework('copy', type, c, s);
end
function nc = wavepaste(type, c, s, n, x)
error(nargchk(5, 5, nargin))
nc = wavework('paste', type, c, s, n, x);
```

Example 8

```
function w = wavedisplay(c, s, scale, border)
error(nargchk(2, 4, nargin));
if (ndims(c) ~= 2) || (size(c, 1) ~= 1)
error('C must be a row vector.');
end
if (ndims(s) ~= 2) || ~isreal(s) || ~isnumeric(s) || ...
((size(s, 2) ~= 2) && (size(s, 2) ~= 3))
error('S must be a real, numeric two- or three-column array.');
end
elements = prod(s, 2);
if (length(c) < elements(end)) || ...
~(elements(1) + 3 * sum(elements(2:end - 1)) >= elements(end))
error(['[C S] must be a standard wavelet ' ...
'decomposition structure.']);
end
if (nargin > 2) && (~isreal(scale) || ~isnumeric(scale))
error('SCALE must be a real, numeric scalar.');
end
if (nargin > 3) && (~ischar(border))
error('BORDER must be character string.');
end
if nargin == 2
scale = 1; % Default scale.
end
if nargin < 4
border = 'absorb'; % Default border.
end
% Scale coefficients and determine pad fill.
absflag = scale < 0;
scale = abs(scale);
if scale == 0
scale = 1;
end
[cd, w] = wavecut('a', c, s); w = mat2gray(w);
cdx = max(abs(cd(:))) / scale;
if absflag
cd = mat2gray(abs(cd), [0, cdx]); fill = 0;
else
cd = mat2gray(cd, [-cdx, cdx]); fill = 0.5;
end
% Build gray image one decomposition at a time.
for i = size(s, 1) - 2:-1:1
ws = size(w);
h = wavecopy('h', cd, s, i);
pad = ws - size(h); frontporch = round(pad / 2);
h = padarray(h, frontporch, fill, 'pre');
h = padarray(h, pad - frontporch, fill, 'post');
v = wavecopy('v', cd, s, i);
pad = ws - size(v); frontporch = round(pad / 2);
v = padarray(v, frontporch, fill, 'pre');
v = padarray(v, pad - frontporch, fill, 'post');
d = wavecopy('d', cd, s, i);
pad = ws - size(d); frontporch = round(pad / 2);
d = padarray(d, frontporch, fill, 'pre');
d = padarray(d, pad - frontporch, fill, 'post');
% Add 1 pixel white border and concatenate coefficients.
switch lower(border)
case 'append'
w = padarray(w, [1 1], 1, 'post');
h = padarray(h, [1 0], 1, 'post');
v = padarray(v, [0 1], 1, 'post');
case 'absorb'
w(:, end, :) = 1; w(end, :, :) = 1;
h(end, :, :) = 1; v(:, end, :) = 1;
otherwise
error('Unrecognized BORDER parameter.');
end
w = [w h; v d];
end
% Display result. If the reconstruction is an extended 2-D array
% with 2 or more pages, display as a time sequence.
if nargout == 0
if size(s, 2) == 2
imshow(w);
else
implay(w);
end
end
```

Example 9

```
function [varargout] = waveback(c, s, varargin)
error(nargchk(3, 5, nargin));
error(nargchk(1, 2, nargout));
if (ndims(c) ~= 2) || (size(c, 1) ~= 1)
error('C must be a row vector.');
end
if (ndims(s) ~= 2) || ~isreal(s) || ~isnumeric(s) || ...
((size(s, 2) ~= 2) && (size(s, 2) ~= 3))
error('S must be a real, numeric two- or three-column array.');
end
elements = prod(s, 2);
if (length(c) < elements(end)) || ...
~(elements(1) + 3 * sum(elements(2:end - 1)) >= elements(end))
error(['[C S] must be a standard wavelet ' ...
'decomposition structure.']);
end
% Maximum levels in [C, S].
nmax = size(s, 1) - 2;
% Get third input parameter and init check flags.
wname = varargin{1}; filterchk = 0; nchk = 0;
switch nargin
case 3
if ischar(wname)
[lp, hp] = wavefilter(wname, 'r'); n = nmax;
else
error('Undefined filter.');
end
if nargout ~= 1
error('Wrong number of output arguments.');
end
case 4
if ischar(wname)
[lp, hp] = wavefilter(wname, 'r');
n = varargin{2}; nchk = 1;
else
lp = varargin{1}; hp = varargin{2};
filterchk = 1; n = nmax;
if nargout ~= 1
error('Wrong number of output arguments.');
end
end
case 5
lp = varargin{1}; hp = varargin{2}; filterchk = 1;
n = varargin{3}; nchk = 1;
otherwise
error('Improper number of input arguments.');
end
fl = length(lp);
if filterchk % Check filters.
if (ndims(lp) ~= 2) || ~isreal(lp) || ~isnumeric(lp) ...
|| (ndims(hp) ~= 2) || ~isreal(hp) || ~isnumeric(hp) ...
|| (fl ~= length(hp)) || rem(fl, 2) ~= 0
error(['LP and HP must be even and equal length real, ' ...
'numeric filter vectors.']);
end
end
if nchk && (~isnumeric(n) || ~isreal(n)) % Check scale N.
error('N must be a real numeric.');
end
if (n > nmax) || (n < 1)
error('Invalid number (N) of reconstructions requested.');
end
if (n ~= nmax) && (nargout ~= 2)
error('Not enough output arguments.');
end
nc = c; ns = s; nnmax = nmax; % Init decomposition.
for i = 1:n
% Compute a new approximation.
a = symconvup(wavecopy('a', nc, ns), lp, lp, fl, ns(3, :)) + ...
symconvup(wavecopy('h', nc, ns, nnmax), ...
hp, lp, fl, ns(3, :)) + ...
symconvup(wavecopy('v', nc, ns, nnmax), ...
lp, hp, fl, ns(3, :)) + ...
symconvup(wavecopy('d', nc, ns, nnmax), ...
hp, hp, fl, ns(3, :));
% Update decomposition.
nc = nc(4 * prod(ns(1, :)) + 1:end); nc = [a(:)' nc];
ns = ns(3:end, :); ns = [ns(1, :); ns];
nnmax = size(ns, 1) - 2;
end
% For complete reconstructions, reformat output as 2-D.
if nargout == 1
a = nc; nc = repmat(0, ns(1, :)); nc(:) = a;
end
varargout{1} = nc;
if nargout == 2
varargout{2} = ns;
end
```

```
function w = symconvup(x, f1, f2, fln, keep)
% Upsample rows and convolve columns with f1; upsample columns and
% convolve rows with f2; then extract center assuming symmetrical
% extension.
% Process each "page" (i.e., 3rd index) of an extended 2-D array
% separately; if 'x' is 2-D, size(x, 3) = 1.
% Preallocate w.
zi = fln - 1:fln + keep(1) - 2;
zj = fln - 1:fln + keep(2) - 2;
w = zeros(numel(zi), numel(zj), size(x, 3));
for i = 1:size(x, 3)
y = zeros([2 1] .* size(x(:, :, i)));
y(1:2:end, :) = x(:, :, i);
y = conv2(y, f1');
z = zeros([1 2] .* size(y)); z(:, 1:2:end) = y;
z = conv2(z, f2);
z = z(zi, zj);
w(:, :, i) = z;
end
```