htswanalysis/MACS/lib/gsl/gsl-1.11/wavelet/dwt.c

   1 /* wavelet/dwt.c
   2  *
   3  * Copyright (C) 2004 Ivo Alxneit
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* function dwt_step is based on the public domain function pwt.c
  21  * available from http://www.numerical-recipes.com
  22  */
  23
  24 #include <config.h>
  25 #include <gsl/gsl_errno.h>
  26 #include <gsl/gsl_wavelet.h>
  27 #include <gsl/gsl_wavelet2d.h>
  28
  29 #define ELEMENT(a,stride,i) ((a)[(stride)*(i)])
  30
  31 static int binary_logn (const size_t n);
  32 static void dwt_step (const gsl_wavelet * w, double *a, size_t stride, size_t n, gsl_wavelet_direction dir, gsl_wavelet_workspace * work);
  33
  34 static int
  35 binary_logn (const size_t n)
  36 {
  37   size_t ntest;
  38   size_t logn = 0;
  39   size_t k = 1;
  40
  41   while (k < n)
  42     {
  43       k *= 2;
  44       logn++;
  45     }
  46
  47   ntest = (1 << logn);
  48
  49   if (n != ntest)
  50     {
  51       return -1;                /* n is not a power of 2 */
  52     }
  53
  54   return logn;
  55 }
  56
  57 static void
  58 dwt_step (const gsl_wavelet * w, double *a, size_t stride, size_t n,
  59           gsl_wavelet_direction dir, gsl_wavelet_workspace * work)
  60 {
  61   double ai, ai1;
  62   size_t i, ii;
  63   size_t jf;
  64   size_t k;
  65   size_t n1, ni, nh, nmod;
  66
  67   for (i = 0; i < work->n; i++)
  68     {
  69       work->scratch[i] = 0.0;
  70     }
  71
  72   nmod = w->nc * n;
  73   nmod -= w->offset;            /* center support */
  74
  75   n1 = n - 1;
  76   nh = n >> 1;
  77
  78   if (dir == gsl_wavelet_forward)
  79     {
  80       for (ii = 0, i = 0; i < n; i += 2, ii++)
  81         {
  82           double h = 0, g = 0;
  83
  84           ni = i + nmod;
  85
  86           for (k = 0; k < w->nc; k++)
  87             {
  88               jf = n1 & (ni + k);
  89               h += w->h1[k] * ELEMENT (a, stride, jf);
  90               g += w->g1[k] * ELEMENT (a, stride, jf);
  91             }
  92
  93           work->scratch[ii] += h;
  94           work->scratch[ii + nh] += g;
  95         }
  96     }
  97   else
  98     {
  99       for (ii = 0, i = 0; i < n; i += 2, ii++)
 100         {
 101           ai = ELEMENT (a, stride, ii);
 102           ai1 = ELEMENT (a, stride, ii + nh);
 103           ni = i + nmod;
 104           for (k = 0; k < w->nc; k++)
 105             {
 106               jf = (n1 & (ni + k));
 107               work->scratch[jf] += (w->h2[k] * ai + w->g2[k] * ai1);
 108             }
 109         }
 110     }
 111
 112   for (i = 0; i < n; i++)
 113     {
 114       ELEMENT (a, stride, i) = work->scratch[i];
 115     }
 116 }
 117
 118 int
 119 gsl_wavelet_transform (const gsl_wavelet * w,
 120                        double *data, size_t stride, size_t n,
 121                        gsl_wavelet_direction dir,
 122                        gsl_wavelet_workspace * work)
 123 {
 124   size_t i;
 125
 126   if (work->n < n)
 127     {
 128       GSL_ERROR ("not enough workspace provided", GSL_EINVAL);
 129     }
 130
 131   if (binary_logn (n) == -1)
 132     {
 133       GSL_ERROR ("n is not a power of 2", GSL_EINVAL);
 134     }
 135
 136   if (n < 2)
 137     {
 138       return GSL_SUCCESS;
 139     }
 140
 141   if (dir == gsl_wavelet_forward)
 142     {
 143       for (i = n; i >= 2; i >>= 1)
 144         {
 145           dwt_step (w, data, stride, i, dir, work);
 146         }
 147     }
 148   else
 149     {
 150       for (i = 2; i <= n; i <<= 1)
 151         {
 152           dwt_step (w, data, stride, i, dir, work);
 153         }
 154     }
 155
 156   return GSL_SUCCESS;
 157 }
 158
 159 int
 160 gsl_wavelet_transform_forward (const gsl_wavelet * w,
 161                                double *data, size_t stride, size_t n,
 162                                gsl_wavelet_workspace * work)
 163 {
 164   return gsl_wavelet_transform (w, data, stride, n, gsl_wavelet_forward, work);
 165 }
 166
 167 int
 168 gsl_wavelet_transform_inverse (const gsl_wavelet * w,
 169                                double *data, size_t stride, size_t n,
 170                                gsl_wavelet_workspace * work)
 171 {
 172   return gsl_wavelet_transform (w, data, stride, n, gsl_wavelet_backward, work);
 173 }
 174
 175
 176 /* Leaving this out for now BJG */
 177 #if 0
 178 int
 179 gsl_dwt_vector (const gsl_wavelet * w, gsl_vector *v, gsl_wavelet_direction
 180                 dir, gsl_wavelet_workspace * work)
 181 {
 182   return gsl_dwt (w, v->data, v->stride, v->size, dir, work);
 183 }
 184 #endif
 185
 186 int
 187 gsl_wavelet2d_transform (const gsl_wavelet * w,
 188                          double *data, size_t tda, size_t size1,
 189                          size_t size2, gsl_wavelet_direction dir,
 190                          gsl_wavelet_workspace * work)
 191 {
 192   size_t i;
 193
 194   if (size1 != size2)
 195     {
 196       GSL_ERROR ("2d dwt works only with square matrix", GSL_EINVAL);
 197     }
 198
 199   if (work->n < size1)
 200     {
 201       GSL_ERROR ("not enough workspace provided", GSL_EINVAL);
 202     }
 203
 204   if (binary_logn (size1) == -1)
 205     {
 206       GSL_ERROR ("n is not a power of 2", GSL_EINVAL);
 207     }
 208
 209   if (size1 < 2)
 210     {
 211       return GSL_SUCCESS;
 212     }
 213
 214   if (dir == gsl_wavelet_forward)
 215     {
 216       for (i = 0; i < size1; i++)       /* for every row j */
 217         {
 218           gsl_wavelet_transform (w, &ELEMENT(data, tda, i), 1, size1, dir, work);
 219         }
 220       for (i = 0; i < size2; i++)       /* for every column j */
 221         {
 222           gsl_wavelet_transform (w, &ELEMENT(data, 1, i), tda, size2, dir, work);
 223         }
 224     }
 225   else
 226     {
 227       for (i = 0; i < size2; i++)       /* for every column j */
 228         {
 229           gsl_wavelet_transform (w, &ELEMENT(data, 1, i), tda, size2, dir, work);
 230         }
 231       for (i = 0; i < size1; i++)       /* for every row j */
 232         {
 233           gsl_wavelet_transform (w, &ELEMENT(data, tda, i), 1, size1, dir, work);
 234         }
 235     }
 236
 237   return GSL_SUCCESS;
 238 }
 239
 240 int
 241 gsl_wavelet2d_nstransform (const gsl_wavelet * w,
 242                            double *data, size_t tda, size_t size1,
 243                            size_t size2, gsl_wavelet_direction dir,
 244                            gsl_wavelet_workspace * work)
 245 {
 246   size_t i, j;
 247
 248   if (size1 != size2)
 249     {
 250       GSL_ERROR ("2d dwt works only with square matrix", GSL_EINVAL);
 251     }
 252
 253   if (work->n < size1)
 254     {
 255       GSL_ERROR ("not enough workspace provided", GSL_EINVAL);
 256     }
 257
 258   if (binary_logn (size1) == -1)
 259     {
 260       GSL_ERROR ("n is not a power of 2", GSL_EINVAL);
 261     }
 262
 263   if (size1 < 2)
 264     {
 265       return GSL_SUCCESS;
 266     }
 267
 268   if (dir == gsl_wavelet_forward)
 269     {
 270       for (i = size1; i >= 2; i >>= 1)
 271         {
 272           for (j = 0; j < i; j++)       /* for every row j */
 273             {
 274               dwt_step (w, &ELEMENT(data, tda, j), 1, i, dir, work);
 275             }
 276           for (j = 0; j < i; j++)       /* for every column j */
 277             {
 278               dwt_step (w, &ELEMENT(data, 1, j), tda, i, dir, work);
 279             }
 280         }
 281     }
 282   else
 283     {
 284       for (i = 2; i <= size1; i <<= 1)
 285         {
 286           for (j = 0; j < i; j++)       /* for every column j */
 287             {
 288               dwt_step (w, &ELEMENT(data, 1, j), tda, i, dir, work);
 289             }
 290           for (j = 0; j < i; j++)       /* for every row j */
 291             {
 292               dwt_step (w, &ELEMENT(data, tda, j), 1, i, dir, work);
 293             }
 294         }
 295     }
 296
 297   return GSL_SUCCESS;
 298 }
 299
 300
 301 int
 302 gsl_wavelet2d_transform_forward (const gsl_wavelet * w,
 303                                  double *data, size_t tda, size_t size1,
 304                                  size_t size2, gsl_wavelet_workspace * work)
 305 {
 306   return gsl_wavelet2d_transform (w, data, tda, size1, size2, gsl_wavelet_forward, work);
 307 }
 308
 309 int
 310 gsl_wavelet2d_transform_inverse (const gsl_wavelet * w,
 311                                   double *data, size_t tda, size_t size1,
 312                                   size_t size2, gsl_wavelet_workspace * work)
 313 {
 314   return gsl_wavelet2d_transform (w, data, tda, size1, size2, gsl_wavelet_backward, work);
 315 }
 316
 317 int
 318 gsl_wavelet2d_nstransform_forward (const gsl_wavelet * w,
 319                                  double *data, size_t tda, size_t size1,
 320                                  size_t size2, gsl_wavelet_workspace * work)
 321 {
 322   return gsl_wavelet2d_nstransform (w, data, tda, size1, size2, gsl_wavelet_forward, work);
 323 }
 324
 325 int
 326 gsl_wavelet2d_nstransform_inverse (const gsl_wavelet * w,
 327                                   double *data, size_t tda, size_t size1,
 328                                   size_t size2, gsl_wavelet_workspace * work)
 329 {
 330   return gsl_wavelet2d_nstransform (w, data, tda, size1, size2, gsl_wavelet_backward, work);
 331 }
 332
 333
 334 int
 335 gsl_wavelet2d_transform_matrix (const gsl_wavelet * w,
 336                                 gsl_matrix * a,
 337                                 gsl_wavelet_direction dir,
 338                                 gsl_wavelet_workspace * work)
 339 {
 340   return gsl_wavelet2d_transform (w, a->data,
 341                                   a->tda, a->size1, a->size2,
 342                                   dir, work);
 343 }
 344
 345 int
 346 gsl_wavelet2d_transform_matrix_forward (const gsl_wavelet * w,
 347                                         gsl_matrix * a,
 348                                         gsl_wavelet_workspace * work)
 349 {
 350   return gsl_wavelet2d_transform (w, a->data,
 351                                   a->tda, a->size1, a->size2,
 352                                   gsl_wavelet_forward, work);
 353 }
 354
 355 int
 356 gsl_wavelet2d_transform_matrix_inverse (const gsl_wavelet * w,
 357                                         gsl_matrix * a,
 358                                         gsl_wavelet_workspace * work)
 359 {
 360   return gsl_wavelet2d_transform (w, a->data,
 361                                   a->tda, a->size1, a->size2,
 362                                   gsl_wavelet_backward, work);
 363 }
 364
 365 int
 366 gsl_wavelet2d_nstransform_matrix (const gsl_wavelet * w,
 367                                 gsl_matrix * a,
 368                                 gsl_wavelet_direction dir,
 369                                 gsl_wavelet_workspace * work)
 370 {
 371   return gsl_wavelet2d_nstransform (w, a->data,
 372                                     a->tda, a->size1, a->size2,
 373                                     dir, work);
 374 }
 375
 376 int
 377 gsl_wavelet2d_nstransform_matrix_forward (const gsl_wavelet * w,
 378                                         gsl_matrix * a,
 379                                         gsl_wavelet_workspace * work)
 380 {
 381   return gsl_wavelet2d_nstransform (w, a->data,
 382                                     a->tda, a->size1, a->size2,
 383                                     gsl_wavelet_forward, work);
 384 }
 385
 386 int
 387 gsl_wavelet2d_nstransform_matrix_inverse (const gsl_wavelet * w,
 388                                         gsl_matrix * a,
 389                                         gsl_wavelet_workspace * work)
 390 {
 391   return gsl_wavelet2d_nstransform (w, a->data,
 392                                     a->tda, a->size1, a->size2,
 393                                     gsl_wavelet_backward, work);
 394 }
 395
 396