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