tarkin/dwt.c

   1 // $Id: dwt.c,v 1.1 2001/02/13 01:06:24 giles Exp $
   2 //
   3 // $Log: dwt.c,v $
   4 // Revision 1.1  2001/02/13 01:06:24  giles
   5 // Initial revision
   6 //
   7
   8 #include <math.h>
   9 #include "tarkin.h"
  10
  11 void daub4(float a[], unsigned long n, int isign)
  12 {
  13   float *wksp;
  14   unsigned long nh, nh1, i, j;
  15
  16   if(n < 4) return;
  17   wksp = (float *)malloc(sizeof(float) * n);
  18   nh1 = (nh=n>>1)+1;
  19   if(isign>=0) {
  20     for(i=1,j=1;j<=n-3;j+=2,i++) {
  21       wksp[i]=C0*a[j]+C1*a[j+1]+C2*a[j+2]+C3*a[j+3];
  22       wksp[i+nh]=C3*a[j]-C2*a[j+1]+C1*a[j+2]-C0*a[j+3];
  23     }
  24     wksp[i]=C0*a[n-1]+C1*a[n]+C2*a[1]+C3*a[2];
  25     wksp[i+nh]=C3*a[n-1]-C2*a[n]+C1*a[1]-C0*a[2];
  26   } else {
  27     wksp[1]=C2*a[nh]+C1*a[n]+C0*a[1]+C3*a[nh1];
  28     wksp[2]=C3*a[nh]-C0*a[n]+C1*a[1]-C2*a[nh1];
  29     for(i=1,j=3;i<nh;i++) {
  30       wksp[j++]=C2*a[i]+C1*a[i+nh]+C0*a[i+1]+C3*a[i+nh1];
  31       wksp[j++]=C3*a[i]-C0*a[i+nh]+C1*a[i+1]-C2*a[i+nh1];
  32     }
  33   }
  34   for(i=1;i<=n;i++) a[i]=wksp[i];
  35   free(wksp);
  36 }
  37
  38 void wtn(float a[], unsigned long nn[], int ndim, int isign)
  39 {
  40   unsigned long i1,i2,i3,k,n,nnew,nprev=1,nt,ntot=1;
  41   int idim;
  42   float *wksp;
  43
  44   for(idim=1;idim<=ndim;idim++) ntot *= nn[idim];
  45   wksp=(float *)malloc(sizeof(float)*ntot);
  46
  47   for(idim=1;idim<=ndim;idim++) {
  48     n = nn[idim];
  49     nnew = n * nprev;
  50     if(n > 4) {
  51       for(i2=0;i2<ntot;i2+=nnew) for(i1=1;i1<=nprev;i1++) {
  52         for(i3=i1+i2,k=1;k<=n;k++,i3+=nprev) wksp[k] = a[i3];
  53         if(isign >= 0) for(nt=n;nt>=4;nt>>=1) daub4(wksp,nt,isign);
  54         else for(nt=4;nt<=n;nt<<=1) daub4(wksp,nt,isign);
  55         for(i3=i1+i2,k=1;k<=n;k++,i3+=nprev) a[i3]=wksp[k];
  56       }
  57     }
  58     nprev = nnew;
  59   }
  60   free(wksp);
  61 }
  62
  63 // Monodimensional DWT/iDWT transform
  64 // if isign is positive - Forward DWT
  65 // if isign is negative - Inverse DWT
  66 //
  67 // Originally snagged from Numeric Recipes in C (2nd Ed)
  68 // Modified to get around braindeadedness with arrays starting at [1].
  69
  70 void zdaub4(float a[], unsigned long n, int isign)
  71 {
  72   float *wksp;
  73   unsigned long nh, nh1, i, j;
  74
  75   if(n < 4) return;
  76   wksp = (float *)malloc(sizeof(float) * n);
  77
  78   nh  = n>>1;
  79   nh1 = nh+1;
  80
  81   if(isign>=0) {
  82     for(i=0,j=0;j<n-3;j+=2,i++) {
  83       wksp[i]=C0*a[j]+C1*a[j+1]+C2*a[j+2]+C3*a[j+3];
  84       wksp[i+nh]=C3*a[j]-C2*a[j+1]+C1*a[j+2]-C0*a[j+3];
  85     }
  86     wksp[i]=C0*a[n-2]+C1*a[n-1]+C2*a[0]+C3*a[1];
  87     wksp[i+nh]=C3*a[n-2]-C2*a[n-1]+C1*a[0]-C0*a[1];
  88   } else {
  89     wksp[0]=C2*a[nh-1]+C1*a[n-1]+C0*a[0]+C3*a[nh];
  90     wksp[1]=C3*a[nh-1]-C0*a[n-1]+C1*a[0]-C2*a[nh];
  91     for(i=0,j=2;i<nh-1;i++,j+=2) {
  92       wksp[j]   = C2*a[i]+C1*a[i+nh]+C0*a[i+1]+C3*a[i+nh1];
  93       wksp[j+1] = C3*a[i]-C0*a[i+nh]+C1*a[i+1]-C2*a[i+nh1];
  94     }
  95   }
  96   for(i=0;i<n;i++) a[i]=wksp[i];
  97   free(wksp);
  98 }
  99
 100 // Multidimensional DWT/iDWT transform
 101 // a[]   - array of floating point values to be converted
 102 // dim[] - dimensions of the a[] array
 103 // ndim  - number of dimensions
 104 // isign - if <0, perform iDWT, otherwise perform DWT
 105 //
 106 // Dimensions *must* be powers of 2. Undefined but usually
 107 // bad things happen if any of the dimensions are not.
 108 //
 109 // This routine will also work with more than 3 dimensions
 110 // should the need arise to compress 3d motion images.
 111
 112 void dwt(float a[], unsigned long dim[], int ndim, int isign)
 113 {
 114   unsigned long i1,i2,i3,k,n,nnew,nprev=1,nt,ntot=1;
 115   int idim;
 116   float *wksp;
 117
 118   // Determine array extents and allocate workspace memory
 119   for(idim=0;idim<ndim;idim++) ntot *= dim[idim];
 120   wksp=(float *)malloc(sizeof(float)*ntot);
 121
 122   // Perform monodimensional DWT/iDWT across float array,
 123   // One dimension at a time...
 124   for(idim=0;idim<ndim;idim++) {
 125     n = dim[idim];
 126
 127     // 'nnew'  - length for 'line' progression
 128     // 'nprev' - length for individual element progression
 129     nnew = n * nprev;
 130
 131     // No point in transforming if there's 4 or less samples...
 132     if(n > 4) {
 133
 134       // Cycle through the entire array
 135       for(i2=0;i2<ntot;i2+=nnew) {
 136         for(i1=0;i1<nprev;i1++) {
 137           // Copy the appropriate elements into a monodimensional workspace...
 138           for(i3=i1+i2,k=0;k<n;k++,i3+=nprev) wksp[k] = a[i3];
 139
 140           // And perform the DWT/iDWT in pyramid fashion
 141           if(isign >= 0) {
 142             for(nt=n;nt>=4;nt>>=1) zdaub4(wksp,nt,isign);
 143           } else {
 144             for(nt=4;nt<=n;nt<<=1) zdaub4(wksp,nt,isign);
 145           }
 146
 147           // Copy the results back from the conversion
 148           for(i3=i1+i2,k=0;k<n;k++,i3+=nprev) a[i3]=wksp[k];
 149         }
 150       }
 151     }
 152     nprev = nnew;
 153   }
 154   free(wksp);
 155 }
 156
 157
 158 #ifdef DWTMAIN
 159 int main() {
 160   float tf[64*64*16];
 161   int x, y, z, dim[3];
 162   ulong l;
 163
 164   for(l=0;l<64*64*16;l++) tf[l] = l;
 165
 166   dim[0] = 64;
 167   dim[1] = 64;
 168   dim[2] = 16;
 169
 170   dwt(tf, dim, 3, 1);
 171   for(l=0;l<64*64*16;l++) tf[l] = floor(tf[l]/2)*2;
 172   dwt(tf, dim, 3, -1);
 173
 174   for(z=0;z<16;z++) for(y=0;y<64;y++) for(x=0;x<64;x++) {
 175     l = x+y*64+z*64*64;
 176     printf("%ld %f\n", l, tf[l]);
 177   }
 178 }
 179
 180 #endif