1 // $Id: dwt.c,v 1.1 2001/02/13 01:06:24 giles Exp $
4 // Revision 1.1 2001/02/13 01:06:24 giles
11 void daub4(float a
[], unsigned long n
, int isign
)
14 unsigned long nh
, nh1
, i
, j
;
17 wksp
= (float *)malloc(sizeof(float) * n
);
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];
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];
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
];
34 for(i
=1;i
<=n
;i
++) a
[i
]=wksp
[i
];
38 void wtn(float a
[], unsigned long nn
[], int ndim
, int isign
)
40 unsigned long i1
,i2
,i3
,k
,n
,nnew
,nprev
=1,nt
,ntot
=1;
44 for(idim
=1;idim
<=ndim
;idim
++) ntot
*= nn
[idim
];
45 wksp
=(float *)malloc(sizeof(float)*ntot
);
47 for(idim
=1;idim
<=ndim
;idim
++) {
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
];
63 // Monodimensional DWT/iDWT transform
64 // if isign is positive - Forward DWT
65 // if isign is negative - Inverse DWT
67 // Originally snagged from Numeric Recipes in C (2nd Ed)
68 // Modified to get around braindeadedness with arrays starting at [1].
70 void zdaub4(float a
[], unsigned long n
, int isign
)
73 unsigned long nh
, nh1
, i
, j
;
76 wksp
= (float *)malloc(sizeof(float) * n
);
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];
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];
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
];
96 for(i
=0;i
<n
;i
++) a
[i
]=wksp
[i
];
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
106 // Dimensions *must* be powers of 2. Undefined but usually
107 // bad things happen if any of the dimensions are not.
109 // This routine will also work with more than 3 dimensions
110 // should the need arise to compress 3d motion images.
112 void dwt(float a
[], unsigned long dim
[], int ndim
, int isign
)
114 unsigned long i1
,i2
,i3
,k
,n
,nnew
,nprev
=1,nt
,ntot
=1;
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
);
122 // Perform monodimensional DWT/iDWT across float array,
123 // One dimension at a time...
124 for(idim
=0;idim
<ndim
;idim
++) {
127 // 'nnew' - length for 'line' progression
128 // 'nprev' - length for individual element progression
131 // No point in transforming if there's 4 or less samples...
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
];
140 // And perform the DWT/iDWT in pyramid fashion
142 for(nt
=n
;nt
>=4;nt
>>=1) zdaub4(wksp
,nt
,isign
);
144 for(nt
=4;nt
<=n
;nt
<<=1) zdaub4(wksp
,nt
,isign
);
147 // Copy the results back from the conversion
148 for(i3
=i1
+i2
,k
=0;k
<n
;k
++,i3
+=nprev
) a
[i3
]=wksp
[k
];
164 for(l
=0;l
<64*64*16;l
++) tf
[l
] = l
;
171 for(l
=0;l
<64*64*16;l
++) tf
[l
] = floor(tf
[l
]/2)*2;
174 for(z
=0;z
<16;z
++) for(y
=0;y
<64;y
++) for(x
=0;x
<64;x
++) {
176 printf("%ld %f\n", l
, tf
[l
]);