00001
00002
00003 #include <stdlib.h>
00004 #include <math.h>
00005 #include "WaveletTransform.h"
00006
00007
00008 static float haar_d[3] = {0.0, 0.7071067811865475, 0.7071067811865475};
00009 static float haar_r[3];
00010
00011
00012 static float daub2_d[5] = {0.0, 4.8296291e-01, 8.3651630e-01, 2.2414387e-01,
00013 -1.2940952e-01};
00014 static float daub2_r[5];
00015 static float daub3_d[7] = {0.0, 3.3267055e-01, 8.0689151e-01, 4.5987750e-01,
00016 -1.3501102e-01,-8.5441274e-02, 3.5226292e-02};
00017 static float daub3_r[5];
00018
00019 static float daub4_d[9] = {0.0, 2.3037781e-01, 7.1484657e-01, 6.3088077e-01,
00020 -2.7983769e-02,-1.8703481e-01, 3.0841382e-02,
00021 3.2883012e-02,-1.0597402e-02};
00022 static float daub4_r[9];
00023 static float daub5_d[11]= {0.0, 1.6010240e-01, 6.0382927e-01, 7.2430853e-01,
00024 1.3842815e-01,-2.4229489e-01,-3.2244870e-02,
00025 7.7571494e-02,-6.2414902e-03,-1.2580752e-02,
00026 3.3357253e-03};
00027 static float daub5_r[11];
00028 static float daub6_d[13]= {0.0, 1.1154074e-01, 4.9462389e-01, 7.5113391e-01,
00029 3.1525035e-01,-2.2626469e-01,-1.2976687e-01,
00030 9.7501606e-02, 2.7522866e-02,-3.1582039e-02,
00031 5.5384220e-04, 4.7772575e-03,-1.0773011e-03};
00032 static float daub6_r[13];
00033 static float daub7_d[15]= {0.0, 7.7852054e-02, 3.9653932e-01, 7.2913209e-01,
00034 4.6978229e-01,-1.4390600e-01,-2.2403618e-01,
00035 7.1309219e-02, 8.0612609e-02,-3.8029937e-02,
00036 -1.6574542e-02, 1.2550999e-02, 4.2957797e-04,
00037 -1.8016407e-03, 3.5371380e-04};
00038 static float daub7_r[15];
00039 static float daub8_d[17]= {0.0, 5.4415842e-02, 3.1287159e-01, 6.7563074e-01,
00040 5.8535468e-01,-1.5829105e-02,-2.8401554e-01,
00041 4.7248457e-04, 1.2874743e-01,-1.7369301e-02,
00042 -4.4088254e-02, 1.3981028e-02, 8.7460940e-03,
00043 -4.8703530e-03,-3.9174037e-04, 6.7544941e-04,
00044 -1.1747678e-04};
00045 static float daub8_r[17];
00046 static float daub9_d[19]= {0.0, 3.8077947e-02, 2.4383467e-01, 6.0482312e-01,
00047 6.5728808e-01, 1.3319739e-01,-2.9327378e-01,
00048 -9.6840783e-02, 1.4854075e-01, 3.0725681e-02,
00049 -6.7632829e-02, 2.5094711e-04, 2.2361662e-02,
00050 -4.7232048e-03,-4.2815037e-03, 1.8476469e-03,
00051 2.3038576e-04,-2.5196319e-04, 3.9347320e-05};
00052 static float daub9_r[19];
00053
00054
00055 static float bior_1_3_lo_d[7] = { 0.0, -0.08838834764832, 0.08838834764832,
00056 0.70710678118655, 0.70710678118655,
00057 0.08838834764832, -0.08838834764832};
00058 static float bior_1_3_hi_d[7] = { 0.0, 0.0, 0.0, -0.70710678118655,
00059 0.70710678118655, 0.0, 0.0};
00060 static float bior_1_3_lo_r[7] = { 0.0, 0.0, 0.0, 0.70710678118655,
00061 0.70710678118655, 0.0, 0.0};
00062 static float bior_1_3_hi_r[7] = { 0.0, -0.08838834764832, -0.08838834764832,
00063 0.70710678118655, -0.70710678118655,
00064 0.08838834764832, 0.08838834764832};
00065
00066 static float bior_1_5_lo_d[11] = { 0.0, 1.6572815e-02, -1.6572815e-02, -1.2153398e-01,
00067 1.2153398e-01, 7.0710678e-01, 7.0710678e-01,
00068 1.2153398e-01, -1.2153398e-01, -1.6572815e-02,
00069 1.6572815e-02};
00070 static float bior_1_5_hi_d[11] = { 0.0, -0.0000000e+00, 0.0000000e+00, -0.0000000e+00,
00071 0.0000000e+00, -7.0710678e-01, 7.0710678e-01,
00072 -0.0000000e+00, 0.0000000e+00, -0.0000000e+00,
00073 0.0000000e+00};
00074 static float bior_1_5_lo_r[11] = { 0.0, 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,
00075 0.0000000e+00, 7.0710678e-01, 7.0710678e-01,
00076 0.0000000e+00, 0.0000000e+00, 0.0000000e+00,
00077 0.0000000e+00};
00078 static float bior_1_5_hi_r[11] = { 0.0, 1.6572815e-02, 1.6572815e-02, -1.2153398e-01,
00079 -1.2153398e-01, 7.0710678e-01, -7.0710678e-01,
00080 1.2153398e-01, 1.2153398e-01, -1.6572815e-02,
00081 -1.6572815e-02};
00082
00083 static float bior_2_2_lo_d[7] = {0.0, 0.0000000e+00, -1.7677670e-01, 3.5355339e-01,
00084 1.0606602e+00, 3.5355339e-01, -1.7677670e-01};
00085 static float bior_2_2_hi_d[7] = {0.0, -0.0000000e+00, 3.5355339e-01, -7.0710678e-01,
00086 3.5355339e-01, -0.0000000e+00, 0.0000000e+00};
00087 static float bior_2_2_lo_r[7] = {0.0, 0.0000000e+00, 3.5355339e-01, 7.0710678e-01,
00088 3.5355339e-01, 0.0000000e+00, 0.0000000e+00};
00089 static float bior_2_2_hi_r[7] = {0.0, 0.0000000e+00, 1.7677670e-01, 3.5355339e-01,
00090 -1.0606602e+00, 3.5355339e-01, 1.7677670e-01};
00091
00092
00093
00094 static float bior_2_8_lo_d[19] = {0.0,
00095 0.0 , 0.00151054305063, -0.00302108610126, -0.01294751186255,
00096 0.02891610982635, 0.05299848189069, -0.13491307360774, -0.16382918343409,
00097 0.46257144047592, 0.95164212189718, 0.46257144047592, -0.16382918343409,
00098 -0.13491307360774, 0.05299848189069, 0.02891610982635, -0.01294751186255,
00099 -0.00302108610126, 0.00151054305063};
00100 static float bior_2_8_hi_d[19] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
00101 0.35355339059327, -0.70710678118655, 0.35355339059327, 0.0, 0.0,
00102 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
00103 static float bior_2_8_lo_r[19] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
00104 0.35355339059327, 0.70710678118655, 0.35355339059327, 0.0, 0.0,
00105 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
00106 static float bior_2_8_hi_r[19] = {0.0,
00107 0.0 , -0.00151054305063, -0.00302108610126, 0.01294751186255,
00108 0.02891610982635, -0.05299848189069, -0.13491307360774, 0.16382918343409,
00109 0.46257144047592, -0.95164212189718, 0.46257144047592, 0.16382918343409,
00110 -0.13491307360774, -0.05299848189069, 0.02891610982635, 0.01294751186255,
00111 -0.00302108610126, -0.00151054305063 };
00112
00113
00114 static float bior_3_1_lo_d[5] ={0.0, -0.35355339059327, 1.06066017177982,
00115 1.06066017177982, -0.35355339059327};
00116 static float bior_3_1_hi_d[5] ={0.0, -0.17677669529664, 0.53033008588991,
00117 -0.53033008588991, 0.17677669529664};
00118 static float bior_3_1_lo_r[5] ={0.0, 0.17677669529664, 0.53033008588991,
00119 0.53033008588991, 0.17677669529664};
00120 static float bior_3_1_hi_r[5] ={0.0, -0.35355339059327, -1.06066017177982,
00121 1.06066017177982, 0.35355339059327};
00122
00123 void WaveletTransform::pwtset(int wavelet_type, wavefilt &wfilt)
00124 {
00125 float sig;
00126 int k;
00127
00128 switch(wavelet_type) {
00129
00130 case HAAR:
00131 wfilt.ncof = 2;
00132 wfilt.cc = haar_d;
00133 wfilt.cr = haar_r;
00134 wfilt.ioff = -1;
00135 wfilt.joff = -1;
00136 break;
00137 case DAUB_2:
00138 wfilt.ncof = 4;
00139 wfilt.cc = daub2_d;
00140 wfilt.cr = daub2_r;
00141 wfilt.ioff = -2;
00142 wfilt.joff = -2;
00143 break;
00144 case DAUB_3:
00145 wfilt.ncof = 6;
00146 wfilt.cc = daub3_d;
00147 wfilt.cr = daub3_r;
00148 wfilt.ioff = -2;
00149 wfilt.joff = -4;
00150 break;
00151 case DAUB_4:
00152 wfilt.ncof = 8;
00153 wfilt.cc = daub4_d;
00154 wfilt.cr = daub4_r;
00155 wfilt.ioff = -2;
00156 wfilt.joff = -6;
00157 break;
00158 case DAUB_5:
00159 wfilt.ncof = 10;
00160 wfilt.cc = daub5_d;
00161 wfilt.cr = daub5_r;
00162 wfilt.ioff = -3;
00163 wfilt.joff = -7;
00164 break;
00165 case DAUB_6:
00166 wfilt.ncof = 12;
00167 wfilt.cc = daub6_d;
00168 wfilt.cr = daub6_r;
00169 wfilt.ioff = -3;
00170 wfilt.joff = -9;
00171 break;
00172 case DAUB_7:
00173 wfilt.ncof = 14;
00174 wfilt.cc = daub7_d;
00175 wfilt.cr = daub7_r;
00176 wfilt.ioff = -3;
00177 wfilt.joff = -11;
00178 break;
00179 case DAUB_8:
00180 wfilt.ncof = 16;
00181 wfilt.cc = daub8_d;
00182 wfilt.cr = daub8_r;
00183 wfilt.ioff = -3;
00184 wfilt.joff = -13;
00185 break;
00186 case DAUB_9:
00187 wfilt.ncof = 18;
00188 wfilt.cc = daub9_d;
00189 wfilt.cr = daub9_r;
00190 wfilt.ioff = -4;
00191 wfilt.joff = -14;
00192 break;
00193 case BIOR_1_3:
00194 wfilt.ncof = 6;
00195 wfilt.cc = bior_1_3_lo_d;
00196 wfilt.cr = bior_1_3_hi_d;
00197 break;
00198 case BIOR_1_3_REC:
00199 wfilt.ncof = 6;
00200 wfilt.cc = bior_1_3_lo_r;
00201 wfilt.cr = bior_1_3_hi_r;
00202 break;
00203 case BIOR_1_5:
00204 wfilt.ncof = 10;
00205 wfilt.cc = bior_1_5_lo_d;
00206 wfilt.cr = bior_1_5_hi_d;
00207 break;
00208 case BIOR_1_5_REC:
00209 wfilt.ncof = 10;
00210 wfilt.cc = bior_1_5_lo_r;
00211 wfilt.cr = bior_1_5_hi_r;
00212 break;
00213 case BIOR_2_2:
00214 wfilt.ncof = 6;
00215 wfilt.cc = bior_2_2_lo_d;
00216 wfilt.cr = bior_2_2_hi_d;
00217 break;
00218 case BIOR_2_2_REC:
00219 wfilt.ncof = 6;
00220 wfilt.cc = bior_2_2_lo_r;
00221 wfilt.cr = bior_2_2_hi_r;
00222 break;
00223 case BIOR_2_8:
00224 wfilt.ncof = 18;
00225 wfilt.cc = bior_2_8_lo_d;
00226 wfilt.cr = bior_2_8_hi_d;
00227 break;
00228 case BIOR_2_8_REC:
00229 wfilt.ncof = 18;
00230 wfilt.cc = bior_2_8_lo_r;
00231 wfilt.cr = bior_2_8_hi_r;
00232 break;
00233 case BIOR_3_1:
00234 wfilt.ncof = 4;
00235 wfilt.cc = bior_3_1_lo_d;
00236 wfilt.cr = bior_3_1_hi_d;
00237 break;
00238 case BIOR_3_1_REC:
00239 wfilt.ncof = 4;
00240 wfilt.cc = bior_3_1_lo_r;
00241 wfilt.cr = bior_3_1_hi_r;
00242 break;
00243 default:
00244 cout << "Unimplemented wavelet filter" << endl;
00245 }
00246
00247 if (wavelet_type < 15) {
00248 sig = -1.0;
00249 for (k=1;k<=wfilt.ncof;k++) {
00250 wfilt.cr[wfilt.ncof+1-k] = sig*wfilt.cc[k];
00251 sig = -sig;
00252 }
00253 }
00254 else {
00255 wfilt.ioff = wfilt.joff = -(wfilt.ncof/2)-1;
00256 }
00257 }
00258
00259 void WaveletTransform::pwt(float *a, int n, int isign)
00260 {
00261 float ai, ai1, *wksp;
00262 int i, ii, j, jf, jr, k, ni, nj, nh, nmod;
00263
00264 if (n < 2) return;
00265 wksp = (float *) calloc(1,(n+1)*sizeof(float));
00266 nmod = wfilt_dec.ncof*n;
00267
00268 nh = n >> 1;
00269
00270 if (isign >= 0)
00271 {
00272 for (ii=1, i=1; i<=n; i+=2, ii++)
00273 {
00274 ni = i + nmod + wfilt_dec.ioff;
00275 nj = i + nmod + wfilt_dec.joff;
00276 for (k=1; k<=wfilt_dec.ncof; k++)
00277 {
00278
00279
00280 jf = ni + k;
00281 jr = nj + k;
00282 if (jf > n-1) jf = jf % n;
00283 if (jr > n-1) jr = jr % n;
00284
00285 if (wavelet_type < 15) {
00286 wksp[ii] += wfilt_dec.cc[k]*a[jf+1];
00287 wksp[ii+nh] += wfilt_dec.cr[k]*a[jr+1];
00288 }
00289 else {
00290 wksp[ii] += wfilt_dec.cc[wfilt_dec.ncof+1-k]*a[jf+1];
00291 wksp[ii+nh] += wfilt_dec.cr[wfilt_dec.ncof+1-k]*a[jr+1];
00292 }
00293 }
00294 }
00295 }
00296 else
00297 {
00298 for (ii=1, i=1; i<=n; i+=2, ii++)
00299 {
00300 ai = a[ii];
00301 ai1 = a[ii+nh];
00302 ni = i + nmod + wfilt_rec.ioff;
00303 nj = i + nmod + wfilt_rec.joff;
00304 for (k=1; k<=wfilt_rec.ncof; k++)
00305 {
00306
00307
00308 jf = ni + k;
00309 jr = nj + k;
00310 if (jf > n-1) jf = jf % n;
00311 if (jr > n-1) jr = jr % n;
00312 wksp[jf+1] += wfilt_rec.cc[k]*ai;
00313 wksp[jr+1] += wfilt_rec.cr[k]*ai1;
00314 }
00315 }
00316 }
00317 for (j=1; j<=n; j++) a[j] = wksp[j];
00318 free(wksp);
00319 }
00320
00321 void WaveletTransform::fwt(int levels)
00322 {
00323 int i,j;
00324 float *wksp1, *wksp2;
00325 int dim1, dim2;
00326 float *a;
00327 a = image->getPixels();
00328
00329 dim1 = width;
00330 dim2 = height;
00331
00332 wksp1 = (float *) calloc(1,(dim1+1) * sizeof(float));
00333 wksp2 = (float *) calloc(1,(dim2+1) * sizeof(float));
00334
00335 while (levels) {
00336 for (j=1; j<=dim2; j++) {
00337 for (i=1; i<=dim1; i++)
00338 wksp1[i] = a[(j-1)*width+(i-1)];
00339 pwt(wksp1, dim1, 1);
00340 for (i=1; i<=dim1; i++)
00341 a[(j-1)*width+(i-1)] = wksp1[i];
00342 }
00343
00344 for (i=1; i<=dim1; i++) {
00345 for (j=1; j<=dim2; j++)
00346 wksp2[j] = a[(j-1)*width+(i-1)];
00347 pwt(wksp2, dim2, 1);
00348 for (j=1; j<=dim2; j++)
00349 a[(j-1)*width+(i-1)] = wksp2[j];
00350 }
00351
00352 dim1 /= 2;
00353 dim2 /= 2;
00354 levels--;
00355 }
00356
00357 free(wksp1);
00358 free(wksp2);
00359
00360 }
00361
00362 void WaveletTransform::ifwt(int levels)
00363 {
00364 int i,j;
00365 float *wksp1, *wksp2;
00366 int dim1, dim2;
00367 float *a;
00368 a = image->getPixels();
00369
00370 dim1 = width;
00371 dim2 = height;
00372
00373 wksp1 = (float *) calloc(1,(dim1+1) * sizeof(float));
00374 wksp2 = (float *) calloc(1,(dim2+1) * sizeof(float));
00375
00376 dim1 = dim1>>(levels-1);
00377 dim2 = dim2>>(levels-1);
00378
00379 while (levels) {
00380 for (j=1; j<=dim2; j++) {
00381 for (i=1; i<=dim1; i++)
00382 wksp1[i] = a[(j-1)*width+(i-1)];
00383 pwt(wksp1, dim1, -1);
00384 for (i=1; i<=dim1; i++)
00385 a[(j-1)*width+(i-1)] = wksp1[i];
00386 }
00387
00388 for (i=1; i<=dim1; i++) {
00389 for (j=1; j<=dim2; j++)
00390 wksp2[j] = a[(j-1)*width+(i-1)];
00391 pwt(wksp2, dim2, -1);
00392 for (j=1; j<=dim2; j++)
00393 a[(j-1)*width+(i-1)] = wksp2[j];
00394 }
00395
00396 dim1 *= 2;
00397 dim2 *= 2;
00398 levels--;
00399 }
00400
00401 free(wksp1);
00402 free(wksp2);
00403 }
00404
00405 void fwt(FloatImage &input, int levels, WAVELET_TYPE wave_type = HAAR)
00406 {
00407 WaveletTransform a(input, wave_type);
00408 a.fwt(levels);
00409 }
00410
00411 void ifwt(FloatImage &input, int levels, WAVELET_TYPE wave_type = HAAR)
00412 {
00413 WaveletTransform a(input, wave_type);
00414 a.ifwt(levels);
00415 }
00416