/* * * Two-dimensional best basis and custom basis development. * * Copyright (C) 1991--94 Wickerhauser Consulting. All Rights Reserved. * May be freely copied for personal, noncommercial use by owners of * ``Adapted Wavelet Analysis from Theory to Software'' ISBN 1-56881-041-5 * by Mladen Victor Wickerhauser [AK Peters, Ltd., Wellesley, Mass., 1994] */ #include #include "real.h" #include "qf.h" #include "hedge.h" #include "infocost.h" #include "cd.h" #include "cdp2d.h" #include "wpt2.h" /******************************************************************* * bbwp2() * * Compute and write the best basis of 2-dimensional isotropic * wavelet packets to an output hedge, using the space-saving * periodic discrete wavelet packet analysis algorithm. * * Calling sequence and basic algorithm: * * bbwp2( GRAPH, KD, IN, X,Y, S, L, WK, H, G ): * Let XY = X*Y * Let COST = infocost( IN, 0, XY-1 ) * If Scontents[0]' is replaced with its coefficients * in the 2-dimensional wavelet packet best basis defined by the cost * function `infocost()', and the other pointers in `graph->contents[]' * are set to their corresponding wavelet packet subspaces. * * Assumptions: * 1. x>0; y>0 * 2. L >= s >= 0 * 3. in != NULL * 4. wk != NULL * 5. kd != NULL * 6. graph != NULL * * External functions called: * scdpe2(), infocost(), assert() * */ extern real bbwp2( hedge *graph, /* Preallocated struct for the best basis */ real *kd, /* Array to hold one branch of 2-D DWPAP */ const real *in, /* Input array, will not be changed */ int x, /* Length of one column of `in[]' */ int y, /* Length of one row of `in[]' */ int s, /* Current level in the DWPA tree */ int L, /* Maximum depth of the decomposition */ real *wk, /* Scratch array for transposition */ pqf *H, /* Low-pass periodic quadrature filter */ pqf *G) /* High-pass periodic quadrature filter */ { int block, i, j, xy; real *out, *k[4]; real cost, kcost; assert(x>0); assert(y>0); assert(s<=L); assert(s>=0); assert(in); assert(wk); assert(kd); assert(graph); xy = x*y; cost = infocost( in, 0, xy-1 ); if( sblocks; for( j=0; j<4; j++ ) k[j] = kd + j*xy/4; kd = k[3] + xy/4; scdpe2( k[0], k[1], k[2], k[3], in, x,y, wk, H, G); kcost = 0; for( j=0; j<4; j++ ) kcost += bbwp2( graph, kd, k[j], x/2, y/2, s+1, L, wk, H, G); if( kcostlevels[block] = s; out = (real *)graph->contents[block]; for( i=0; iblocks = block+1; graph->contents[graph->blocks] = out + xy; } } else { graph->levels[graph->blocks] = L; out = graph->contents[graph->blocks]; for( i=0; iblocks += 1; graph->contents[graph->blocks] = out + xy; } return(cost); } /******************************************************************* * cbwp2() * * Compute and write a custom basis of 2-dimensional isotropic * wavelet packets to an output hedge, doing all development * into periodic discrete wavelet packets in place. * * Calling sequence and basic algorithm: * * cbwp2( GRAPH, LEVEL, IX, IY, WORK, H, G ): * If LEVEL < GRAPH.LEVELS[GRAPH.BLOCKS] then * scdpi2( GRAPH.CONTENTS[GRAPH.BLOCKS], IX,IY, WORK, H, G ) * For K = 0 to 3 * cbwp2( GRAPH, LEVEL+1, IX/2, IY/2, WORK, H, G ) * Else * Let GRAPH.CONTENTS[GRAPH.BLOCKS+1] = * GRAPH.CONTENTS[GRAPH.BLOCKS] + IX*IY * GRAPH.BLOCKS += 1 * Return * * * Inputs: * * (hedge *)graph This preallocated HEDGE must contain enough * space for the desired graph basis, up to * `1<<(2*L)' elements in its * `contents' and `levels' arrays. * * (int)level This small integer is the current level. * * (int)ix These positive integers are the numbers of * (int)iy rows and columns of the input. * * (real *)work This is a scratch array for transposition. * * (pqf *)H Use these low-pass and high-pass periodized * (pqf *)G filters for convolution-decimation. * * Output: * The array at `graph->contents[0]' is replaced with its coefficients * in the 2-dimensional wavelet packet basis defined by the encounter * list `graph->levels[]', and the other pointers in `graph->contents[]' * are set to their corresponding wavelet packet subspaces. * * Assumptions: * 1. ix>0; iy>0 * 2. level >= 0 * 3. graph != NULL * 4. work != NULL * * External functions called: * assert(), scdpi2() * */ extern void cbwp2( hedge *graph, /* Preallocated struct for the best basis */ int level, /* Current level in the DWPA tree */ int ix, /* Length of one column of `in[]' */ int iy, /* Length of one row of `in[]' */ real *work, /* Scratch array for transposition */ pqf *H, /* Low-pass periodic quadrature filter */ pqf *G) /* High-pass periodic quadrature filter */ { real *out; int k; assert(ix>0); assert(iy>0); assert(level>=0); assert(work); assert(graph); if( levellevels[graph->blocks] ) { out = (real *)graph->contents[graph->blocks]; scdpi2( out, ix,iy, work, H, G ); for( k=0; k<4; k++ ) cbwp2( graph, level+1, ix/2, iy/2, work, H, G ); } else { graph->contents[graph->blocks+1] = (real *)graph->contents[graph->blocks] + ix*iy; graph->blocks += 1; } return; }