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) * sizeof(pi_iterator_t));
  48
  49         for (pino = 0; pino < tcp->numpocs + 1; pino++) {       /* change */
  50                 p = tileno % cp->tw;
  51                 q = tileno / cp->tw;
  52
  53                 pi[pino].tx0 = int_max(cp->tx0 + p * cp->tdx, img->x0);
  54                 pi[pino].ty0 = int_max(cp->ty0 + q * cp->tdy, img->y0);
  55                 pi[pino].tx1 = int_min(cp->tx0 + (p + 1) * cp->tdx, img->x1);
  56                 pi[pino].ty1 = int_min(cp->ty0 + (q + 1) * cp->tdy, img->y1);
  57                 pi[pino].numcomps = img->numcomps;
  58                 pi[pino].comps = (pi_comp_t *) malloc(img->numcomps * sizeof(pi_comp_t));
  59
  60                 for (compno = 0; compno < pi->numcomps; compno++) {
  61                         int tcx0, tcy0, tcx1, tcy1;
  62                         pi_comp_t *comp = &pi[pino].comps[compno];
  63                         tccp = &tcp->tccps[compno];
  64                         comp->dx = img->comps[compno].dx;
  65                         comp->dy = img->comps[compno].dy;
  66                         comp->numresolutions = tccp->numresolutions;
  67                         comp->resolutions = (pi_resolution_t *) malloc(comp->numresolutions * sizeof(pi_resolution_t));
  68                         tcx0 = int_ceildiv(pi->tx0, comp->dx);
  69                         tcy0 = int_ceildiv(pi->ty0, comp->dy);
  70                         tcx1 = int_ceildiv(pi->tx1, comp->dx);
  71                         tcy1 = int_ceildiv(pi->ty1, comp->dy);
  72                         if (comp->numresolutions > maxres) {
  73                                 maxres = comp->numresolutions;
  74                         }
  75                         for (resno = 0; resno < comp->numresolutions; resno++) {
  76                                 int levelno;
  77                                 int rx0, ry0, rx1, ry1;
  78                                 int px0, py0, px1, py1;
  79                                 pi_resolution_t *res = &comp->resolutions[resno];
  80                                 if (tccp->csty & J2K_CCP_CSTY_PRT) {
  81                                         res->pdx = tccp->prcw[resno];
  82                                         res->pdy = tccp->prch[resno];
  83                                 } else {
  84                                         res->pdx = 15;
  85                                         res->pdy = 15;
  86                                 }
  87                                 levelno = comp->numresolutions - 1 - resno;
  88                                 rx0 = int_ceildivpow2(tcx0, levelno);
  89                                 ry0 = int_ceildivpow2(tcy0, levelno);
  90                                 rx1 = int_ceildivpow2(tcx1, levelno);
  91                                 ry1 = int_ceildivpow2(tcy1, levelno);
  92                                 px0 = int_floordivpow2(rx0, res->pdx) << res->pdx;
  93                                 py0 = int_floordivpow2(ry0, res->pdy) << res->pdy;
  94                                 px1 = int_ceildivpow2(rx1, res->pdx) << res->pdx;
  95                                 py1 = int_ceildivpow2(ry1, res->pdy) << res->pdy;
  96                                 res->pw = (px1 - px0) >> res->pdx;
  97                                 res->ph = (py1 - py0) >> res->pdy;
  98                         }
  99                 }
 100
 101                 tccp = &tcp->tccps[0];
 102                 pi[pino].step_p=1;
 103                 pi[pino].step_c=100*pi[pino].step_p;
 104                 pi[pino].step_r=img->numcomps*pi[pino].step_c;
 105                 pi[pino].step_l=maxres*pi[pino].step_r;
 106
 107                 if (pino==0)
 108                   {
 109                     pi[pino].include = (short int*)malloc(img->numcomps*maxres*tcp->numlayers*100*sizeof(short int));
 110                     for (i=0 ; i<img->numcomps*maxres*tcp->numlayers*100; i++)
 111                       pi[pino].include[i]=0;
 112                   }
 113                 /* pi[pino].include=(short int*)calloc(img->numcomps*maxres*tcp->numlayers*1000,sizeof(short int));*/
 114                 else
 115                   pi[pino].include=pi[pino-1].include;
 116
 117                 if (tcp->POC == 0) {
 118                         pi[pino].first = 1;
 119                         pi[pino].poc.resno0 = 0;
 120                         pi[pino].poc.compno0 = 0;
 121                         pi[pino].poc.layno1 = tcp->numlayers;
 122                         pi[pino].poc.resno1 = maxres;
 123                         pi[pino].poc.compno1 = img->numcomps;
 124                         pi[pino].poc.prg = tcp->prg;
 125                 } else {
 126                         pi[pino].first = 1;
 127                         pi[pino].poc.resno0 = tcp->pocs[pino].resno0;
 128                         pi[pino].poc.compno0 = tcp->pocs[pino].compno0;
 129                         pi[pino].poc.layno1 = tcp->pocs[pino].layno1;
 130                         pi[pino].poc.resno1 = tcp->pocs[pino].resno1;
 131                         pi[pino].poc.compno1 = tcp->pocs[pino].compno1;
 132                         pi[pino].poc.prg = tcp->pocs[pino].prg;
 133                 }
 134         }
 135         return pi;
 136 }
 137
 138 /* <summary>
 139  * Get next packet in layer-resolution-component-precinct order.
 140  *
 141  * pi: packet iterator to modify
 142  * </summary> */
 143 int pi_next_lrcp(pi_iterator_t * pi)
 144 {
 145         pi_comp_t *comp;
 146         pi_resolution_t *res;
 147
 148         if (!pi->first) {
 149                 comp = &pi->comps[pi->compno];
 150                 res = &comp->resolutions[pi->resno];
 151                 goto skip;
 152         } else {
 153                 pi->first = 0;
 154         }
 155         for (pi->layno = 0; pi->layno < pi->poc.layno1; pi->layno++) {
 156           for (pi->resno = pi->poc.resno0; pi->resno < pi->poc.resno1; pi->resno++) {
 157             for (pi->compno = pi->poc.compno0; pi->compno < pi->poc.compno1; pi->compno++) {
 158               comp = &pi->comps[pi->compno];
 159               if (pi->resno >= comp->numresolutions) {
 160                 continue;
 161               }
 162               res = &comp->resolutions[pi->resno];
 163               for (pi->precno = 0; pi->precno < res->pw * res->ph; pi->precno++)
 164                 {
 165                   if (!pi->include[pi->layno*pi->step_l+pi->resno*pi->step_r+pi->compno*pi->step_c+pi->precno*pi->step_p]){
 166                     pi->include[pi->layno*pi->step_l+pi->resno*pi->step_r+pi->compno*pi->step_c+pi->precno*pi->step_p] = 1;
 167                     return 1;
 168                   }
 169                 skip:;
 170                 }
 171             }
 172           }
 173         }
 174         return 0;
 175 }
 176
 177 /* <summary>
 178  * Get next packet in resolution-layer-component-precinct order.
 179  *
 180  * pi: packet iterator to modify
 181  * </summary> */
 182 int pi_next_rlcp(pi_iterator_t * pi)
 183 {
 184         pi_comp_t *comp;
 185         pi_resolution_t *res;
 186         if (!pi->first) {
 187                 comp = &pi->comps[pi->compno];
 188                 res = &comp->resolutions[pi->resno];
 189                 goto skip;
 190         } else {
 191                 pi->first = 0;
 192         }
 193         for (pi->resno = pi->poc.resno0; pi->resno < pi->poc.resno1; pi->resno++) {
 194                 for (pi->layno = 0; pi->layno < pi->poc.layno1; pi->layno++) {
 195                         for (pi->compno = pi->poc.compno0; pi->compno < pi->poc.compno1;
 196                              pi->compno++) {
 197                           comp = &pi->comps[pi->compno];
 198                           if (pi->resno >= comp->numresolutions) {
 199                             continue;
 200                           }
 201                           res = &comp->resolutions[pi->resno];
 202                           for (pi->precno = 0; pi->precno < res->pw * res->ph; pi->precno++) {
 203                             if (!pi->include[pi->layno*pi->step_l+pi->resno*pi->step_r+pi->compno*pi->step_c+pi->precno*pi->step_p]){
 204                               pi->include[pi->layno*pi->step_l+pi->resno*pi->step_r+pi->compno*pi->step_c+pi->precno*pi->step_p] = 1;
 205                               return 1;
 206                             }
 207                           skip:;
 208                           }
 209                         }
 210                 }
 211         }
 212         return 0;
 213 }
 214
 215 /* <summary>
 216  * Get next packet in resolution-precinct-component-layer order.
 217  *
 218  * pi: packet iterator to modify
 219  * </summary> */
 220 int pi_next_rpcl(pi_iterator_t * pi)
 221 {
 222         pi_comp_t *comp;
 223         pi_resolution_t *res;
 224         if (!pi->first) {
 225                 goto skip;
 226         } else {
 227                 int compno, resno;
 228                 pi->first = 0;
 229                 pi->dx = 0;
 230                 pi->dy = 0;
 231                 for (compno = 0; compno < pi->numcomps; compno++) {
 232                         comp = &pi->comps[compno];
 233                         for (resno = 0; resno < comp->numresolutions; resno++) {
 234                                 int dx, dy;
 235                                 res = &comp->resolutions[resno];
 236                                 dx = comp->dx * (1 << (res->pdx + comp->numresolutions - 1 - resno));
 237                                 dy = comp->dy * (1 << (res->pdy + comp->numresolutions - 1 - resno));
 238                                 pi->dx = !pi->dx ? dx : int_min(pi->dx, dx);
 239                                 pi->dy = !pi->dy ? dy : int_min(pi->dy, dy);
 240                         }
 241                 }
 242         }
 243         for (pi->resno = pi->poc.resno0; pi->resno < pi->poc.resno1; pi->resno++) {
 244                 for (pi->y = pi->ty0; pi->y < pi->ty1; pi->y += pi->dy - (pi->y % pi->dy)) {
 245                         for (pi->x = pi->tx0; pi->x < pi->tx1; pi->x += pi->dx - (pi->x % pi->dx)) {
 246                                 for (pi->compno = pi->poc.compno0; pi->compno < pi->poc.compno1; pi->compno++) {
 247                                         int levelno;
 248                                         int trx0, try0;
 249                                         int rpx, rpy;
 250                                         int prci, prcj;
 251                                         comp = &pi->comps[pi->compno];
 252                                         if (pi->resno >= comp->numresolutions) {
 253                                                 continue;
 254                                         }
 255                                         res = &comp->resolutions[pi->resno];
 256                                         levelno = comp->numresolutions - 1 - pi->resno;
 257                                         trx0 = int_ceildiv(pi->tx0, comp->dx << levelno);
 258                                         try0 = int_ceildiv(pi->ty0, comp->dy << levelno);
 259                                         rpx = res->pdx + levelno;
 260                                         rpy = res->pdy + levelno;
 261                                         if (!(pi->x % (comp->dx << rpx) == 0 || (pi->x == pi->tx0 && (trx0 << levelno) % (1 << rpx)))) {
 262                                                 continue;
 263                                         }
 264                                         if (!(pi->y % (comp->dy << rpy) == 0 || (pi->y == pi->ty0 && (try0 << levelno) % (1 << rpx)))) {
 265                                                 continue;
 266                                         }
 267                                         prci = int_floordivpow2(int_ceildiv(pi->x, comp->dx << levelno),
 268                                                                 res->pdx) - int_floordivpow2(trx0, res->pdx);
 269                                         prcj = int_floordivpow2(int_ceildiv(pi->y, comp->dy << levelno),
 270                                                                 res->pdy) - int_floordivpow2(try0, res->pdy);
 271                                         pi->precno = prci + prcj * res->pw;
 272                                         for (pi->layno = 0; pi->layno < pi->poc.layno1; pi->layno++) {
 273                                           if (!pi->include[pi->layno*pi->step_l+pi->resno*pi->step_r+pi->compno*pi->step_c+pi->precno*pi->step_p]){
 274                                             pi->include[pi->layno*pi->step_l+pi->resno*pi->step_r+pi->compno*pi->step_c+pi->precno*pi->step_p] = 1;
 275                                             return 1;
 276                                                 }
 277                                         skip:;
 278                                         }
 279                                 }
 280                         }
 281                 }
 282         }
 283         return 0;
 284 }
 285
 286 /* <summary>
 287  * Get next packet in precinct-component-resolution-layer order.
 288  *
 289  * pi: packet iterator to modify
 290  * </summary> */
 291 int pi_next_pcrl(pi_iterator_t * pi)
 292 {
 293         pi_comp_t *comp;
 294         pi_resolution_t *res;
 295         if (!pi->first) {
 296                 comp = &pi->comps[pi->compno];
 297                 goto skip;
 298         } else {
 299                 int compno, resno;
 300                 pi->first = 0;
 301                 pi->dx = 0;
 302                 pi->dy = 0;
 303                 for (compno = 0; compno < pi->numcomps; compno++) {
 304                         comp = &pi->comps[compno];
 305                         for (resno = 0; resno < comp->numresolutions; resno++) {
 306                                 int dx, dy;
 307                                 res = &comp->resolutions[resno];
 308                                 dx = comp->dx * (1 << (res->pdx + comp->numresolutions - 1 - resno));
 309                                 dy = comp->dy * (1 << (res->pdy + comp->numresolutions - 1 - resno));
 310                                 pi->dx = !pi->dx ? dx : int_min(pi->dx, dx);
 311                                 pi->dy = !pi->dy ? dy : int_min(pi->dy, dy);
 312                         }
 313                 }
 314         }
 315         for (pi->y = pi->ty0; pi->y < pi->ty1; pi->y += pi->dy - (pi->y % pi->dy)) {
 316                 for (pi->x = pi->tx0; pi->x < pi->tx1; pi->x += pi->dx - (pi->x % pi->dx)) {
 317                         for (pi->compno = pi->poc.compno0; pi->compno < pi->poc.compno1; pi->compno++) {
 318                                 comp = &pi->comps[pi->compno];
 319                                 for (pi->resno = pi->poc.resno0; pi->resno < int_min(pi->poc.resno1, comp->numresolutions); pi->resno++) {
 320                                         int levelno;
 321                                         int trx0, try0;
 322                                         int rpx, rpy;
 323                                         int prci, prcj;
 324                                         res = &comp->resolutions[pi->resno];
 325                                         levelno = comp->numresolutions - 1 - pi->resno;
 326                                         trx0 = int_ceildiv(pi->tx0, comp->dx << levelno);
 327                                         try0 = int_ceildiv(pi->ty0, comp->dy << levelno);
 328                                         rpx = res->pdx + levelno;
 329                                         rpy = res->pdy + levelno;
 330                                         if (!(pi->x % (comp->dx << rpx) == 0 || (pi->x == pi->tx0 && (trx0 << levelno) % (1 << rpx)))) {
 331                                                 continue;
 332                                         }
 333                                         if (!(pi->y % (comp->dy << rpy) == 0 || (pi->y == pi->ty0 && (try0 << levelno) % (1 << rpx)))) {
 334                                                 continue;
 335                                         }
 336                                         prci = int_floordivpow2(int_ceildiv(pi->x, comp->dx << levelno), res->pdx) - int_floordivpow2(trx0, res->pdx);
 337                                         prcj = int_floordivpow2(int_ceildiv(pi->y, comp->dy << levelno), res->pdy) - int_floordivpow2(try0, res->pdy);
 338                                         pi->precno = prci + prcj * res->pw;
 339                                         for (pi->layno = 0; pi->layno < pi->poc.layno1; pi->layno++) {
 340                                           if (! pi->include[pi->layno*pi->step_l+pi->resno*pi->step_r+pi->compno*pi->step_c+pi->precno*pi->step_p]){
 341                                             pi->include[pi->layno*pi->step_l+pi->resno*pi->step_r+pi->compno*pi->step_c+pi->precno*pi->step_p] = 1;
 342                                                         return 1;
 343                                                 }
 344                                         skip:;
 345                                         }
 346                                 }
 347                         }
 348                 }
 349         }
 350         return 0;
 351 }
 352
 353 /* <summary>
 354  * Get next packet in component-precinct-resolution-layer order.
 355  *
 356  * pi: packet iterator to modify
 357  * </summary> */
 358 int pi_next_cprl(pi_iterator_t * pi)
 359 {
 360         pi_comp_t *comp;
 361         pi_resolution_t *res;
 362         if (!pi->first) {
 363                 comp = &pi->comps[pi->compno];
 364                 goto skip;
 365         } else {
 366                 pi->first = 0;
 367         }
 368         for (pi->compno = pi->poc.compno0; pi->compno < pi->poc.compno1;
 369                          pi->compno++) {
 370                 int resno;
 371                 comp = &pi->comps[pi->compno];
 372                 pi->dx = 0;
 373                 pi->dy = 0;
 374                 for (resno = 0; resno < comp->numresolutions; resno++) {
 375                         int dx, dy;
 376                         res = &comp->resolutions[resno];
 377                         dx = comp->dx * (1 << (res->pdx + comp->numresolutions - 1 - resno));
 378                         dy = comp->dy * (1 << (res->pdy + comp->numresolutions - 1 - resno));
 379                         pi->dx = !pi->dx ? dx : int_min(pi->dx, dx);
 380                         pi->dy = !pi->dy ? dy : int_min(pi->dy, dy);
 381                 }
 382                 for (pi->y = pi->ty0; pi->y < pi->ty1;
 383                                  pi->y += pi->dy - (pi->y % pi->dy)) {
 384                         for (pi->x = pi->tx0; pi->x < pi->tx1;
 385                                          pi->x += pi->dx - (pi->x % pi->dx)) {
 386                                 for (pi->resno = pi->poc.resno0;
 387                                                  pi->resno < int_min(pi->poc.resno1, comp->numresolutions);
 388                                                  pi->resno++) {
 389                                         int levelno;
 390                                         int trx0, try0;
 391                                         int rpx, rpy;
 392                                         int prci, prcj;
 393                                         res = &comp->resolutions[pi->resno];
 394                                         levelno = comp->numresolutions - 1 - pi->resno;
 395                                         trx0 = int_ceildiv(pi->tx0, comp->dx << levelno);
 396                                         try0 = int_ceildiv(pi->ty0, comp->dy << levelno);
 397                                         rpx = res->pdx + levelno;
 398                                         rpy = res->pdy + levelno;
 399                                         if (!(pi->x % (comp->dx << rpx) == 0 || (pi->x == pi->tx0 && (trx0 << levelno) % (1 << rpx)))) {
 400                                                 continue;
 401                                         }
 402                                         if (!(pi->y % (comp->dy << rpy) == 0 || (pi->y == pi->ty0 && (try0 << levelno) % (1 << rpx)))) {
 403                                                 continue;
 404                                         }
 405                                         prci = int_floordivpow2(int_ceildiv(pi->x, comp->dx << levelno), res->pdx) - int_floordivpow2(trx0, res->pdx);
 406                                         prcj = int_floordivpow2(int_ceildiv(pi->y, comp->dy << levelno), res->pdy) - int_floordivpow2(try0, res->pdy);
 407                                         pi->precno = prci + prcj * res->pw;
 408                                         for (pi->layno = 0; pi->layno < pi->poc.layno1; pi->layno++) {
 409                                           if (! pi->include[pi->layno*pi->step_l+pi->resno*pi->step_r+pi->compno*pi->step_c+pi->precno*pi->step_p]){
 410                                             pi->include[pi->layno*pi->step_l+pi->resno*pi->step_r+pi->compno*pi->step_c+pi->precno*pi->step_p] = 1;
 411                                             return 1;
 412                                                 }
 413                                         skip:;
 414                                         }
 415                                 }
 416                         }
 417                 }
 418         }
 419         return 0;
 420 }
 421
 422 /* <summary>
 423  * Get next packet.
 424   *
 425   * pi: packet iterator to modify
 426   * </summary> */
 427 int pi_next(pi_iterator_t * pi)
 428 {
 429         switch (pi->poc.prg) {
 430         case 0:
 431                 return pi_next_lrcp(pi);
 432         case 1:
 433                 return pi_next_rlcp(pi);
 434         case 2:
 435                 return pi_next_rpcl(pi);
 436         case 3:
 437                 return pi_next_pcrl(pi);
 438         case 4:
 439                 return pi_next_cprl(pi);
 440         }
 441         return 0;
 442 }