Toolbox Wavelets: perform_spiht_coding(x,options) - File Exchange

Code covered by the BSD License

Highlights from
Toolbox Wavelets

Toolox Wavelets - A Toolbox f...
apply_tensorial(x,func_name,f... apply_tensorial - Multimensional application of a 1D function.
clamp(x,a,b)
compute_best_threshold(type, ... compute_best_threshold - oracle estimation of the threshold for denoising
compute_gabor_features(M,opti... compute_gabor_features - compute a set of gabor filtered features
compute_haar_matrix(w,Jmin, o... compute_haar_matrix - compute the haar matrix for a 2D haar transform.
compute_quadrant_selection(j,... compute_quadrant_selection - compute the indices for selecting subband of a wavelet transform.
compute_spline_filter(s,N) compute_spline_filter - compute the spline filter of a wavelete transform
compute_symmetric_conditional... compute_symmetric_conditional_histogram - compute a symmetric histogram
convert_wavelet2vect(x, Jmin,... convert_wavelet2vect - turn wavelet coefficients into a vector
convert_wavelets2list(x, Jmin... convert_wavelets2list - convert a wavelet transform into a cell array.
crop(M,n,c) crop - crop an image to reduce its size
cwpt2_interface(in, direction... function out = cwpt2_interface(in, direction, filter_type, trans_type, btree)
get_lifting_param(type) get_lifting_param - get the parameters of the lifting
getoptions(options, name, v, ... getoptions - retrieve options parameter
histo(mtx, nbins, binCtr)
imageplot(M,str, a,b,c) imageplot - diplay an image and a title
imden=den2d_block(x,qmf,thd,L... wc=FWT2_PO(x,scale,qmf);
jp2_codec(in_data, compressed... this function calls via the compiled function jp2_codec.mexglx
jp2_degrade(Isrc, comp_ratio,... Original jp2k compression/decompression function
jp2codec_degrade(Isrc, comp_r... % Attention: a) no rounding operation; it should be performed beforehand
keep_biggest(x,n) keep_biggest - keep only the n biggest coef,
load_image(type, n, options) load_image - load benchmark images.
load_signal(name, n, options) load_signal - load a 1D signal
perform_79_transform_slow(x, ... perform_79_transform_slow - biorthogonal 7/9 wavelet transform
perform_atrou_transform(x,Jmi... perform_atrou_transform - compute the "a trou" wavelet transform,
perform_blsgsm_denoising(x, o... perform_blsgsm_denoising - denoise an image using BLS-GSM
perform_circle_haar_transform... perform_circle_haar_transform - perform a circle-valued Haar transform.
perform_cpx_dualtree_transfor... perform_cpx_dualtree_transform - perform a complex valued dual-tree wavelet transform
perform_curvelet_transform(x,... perform_curvelet_transform - a wrapper to curvlab
perform_haar_transform_1d(x,d... perform_haar_transform - perform the 1D haar transform along columns
perform_haar_transform_slow(x... perform_haar_transform - compute a wavelet haar transform
perform_histogram_equalizatio... perform_histogram_equalization - perform histogram equalization
perform_histogram_matching_wa... perform_histogram_matching_wavelet - match the histogram of a wavelet transform
perform_image_extension(M,n, ... perform_image_extension - extend the size by symetry
perform_jbig_coding(x) perform_jbig_coding - perform binary image coding
perform_jp2_rescaling(MW,Jmin...
perform_jp2k_compression(x,op... perform_jp2k_compression - compress/decompress an image using JPEG2000
perform_jp2k_degradation(MW,J... perform_jp2k_degradation - perform coding/decoding of wavelet coefficients.
perform_lifting_transform(x, ... perform_lifting_transform - peform fast lifting transform
perform_lifting_transform_byn... perform_lifting_transform_byname - wavelet transform via lifting
perform_lifting_transform_slo... perform_lifting_transform_slow - perform a lifting transform without relying on Mex file
perform_pyramid_transform(x,g... perform_pyramid_transform - fast pyramidal transform (1D or 2D).
perform_pyramid_transform_non... perform_pyramid_transform_nonframe - fast pyramidal transform (1D or 2D).
perform_pyramid_transform_sim... perform_pyramid_transform_simoncelli - Laplacian pyramidal transform
perform_pyramid_transform_ti(... perform_pyramid_transform_ti - gaussian pyramidal transform (1D or 2D).
perform_quantization(x, T, di... perform_quantization - perform a nearly uniform quantization of the signal.
perform_quicunx_wavelet_trans... perform_quicunx_wavelet_transform_ti - translation invariant quincunx wavelets
perform_quincunx_wavelet_tran... perform_quincunx_wavelet_transform - compute quincunx transform
perform_real_dualtree_transfo... perform_real_dualtree_transform - perform a real valued dual-tree wavelet transform
perform_segmentation(E,option... perform_segmentation - perform image segmentation
perform_segmentation_coloriza... perform_segmentation_colorization - return a color image
perform_shannon_estimation_wa... evaluate_nbr_bits_wavelets - evaluate the number of bits to code a wavelet transform
perform_spaciowav_matching(M_... perform_spaciowav_matching - match both the spacial and wavelet histogram
perform_spiht_coding(x,option... perform_spiht_coding - SPIHT coding of wavelet coefficients
perform_steerable_matching(M1... perform_steerable_matching - match multiscale histograms
perform_steerable_transform(x... perform_steerable_transform - steerable pyramidal transform
perform_thresholding(x, t, ty... perform_thresholding - perform hard or soft thresholding
perform_waveatoms_transform(x... perform_waveatoms_transform - interface to WaveAtom transform
perform_wavelet_arithmetic_co... perform_wavelet_arithmetic_coding - arithmetic code a wavelet transformed image.
perform_wavelet_jp2k_coding(x... perform_wavelet_jp2k_coding - code a wavelet transformed image using JPEG2000.
perform_wavelet_matching(M1,M... perform_wavelet_matching - match multiscale histograms
perform_wavelet_transform(x, ... perform_wavelet_transform - wrapper to several wavelet transform.
perform_wavelet_transform_hyp... fwt_hyperbolic - multidimensional hyperbolic (i.e. fully tensorial) wavelet transform
perform_wavelet_transform_irr... perform_irregular_lifting - wavelet for irregular sampling
perform_wavelet_transform_iso... perform_wavelet_transform_isotropic - multidimensional isotropic wavelet transform (aka classical ND transform).
plot_quincunx_wavelet(MW, Jmi... plot_wavelet - plot 2D wavelet transform stored using Mallat's ordering.
plot_wavelet(M, Jmin, options... plot_wavelet - plot 2D wavelet transform stored using Mallat's ordering.
psnr(x,y, vmax) psnr - compute the Peack Signal to Noise Ratio
publish_html(filename, output... publish_html - publish a file to HTML format
reorder_coefs(x, Jmin ,dir) reorder_coefs - reorder the wavelet coefficients.
rescale(x,a,b)
compile_jp2k.m
compile_mex.m
perform_79_transform.m
perform_haar_transform.m
perform_lifting_transform.m
test_atrou.m
test_block_thresholding.m
test_circle_haar.m
test_classical.m
test_compression.m
test_denoising.m
test_denoising_blsgsm.m
test_dualtree.m
test_haar_transform_1d.m
test_j2k.m
test_laplacian_arithmetic_cod...
test_laplacian_do.m
test_laplacian_pyramid.m
test_lifting.m
test_lifting_irregular.m
test_mex.m
test_openjpeg.m
test_pyramidal.m
test_quincunx.m
test_quincunx_ti.m
test_reorder_coefs.m
test_segmentation_gabor.m
test_statistics.m
test_statistics_batch.m
test_steerable.m
test_texture_synthesis_wavele...
test_wavelet.m
test_wavelets_ti.m
View all files

from Toolbox Wavelets by Gabriel Peyre
Wavelet transform and coding functions, including other more exotic transforms (laplacian, steerable

perform_spiht_coding(x,options)

function [y,nbr_bits] = perform_spiht_coding(x,options)

% perform_spiht_coding - SPIHT coding of wavelet coefficients
%
% Coding : 
%   options.Jmin = ??;      % minimum scale of the transform
%   options.nb_bits = ??;   % target number of bits
%   [stream,nbr_bits] = perform_spiht_coding(MW,options);
% Decoding : 
%   MW = perform_spiht_coding(stream);
%
% This is a simple wrapper of the code of Jing Tian
% <scuteejtian at hotmail.com>

options.null = 0;

if isfield(options, 'arithmetic_coding')
    arithmetic_coding = options.arithmetic_coding;
else
    arithmetic_coding = 1;
end

if size(x,1)>1 && size(x,2)>1
    if isfield(options, 'Jmin')
        Jmin = options.Jmin;
    else
        warning('You should provide options.Jmin');
        Jmin = 4;
    end
    if isfield(options, 'nb_bits')
        nb_bits = options.nb_bits;
    else
        nb_bits = floor( prod(size(x))*0.1 );
    end
    %-----------   Coding   ----------------
    Jmax = log2(size(x,1))-1;
    level = Jmax-Jmin+1;
    y = func_SPIHT_Enc(x, nb_bits, level); y = y(:);
    if arithmetic_coding
        % remove trailing 2
        I = find(y==2); y(I) = [];
        % the 3 first entry are [size nb_bitsplanes level]
        [z,nbr_bits] = perform_arithmetic_coding(y(4:end), +1, options);
        y = [y(1:3); z];
        nbr_bits = nbr_bits + 16; % approx 16 bits
    else
        nbr_bits = nb_bits;
    end
else
    %-----------   Decoding   ----------------
    if arithmetic_coding
        x = x(:);
        % the 3 first entry are [size nb_bitsplanes level]
        z = perform_arithmetic_coding(x(4:end), -1, options);
        x = [x(1:3); z];
        nbr_bits = -1;
    end
    y = func_SPIHT_Dec(x(:)');
end

function out = func_SPIHT_Enc(m, max_bits, level)
% Matlab implementation of SPIHT (without Arithmatic coding stage)
%
% Encoder
%
% input:    m : input image in wavelet domain
%           max_bits : maximum bits can be used
%           level : wavelet decomposition level
%
% output:   out : bit stream
%
% Jing Tian
% Contact me : scuteejtian@hotmail.com
% This program is part of my undergraduate project in GuangZhou, P. R. China.
% April - July 1999


%-----------   Initialization  -----------------
bitctr = 0;
out = 2*ones(1,max_bits - 14);
n_max = floor(log2(abs(max(max(m)'))));
Bits_Header = 0;
Bits_LSP = 0;
Bits_LIP = 0;
Bits_LIS = 0;

%-----------   output bit stream header   ----------------
% image size, number of bit plane, wavelet decomposition level should be
% written as bit stream header.
out(1,[1 2 3]) = [size(m,1) n_max level]; bitctr = bitctr + 24;
index = 4;
Bits_Header = Bits_Header + 24;

%-----------   Initialize LIP, LSP, LIS   ----------------
temp = [];
bandsize = 2.^(log2(size(m, 1)) - level + 1);
temp1 = 1 : bandsize;
for i = 1 : bandsize
    temp = [temp; temp1];
end
LIP(:, 1) = temp(:);
temp = temp';
LIP(:, 2) = temp(:);
LIS(:, 1) = LIP(:, 1);
LIS(:, 2) = LIP(:, 2);
LIS(:, 3) = zeros(length(LIP(:, 1)), 1);
pstart = 1;
pend = bandsize / 2;
for i = 1 : bandsize / 2
    LIS(pstart : pend, :) = [];
    pdel = pend - pstart + 1;
    pstart = pstart + bandsize - pdel;
    pend = pend + bandsize - pdel;
end
LSP = [];

n = n_max;

%-----------   coding   ----------------
while(bitctr < max_bits)
        
    % Sorting Pass
    LIPtemp = LIP; temp = 0;
    for i = 1:size(LIPtemp,1)
        temp = temp+1;
        if (bitctr + 1) >= max_bits
            if (bitctr < max_bits)
                out(length(out))=[];
            end
            return
        end
        if abs(m(LIPtemp(i,1),LIPtemp(i,2))) >= 2^n % 1: positive; 0: negative
            out(index) = 1; bitctr = bitctr + 1;
            index = index +1; Bits_LIP = Bits_LIP + 1;
            sgn = m(LIPtemp(i,1),LIPtemp(i,2))>=0;
            out(index) = sgn; bitctr = bitctr + 1;
            index = index +1; Bits_LIP = Bits_LIP + 1;
            LSP = [LSP; LIPtemp(i,:)];
            LIP(temp,:) = []; temp = temp - 1;
        else
            out(index) = 0; bitctr = bitctr + 1;
            index = index +1;
            Bits_LIP = Bits_LIP + 1;
        end
    end
    
    LIStemp = LIS; temp = 0; i = 1;
    while ( i <= size(LIStemp,1))
        temp = temp + 1;
        if LIStemp(i,3) == 0
            if bitctr >= max_bits
                return
            end
            max_d = func_MyDescendant(LIStemp(i,1),LIStemp(i,2),LIStemp(i,3),m);
            if max_d >= 2^n
                out(index) = 1; bitctr = bitctr + 1;
                index = index +1; Bits_LIS = Bits_LIS + 1;
                x = LIStemp(i,1); y = LIStemp(i,2);
                
                if (bitctr + 1) >= max_bits
                    if (bitctr < max_bits)
                        out(length(out))=[];
                    end
                    return
                end
                if abs(m(2*x-1,2*y-1)) >= 2^n
                    LSP = [LSP; 2*x-1 2*y-1];
                    out(index) = 1; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    sgn = m(2*x-1,2*y-1)>=0;
                    out(index) = sgn; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                else
                    out(index) = 0; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    LIP = [LIP; 2*x-1 2*y-1];
                end
                
                if (bitctr + 1) >= max_bits
                    if (bitctr < max_bits)
                        out(length(out))=[];
                    end
                    return
                end
                if abs(m(2*x-1,2*y)) >= 2^n
                    LSP = [LSP; 2*x-1 2*y];
                    out(index) = 1; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    sgn = m(2*x-1,2*y)>=0;
                    out(index) = sgn; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                else
                    out(index) = 0; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    LIP = [LIP; 2*x-1 2*y];
                end
                
                if (bitctr + 1) >= max_bits
                    if (bitctr < max_bits)
                        out(length(out))=[];
                    end
                    return
                end
                if abs(m(2*x,2*y-1)) >= 2^n
                    LSP = [LSP; 2*x 2*y-1];
                    out(index) = 1; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    sgn = m(2*x,2*y-1)>=0;
                    out(index) = sgn; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                else
                    out(index) = 0; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    LIP = [LIP; 2*x 2*y-1];
                end
                
                if (bitctr + 1) >= max_bits
                    if (bitctr < max_bits)
                        out(length(out))=[];
                    end
                    return
                end
                if abs(m(2*x,2*y)) >= 2^n
                    LSP = [LSP; 2*x 2*y];
                    out(index) = 1; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    sgn = m(2*x,2*y)>=0;
                    out(index) = sgn; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                else
                    out(index) = 0; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    LIP = [LIP; 2*x 2*y];
                end
                
                if ((2*(2*x)-1) < size(m) & (2*(2*y)-1) < size(m))
                    LIS = [LIS; LIStemp(i,1) LIStemp(i,2) 1];
                    LIStemp = [LIStemp; LIStemp(i,1) LIStemp(i,2) 1];
                end
                LIS(temp,:) = []; temp = temp-1;
                
            else
                out(index) = 0; bitctr = bitctr + 1;
                index = index +1; Bits_LIS = Bits_LIS + 1;
            end
        else
            if bitctr >= max_bits
                return
            end
            max_d = func_MyDescendant(LIStemp(i,1),LIStemp(i,2),LIStemp(i,3),m);
            if max_d >= 2^n
                out(index) = 1; bitctr = bitctr + 1;
                index = index +1;
                x = LIStemp(i,1); y = LIStemp(i,2);
                LIS = [LIS; 2*x-1 2*y-1 0; 2*x-1 2*y 0; 2*x 2*y-1 0; 2*x 2*y 0];
                LIStemp = [LIStemp; 2*x-1 2*y-1 0; 2*x-1 2*y 0; 2*x 2*y-1 0; 2*x 2*y 0];
                LIS(temp,:) = []; temp = temp - 1;
            else
                out(index) = 0; bitctr = bitctr + 1;
                index = index +1; Bits_LIS = Bits_LIS + 1;
            end
        end
        i = i+1;
    end
    
    % Refinement Pass
    temp = 1;
    value = floor(abs(2^(n_max-n+1)*m(LSP(temp,1),LSP(temp,2))));
    while (value >= 2^(n_max+2) & (temp <= size(LSP,1)))
        if bitctr >= max_bits
            return
        end
        s = bitget(value,n_max+2);
        out(index) = s; bitctr = bitctr + 1;
        index = index +1; Bits_LSP = Bits_LSP + 1;
        temp = temp + 1;
        if temp <= size(LSP,1)
            value = floor(abs(2^(n_max-n+1)*m(LSP(temp,1),LSP(temp,2))));
        end
    end
    
    n = n - 1;
end




function m = func_SPIHT_Dec(in)
% Matlab implementation of SPIHT (without Arithmatic coding stage)
%
% Decoder
%
% input:    in : bit stream
%
% output:   m : reconstructed image in wavelet domain
%
% Jing Tian
% Contact me : scuteejtian@hotmail.com
% This program is part of my undergraduate project in GuangZhou, P. R. China.
% April - July 1999

%-----------   Initialization  -----------------
% image size, number of bit plane, wavelet decomposition level should be
% written as bit stream header.
m = zeros(in(1,1));
n_max = in(1,2);
level = in(1,3);
ctr = 4;
 
%-----------   Initialize LIP, LSP, LIS   ----------------
temp = [];
bandsize = 2.^(log2(in(1,1)) - level + 1);
temp1 = 1 : bandsize;
for i = 1 : bandsize
    temp = [temp; temp1];
end
LIP(:, 1) = temp(:);
temp = temp';
LIP(:, 2) = temp(:);

LIS(:, 1) = LIP(:, 1);
LIS(:, 2) = LIP(:, 2);
LIS(:, 3) = zeros(length(LIP(:, 1)), 1);
pstart = 1;
pend = bandsize / 2;
for i = 1 : bandsize / 2
    LIS(pstart : pend, :) = [];
    pdel = pend - pstart + 1;
    pstart = pstart + bandsize - pdel;
    pend = pend + bandsize - pdel;
end
LSP = [];

%-----------   coding   ----------------
n = n_max;
while (ctr <= size(in,2))
    
    %Sorting Pass
    LIPtemp = LIP; temp = 0;
    for i = 1:size(LIPtemp,1)
        temp = temp+1;
        if ctr > size(in,2)
            return
        end
        if in(1,ctr) == 1
            ctr = ctr + 1;
            if in(1,ctr) > 0
                m(LIPtemp(i,1),LIPtemp(i,2)) = 2^n + 2^(n-1);  
            else
                m(LIPtemp(i,1),LIPtemp(i,2)) = -2^n  - 2^(n-1); 
            end
            LSP = [LSP; LIPtemp(i,:)];
            LIP(temp,:) = []; temp = temp - 1;
        end
        ctr = ctr + 1;
    end
    
    LIStemp = LIS; temp = 0; i = 1;
    while ( i <= size(LIStemp,1))
        temp = temp + 1;
        if ctr > size(in,2)
            return
        end
        if LIStemp(i,3) == 0
            if in(1,ctr) == 1 
                ctr = ctr + 1;
                x = LIStemp(i,1); y = LIStemp(i,2);
                
                if ctr > size(in,2)
                    return
                end
                if in(1,ctr) == 1
                    LSP = [LSP; 2*x-1 2*y-1];
                    ctr = ctr + 1;
                    if in(1,ctr) == 1
                        m(2*x-1,2*y-1) = 2^n + 2^(n-1); 
                    else
                        m(2*x-1,2*y-1) = -2^n  - 2^(n-1); 
                    end
                    ctr = ctr + 1;
                else
                    LIP = [LIP; 2*x-1 2*y-1];
                    ctr = ctr + 1;
                end
                
                if ctr > size(in,2)
                    return
                end
                if in(1,ctr) == 1
                    ctr = ctr + 1;
                    LSP = [LSP; 2*x-1 2*y];
                    if in(1,ctr) == 1;
                        m(2*x-1,2*y) = 2^n + 2^(n-1); 
                    else
                        m(2*x-1,2*y) = -2^n  - 2^(n-1); 
                    end
                    ctr = ctr + 1;
                else
                    LIP = [LIP; 2*x-1 2*y];
                    ctr = ctr + 1;
                end
                
                if ctr > size(in,2)
                    return
                end
                if in(1,ctr) == 1
                    ctr = ctr + 1;
                    LSP = [LSP; 2*x 2*y-1];
                    if in(1,ctr) == 1
                        m(2*x,2*y-1) = 2^n + 2^(n-1); 
                    else
                        m(2*x,2*y-1) = -2^n  - 2^(n-1);
                    end
                    ctr = ctr + 1;
                else
                    LIP = [LIP; 2*x 2*y-1];
                    ctr = ctr + 1;
                end
                
                if ctr > size(in,2)
                    return
                end
                if in(1,ctr) == 1
                    ctr = ctr + 1;
                    LSP = [LSP; 2*x 2*y];
                    if in(1,ctr) == 1
                        m(2*x,2*y) = 2^n + 2^(n-1); 
                    else
                        m(2*x,2*y) = -2^n  - 2^(n-1); 
                    end
                    ctr = ctr + 1;
                else
                    LIP = [LIP; 2*x 2*y];
                    ctr = ctr + 1;
                end
                
                if ((2*(2*x)-1) < size(m) & (2*(2*y)-1) < size(m))
                    LIS = [LIS; LIStemp(i,1) LIStemp(i,2) 1];
                    LIStemp = [LIStemp; LIStemp(i,1) LIStemp(i,2) 1];
                end
                LIS(temp,:) = []; temp = temp-1;
                
            else
                ctr = ctr + 1;
            end
        else
            if in(1,ctr) == 1
                x = LIStemp(i,1); y = LIStemp(i,2);
                LIS = [LIS; 2*x-1 2*y-1 0; 2*x-1 2*y 0; 2*x 2*y-1 0; 2*x 2*y 0];
                LIStemp = [LIStemp; 2*x-1 2*y-1 0; 2*x-1 2*y 0; 2*x 2*y-1 0; 2*x 2*y 0];
                LIS(temp,:) = []; temp = temp - 1;
            end
            ctr = ctr + 1;
        end
        i = i+1;
    end
    
    % Refinement Pass
    temp = 1;
    value = m(LSP(temp,1), LSP(temp,2));
    while (abs(value) >= 2^(n+1) & (temp <= size(LSP,1)))
        if ctr > size(in,2)
            return
        end

        value = value + ((-1)^(in(1,ctr) + 1)) * (2^(n-1))*sign(m(LSP(temp,1),LSP(temp,2))); 
        m(LSP(temp,1),LSP(temp,2)) = value;
        ctr = ctr + 1;
        temp = temp + 1;    
        if temp <= size(LSP,1)
            value = m(LSP(temp,1),LSP(temp,2));
        end
    end
    
    n = n-1;
end




function value = func_MyDescendant(i, j, type, m)
% Matlab implementation of SPIHT (without Arithmatic coding stage)
%
% Find the descendant with largest absolute value of pixel (i,j)
%
% input:    i : row coordinate
%           j : column coordinate
%           type : type of descendant
%           m : whole image
%
% output:   value : largest absolute value
%
% Jing Tian
% Contact me : scuteejtian@hotmail.com
% This program is part of my undergraduate project in GuangZhou, P. R. China.
% April - July 1999

s = size(m,1);

S = [];

index = 0; a = 0; b = 0;

while ((2*i-1)<s & (2*j-1)<s)
    a = i-1; b = j-1;

    mind = [2*(a+1)-1:2*(a+2^index)];
    nind = [2*(b+1)-1:2*(b+2^index)];
    
    
    chk = mind <= s;
    len = sum(chk);
    if len < length(mind)
        mind(len+1:length(mind)) = [];
    end
    
    
    chk = nind <= s;
    len = sum(chk);
    if len < length(nind)
        nind(len+1:length(nind)) = [];
    end
    
    S = [S reshape(m(mind,nind),1,[])];
    
    index = index + 1;
    i = 2*a+1; j = 2*b+1;
end

if type == 1
    S(:,1:4) = [];; 
end

value = max(abs(S));

Highlights from Toolbox Wavelets

Highlights from
Toolbox Wavelets