Main Page   Modules   Class Hierarchy   Compound List   File List   Compound Members   Related Pages   Examples  

WaveletTransform.cpp

00001 /* Copyright (c) 2001 C. Grigorescu */
00002 
00003 #include <stdlib.h>
00004 #include <math.h>
00005 #include "WaveletTransform.h"
00006 
00007 /* Haar */
00008 static float haar_d[3] = {0.0,  0.7071067811865475,  0.7071067811865475};
00009 static float haar_r[3];
00010 
00011 /* Daubechies */
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 /* Biorthogonal 1.3 */
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 /* Biorthogonal 1.5 */
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 /* Biorthogonal 2.2 */
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 /* Biorthogonal 2.8 */
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 /* Biorthogonal 3.1 */
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   //  n1 = n-1;
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 //            jf = (n1 & (ni + k));
00279 //            jr = (n1 & (nj + k));
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 //            jf = (n1 & (ni + k));
00307 //            jr = (n1 & (nj + k));
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