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