libopenjpeg/dwt.c

   1 /*
   2  * Copyright (c) 2001-2002, David Janssens
   3  * Copyright (c) 2002-2004, Yannick Verschueren
   4  * Copyright (c) 2002-2004, Communications and remote sensing Laboratory, Universite catholique de Louvain, Belgium
   5  * All rights reserved.
   6  *
   7  * Redistribution and use in source and binary forms, with or without
   8  * modification, are permitted provided that the following conditions
   9  * are met:
  10  * 1. Redistributions of source code must retain the above copyright
  11  *    notice, this list of conditions and the following disclaimer.
  12  * 2. Redistributions in binary form must reproduce the above copyright
  13  *    notice, this list of conditions and the following disclaimer in the
  14  *    documentation and/or other materials provided with the distribution.
  15  *
  16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS `AS IS'
  17  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  18  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  19  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  20  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  21  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  22  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  23  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  24  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  25  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  26  * POSSIBILITY OF SUCH DAMAGE.
  27  */
  28
  29 #include "dwt.h"
  30 #include "int.h"
  31 #include "fix.h"
  32 #include "tcd.h"
  33 #include <stdlib.h>
  34 #include <stdio.h>
  35 //#include <math.h>
  36
  37 #define S(i) a[x*(i)*2]
  38 #define D(i) a[x*(1+(i)*2)]
  39 #define S_(i) ((i)<0?S(0):((i)>=sn?S(sn-1):S(i)))
  40 #define D_(i) ((i)<0?D(0):((i)>=dn?D(dn-1):D(i)))
  41 /* new */
  42 #define SS_(i) ((i)<0?S(0):((i)>=dn?S(dn-1):S(i)))
  43 #define DD_(i) ((i)<0?D(0):((i)>=sn?D(sn-1):D(i)))
  44
  45 /* <summary>                                                              */
  46 /* This table contains the norms of the 5-3 wavelets for different bands. */
  47 /* </summary>                                                             */
  48 double dwt_norms[4][10] = {
  49     {1.000, 1.500, 2.750, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
  50     {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
  51     {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
  52     {.7186, .9218, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93}
  53 };
  54
  55 /* <summary>                                                              */
  56 /* This table contains the norms of the 9-7 wavelets for different bands. */
  57 /* </summary>                                                             */
  58 double dwt_norms_real[4][10] = {
  59     {1.000, 1.965, 4.177, 8.403, 16.90, 33.84, 67.69, 135.3, 270.6, 540.9},
  60     {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
  61     {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
  62     {2.080, 3.865, 8.307, 17.18, 34.71, 69.59, 139.3, 278.6, 557.2}
  63 };
  64
  65 /* Add Patrick */
  66 static int *b = NULL;
  67 static int lastSizeOfB = 0;
  68
  69 /* <summary>       */
  70 /* Claning memory. */
  71 /* </summary>      */
  72
  73 void dwt_clean()
  74 {
  75     if (b != NULL) {
  76         free(b);
  77     }
  78     b = NULL;
  79     lastSizeOfB = 0;
  80 }
  81
  82 /* \ Add Patrick */
  83
  84 /* <summary>               */
  85 /* Forward lazy transform. */
  86 /* </summary>              */
  87 void dwt_deinterleave(int *a, int n, int x, int res, int cas)
  88 {
  89     int dn, sn, i;
  90     sn = res;
  91     dn = n - res;
  92     if (lastSizeOfB != n) {
  93         if (b != NULL)
  94             free(b);
  95         b = (int *) malloc(n * sizeof(int));
  96         lastSizeOfB = n;
  97     }
  98
  99     if (cas) {
 100         for (i = 0; i < sn; i++)
 101             b[i] = a[(2 * i + 1) * x];
 102         for (i = 0; i < dn; i++)
 103             b[sn + i] = a[2 * i * x];
 104     } else {
 105         for (i = 0; i < sn; i++)
 106             b[i] = a[2 * i * x];
 107         for (i = 0; i < dn; i++)
 108             b[sn + i] = a[(2 * i + 1) * x];
 109     }
 110     for (i = 0; i < n; i++)
 111         a[i * x] = b[i];
 112 }
 113
 114 /* <summary>               */
 115 /* Inverse lazy transform. */
 116 /* </summary>              */
 117 void dwt_interleave(int *a, int n, int x, int res, int cas)
 118 {
 119     int dn, sn, i;
 120     sn = res;
 121     dn = n - res;
 122
 123     if (lastSizeOfB != n) {
 124         if (b != NULL)
 125             free(b);
 126         b = (int *) malloc(n * sizeof(int));
 127         lastSizeOfB = n;
 128     }
 129
 130     if (cas) {
 131         for (i = 0; i < sn; i++)
 132             b[2 * i + 1] = a[i * x];
 133         for (i = 0; i < dn; i++)
 134             b[2 * i] = a[(sn + i) * x];
 135     } else {
 136         for (i = 0; i < sn; i++)
 137             b[2 * i] = a[i * x];
 138         for (i = 0; i < dn; i++)
 139             b[2 * i + 1] = a[(sn + i) * x];
 140     }
 141     for (i = 0; i < n; i++)
 142         a[i * x] = b[i];
 143 }
 144
 145 /* <summary>                            */
 146 /* Forward 5-3 wavelet tranform in 1-D. */
 147 /* </summary>                           */
 148 void dwt_encode_1(int *a, int n, int x, int res, int cas)
 149 {
 150     int dn, sn, i = 0;
 151     sn = res;
 152     dn = n - res;
 153
 154     if (cas) {
 155         if (!sn && dn == 1)     /* NEW :  CASE ONE ELEMENT */
 156             S(i) *= 2;
 157         else {
 158             for (i = 0; i < dn; i++)
 159                 S(i) -= (DD_(i) + DD_(i - 1)) >> 1;
 160             for (i = 0; i < sn; i++)
 161                 D(i) += (SS_(i) + SS_(i + 1) + 2) >> 2;
 162         }
 163     } else {
 164         if ((dn > 0) || (sn > 1)) {     /* NEW :  CASE ONE ELEMENT */
 165             for (i = 0; i < dn; i++)
 166                 D(i) -= (S_(i) + S_(i + 1)) >> 1;
 167             for (i = 0; i < sn; i++)
 168                 S(i) += (D_(i - 1) + D_(i) + 2) >> 2;
 169         }
 170     }
 171     dwt_deinterleave(a, n, x, res, cas);
 172 }
 173
 174 /* <summary>                            */
 175 /* Inverse 5-3 wavelet tranform in 1-D. */
 176 /* </summary>                           */
 177 void dwt_decode_1(int *a, int n, int x, int res, int cas)
 178 {
 179     int dn, sn, i = 0;
 180     sn = res;
 181     dn = n - res;
 182
 183     dwt_interleave(a, n, x, res, cas);
 184     if (cas) {
 185         if (!sn && dn == 1)     /* NEW :  CASE ONE ELEMENT */
 186             S(i) /= 2;
 187         else {
 188             for (i = 0; i < sn; i++)
 189                 D(i) -= (SS_(i) + SS_(i + 1) + 2) >> 2;
 190             for (i = 0; i < dn; i++)
 191                 S(i) += (DD_(i) + DD_(i - 1)) >> 1;
 192         }
 193     } else {
 194         if ((dn > 0) || (sn > 1)) {     /* NEW :  CASE ONE ELEMENT */
 195             for (i = 0; i < sn; i++)
 196                 S(i) -= (D_(i - 1) + D_(i) + 2) >> 2;
 197             for (i = 0; i < dn; i++)
 198                 D(i) += (S_(i) + S_(i + 1)) >> 1;
 199         }
 200     }
 201 }
 202
 203 /* <summary>                            */
 204 /* Forward 5-3 wavelet tranform in 2-D. */
 205 /* </summary>                           */
 206 void dwt_encode(int *a, int w, int h, tcd_tilecomp_t * tilec, int l)
 207 {
 208     int i, j;
 209     int rw;                     /* width of the resolution level computed                                                           */
 210     int rh;                     /* heigth of the resolution level computed                                                          */
 211     int rw1;                    /* width of the resolution level once lower than computed one                                       */
 212     int rh1;                    /* height of the resolution level once lower than computed one                                      */
 213
 214     for (i = 0; i < l; i++) {
 215         int cas_col = 0;        /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 216         int cas_row = 0;        /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 217         rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
 218         rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
 219         rw1 =
 220             tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i -
 221                                                                   1].x0;
 222         rh1 =
 223             tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i -
 224                                                                   1].y0;
 225
 226         cas_row = tilec->resolutions[l - i].x0 % 2;
 227         cas_col = tilec->resolutions[l - i].y0 % 2;
 228
 229         for (j = 0; j < rw; j++)
 230             dwt_encode_1(a + j, rh, w, rh1, cas_col);
 231         for (j = 0; j < rh; j++)
 232             dwt_encode_1(a + j * w, rw, 1, rw1, cas_row);
 233     }
 234
 235     dwt_clean();
 236 }
 237
 238 /* <summary>                            */
 239 /* Inverse 5-3 wavelet tranform in 2-D. */
 240 /* </summary>                           */
 241 void dwt_decode(int *a, int w, int h, tcd_tilecomp_t * tilec, int l,
 242                 int stop)
 243 {
 244     int i, j;
 245     int rw;                     /* width of the resolution level computed                                                           */
 246     int rh;                     /* heigth of the resolution level computed                                                          */
 247     int rw1;                    /* width of the resolution level once lower than computed one                                       */
 248     int rh1;                    /* height of the resolution level once lower than computed one                                      */
 249
 250     for (i = l - 1; i >= stop; i--) {
 251         int cas_col = 0;        /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 252         int cas_row = 0;        /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 253
 254         rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
 255         rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
 256         rw1 =
 257             tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i -
 258                                                                   1].x0;
 259         rh1 =
 260             tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i -
 261                                                                   1].y0;
 262
 263         cas_row = tilec->resolutions[l - i].x0 % 2;
 264         cas_col = tilec->resolutions[l - i].y0 % 2;
 265
 266         for (j = 0; j < rh; j++)
 267             dwt_decode_1(a + j * w, rw, 1, rw1, cas_row);
 268         for (j = 0; j < rw; j++)
 269             dwt_decode_1(a + j, rh, w, rh1, cas_col);
 270     }
 271     dwt_clean();
 272 }
 273
 274 /* <summary>                          */
 275 /* Get gain of 5-3 wavelet transform. */
 276 /* </summary>                         */
 277 int dwt_getgain(int orient)
 278 {
 279     if (orient == 0)
 280         return 0;
 281     if (orient == 1 || orient == 2)
 282         return 1;
 283     return 2;
 284 }
 285
 286 /* <summary>                */
 287 /* Get norm of 5-3 wavelet. */
 288 /* </summary>               */
 289 double dwt_getnorm(int level, int orient)
 290 {
 291     return dwt_norms[orient][level];
 292 }
 293
 294 /* <summary>                             */
 295 /* Forward 9-7 wavelet transform in 1-D. */
 296 /* </summary>                            */
 297 void dwt_encode_1_real(int *a, int n, int x, int res, int cas)
 298 {
 299     int dn, sn, i = 0;
 300     dn = n - res;
 301     sn = res;
 302
 303     if (cas) {
 304         if ((sn > 0) || (dn > 1)) {     /* NEW :  CASE ONE ELEMENT */
 305             for (i = 0; i < dn; i++)
 306                 S(i) -= fix_mul(DD_(i) + DD_(i - 1), 12993);
 307             for (i = 0; i < sn; i++)
 308                 D(i) -= fix_mul(SS_(i) + SS_(i + 1), 434);
 309             for (i = 0; i < dn; i++)
 310                 S(i) += fix_mul(DD_(i) + DD_(i - 1), 7233);
 311             for (i = 0; i < sn; i++)
 312                 D(i) += fix_mul(SS_(i) + SS_(i + 1), 3633);
 313             for (i = 0; i < dn; i++)
 314                 S(i) = fix_mul(S(i), 5038);     /*5038 */
 315             for (i = 0; i < sn; i++)
 316                 D(i) = fix_mul(D(i), 6659);     /*6660 */
 317         }
 318     } else {
 319         if ((dn > 0) || (sn > 1)) {     /* NEW :  CASE ONE ELEMENT */
 320             for (i = 0; i < dn; i++)
 321                 D(i) -= fix_mul(S_(i) + S_(i + 1), 12993);
 322             for (i = 0; i < sn; i++)
 323                 S(i) -= fix_mul(D_(i - 1) + D_(i), 434);
 324             for (i = 0; i < dn; i++)
 325                 D(i) += fix_mul(S_(i) + S_(i + 1), 7233);
 326             for (i = 0; i < sn; i++)
 327                 S(i) += fix_mul(D_(i - 1) + D_(i), 3633);
 328             for (i = 0; i < dn; i++)
 329                 D(i) = fix_mul(D(i), 5038);     /*5038 */
 330             for (i = 0; i < sn; i++)
 331                 S(i) = fix_mul(S(i), 6659);     /*6660 */
 332         }
 333     }
 334     dwt_deinterleave(a, n, x, res, cas);
 335 }
 336
 337 /* <summary>                             */
 338 /* Inverse 9-7 wavelet transform in 1-D. */
 339 /* </summary>                            */
 340 void dwt_decode_1_real(int *a, int n, int x, int res, int cas)
 341 {
 342     int dn, sn, i = 0;
 343     dn = n - res;
 344     sn = res;
 345     dwt_interleave(a, n, x, res, cas);
 346     if (cas) {
 347         if ((sn > 0) || (dn > 1)) {     /* NEW :  CASE ONE ELEMENT */
 348             for (i = 0; i < sn; i++)
 349                 D(i) = fix_mul(D(i), 10078);    /* 10076 */
 350             for (i = 0; i < dn; i++)
 351                 S(i) = fix_mul(S(i), 13318);    /* 13320 */
 352             for (i = 0; i < sn; i++)
 353                 D(i) -= fix_mul(SS_(i) + SS_(i + 1), 3633);
 354             for (i = 0; i < dn; i++)
 355                 S(i) -= fix_mul(DD_(i) + DD_(i - 1), 7233);
 356             for (i = 0; i < sn; i++)
 357                 D(i) += fix_mul(SS_(i) + SS_(i + 1), 434);
 358             for (i = 0; i < dn; i++)
 359                 S(i) += fix_mul(DD_(i) + DD_(i - 1), 12993);
 360         }
 361     } else {
 362         if ((dn > 0) || (sn > 1)) {     /* NEW :  CASE ONE ELEMENT */
 363             for (i = 0; i < sn; i++)
 364                 S(i) = fix_mul(S(i), 10078);    /* 10076 */
 365             for (i = 0; i < dn; i++)
 366                 D(i) = fix_mul(D(i), 13318);    /* 13320 */
 367             for (i = 0; i < sn; i++)
 368                 S(i) -= fix_mul(D_(i - 1) + D_(i), 3633);
 369             for (i = 0; i < dn; i++)
 370                 D(i) -= fix_mul(S_(i) + S_(i + 1), 7233);
 371             for (i = 0; i < sn; i++)
 372                 S(i) += fix_mul(D_(i - 1) + D_(i), 434);
 373             for (i = 0; i < dn; i++)
 374                 D(i) += fix_mul(S_(i) + S_(i + 1), 12993);
 375         }
 376     }
 377 }
 378
 379 /* <summary>                             */
 380 /* Forward 9-7 wavelet transform in 2-D. */
 381 /* </summary>                            */
 382
 383 void dwt_encode_real(int *a, int w, int h, tcd_tilecomp_t * tilec, int l)
 384 {
 385     int i, j;
 386     int rw;                     /* width of the resolution level computed                                                     */
 387     int rh;                     /* heigth of the resolution level computed                                                    */
 388     int rw1;                    /* width of the resolution level once lower than computed one                                 */
 389     int rh1;                    /* height of the resolution level once lower than computed one                                */
 390
 391     for (i = 0; i < l; i++) {
 392         int cas_col = 0;        /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 393         int cas_row = 0;        /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 394         rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
 395         rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
 396         rw1 =
 397             tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i -
 398                                                                   1].x0;
 399         rh1 =
 400             tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i -
 401                                                                   1].y0;
 402
 403         cas_row = tilec->resolutions[l - i].x0 % 2;
 404         cas_col = tilec->resolutions[l - i].y0 % 2;
 405
 406         for (j = 0; j < rw; j++)
 407             dwt_encode_1_real(a + j, rh, w, rh1, cas_col);
 408         for (j = 0; j < rh; j++)
 409             dwt_encode_1_real(a + j * w, rw, 1, rw1, cas_row);
 410     }
 411 }
 412
 413 /* <summary>                             */
 414 /* Inverse 9-7 wavelet transform in 2-D. */
 415 /* </summary>                            */
 416 void dwt_decode_real(int *a, int w, int h, tcd_tilecomp_t * tilec, int l,
 417                      int stop)
 418 {
 419     int i, j;
 420     int rw;                     /* width of the resolution level computed                                                           */
 421     int rh;                     /* heigth of the resolution level computed                                                          */
 422     int rw1;                    /* width of the resolution level once lower than computed one                                       */
 423     int rh1;                    /* height of the resolution level once lower than computed one                                      */
 424
 425     for (i = l - 1; i >= stop; i--) {
 426         int cas_col = 0;        /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 427         int cas_row = 0;        /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 428
 429         rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
 430         rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
 431         rw1 =
 432             tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i -
 433                                                                   1].x0;
 434         rh1 =
 435             tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i -
 436                                                                   1].y0;
 437
 438         cas_row = tilec->resolutions[l - i].x0 % 2;
 439         cas_col = tilec->resolutions[l - i].y0 % 2;
 440
 441         for (j = 0; j < rh; j++)
 442             dwt_decode_1_real(a + j * w, rw, 1, rw1, cas_row);
 443         for (j = 0; j < rw; j++)
 444             dwt_decode_1_real(a + j, rh, w, rh1, cas_col);
 445     }
 446 }
 447
 448 /* <summary>                          */
 449 /* Get gain of 9-7 wavelet transform. */
 450 /* </summary>                         */
 451 int dwt_getgain_real(int orient)
 452 {
 453     return 0;
 454 }
 455
 456 /* <summary>                */
 457 /* Get norm of 9-7 wavelet. */
 458 /* </summary>               */
 459 double dwt_getnorm_real(int level, int orient)
 460 {
 461     return dwt_norms_real[orient][level];
 462 }