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 = tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
 220                 rh1 = tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
 221
 222                 cas_row = tilec->resolutions[l - i].x0 % 2;
 223                 cas_col = tilec->resolutions[l - i].y0 % 2;
 224
 225                 for (j = 0; j < rw; j++)
 226                         dwt_encode_1(a + j, rh, w, rh1, cas_col);
 227                 for (j = 0; j < rh; j++)
 228                         dwt_encode_1(a + j * w, rw, 1, rw1, cas_row);
 229         }
 230
 231         dwt_clean();
 232 }
 233
 234 /* <summary>                            */
 235 /* Inverse 5-3 wavelet tranform in 2-D. */
 236 /* </summary>                           */
 237 void dwt_decode(int *a, int w, int h, tcd_tilecomp_t * tilec, int l, int stop)
 238 {
 239         int i, j;
 240         int rw;            /* width of the resolution level computed                                                           */
 241         int rh;            /* heigth of the resolution level computed                                                          */
 242         int rw1;           /* width of the resolution level once lower than computed one                                       */
 243         int rh1;           /* height of the resolution level once lower than computed one                                      */
 244
 245         for (i = l - 1; i >= stop; i--) {
 246           int cas_col = 0; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 247           int cas_row = 0; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 248
 249                 rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
 250                 rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
 251                 rw1 = tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
 252                 rh1 = tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
 253
 254                 cas_row = tilec->resolutions[l - i].x0 % 2;
 255                 cas_col = tilec->resolutions[l - i].y0 % 2;
 256
 257                 for (j = 0; j < rh; j++)
 258                         dwt_decode_1(a + j * w, rw, 1, rw1, cas_row);
 259                 for (j = 0; j < rw; j++)
 260                         dwt_decode_1(a + j, rh, w, rh1, cas_col);
 261         }
 262         dwt_clean();
 263 }
 264
 265 /* <summary>                          */
 266 /* Get gain of 5-3 wavelet transform. */
 267 /* </summary>                         */
 268 int dwt_getgain(int orient)
 269 {
 270         if (orient == 0)
 271                 return 0;
 272         if (orient == 1 || orient == 2)
 273                 return 1;
 274         return 2;
 275 }
 276
 277 /* <summary>                */
 278 /* Get norm of 5-3 wavelet. */
 279 /* </summary>               */
 280 double dwt_getnorm(int level, int orient)
 281 {
 282         return dwt_norms[orient][level];
 283 }
 284
 285 /* <summary>                             */
 286 /* Forward 9-7 wavelet transform in 1-D. */
 287 /* </summary>                            */
 288 void dwt_encode_1_real(int *a, int n, int x, int res, int cas)
 289 {
 290         int dn, sn, i = 0;
 291         dn = n - res;
 292         sn = res;
 293
 294         if (cas) {
 295                 if ((sn > 0) || (dn > 1)) {  /* NEW :  CASE ONE ELEMENT */
 296                         for (i = 0; i < dn; i++)
 297                                 S(i) -= fix_mul(DD_(i) + DD_(i - 1), 12993);
 298                         for (i = 0; i < sn; i++)
 299                                 D(i) -= fix_mul(SS_(i) + SS_(i + 1), 434);
 300                         for (i = 0; i < dn; i++)
 301                                 S(i) += fix_mul(DD_(i) + DD_(i - 1), 7233);
 302                         for (i = 0; i < sn; i++)
 303                                 D(i) += fix_mul(SS_(i) + SS_(i + 1), 3633);
 304                         for (i = 0; i < dn; i++)
 305                           S(i) = fix_mul(S(i), 5038); /*5038*/
 306                         for (i = 0; i < sn; i++)
 307                           D(i) = fix_mul(D(i), 6659); /*6660*/
 308                 }
 309         } else {
 310           if ((dn > 0) || (sn > 1)) {  /* NEW :  CASE ONE ELEMENT */
 311                         for (i = 0; i < dn; i++)
 312                                 D(i) -= fix_mul(S_(i) + S_(i + 1), 12993);
 313                         for (i = 0; i < sn; i++)
 314                                 S(i) -= fix_mul(D_(i - 1) + D_(i), 434);
 315                         for (i = 0; i < dn; i++)
 316                                 D(i) += fix_mul(S_(i) + S_(i + 1), 7233);
 317                         for (i = 0; i < sn; i++)
 318                                 S(i) += fix_mul(D_(i - 1) + D_(i), 3633);
 319                         for (i = 0; i < dn; i++)
 320                                 D(i) = fix_mul(D(i), 5038); /*5038*/
 321                         for (i = 0; i < sn; i++)
 322                                 S(i) = fix_mul(S(i), 6659); /*6660*/
 323                 }
 324         }
 325         dwt_deinterleave(a, n, x, res, cas);
 326 }
 327
 328 /* <summary>                             */
 329 /* Inverse 9-7 wavelet transform in 1-D. */
 330 /* </summary>                            */
 331 void dwt_decode_1_real(int *a, int n, int x, int res, int cas)
 332 {
 333         int dn, sn, i = 0;
 334         dn = n - res;
 335         sn = res;
 336         dwt_interleave(a, n, x, res, cas);
 337         if (cas) {
 338                 if ((sn > 0) || (dn > 1)) {  /* NEW :  CASE ONE ELEMENT */
 339                         for (i = 0; i < sn; i++)
 340                           D(i) = fix_mul(D(i), 10078); /* 10076 */
 341                         for (i = 0; i < dn; i++)
 342                                 S(i) = fix_mul(S(i), 13318); /* 13320*/
 343                         for (i = 0; i < sn; i++)
 344                                 D(i) -= fix_mul(SS_(i) + SS_(i + 1), 3633);
 345                         for (i = 0; i < dn; i++)
 346                                 S(i) -= fix_mul(DD_(i) + DD_(i - 1), 7233);
 347                         for (i = 0; i < sn; i++)
 348                                 D(i) += fix_mul(SS_(i) + SS_(i + 1), 434);
 349                         for (i = 0; i < dn; i++)
 350                                 S(i) += fix_mul(DD_(i) + DD_(i - 1), 12993);
 351                 }
 352         } else {
 353                 if ((dn > 0) || (sn > 1)) {  /* NEW :  CASE ONE ELEMENT */
 354                         for (i = 0; i < sn; i++)
 355                                 S(i) = fix_mul(S(i), 10078); /* 10076 */
 356                         for (i = 0; i < dn; i++)
 357                           D(i) = fix_mul(D(i), 13318); /* 13320*/
 358                         for (i = 0; i < sn; i++)
 359                                 S(i) -= fix_mul(D_(i - 1) + D_(i), 3633);
 360                         for (i = 0; i < dn; i++)
 361                                 D(i) -= fix_mul(S_(i) + S_(i + 1), 7233);
 362                         for (i = 0; i < sn; i++)
 363                                 S(i) += fix_mul(D_(i - 1) + D_(i), 434);
 364                         for (i = 0; i < dn; i++)
 365                                 D(i) += fix_mul(S_(i) + S_(i + 1), 12993);
 366                 }
 367         }
 368 }
 369
 370 /* <summary>                             */
 371 /* Forward 9-7 wavelet transform in 2-D. */
 372 /* </summary>                            */
 373
 374 void dwt_encode_real(int *a, int w, int h, tcd_tilecomp_t * tilec, int l)
 375 {
 376         int i, j;
 377         int rw;            /* width of the resolution level computed                                                     */
 378         int rh;            /* heigth of the resolution level computed                                                    */
 379         int rw1;           /* width of the resolution level once lower than computed one                                 */
 380         int rh1;           /* height of the resolution level once lower than computed one                                */
 381
 382         for (i = 0; i < l; i++) {
 383                 int cas_col = 0;  /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 384                 int cas_row = 0;  /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 385                 rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
 386                 rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
 387                 rw1 = tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
 388                 rh1 = tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
 389
 390                 cas_row = tilec->resolutions[l - i].x0 % 2;
 391                 cas_col = tilec->resolutions[l - i].y0 % 2;
 392
 393                 for (j = 0; j < rw; j++)
 394                         dwt_encode_1_real(a + j, rh, w, rh1, cas_col);
 395                 for (j = 0; j < rh; j++)
 396                         dwt_encode_1_real(a + j * w, rw, 1, rw1, cas_row);
 397         }
 398 }
 399
 400 /* <summary>                             */
 401 /* Inverse 9-7 wavelet transform in 2-D. */
 402 /* </summary>                            */
 403 void dwt_decode_real(int *a, int w, int h, tcd_tilecomp_t * tilec, int l, int stop)
 404 {
 405         int i, j;
 406         int rw;            /* width of the resolution level computed                                                           */
 407         int rh;            /* heigth of the resolution level computed                                                          */
 408         int rw1;           /* width of the resolution level once lower than computed one                                       */
 409         int rh1;           /* height of the resolution level once lower than computed one                                      */
 410
 411         for (i = l - 1; i >= stop; i--) {
 412           int cas_col = 0; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 413           int cas_row = 0; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 414
 415                 rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
 416                 rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
 417                 rw1 = tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
 418                 rh1 = tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
 419
 420                 cas_row = tilec->resolutions[l - i].x0 % 2;
 421                 cas_col = tilec->resolutions[l - i].y0 % 2;
 422
 423                 for (j = 0; j < rh; j++)
 424                         dwt_decode_1_real(a + j * w, rw, 1, rw1, cas_row);
 425                 for (j = 0; j < rw; j++)
 426                         dwt_decode_1_real(a + j, rh, w, rh1, cas_col);
 427         }
 428 }
 429
 430 /* <summary>                          */
 431 /* Get gain of 9-7 wavelet transform. */
 432 /* </summary>                         */
 433 int dwt_getgain_real(int orient)
 434 {
 435         return 0;
 436 }
 437
 438 /* <summary>                */
 439 /* Get norm of 9-7 wavelet. */
 440 /* </summary>               */
 441 double dwt_getnorm_real(int level, int orient)
 442 {
 443         return dwt_norms_real[orient][level];
 444 }