/* * * Separable 2-D convolution-decimation routines. * * 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 "cd.h" #include "cdp2d.h" /********************************************************************** * sacdpd2() * * [S]eparable [A]djoint [C]onvolution-[D]ecimation; * [P]eriodic, [D]isjoint, [2]-dimensional. * * This function inverts the separable 2-dimensional periodic * convolution-decimation performed by `scdpd2()' or `scdpi2()', * using the quadrature mirror filters applied by the function * `acdpo()'. It superposes the results into the array `out[]'. * It takes 4 input arrays as arguments. It assumes that the * input and output are disjoint. * * Note: if any of the 4 output arrays are NULL, then they are * not included in the reconstruction. This provides a means of * cheaply including totally zero descendents. * * Calling sequence and basic algorithm: * * sacdpd2( OUT, IN0,IN1,IN2,IN3, IX,IY, WORK, HQF,GQF ): * Let OY = 2*IY * For I=0 to IX-1 * acdpe( WORK+I, IX, IN0+I*IY, IY, HQF ) * acdpo( WORK+I, IX, IN1+I*IY, IY, GQF ) * For I=0 to OY-1 * acdpo( OUT+I, OY, WORK+I*IX, IX, HQF ) * For I=0 to IX-1 * acdpe( WORK+I, IX, IN2+I*IY, IY, HQF ) * acdpo( WORK+I, IX, IN3+I*IY, IY, GQF ) * For I=0 to OY-1 * acdpo( OUT+I, OY, WORK+I*IX, IX, GQF ) * * * Inputs: * (real *)out This array must be preallocated with at * least `4*iy*ix' elements. * It will be overwritten even before all * the arrays `in0[]...in3[]' have been read. * * (const real *)in0 These arrays must be preallocated and defined * (const real *)in1 with at least `iy*ix' locations each, * (const real *)in2 or else be NULL. They are not changed by * (const real *)in3 this function. * * (int)ix These positive integers are the number of rows * (int)iy and columns in the arrays `in0[]...in3[]' * * (real *)work This temporary array must be preallocated * with at least `2*ix*iy' elements. * It will be overwritten even before all the * arrays `in0[]...in3[],' have been read. * * (const pqf *)HQF These point to structs containing filter * (const pqf *)GQF coefficients for adjoint convolution- * decimation, used in reconstruction. * * Outputs: * (real *)out This array is filled by side effect. * * (real *)work This array is trashed by side effect. * * External functions called: * acdpo(), acdpe(), assert * * Assumptions: * 1. ix>0 * 2. iy>0 * 3. HQF != NULL * 4. GQF != NULL * 5. `out[0,...,4*ix*iy-1]' is disjoint from each of * `in0[0,...,ix*iy-1]',...,`in3[0,...,ix*iy-1]'. * 6. work!=NULL * 7. out!=NULL * 8. in0!=NULL || in1!=NULL || in2!=NULL || in3!=NULL */ extern void sacdpd2( real *out, /* Pointer to a preallocated output array */ const real *in0, /* Pointer to the (HxHy) input array or NULL */ const real *in1, /* Pointer to the (HxGy) input array or NULL */ const real *in2, /* Pointer to the (GxHy) input array or NULL */ const real *in3, /* Pointer to the (GxGy) input array or NULL */ int ix, /* Rows in each of `in0[]...in3[]' */ int iy, /* Columns in each of `in0[]...in3[]' */ real *work, /* Scratch array of size `2*ix*iy' */ const pqf *HQF, /* 1-D Adjoint low-pass quadrature filter */ const pqf *GQF) /* 1-D Adjoint high-pass quadrature filter */ { int i, oy; real *wptr; assert(ix>0); assert(iy>0); assert(HQF); assert(GQF); assert(work); assert(out); assert( in0 || in1 || in2 || in3 ); oy = 2*iy; if( in0 || in1 ) { /* Y-direction 1: Loop over the rows of `in0[]' and/or `in1[]': */ if( in0 ) { for( i=0; i0 * 2. iy>0 * 3. HQF != NULL * 4. GQF != NULL * 5. `out[0,...,4*ix*iy-1]' is disjoint from each of * `in0[0,...,ix*iy-1]',...,`in3[0,...,ix*iy-1]'. * 6. work!=NULL * 7. out!=NULL * 8. in0!=NULL || in1!=NULL || in2!=NULL || in3!=NULL */ extern void sacdpe2( real *out, /* Pointer to a preallocated output array */ const real *in0, /* Pointer to the (HxHy) input array or NULL */ const real *in1, /* Pointer to the (HxGy) input array or NULL */ const real *in2, /* Pointer to the (GxHy) input array or NULL */ const real *in3, /* Pointer to the (GxGy) input array or NULL */ int ix, /* Rows in each of `in0[]...in3[]' */ int iy, /* Columns in each of `in0[]...in3[]' */ real *work, /* Scratch array of size `2*ix*iy' */ const pqf *HQF, /* 1-D Adjoint low-pass quadrature filter */ const pqf *GQF) /* 1-D Adjoint high-pass quadrature filter */ { int i, oy; real *wptr; assert(ix>0); assert(iy>0); assert(HQF); assert(GQF); assert(work); assert(out); assert( in0 || in1 || in2 || in3 ); oy = 2*iy; if( in0 || in1 ) { /* Y-direction 1: Loop over the rows of `in0[]' and/or `in1[]': */ if( in0 ) { for( i=0; i0); assert(iy>0); assert(work); assert(HQF); assert(GQF); oy = 2*iy; n = iy*ix; work1 = work; work2 = work1 + 2*n; in0 = data; in1 = in0 + n; in2 = in1 + n; in3 = in2 + n; for( i=0; i0); assert( !(ix&1) ); assert(iy>0); assert( !(iy&1) ); oy = iy/2; /* First step (`out0[]' and `out2[]'): */ if( out0 || out2 ) { /* Y-direction 1: Low-pass filter the rows of `in[]': */ iptr = in; for( i=0; i0); assert( !(ix&1) ); assert(iy>0); assert( !(iy&1) ); oy = iy/2; /* First step (`out0[]' and `out2[]'): */ if( out0 || out2 ) { /* Y-direction 1: Low-pass filter the rows of `in[]': */ iptr = in; for( i=0; i0); assert(ix>0); assert( !(iy&1) ); assert( !(ix&1) ); assert( data ); assert( work ); assert(HQF); assert(GQF); oy = iy/2; n = oy*ix/2; work1 = work; work2 = work + oy*ix; for( i=0; i