src/gdcmopenjpeg/libopenjpeg/pi.c

   1 /*
   2  * Copyright (c) 2001-2002, David Janssens
   3  * Copyright (c) 2003-2004, Yannick Verschueren
   4  * Copyright (c) 2003-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 "pi.h"
  30 #include "int.h"
  31 #include <stdlib.h>
  32 #include <stdio.h>
  33
  34 /* <summary>
  35  * Create a packet iterator.
  36  * </summary> */
  37 pi_iterator_t *pi_create(j2k_image_t * img, j2k_cp_t * cp, int tileno)
  38 {
  39   int p, q, i;
  40   int compno, resno, pino;
  41   int maxres = 0;
  42   pi_iterator_t *pi;
  43   j2k_tcp_t *tcp;
  44   j2k_tccp_t *tccp;
  45
  46   tcp = &cp->tcps[tileno];
  47   pi = (pi_iterator_t *) malloc((tcp->numpocs + 1) *
  48             sizeof(pi_iterator_t));
  49
  50   for (pino = 0; pino < tcp->numpocs + 1; pino++) {   /* change */
  51     p = tileno % cp->tw;
  52     q = tileno / cp->tw;
  53
  54     pi[pino].tx0 = int_max(cp->tx0 + p * cp->tdx, img->x0);
  55     pi[pino].ty0 = int_max(cp->ty0 + q * cp->tdy, img->y0);
  56     pi[pino].tx1 = int_min(cp->tx0 + (p + 1) * cp->tdx, img->x1);
  57     pi[pino].ty1 = int_min(cp->ty0 + (q + 1) * cp->tdy, img->y1);
  58     pi[pino].numcomps = img->numcomps;
  59     pi[pino].comps =
  60       (pi_comp_t *) malloc(img->numcomps * sizeof(pi_comp_t));
  61
  62     for (compno = 0; compno < pi->numcomps; compno++) {
  63       int tcx0, tcy0, tcx1, tcy1;
  64       pi_comp_t *comp = &pi[pino].comps[compno];
  65       tccp = &tcp->tccps[compno];
  66       comp->dx = img->comps[compno].dx;
  67       comp->dy = img->comps[compno].dy;
  68       comp->numresolutions = tccp->numresolutions;
  69       comp->resolutions =
  70    (pi_resolution_t *) malloc(comp->numresolutions *
  71                sizeof(pi_resolution_t));
  72       tcx0 = int_ceildiv(pi->tx0, comp->dx);
  73       tcy0 = int_ceildiv(pi->ty0, comp->dy);
  74       tcx1 = int_ceildiv(pi->tx1, comp->dx);
  75       tcy1 = int_ceildiv(pi->ty1, comp->dy);
  76       if (comp->numresolutions > maxres) {
  77    maxres = comp->numresolutions;
  78       }
  79       for (resno = 0; resno < comp->numresolutions; resno++) {
  80    int levelno;
  81    int rx0, ry0, rx1, ry1;
  82    int px0, py0, px1, py1;
  83    pi_resolution_t *res = &comp->resolutions[resno];
  84    if (tccp->csty & J2K_CCP_CSTY_PRT) {
  85      res->pdx = tccp->prcw[resno];
  86      res->pdy = tccp->prch[resno];
  87    } else {
  88      res->pdx = 15;
  89      res->pdy = 15;
  90    }
  91    levelno = comp->numresolutions - 1 - resno;
  92    rx0 = int_ceildivpow2(tcx0, levelno);
  93    ry0 = int_ceildivpow2(tcy0, levelno);
  94    rx1 = int_ceildivpow2(tcx1, levelno);
  95    ry1 = int_ceildivpow2(tcy1, levelno);
  96    px0 = int_floordivpow2(rx0, res->pdx) << res->pdx;
  97    py0 = int_floordivpow2(ry0, res->pdy) << res->pdy;
  98    px1 = int_ceildivpow2(rx1, res->pdx) << res->pdx;
  99    py1 = int_ceildivpow2(ry1, res->pdy) << res->pdy;
 100    res->pw = (rx0==rx1)?0:((px1 - px0) >> res->pdx); /*Mod Antonin : sizebug1*/
 101    res->ph = (ry0==ry1)?0:((py1 - py0) >> res->pdy); /*Mod Antonin : sizebug1*/
 102       }
 103     }
 104
 105     tccp = &tcp->tccps[0];
 106     pi[pino].step_p = 1;
 107     pi[pino].step_c = 100 * pi[pino].step_p;
 108     pi[pino].step_r = img->numcomps * pi[pino].step_c;
 109     pi[pino].step_l = maxres * pi[pino].step_r;
 110
 111     if (pino == 0) {
 112       pi[pino].include =
 113    (short int *) malloc(img->numcomps * maxres *
 114               tcp->numlayers * 100 * sizeof(short int));
 115       for (i = 0; i < img->numcomps * maxres * tcp->numlayers * 100; i++)
 116    pi[pino].include[i] = 0;
 117     }
 118     /* pi[pino].include=(short int*)calloc(img->numcomps*maxres*tcp->numlayers*1000,sizeof(short int)); */
 119     else
 120       pi[pino].include = pi[pino - 1].include;
 121
 122     if (tcp->POC == 0) {
 123       pi[pino].first = 1;
 124       pi[pino].poc.resno0 = 0;
 125       pi[pino].poc.compno0 = 0;
 126       pi[pino].poc.layno1 = tcp->numlayers;
 127       pi[pino].poc.resno1 = maxres;
 128       pi[pino].poc.compno1 = img->numcomps;
 129       pi[pino].poc.prg = tcp->prg;
 130     } else {
 131       pi[pino].first = 1;
 132       pi[pino].poc.resno0 = tcp->pocs[pino].resno0;
 133       pi[pino].poc.compno0 = tcp->pocs[pino].compno0;
 134       pi[pino].poc.layno1 = tcp->pocs[pino].layno1;
 135       pi[pino].poc.resno1 = tcp->pocs[pino].resno1;
 136       pi[pino].poc.compno1 = tcp->pocs[pino].compno1;
 137       pi[pino].poc.prg = tcp->pocs[pino].prg;
 138     }
 139   }
 140   return pi;
 141 }
 142
 143 /* <summary>
 144  * Get next packet in layer-resolution-component-precinct order.
 145  *
 146  * pi: packet iterator to modify
 147  * </summary> */
 148 int pi_next_lrcp(pi_iterator_t * pi)
 149 {
 150   pi_comp_t *comp;
 151   pi_resolution_t *res;
 152
 153   if (!pi->first) {
 154     comp = &pi->comps[pi->compno];
 155     res = &comp->resolutions[pi->resno];
 156     goto skip;
 157   } else {
 158     pi->first = 0;
 159   }
 160   for (pi->layno = 0; pi->layno < pi->poc.layno1; pi->layno++) {
 161     for (pi->resno = pi->poc.resno0; pi->resno < pi->poc.resno1;
 162     pi->resno++) {
 163       for (pi->compno = pi->poc.compno0;
 164       pi->compno < pi->poc.compno1; pi->compno++) {
 165    comp = &pi->comps[pi->compno];
 166    if (pi->resno >= comp->numresolutions) {
 167      continue;
 168    }
 169    res = &comp->resolutions[pi->resno];
 170    for (pi->precno = 0; pi->precno < res->pw * res->ph; pi->precno++) {
 171      if (!pi->
 172          include[pi->layno * pi->step_l +
 173             pi->resno * pi->step_r +
 174             pi->compno * pi->step_c + pi->precno * pi->step_p]) {
 175        pi->include[pi->layno * pi->step_l +
 176          pi->resno * pi->step_r +
 177          pi->compno * pi->step_c +
 178          pi->precno * pi->step_p] = 1;
 179        return 1;
 180      }
 181    skip:;
 182    }
 183       }
 184     }
 185   }
 186   return 0;
 187 }
 188
 189 /* <summary>
 190  * Get next packet in resolution-layer-component-precinct order.
 191  *
 192  * pi: packet iterator to modify
 193  * </summary> */
 194 int pi_next_rlcp(pi_iterator_t * pi)
 195 {
 196   pi_comp_t *comp;
 197   pi_resolution_t *res;
 198   if (!pi->first) {
 199     comp = &pi->comps[pi->compno];
 200     res = &comp->resolutions[pi->resno];
 201     goto skip;
 202   } else {
 203     pi->first = 0;
 204   }
 205   for (pi->resno = pi->poc.resno0; pi->resno < pi->poc.resno1; pi->resno++) {
 206     for (pi->layno = 0; pi->layno < pi->poc.layno1; pi->layno++) {
 207       for (pi->compno = pi->poc.compno0;
 208       pi->compno < pi->poc.compno1; pi->compno++) {
 209    comp = &pi->comps[pi->compno];
 210    if (pi->resno >= comp->numresolutions) {
 211      continue;
 212    }
 213    res = &comp->resolutions[pi->resno];
 214    for (pi->precno = 0; pi->precno < res->pw * res->ph; pi->precno++) {
 215      if (!pi->
 216          include[pi->layno * pi->step_l +
 217             pi->resno * pi->step_r +
 218             pi->compno * pi->step_c + pi->precno * pi->step_p]) {
 219        pi->include[pi->layno * pi->step_l +
 220          pi->resno * pi->step_r +
 221          pi->compno * pi->step_c +
 222          pi->precno * pi->step_p] = 1;
 223        return 1;
 224      }
 225    skip:;
 226    }
 227       }
 228     }
 229   }
 230   return 0;
 231 }
 232
 233 /* <summary>
 234  * Get next packet in resolution-precinct-component-layer order.
 235  *
 236  * pi: packet iterator to modify
 237  * </summary> */
 238 int pi_next_rpcl(pi_iterator_t * pi)
 239 {
 240   pi_comp_t *comp;
 241   pi_resolution_t *res;
 242   if (!pi->first) {
 243     goto skip;
 244   } else {
 245     int compno, resno;
 246     pi->first = 0;
 247     pi->dx = 0;
 248     pi->dy = 0;
 249     for (compno = 0; compno < pi->numcomps; compno++) {
 250       comp = &pi->comps[compno];
 251       for (resno = 0; resno < comp->numresolutions; resno++) {
 252    int dx, dy;
 253    res = &comp->resolutions[resno];
 254    dx = comp->dx *
 255      (1 << (res->pdx + comp->numresolutions - 1 - resno));
 256    dy = comp->dy *
 257      (1 << (res->pdy + comp->numresolutions - 1 - resno));
 258    pi->dx = !pi->dx ? dx : int_min(pi->dx, dx);
 259    pi->dy = !pi->dy ? dy : int_min(pi->dy, dy);
 260       }
 261     }
 262   }
 263   for (pi->resno = pi->poc.resno0; pi->resno < pi->poc.resno1; pi->resno++) {
 264     for (pi->y = pi->ty0; pi->y < pi->ty1;
 265     pi->y += pi->dy - (pi->y % pi->dy)) {
 266       for (pi->x = pi->tx0; pi->x < pi->tx1;
 267       pi->x += pi->dx - (pi->x % pi->dx)) {
 268    for (pi->compno = pi->poc.compno0;
 269         pi->compno < pi->poc.compno1; pi->compno++) {
 270      int levelno;
 271      int trx0, try0;
 272      int trx1, try1;/* Add antonin pcrl*/
 273      int rpx, rpy;
 274      int prci, prcj;
 275      comp = &pi->comps[pi->compno];
 276      if (pi->resno >= comp->numresolutions) {
 277        continue;
 278      }
 279      res = &comp->resolutions[pi->resno];
 280      levelno = comp->numresolutions - 1 - pi->resno;
 281      trx0 = int_ceildiv(pi->tx0, comp->dx << levelno);
 282      try0 = int_ceildiv(pi->ty0, comp->dy << levelno);
 283      trx1 = int_ceildiv(pi->tx1, comp->dx << levelno);/* Add antonin pcrl*/
 284      try1 = int_ceildiv(pi->ty1, comp->dy << levelno);/* Add antonin pcrl*/
 285      rpx = res->pdx + levelno;
 286      rpy = res->pdy + levelno;
 287      if (!
 288          (pi->x % (comp->dx << rpx) == 0
 289           || (pi->x == pi->tx0 && (trx0 << levelno) % (1 << rpx)))) {
 290        continue;
 291      }
 292      if (!
 293          (pi->y % (comp->dy << rpy) == 0
 294           || (pi->y == pi->ty0 && (try0 << levelno) % (1 << rpx)))) {
 295        continue;
 296      }
 297
 298      /*Add Antonin : sizebug1*/
 299      if ((res->pw==0)||(res->pw==0)) continue;
 300      /*ddA*/
 301
 302      /*Add Antonin : pcrl*/
 303      if ((trx0==trx1)||(try0==try1)) continue;
 304      /*ddA*/
 305
 306      prci =
 307        int_floordivpow2(int_ceildiv
 308               (pi->x, comp->dx << levelno),
 309               res->pdx) - int_floordivpow2(trx0, res->pdx);
 310      prcj =
 311        int_floordivpow2(int_ceildiv
 312               (pi->y, comp->dy << levelno),
 313               res->pdy) - int_floordivpow2(try0, res->pdy);
 314      pi->precno = prci + prcj * res->pw;
 315      for (pi->layno = 0; pi->layno < pi->poc.layno1; pi->layno++) {
 316        if (!pi->
 317       include[pi->layno * pi->step_l +
 318          pi->resno * pi->step_r +
 319          pi->compno * pi->step_c +
 320          pi->precno * pi->step_p]) {
 321          pi->include[pi->layno * pi->step_l +
 322            pi->resno * pi->step_r +
 323            pi->compno * pi->step_c +
 324            pi->precno * pi->step_p] = 1;
 325          return 1;
 326        }
 327      skip:;
 328      }
 329    }
 330       }
 331     }
 332   }
 333   return 0;
 334 }
 335
 336 /* <summary>
 337  * Get next packet in precinct-component-resolution-layer order.
 338  *
 339  * pi: packet iterator to modify
 340  * </summary> */
 341 int pi_next_pcrl(pi_iterator_t * pi)
 342 {
 343   pi_comp_t *comp;
 344   pi_resolution_t *res;
 345   if (!pi->first) {
 346     comp = &pi->comps[pi->compno];
 347     goto skip;
 348   } else {
 349     int compno, resno;
 350     pi->first = 0;
 351     pi->dx = 0;
 352     pi->dy = 0;
 353     for (compno = 0; compno < pi->numcomps; compno++) {
 354       comp = &pi->comps[compno];
 355       for (resno = 0; resno < comp->numresolutions; resno++) {
 356    int dx, dy;
 357    res = &comp->resolutions[resno];
 358    dx = comp->dx *
 359      (1 << (res->pdx + comp->numresolutions - 1 - resno));
 360    dy = comp->dy *
 361      (1 << (res->pdy + comp->numresolutions - 1 - resno));
 362    pi->dx = !pi->dx ? dx : int_min(pi->dx, dx);
 363    pi->dy = !pi->dy ? dy : int_min(pi->dy, dy);
 364       }
 365     }
 366   }
 367   for (pi->y = pi->ty0; pi->y < pi->ty1;
 368   pi->y += pi->dy - (pi->y % pi->dy)) {
 369     for (pi->x = pi->tx0; pi->x < pi->tx1;
 370     pi->x += pi->dx - (pi->x % pi->dx)) {
 371       for (pi->compno = pi->poc.compno0;
 372       pi->compno < pi->poc.compno1; pi->compno++) {
 373    comp = &pi->comps[pi->compno];
 374    for (pi->resno = pi->poc.resno0;
 375    pi->resno < int_min(pi->poc.resno1,
 376      comp->numresolutions); pi->resno++) {
 377      int levelno;
 378      int trx0, try0;
 379      int trx1, try1;/* Add antonin pcrl*/
 380      int rpx, rpy;
 381      int prci, prcj;
 382      res = &comp->resolutions[pi->resno];
 383      levelno = comp->numresolutions - 1 - pi->resno;
 384      trx0 = int_ceildiv(pi->tx0, comp->dx << levelno);
 385      try0 = int_ceildiv(pi->ty0, comp->dy << levelno);
 386      trx1 = int_ceildiv(pi->tx1, comp->dx << levelno);/* Add antonin pcrl*/
 387      try1 = int_ceildiv(pi->ty1, comp->dy << levelno);/* Add antonin pcrl*/
 388      rpx = res->pdx + levelno;
 389      rpy = res->pdy + levelno;
 390      if (!
 391        (pi->x % (comp->dx << rpx) == 0
 392        || (pi->x == pi->tx0 && (trx0 << levelno) % (1 << rpx)))) {
 393        continue;
 394      }
 395      if (!
 396        (pi->y % (comp->dy << rpy) == 0
 397        || (pi->y == pi->ty0 && (try0 << levelno) % (1 << rpx)))) {
 398        continue;
 399      }
 400
 401      /*Add Antonin : sizebug1*/
 402      if ((res->pw==0)||(res->pw==0)) continue;
 403      /*ddA*/
 404
 405      /*Add Antonin : pcrl*/
 406      if ((trx0==trx1)||(try0==try1)) continue;
 407      /*ddA*/
 408
 409      prci =
 410        int_floordivpow2(int_ceildiv
 411        (pi->x, comp->dx << levelno),
 412        res->pdx) - int_floordivpow2(trx0, res->pdx);
 413      prcj =
 414        int_floordivpow2(int_ceildiv
 415        (pi->y, comp->dy << levelno),
 416        res->pdy) - int_floordivpow2(try0, res->pdy);
 417      pi->precno = prci + prcj * res->pw;
 418      for (pi->layno = 0; pi->layno < pi->poc.layno1; pi->layno++) {
 419        if (!pi->
 420          include[pi->layno * pi->step_l +
 421          pi->resno * pi->step_r +
 422          pi->compno * pi->step_c +
 423          pi->precno * pi->step_p]) {
 424          pi->include[pi->layno * pi->step_l +
 425       pi->resno * pi->step_r +
 426       pi->compno * pi->step_c +
 427       pi->precno * pi->step_p] = 1;
 428          return 1;
 429        }
 430 skip:;
 431      }
 432    }
 433       }
 434     }
 435   }
 436   return 0;
 437 }
 438
 439 /* <summary>
 440  * Get next packet in component-precinct-resolution-layer order.
 441  *
 442  * pi: packet iterator to modify
 443  * </summary> */
 444 int pi_next_cprl(pi_iterator_t * pi)
 445 {
 446   pi_comp_t *comp;
 447   pi_resolution_t *res;
 448   if (!pi->first) {
 449     comp = &pi->comps[pi->compno];
 450     goto skip;
 451   } else {
 452     pi->first = 0;
 453   }
 454   for (pi->compno = pi->poc.compno0; pi->compno < pi->poc.compno1;
 455        pi->compno++) {
 456     int resno;
 457     comp = &pi->comps[pi->compno];
 458     pi->dx = 0;
 459     pi->dy = 0;
 460     for (resno = 0; resno < comp->numresolutions; resno++) {
 461       int dx, dy;
 462       res = &comp->resolutions[resno];
 463       dx = comp->dx * (1 << (res->pdx + comp->numresolutions - 1 - resno));
 464       dy = comp->dy * (1 << (res->pdy + comp->numresolutions - 1 - resno));
 465       pi->dx = !pi->dx ? dx : int_min(pi->dx, dx);
 466       pi->dy = !pi->dy ? dy : int_min(pi->dy, dy);
 467     }
 468     for (pi->y = pi->ty0; pi->y < pi->ty1;
 469     pi->y += pi->dy - (pi->y % pi->dy)) {
 470       for (pi->x = pi->tx0; pi->x < pi->tx1;
 471       pi->x += pi->dx - (pi->x % pi->dx)) {
 472    for (pi->resno = pi->poc.resno0;
 473         pi->resno < int_min(pi->poc.resno1,
 474              comp->numresolutions); pi->resno++) {
 475      int levelno;
 476      int trx0, try0;
 477      int trx1, try1;/* Add antonin pcrl*/
 478      int rpx, rpy;
 479      int prci, prcj;
 480      res = &comp->resolutions[pi->resno];
 481      levelno = comp->numresolutions - 1 - pi->resno;
 482      trx0 = int_ceildiv(pi->tx0, comp->dx << levelno);
 483      try0 = int_ceildiv(pi->ty0, comp->dy << levelno);
 484      trx1 = int_ceildiv(pi->tx1, comp->dx << levelno);/* Add antonin pcrl*/
 485      try1 = int_ceildiv(pi->ty1, comp->dy << levelno);/* Add antonin pcrl*/
 486      rpx = res->pdx + levelno;
 487      rpy = res->pdy + levelno;
 488      if (!
 489          (pi->x % (comp->dx << rpx) == 0
 490           || (pi->x == pi->tx0 && (trx0 << levelno) % (1 << rpx)))) {
 491        continue;
 492      }
 493      if (!
 494          (pi->y % (comp->dy << rpy) == 0
 495           || (pi->y == pi->ty0 && (try0 << levelno) % (1 << rpx)))) {
 496        continue;
 497      }
 498
 499      /*Add Antonin : sizebug1*/
 500      if ((res->pw==0)||(res->pw==0)) continue;
 501      /*ddA*/
 502
 503      /*Add Antonin : pcrl*/
 504      if ((trx0==trx1)||(try0==try1)) continue;
 505      /*ddA*/
 506
 507      prci =
 508        int_floordivpow2(int_ceildiv
 509               (pi->x, comp->dx << levelno),
 510               res->pdx) - int_floordivpow2(trx0, res->pdx);
 511      prcj =
 512        int_floordivpow2(int_ceildiv
 513               (pi->y, comp->dy << levelno),
 514               res->pdy) - int_floordivpow2(try0, res->pdy);
 515      pi->precno = prci + prcj * res->pw;
 516      for (pi->layno = 0; pi->layno < pi->poc.layno1; pi->layno++) {
 517        if (!pi->
 518       include[pi->layno * pi->step_l +
 519          pi->resno * pi->step_r +
 520          pi->compno * pi->step_c +
 521          pi->precno * pi->step_p]) {
 522          pi->include[pi->layno * pi->step_l +
 523            pi->resno * pi->step_r +
 524            pi->compno * pi->step_c +
 525            pi->precno * pi->step_p] = 1;
 526          return 1;
 527        }
 528      skip:;
 529      }
 530    }
 531       }
 532     }
 533   }
 534   return 0;
 535 }
 536
 537 /* <summary>
 538  * Get next packet.
 539   *
 540   * pi: packet iterator to modify
 541   * </summary> */
 542 int pi_next(pi_iterator_t * pi)
 543 {
 544   switch (pi->poc.prg) {
 545   case 0:
 546     return pi_next_lrcp(pi);
 547   case 1:
 548     return pi_next_rlcp(pi);
 549   case 2:
 550     return pi_next_rpcl(pi);
 551   case 3:
 552     return pi_next_pcrl(pi);
 553   case 4:
 554     return pi_next_cprl(pi);
 555   }
 556   return 0;
 557 }