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

   1 /* wavelet/bspline.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 /* Coefficients are from A. Cohen, I. Daubechies, and J.-C. Feauveau;
  21  * "Biorthogonal Bases of Compactly Supported Wavelets", Communications
  22  * on Pure and Applied Mathematics, 45 (1992) 485--560 (table 6.1).
  23  *
  24  * Note the following errors in table 1:
  25  *
  26  * N = 2, N~ = 4, m0~
  27  * the second term in z^-1 (45/64 z^-1) should be left out.
  28  *
  29  * N = 3, N~ = 7, m0~
  30  * the term 336z^-3 should read 363z^-3.
  31  */
  32
  33 #include <config.h>
  34 #include <gsl/gsl_errno.h>
  35 #include <gsl/gsl_math.h>
  36 #include <gsl/gsl_wavelet.h>
  37
  38 static const double h1_103[6] = { -0.0883883476483184405501055452631,
  39   0.0883883476483184405501055452631,
  40   M_SQRT1_2,
  41   M_SQRT1_2,
  42   0.0883883476483184405501055452631,
  43   -0.0883883476483184405501055452631
  44 };
  45
  46 static const double g2_103[6] = { -0.0883883476483184405501055452631,
  47   -0.0883883476483184405501055452631,
  48   M_SQRT1_2,
  49   -(M_SQRT1_2),
  50   0.0883883476483184405501055452631,
  51   0.0883883476483184405501055452631
  52 };
  53
  54 static const double h1_105[10] = { 0.0165728151840597076031447897368,
  55   -0.0165728151840597076031447897368,
  56   -0.1215339780164378557563951247368,
  57   0.1215339780164378557563951247368,
  58   M_SQRT1_2,
  59   M_SQRT1_2,
  60   0.1215339780164378557563951247368,
  61   -0.1215339780164378557563951247368,
  62   -0.0165728151840597076031447897368,
  63   0.0165728151840597076031447897368
  64 };
  65
  66 static const double g2_105[10] = { 0.0165728151840597076031447897368,
  67   0.0165728151840597076031447897368,
  68   -0.1215339780164378557563951247368,
  69   -0.1215339780164378557563951247368,
  70   M_SQRT1_2,
  71   -(M_SQRT1_2),
  72   0.1215339780164378557563951247368,
  73   0.1215339780164378557563951247368,
  74   -0.0165728151840597076031447897368,
  75   -0.0165728151840597076031447897368
  76 };
  77
  78 static const double g1_1[10] = { 0.0, 0.0, 0.0, 0.0,
  79   M_SQRT1_2,
  80   -(M_SQRT1_2),
  81   0.0, 0.0, 0.0, 0.0
  82 };
  83
  84 static const double h2_1[10] = { 0.0, 0.0, 0.0, 0.0,
  85   M_SQRT1_2,
  86   M_SQRT1_2,
  87   0.0, 0.0, 0.0, 0.0
  88 };
  89
  90 static const double h1_202[6] = { -0.1767766952966368811002110905262,
  91   0.3535533905932737622004221810524,
  92   1.0606601717798212866012665431573,
  93   0.3535533905932737622004221810524,
  94   -0.1767766952966368811002110905262,
  95   0.0
  96 };
  97
  98 static const double g2_202[6] = { 0.0,
  99   -0.1767766952966368811002110905262,
 100   -0.3535533905932737622004221810524,
 101   1.0606601717798212866012665431573,
 102   -0.3535533905932737622004221810524,
 103   -0.1767766952966368811002110905262
 104 };
 105
 106 static const double h1_204[10] = { 0.0331456303681194152062895794737,
 107   -0.0662912607362388304125791589473,
 108   -0.1767766952966368811002110905262,
 109   0.4198446513295125926130013399998,
 110   0.9943689110435824561886873842099,
 111   0.4198446513295125926130013399998,
 112   -0.1767766952966368811002110905262,
 113   -0.0662912607362388304125791589473,
 114   0.0331456303681194152062895794737,
 115   0.0
 116 };
 117
 118 static const double g2_204[10] = { 0.0,
 119   0.0331456303681194152062895794737,
 120   0.0662912607362388304125791589473,
 121   -0.1767766952966368811002110905262,
 122   -0.4198446513295125926130013399998,
 123   0.9943689110435824561886873842099,
 124   -0.4198446513295125926130013399998,
 125   -0.1767766952966368811002110905262,
 126   0.0662912607362388304125791589473,
 127   0.0331456303681194152062895794737
 128 };
 129
 130 static const double h1_206[14] = { -0.0069053396600248781679769957237,
 131   0.0138106793200497563359539914474,
 132   0.0469563096881691715422435709210,
 133   -0.1077232986963880994204411332894,
 134   -0.1698713556366120029322340948025,
 135   0.4474660099696121052849093228945,
 136   0.9667475524034829435167794013152,
 137   0.4474660099696121052849093228945,
 138   -0.1698713556366120029322340948025,
 139   -0.1077232986963880994204411332894,
 140   0.0469563096881691715422435709210,
 141   0.0138106793200497563359539914474,
 142   -0.0069053396600248781679769957237,
 143   0.0
 144 };
 145
 146 static const double g2_206[14] = { 0.0,
 147   -0.0069053396600248781679769957237,
 148   -0.0138106793200497563359539914474,
 149   0.0469563096881691715422435709210,
 150   0.1077232986963880994204411332894,
 151   -0.1698713556366120029322340948025,
 152   -0.4474660099696121052849093228945,
 153   0.9667475524034829435167794013152,
 154   -0.4474660099696121052849093228945,
 155   -0.1698713556366120029322340948025,
 156   0.1077232986963880994204411332894,
 157   0.0469563096881691715422435709210,
 158   -0.0138106793200497563359539914474,
 159   -0.0069053396600248781679769957237,
 160 };
 161
 162 static const double h1_208[18] = { 0.0015105430506304420992449678146,
 163   -0.0030210861012608841984899356291,
 164   -0.0129475118625466465649568669819,
 165   0.0289161098263541773284036695929,
 166   0.0529984818906909399392234421792,
 167   -0.1349130736077360572068505539514,
 168   -0.1638291834340902345352542235443,
 169   0.4625714404759165262773590010400,
 170   0.9516421218971785225243297231697,
 171   0.4625714404759165262773590010400,
 172   -0.1638291834340902345352542235443,
 173   -0.1349130736077360572068505539514,
 174   0.0529984818906909399392234421792,
 175   0.0289161098263541773284036695929,
 176   -0.0129475118625466465649568669819,
 177   -0.0030210861012608841984899356291,
 178   0.0015105430506304420992449678146,
 179   0.0
 180 };
 181
 182 static const double g2_208[18] = { 0.0,
 183   0.0015105430506304420992449678146,
 184   0.0030210861012608841984899356291,
 185   -0.0129475118625466465649568669819,
 186   -0.0289161098263541773284036695929,
 187   0.0529984818906909399392234421792,
 188   0.1349130736077360572068505539514,
 189   -0.1638291834340902345352542235443,
 190   -0.4625714404759165262773590010400,
 191   0.9516421218971785225243297231697,
 192   -0.4625714404759165262773590010400,
 193   -0.1638291834340902345352542235443,
 194   0.1349130736077360572068505539514,
 195   0.0529984818906909399392234421792,
 196   -0.0289161098263541773284036695929,
 197   -0.0129475118625466465649568669819,
 198   0.0030210861012608841984899356291,
 199   0.0015105430506304420992449678146,
 200 };
 201
 202 static const double h2_2[18] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
 203   0.3535533905932737622004221810524,
 204   0.7071067811865475244008443621048,
 205   0.3535533905932737622004221810524,
 206   0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
 207 };
 208
 209 static const double g1_2[18] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
 210   -0.3535533905932737622004221810524,
 211   0.7071067811865475244008443621048,
 212   -0.3535533905932737622004221810524,
 213   0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
 214 };
 215
 216 static const double h1_301[4] = { -0.3535533905932737622004221810524,
 217   1.0606601717798212866012665431573,
 218   1.0606601717798212866012665431573,
 219   -0.3535533905932737622004221810524
 220 };
 221
 222 static const double g2_301[4] = { 0.3535533905932737622004221810524,
 223   1.0606601717798212866012665431573,
 224   -1.0606601717798212866012665431573,
 225   -0.3535533905932737622004221810524
 226 };
 227
 228 static const double h1_303[8] = { 0.0662912607362388304125791589473,
 229   -0.1988737822087164912377374768420,
 230   -0.1546796083845572709626847042104,
 231   0.9943689110435824561886873842099,
 232   0.9943689110435824561886873842099,
 233   -0.1546796083845572709626847042104,
 234   -0.1988737822087164912377374768420,
 235   0.0662912607362388304125791589473
 236 };
 237
 238 static const double g2_303[8] = { -0.0662912607362388304125791589473,
 239   -0.1988737822087164912377374768420,
 240   0.1546796083845572709626847042104,
 241   0.9943689110435824561886873842099,
 242   -0.9943689110435824561886873842099,
 243   -0.1546796083845572709626847042104,
 244   0.1988737822087164912377374768420,
 245   0.0662912607362388304125791589473
 246 };
 247
 248 static const double h1_305[12] = { -0.0138106793200497563359539914474,
 249   0.0414320379601492690078619743421,
 250   0.0524805814161890740766251675000,
 251   -0.2679271788089652729175074340788,
 252   -0.0718155324642587329469607555263,
 253   0.9667475524034829435167794013152,
 254   0.9667475524034829435167794013152,
 255   -0.0718155324642587329469607555263,
 256   -0.2679271788089652729175074340788,
 257   0.0524805814161890740766251675000,
 258   0.0414320379601492690078619743421,
 259   -0.0138106793200497563359539914474
 260 };
 261
 262 static const double g2_305[12] = { 0.0138106793200497563359539914474,
 263   0.0414320379601492690078619743421,
 264   -0.0524805814161890740766251675000,
 265   -0.2679271788089652729175074340788,
 266   0.0718155324642587329469607555263,
 267   0.9667475524034829435167794013152,
 268   -0.9667475524034829435167794013152,
 269   -0.0718155324642587329469607555263,
 270   0.2679271788089652729175074340788,
 271   0.0524805814161890740766251675000,
 272   -0.0414320379601492690078619743421,
 273   -0.0138106793200497563359539914474
 274 };
 275
 276 static const double h1_307[16] = { 0.0030210861012608841984899356291,
 277   -0.0090632583037826525954698068873,
 278   -0.0168317654213106405344439270765,
 279   0.0746639850740189951912512662623,
 280   0.0313329787073628846871956180962,
 281   -0.3011591259228349991008967259990,
 282   -0.0264992409453454699696117210896,
 283   0.9516421218971785225243297231697,
 284   0.9516421218971785225243297231697,
 285   -0.0264992409453454699696117210896,
 286   -0.3011591259228349991008967259990,
 287   0.0313329787073628846871956180962,
 288   0.0746639850740189951912512662623,
 289   -0.0168317654213106405344439270765,
 290   -0.0090632583037826525954698068873,
 291   0.0030210861012608841984899356291
 292 };
 293
 294 static const double g2_307[16] = { -0.0030210861012608841984899356291,
 295   -0.0090632583037826525954698068873,
 296   0.0168317654213106405344439270765,
 297   0.0746639850740189951912512662623,
 298   -0.0313329787073628846871956180962,
 299   -0.3011591259228349991008967259990,
 300   0.0264992409453454699696117210896,
 301   0.9516421218971785225243297231697,
 302   -0.9516421218971785225243297231697,
 303   -0.0264992409453454699696117210896,
 304   0.3011591259228349991008967259990,
 305   0.0313329787073628846871956180962,
 306   -0.0746639850740189951912512662623,
 307   -0.0168317654213106405344439270765,
 308   0.0090632583037826525954698068873,
 309   0.0030210861012608841984899356291
 310 };
 311
 312 static const double h1_309[20] = { -0.0006797443727836989446602355165,
 313   0.0020392331183510968339807065496,
 314   0.0050603192196119810324706421788,
 315   -0.0206189126411055346546938106687,
 316   -0.0141127879301758447558029850103,
 317   0.0991347824942321571990197448581,
 318   0.0123001362694193142367090236328,
 319   -0.3201919683607785695513833204624,
 320   0.0020500227115698857061181706055,
 321   0.9421257006782067372990864259380,
 322   0.9421257006782067372990864259380,
 323   0.0020500227115698857061181706055,
 324   -0.3201919683607785695513833204624,
 325   0.0123001362694193142367090236328,
 326   0.0991347824942321571990197448581,
 327   -0.0141127879301758447558029850103,
 328   -0.0206189126411055346546938106687,
 329   0.0050603192196119810324706421788,
 330   0.0020392331183510968339807065496,
 331   -0.0006797443727836989446602355165
 332 };
 333
 334 static const double g2_309[20] = { 0.0006797443727836989446602355165,
 335   0.0020392331183510968339807065496,
 336   -0.0050603192196119810324706421788,
 337   -0.0206189126411055346546938106687,
 338   0.0141127879301758447558029850103,
 339   0.0991347824942321571990197448581,
 340   -0.0123001362694193142367090236328,
 341   -0.3201919683607785695513833204624,
 342   -0.0020500227115698857061181706055,
 343   0.9421257006782067372990864259380,
 344   -0.9421257006782067372990864259380,
 345   0.0020500227115698857061181706055,
 346   0.3201919683607785695513833204624,
 347   0.0123001362694193142367090236328,
 348   -0.0991347824942321571990197448581,
 349   -0.0141127879301758447558029850103,
 350   0.0206189126411055346546938106687,
 351   0.0050603192196119810324706421788,
 352   -0.0020392331183510968339807065496,
 353   -0.0006797443727836989446602355165
 354 };
 355
 356 static const double h2_3[20] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
 357   0.1767766952966368811002110905262,
 358   0.5303300858899106433006332715786,
 359   0.5303300858899106433006332715786,
 360   0.1767766952966368811002110905262,
 361   0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
 362 };
 363
 364 static const double g1_3[20] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
 365   -0.1767766952966368811002110905262,
 366   0.5303300858899106433006332715786,
 367   -0.5303300858899106433006332715786,
 368   0.1767766952966368811002110905262,
 369   0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
 370 };
 371
 372 static int
 373 bspline_init (const double **h1, const double **g1,
 374               const double **h2, const double **g2, size_t * nc,
 375               size_t * offset, size_t member)
 376 {
 377   switch (member)
 378     {
 379     case 103:
 380       *nc = 6;
 381       *h1 = h1_103;
 382       *g1 = &g1_1[2];
 383       *h2 = &h2_1[2];
 384       *g2 = g2_103;
 385       break;
 386
 387     case 105:
 388       *nc = 10;
 389       *h1 = h1_105;
 390       *g1 = g1_1;
 391       *h2 = h2_1;
 392       *g2 = g2_105;
 393       break;
 394
 395     case 202:
 396       *nc = 6;
 397       *h1 = h1_202;
 398       *g1 = &g1_2[6];
 399       *h2 = &h2_2[6];
 400       *g2 = g2_202;
 401       break;
 402
 403     case 204:
 404       *nc = 10;
 405       *h1 = h1_204;
 406       *g1 = &g1_2[4];
 407       *h2 = &h2_2[4];
 408       *g2 = g2_204;
 409       break;
 410
 411     case 206:
 412       *nc = 14;
 413       *h1 = h1_206;
 414       *g1 = &g1_2[2];
 415       *h2 = &h2_2[2];
 416       *g2 = g2_206;
 417       break;
 418
 419     case 208:
 420       *nc = 18;
 421       *h1 = h1_208;
 422       *g1 = g1_2;
 423       *h2 = h2_2;
 424       *g2 = g2_208;
 425       break;
 426
 427     case 301:
 428       *nc = 4;
 429       *h1 = h1_301;
 430       *g1 = &g1_3[8];
 431       *h2 = &h2_3[8];
 432       *g2 = g2_301;
 433       break;
 434
 435     case 303:
 436       *nc = 8;
 437       *h1 = h1_303;
 438       *g1 = &g1_3[6];
 439       *h2 = &h2_3[6];
 440       *g2 = g2_303;
 441       break;
 442
 443     case 305:
 444       *nc = 12;
 445       *h1 = h1_305;
 446       *g1 = &g1_3[4];
 447       *h2 = &h2_3[4];
 448       *g2 = g2_305;
 449       break;
 450
 451     case 307:
 452       *nc = 16;
 453       *h1 = h1_307;
 454       *g1 = &g1_3[2];
 455       *h2 = &h2_3[2];
 456       *g2 = g2_307;
 457       break;
 458
 459     case 309:
 460       *nc = 20;
 461       *h1 = h1_309;
 462       *g1 = g1_3;
 463       *h2 = h2_3;
 464       *g2 = g2_309;
 465       break;
 466
 467     default:
 468       return GSL_FAILURE;
 469     }
 470
 471   *offset = 0;
 472
 473   return GSL_SUCCESS;
 474 }
 475
 476 static int
 477 bspline_centered_init (const double **h1, const double **g1,
 478                        const double **h2, const double **g2, size_t * nc,
 479                        size_t * offset, size_t member)
 480 {
 481   switch (member)
 482     {
 483     case 103:
 484       *nc = 6;
 485       *h1 = h1_103;
 486       *g1 = &g1_1[2];
 487       *h2 = &h2_1[2];
 488       *g2 = g2_103;
 489       break;
 490
 491     case 105:
 492       *nc = 10;
 493       *h1 = h1_105;
 494       *g1 = g1_1;
 495       *h2 = h2_1;
 496       *g2 = g2_105;
 497       break;
 498
 499     case 202:
 500       *nc = 6;
 501       *h1 = h1_202;
 502       *g1 = &g1_2[6];
 503       *h2 = &h2_2[6];
 504       *g2 = g2_202;
 505       break;
 506
 507     case 204:
 508       *nc = 10;
 509       *h1 = h1_204;
 510       *g1 = &g1_2[4];
 511       *h2 = &h2_2[4];
 512       *g2 = g2_204;
 513       break;
 514
 515     case 206:
 516       *nc = 14;
 517       *h1 = h1_206;
 518       *g1 = &g1_2[2];
 519       *h2 = &h2_2[2];
 520       *g2 = g2_206;
 521       break;
 522
 523     case 208:
 524       *nc = 18;
 525       *h1 = h1_208;
 526       *g1 = g1_2;
 527       *h2 = h2_2;
 528       *g2 = g2_208;
 529       break;
 530
 531     case 301:
 532       *nc = 4;
 533       *h1 = h1_301;
 534       *g1 = &g1_3[8];
 535       *h2 = &h2_3[8];
 536       *g2 = g2_301;
 537       break;
 538
 539     case 303:
 540       *nc = 8;
 541       *h1 = h1_303;
 542       *g1 = &g1_3[6];
 543       *h2 = &h2_3[6];
 544       *g2 = g2_303;
 545       break;
 546
 547     case 305:
 548       *nc = 12;
 549       *h1 = h1_305;
 550       *g1 = &g1_3[4];
 551       *h2 = &h2_3[4];
 552       *g2 = g2_305;
 553       break;
 554
 555     case 307:
 556       *nc = 16;
 557       *h1 = h1_307;
 558       *g1 = &g1_3[2];
 559       *h2 = &h2_3[2];
 560       *g2 = g2_307;
 561       break;
 562
 563     case 309:
 564       *nc = 20;
 565       *h1 = h1_309;
 566       *g1 = g1_3;
 567       *h2 = h2_3;
 568       *g2 = g2_309;
 569       break;
 570
 571     default:
 572       return GSL_FAILURE;
 573     }
 574
 575   *offset = ((*nc) >> 1);
 576
 577   return GSL_SUCCESS;
 578 }
 579
 580 static const gsl_wavelet_type bspline_type = {
 581   "bspline",
 582   &bspline_init
 583 };
 584
 585 static const gsl_wavelet_type bspline_centered_type = {
 586   "bspline-centered",
 587   &bspline_centered_init
 588 };
 589
 590 const gsl_wavelet_type *gsl_wavelet_bspline = &bspline_type;
 591 const gsl_wavelet_type *gsl_wavelet_bspline_centered = &bspline_centered_type;