libopenjpeg/tcd.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 "tcd.h"
  30 #include "int.h"
  31 #include "t1.h"
  32 #include "t2.h"
  33 #include "dwt.h"
  34 #include "mct.h"
  35 #include <setjmp.h>
  36 #include <float.h>
  37 #include <stdio.h>
  38 #include <time.h>
  39 #include <math.h>
  40 #include <stdlib.h>
  41 #include <string.h>
  42
  43 static tcd_image_t tcd_image;
  44
  45 static j2k_image_t *tcd_img;
  46 static j2k_cp_t *tcd_cp;
  47
  48 static tcd_tile_t *tcd_tile;
  49 static j2k_tcp_t *tcd_tcp;
  50 static int tcd_tileno;
  51
  52 static tcd_tile_t *tile;
  53 static tcd_tilecomp_t *tilec;
  54 static tcd_resolution_t *res;
  55 static tcd_band_t *band;
  56 static tcd_precinct_t *prc;
  57 static tcd_cblk_t *cblk;
  58
  59 extern jmp_buf j2k_error;
  60
  61 void
  62 tcd_dump (tcd_image_t * img, int curtileno)
  63 {
  64   int tileno, compno, resno, bandno, precno, cblkno;
  65   fprintf (stderr, "image {\n");
  66   fprintf (stderr, "  tw=%d, th=%d x0=%d x1=%d y0=%d y1=%d\n", img->tw,
  67            img->th, tcd_img->x0, tcd_img->x1, tcd_img->y0, tcd_img->y1);
  68   for (tileno = 0; tileno < img->th * img->tw; tileno++)
  69     {
  70       tcd_tile_t *tile = &tcd_image.tiles[tileno];
  71       fprintf (stderr, "  tile {\n");
  72       fprintf (stderr, "    x0=%d, y0=%d, x1=%d, y1=%d, numcomps=%d\n",
  73                tile->x0, tile->y0, tile->x1, tile->y1, tile->numcomps);
  74       for (compno = 0; compno < tile->numcomps; compno++)
  75         {
  76           tcd_tilecomp_t *tilec = &tile->comps[compno];
  77           fprintf (stderr, "    tilec {\n");
  78           fprintf (stderr,
  79                    "      x0=%d, y0=%d, x1=%d, y1=%d, numresolutions=%d\n",
  80                    tilec->x0, tilec->y0, tilec->x1, tilec->y1,
  81                    tilec->numresolutions);
  82           for (resno = 0; resno < tilec->numresolutions; resno++)
  83             {
  84               tcd_resolution_t *res = &tilec->resolutions[resno];
  85               fprintf (stderr, "\n   res {\n");
  86               fprintf (stderr,
  87                        "          x0=%d, y0=%d, x1=%d, y1=%d, pw=%d, ph=%d, numbands=%d\n",
  88                        res->x0, res->y0, res->x1, res->y1, res->pw, res->ph,
  89                        res->numbands);
  90               for (bandno = 0; bandno < res->numbands; bandno++)
  91                 {
  92                   tcd_band_t *band = &res->bands[bandno];
  93                   fprintf (stderr, "        band {\n");
  94                   fprintf (stderr,
  95                            "          x0=%d, y0=%d, x1=%d, y1=%d, stepsize=%d, numbps=%d\n",
  96                            band->x0, band->y0, band->x1, band->y1,
  97                            band->stepsize, band->numbps);
  98                   for (precno = 0; precno < res->pw * res->ph; precno++)
  99                     {
 100                       tcd_precinct_t *prec = &band->precincts[precno];
 101                       fprintf (stderr, "          prec {\n");
 102                       fprintf (stderr,
 103                                "            x0=%d, y0=%d, x1=%d, y1=%d, cw=%d, ch=%d\n",
 104                                prec->x0, prec->y0, prec->x1, prec->y1,
 105                                prec->cw, prec->ch);
 106                       for (cblkno = 0; cblkno < prec->cw * prec->ch; cblkno++)
 107                         {
 108                           tcd_cblk_t *cblk = &prec->cblks[cblkno];
 109                           fprintf (stderr, "            cblk {\n");
 110                           fprintf (stderr,
 111                                    "              x0=%d, y0=%d, x1=%d, y1=%d\n",
 112                                    cblk->x0, cblk->y0, cblk->x1, cblk->y1);
 113                           fprintf (stderr, "            }\n");
 114                         }
 115                       fprintf (stderr, "          }\n");
 116                     }
 117                   fprintf (stderr, "        }\n");
 118                 }
 119               fprintf (stderr, "      }\n");
 120             }
 121           fprintf (stderr, "    }\n");
 122         }
 123       fprintf (stderr, "  }\n");
 124     }
 125   fprintf (stderr, "}\n");
 126 }
 127
 128 void
 129 tcd_malloc_encode (j2k_image_t * img, j2k_cp_t * cp, int curtileno)
 130 {
 131   int tileno, compno, resno, bandno, precno, cblkno;
 132   tcd_img = img;
 133   tcd_cp = cp;
 134   tcd_image.tw = cp->tw;
 135   tcd_image.th = cp->th;
 136   tcd_image.tiles = (tcd_tile_t *) malloc (sizeof (tcd_tile_t));
 137
 138   for (tileno = 0; tileno < 1; tileno++)
 139     {
 140       j2k_tcp_t *tcp = &cp->tcps[curtileno];
 141       int j;
 142       /* cfr p59 ISO/IEC FDIS15444-1 : 2000 (18 august 2000) */
 143       int p = curtileno % cp->tw;       /* si numerotation matricielle .. */
 144       int q = curtileno / cp->tw;       /* .. coordonnees de la tile (q,p) q pour ligne et p pour colonne */
 145       /* tcd_tile_t *tile=&tcd_image.tiles[tileno]; */
 146       tile = tcd_image.tiles;
 147       /* 4 borders of the tile rescale on the image if necessary */
 148       tile->x0 = int_max (cp->tx0 + p * cp->tdx, img->x0);
 149       tile->y0 = int_max (cp->ty0 + q * cp->tdy, img->y0);
 150       tile->x1 = int_min (cp->tx0 + (p + 1) * cp->tdx, img->x1);
 151       tile->y1 = int_min (cp->ty0 + (q + 1) * cp->tdy, img->y1);
 152       tile->numcomps = img->numcomps;
 153       /* tile->PPT=img->PPT;  */
 154       /* Modification of the RATE >> */
 155       for (j = 0; j < tcp->numlayers; j++)
 156         {
 157           tcp->rates[j] = tcp->rates[j] ? int_ceildiv (tile->numcomps * (tile->x1 - tile->x0) * (tile->y1 - tile->y0) * img->comps[0].prec, (tcp->rates[j] * 8 * img->comps[0].dx * img->comps[0].dy)) : 0;     //Mod antonin losslessbug
 158           if (tcp->rates[j])
 159             {
 160               if (j && tcp->rates[j] < tcp->rates[j - 1] + 10)
 161                 {
 162                   tcp->rates[j] = tcp->rates[j - 1] + 20;
 163                 }
 164               else
 165                 {
 166                   if (!j && tcp->rates[j] < 30)
 167                     tcp->rates[j] = 30;
 168                 }
 169             }
 170         }
 171       /* << Modification of the RATE */
 172
 173       tile->comps =
 174         (tcd_tilecomp_t *) malloc (img->numcomps * sizeof (tcd_tilecomp_t));
 175       for (compno = 0; compno < tile->numcomps; compno++)
 176         {
 177           j2k_tccp_t *tccp = &tcp->tccps[compno];
 178           /* tcd_tilecomp_t *tilec=&tile->comps[compno]; */
 179           tilec = &tile->comps[compno];
 180           /* border of each tile component (global) */
 181           tilec->x0 = int_ceildiv (tile->x0, img->comps[compno].dx);
 182
 183           tilec->y0 = int_ceildiv (tile->y0, img->comps[compno].dy);
 184           tilec->x1 = int_ceildiv (tile->x1, img->comps[compno].dx);
 185           tilec->y1 = int_ceildiv (tile->y1, img->comps[compno].dy);
 186
 187           tilec->data =
 188             (int *) malloc ((tilec->x1 - tilec->x0) *
 189                             (tilec->y1 - tilec->y0) * sizeof (int));
 190           tilec->numresolutions = tccp->numresolutions;
 191
 192           tilec->resolutions =
 193             (tcd_resolution_t *) malloc (tilec->numresolutions *
 194                                          sizeof (tcd_resolution_t));
 195
 196           for (resno = 0; resno < tilec->numresolutions; resno++)
 197             {
 198               int pdx, pdy;
 199               int levelno = tilec->numresolutions - 1 - resno;
 200               int tlprcxstart, tlprcystart, brprcxend, brprcyend;
 201               int tlcbgxstart, tlcbgystart, brcbgxend, brcbgyend;
 202               int cbgwidthexpn, cbgheightexpn;
 203               int cblkwidthexpn, cblkheightexpn;
 204               /* tcd_resolution_t *res=&tilec->resolutions[resno]; */
 205
 206               res = &tilec->resolutions[resno];
 207
 208               /* border for each resolution level (global) */
 209               res->x0 = int_ceildivpow2 (tilec->x0, levelno);
 210               res->y0 = int_ceildivpow2 (tilec->y0, levelno);
 211               res->x1 = int_ceildivpow2 (tilec->x1, levelno);
 212               res->y1 = int_ceildivpow2 (tilec->y1, levelno);
 213
 214               res->numbands = resno == 0 ? 1 : 3;
 215               /* p. 35, table A-23, ISO/IEC FDIS154444-1 : 2000 (18 august 2000) */
 216               if (tccp->csty & J2K_CCP_CSTY_PRT)
 217                 {
 218                   pdx = tccp->prcw[resno];
 219                   pdy = tccp->prch[resno];
 220                 }
 221               else
 222                 {
 223                   pdx = 15;
 224                   pdy = 15;
 225                 }
 226               /* p. 64, B.6, ISO/IEC FDIS15444-1 : 2000 (18 august 2000)  */
 227               tlprcxstart = int_floordivpow2 (res->x0, pdx) << pdx;
 228               tlprcystart = int_floordivpow2 (res->y0, pdy) << pdy;
 229               brprcxend = int_ceildivpow2 (res->x1, pdx) << pdx;
 230               brprcyend = int_ceildivpow2 (res->y1, pdy) << pdy;
 231
 232               res->pw = (brprcxend - tlprcxstart) >> pdx;
 233               res->ph = (brprcyend - tlprcystart) >> pdy;
 234
 235               if (resno == 0)
 236                 {
 237                   tlcbgxstart = tlprcxstart;
 238                   tlcbgystart = tlprcystart;
 239                   brcbgxend = brprcxend;
 240                   brcbgyend = brprcyend;
 241                   cbgwidthexpn = pdx;
 242                   cbgheightexpn = pdy;
 243                 }
 244               else
 245                 {
 246                   tlcbgxstart = int_ceildivpow2 (tlprcxstart, 1);
 247                   tlcbgystart = int_ceildivpow2 (tlprcystart, 1);
 248                   brcbgxend = int_ceildivpow2 (brprcxend, 1);
 249                   brcbgyend = int_ceildivpow2 (brprcyend, 1);
 250                   cbgwidthexpn = pdx - 1;
 251                   cbgheightexpn = pdy - 1;
 252                 }
 253
 254               cblkwidthexpn = int_min (tccp->cblkw, cbgwidthexpn);
 255               cblkheightexpn = int_min (tccp->cblkh, cbgheightexpn);
 256
 257               for (bandno = 0; bandno < res->numbands; bandno++)
 258                 {
 259                   int x0b, y0b, i;
 260                   int gain, numbps;
 261                   j2k_stepsize_t *ss;
 262                   band = &res->bands[bandno];
 263                   band->bandno = resno == 0 ? 0 : bandno + 1;
 264                   x0b = (band->bandno == 1) || (band->bandno == 3) ? 1 : 0;
 265                   y0b = (band->bandno == 2) || (band->bandno == 3) ? 1 : 0;
 266
 267                   if (band->bandno == 0)
 268                     {
 269                       /* band border (global) */
 270                       band->x0 = int_ceildivpow2 (tilec->x0, levelno);
 271                       band->y0 = int_ceildivpow2 (tilec->y0, levelno);
 272                       band->x1 = int_ceildivpow2 (tilec->x1, levelno);
 273                       band->y1 = int_ceildivpow2 (tilec->y1, levelno);
 274                     }
 275                   else
 276                     {
 277                       /* band border (global) */
 278                       band->x0 =
 279                         int_ceildivpow2 (tilec->x0 -
 280                                          (1 << levelno) * x0b, levelno + 1);
 281                       band->y0 =
 282                         int_ceildivpow2 (tilec->y0 -
 283                                          (1 << levelno) * y0b, levelno + 1);
 284                       band->x1 =
 285                         int_ceildivpow2 (tilec->x1 -
 286                                          (1 << levelno) * x0b, levelno + 1);
 287                       band->y1 =
 288                         int_ceildivpow2 (tilec->y1 -
 289                                          (1 << levelno) * y0b, levelno + 1);
 290
 291                     }
 292
 293                   ss = &tccp->stepsizes[resno ==
 294                                         0 ? 0 : 3 * (resno - 1) + bandno + 1];
 295                   gain =
 296                     tccp->qmfbid ==
 297                     0 ? dwt_getgain_real (band->bandno) : dwt_getgain (band->
 298                                                                        bandno);
 299                   numbps = img->comps[compno].prec + gain;
 300                   band->stepsize =
 301                     (int) floor ((1.0 + ss->mant / 2048.0) *
 302                                  pow (2.0, numbps - ss->expn) * 8192.0);
 303                   band->numbps = ss->expn + tccp->numgbits - 1; /* WHY -1 ? */
 304
 305                   band->precincts =
 306                     (tcd_precinct_t *) malloc (3 * res->pw * res->ph *
 307                                                sizeof (tcd_precinct_t));
 308
 309                   for (i = 0; i < res->pw * res->ph * 3; i++)
 310                     {
 311                       band->precincts[i].imsbtree = NULL;
 312                       band->precincts[i].incltree = NULL;
 313                     }
 314
 315                   for (precno = 0; precno < res->pw * res->ph; precno++)
 316                     {
 317                       int tlcblkxstart, tlcblkystart, brcblkxend, brcblkyend;
 318                       int cbgxstart =
 319                         tlcbgxstart +
 320                         (precno % res->pw) * (1 << cbgwidthexpn);
 321                       int cbgystart =
 322                         tlcbgystart +
 323                         (precno / res->pw) * (1 << cbgheightexpn);
 324                       int cbgxend = cbgxstart + (1 << cbgwidthexpn);
 325                       int cbgyend = cbgystart + (1 << cbgheightexpn);
 326                       /* tcd_precinct_t *prc=&band->precincts[precno]; */
 327                       prc = &band->precincts[precno];
 328                       /* precinct size (global) */
 329                       prc->x0 = int_max (cbgxstart, band->x0);
 330                       prc->y0 = int_max (cbgystart, band->y0);
 331                       prc->x1 = int_min (cbgxend, band->x1);
 332                       prc->y1 = int_min (cbgyend, band->y1);
 333
 334                       tlcblkxstart =
 335                         int_floordivpow2 (prc->x0,
 336                                           cblkwidthexpn) << cblkwidthexpn;
 337                       tlcblkystart =
 338                         int_floordivpow2 (prc->y0,
 339                                           cblkheightexpn) << cblkheightexpn;
 340                       brcblkxend =
 341                         int_ceildivpow2 (prc->x1,
 342                                          cblkwidthexpn) << cblkwidthexpn;
 343                       brcblkyend =
 344                         int_ceildivpow2 (prc->y1,
 345                                          cblkheightexpn) << cblkheightexpn;
 346                       prc->cw = (brcblkxend - tlcblkxstart) >> cblkwidthexpn;
 347                       prc->ch = (brcblkyend - tlcblkystart) >> cblkheightexpn;
 348
 349                       prc->cblks =
 350                         (tcd_cblk_t *) malloc ((prc->cw * prc->ch) *
 351                                                sizeof (tcd_cblk_t));
 352                       prc->incltree = tgt_create (prc->cw, prc->ch);
 353                       prc->imsbtree = tgt_create (prc->cw, prc->ch);
 354
 355                       for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++)
 356                         {
 357                           int cblkxstart =
 358                             tlcblkxstart +
 359                             (cblkno % prc->cw) * (1 << cblkwidthexpn);
 360                           int cblkystart =
 361                             tlcblkystart +
 362                             (cblkno / prc->cw) * (1 << cblkheightexpn);
 363                           int cblkxend = cblkxstart + (1 << cblkwidthexpn);
 364                           int cblkyend = cblkystart + (1 << cblkheightexpn);
 365
 366                           cblk = &prc->cblks[cblkno];
 367                           /* code-block size (global) */
 368                           cblk->x0 = int_max (cblkxstart, prc->x0);
 369                           cblk->y0 = int_max (cblkystart, prc->y0);
 370                           cblk->x1 = int_min (cblkxend, prc->x1);
 371                           cblk->y1 = int_min (cblkyend, prc->y1);
 372                         }
 373                     }
 374                 }
 375             }
 376         }
 377     }
 378   /* tcd_dump(&tcd_image,curtileno); */
 379 }
 380
 381 void
 382 tcd_free_encode (j2k_image_t * img, j2k_cp_t * cp, int curtileno)
 383 {
 384   int tileno, compno, resno, bandno, precno;
 385   tcd_img = img;
 386   tcd_cp = cp;
 387   tcd_image.tw = cp->tw;
 388   tcd_image.th = cp->th;
 389   for (tileno = 0; tileno < 1; tileno++)
 390     {
 391       /* j2k_tcp_t *tcp=&cp->tcps[curtileno]; */
 392       tile = tcd_image.tiles;
 393       for (compno = 0; compno < tile->numcomps; compno++)
 394         {
 395           tilec = &tile->comps[compno];
 396           for (resno = 0; resno < tilec->numresolutions; resno++)
 397             {
 398               res = &tilec->resolutions[resno];
 399               for (bandno = 0; bandno < res->numbands; bandno++)
 400                 {
 401                   band = &res->bands[bandno];
 402                   for (precno = 0; precno < res->pw * res->ph; precno++)
 403                     {
 404                       prc = &band->precincts[precno];
 405
 406                       if (prc->incltree != NULL)
 407                         tgt_destroy (prc->incltree);
 408                       if (prc->imsbtree != NULL)
 409                         tgt_destroy (prc->imsbtree);
 410                       free (prc->cblks);
 411                     }           /* for (precno */
 412                   free (band->precincts);
 413                 }               /* for (bandno */
 414             }                   /* for (resno */
 415           free (tilec->resolutions);
 416         }                       /* for (compno */
 417       free (tile->comps);
 418     }                           /* for (tileno */
 419   free (tcd_image.tiles);
 420 }
 421
 422 void
 423 tcd_init_encode (j2k_image_t * img, j2k_cp_t * cp, int curtileno)
 424 {
 425   int tileno, compno, resno, bandno, precno, cblkno;
 426
 427   for (tileno = 0; tileno < 1; tileno++)
 428     {
 429       j2k_tcp_t *tcp = &cp->tcps[curtileno];
 430       int j;
 431       //              int previous_x0, previous_x1, previous_y0, previous_y1;
 432       /* cfr p59 ISO/IEC FDIS15444-1 : 2000 (18 august 2000) */
 433       int p = curtileno % cp->tw;
 434       int q = curtileno / cp->tw;
 435       tile = tcd_image.tiles;
 436
 437       /* 4 borders of the tile rescale on the image if necessary */
 438       tile->x0 = int_max (cp->tx0 + p * cp->tdx, img->x0);
 439       tile->y0 = int_max (cp->ty0 + q * cp->tdy, img->y0);
 440       tile->x1 = int_min (cp->tx0 + (p + 1) * cp->tdx, img->x1);
 441       tile->y1 = int_min (cp->ty0 + (q + 1) * cp->tdy, img->y1);
 442
 443       tile->numcomps = img->numcomps;
 444       /* tile->PPT=img->PPT; */
 445       /* Modification of the RATE >> */
 446       for (j = 0; j < tcp->numlayers; j++)
 447         {
 448           tcp->rates[j] = tcp->rates[j] ? int_ceildiv (tile->numcomps * (tile->x1 - tile->x0) * (tile->y1 - tile->y0) * img->comps[0].prec, (tcp->rates[j] * 8 * img->comps[0].dx * img->comps[0].dy)) : 0;     //Mod antonin losslessbug
 449           if (tcp->rates[j])
 450             {
 451               if (j && tcp->rates[j] < tcp->rates[j - 1] + 10)
 452                 {
 453                   tcp->rates[j] = tcp->rates[j - 1] + 20;
 454                 }
 455               else
 456                 {
 457                   if (!j && tcp->rates[j] < 30)
 458                     tcp->rates[j] = 30;
 459                 }
 460             }
 461         }
 462       /* << Modification of the RATE */
 463       /* tile->comps=(tcd_tilecomp_t*)realloc(tile->comps,img->numcomps*sizeof(tcd_tilecomp_t)); */
 464       for (compno = 0; compno < tile->numcomps; compno++)
 465         {
 466           j2k_tccp_t *tccp = &tcp->tccps[compno];
 467           /* int realloc_op; */
 468
 469           tilec = &tile->comps[compno];
 470           /* border of each tile component (global) */
 471           tilec->x0 = int_ceildiv (tile->x0, img->comps[compno].dx);
 472           tilec->y0 = int_ceildiv (tile->y0, img->comps[compno].dy);
 473           tilec->x1 = int_ceildiv (tile->x1, img->comps[compno].dx);
 474           tilec->y1 = int_ceildiv (tile->y1, img->comps[compno].dy);
 475
 476           tilec->data =
 477             (int *) malloc ((tilec->x1 - tilec->x0) *
 478                             (tilec->y1 - tilec->y0) * sizeof (int));
 479           tilec->numresolutions = tccp->numresolutions;
 480           /* tilec->resolutions=(tcd_resolution_t*)realloc(tilec->resolutions,tilec->numresolutions*sizeof(tcd_resolution_t)); */
 481           for (resno = 0; resno < tilec->numresolutions; resno++)
 482             {
 483               int pdx, pdy;
 484               int levelno = tilec->numresolutions - 1 - resno;
 485               int tlprcxstart, tlprcystart, brprcxend, brprcyend;
 486               int tlcbgxstart, tlcbgystart, brcbgxend, brcbgyend;
 487               int cbgwidthexpn, cbgheightexpn;
 488               int cblkwidthexpn, cblkheightexpn;
 489
 490               res = &tilec->resolutions[resno];
 491               /* border for each resolution level (global) */
 492               res->x0 = int_ceildivpow2 (tilec->x0, levelno);
 493               res->y0 = int_ceildivpow2 (tilec->y0, levelno);
 494               res->x1 = int_ceildivpow2 (tilec->x1, levelno);
 495               res->y1 = int_ceildivpow2 (tilec->y1, levelno);
 496
 497               res->numbands = resno == 0 ? 1 : 3;
 498               /* p. 35, table A-23, ISO/IEC FDIS154444-1 : 2000 (18 august 2000) */
 499               if (tccp->csty & J2K_CCP_CSTY_PRT)
 500                 {
 501                   pdx = tccp->prcw[resno];
 502                   pdy = tccp->prch[resno];
 503                 }
 504               else
 505                 {
 506                   pdx = 15;
 507                   pdy = 15;
 508                 }
 509               /* p. 64, B.6, ISO/IEC FDIS15444-1 : 2000 (18 august 2000)  */
 510               tlprcxstart = int_floordivpow2 (res->x0, pdx) << pdx;
 511               tlprcystart = int_floordivpow2 (res->y0, pdy) << pdy;
 512               brprcxend = int_ceildivpow2 (res->x1, pdx) << pdx;
 513               brprcyend = int_ceildivpow2 (res->y1, pdy) << pdy;
 514
 515               res->pw = (brprcxend - tlprcxstart) >> pdx;
 516               res->ph = (brprcyend - tlprcystart) >> pdy;
 517
 518               if (resno == 0)
 519                 {
 520                   tlcbgxstart = tlprcxstart;
 521                   tlcbgystart = tlprcystart;
 522                   brcbgxend = brprcxend;
 523                   brcbgyend = brprcyend;
 524                   cbgwidthexpn = pdx;
 525                   cbgheightexpn = pdy;
 526                 }
 527               else
 528                 {
 529                   tlcbgxstart = int_ceildivpow2 (tlprcxstart, 1);
 530                   tlcbgystart = int_ceildivpow2 (tlprcystart, 1);
 531                   brcbgxend = int_ceildivpow2 (brprcxend, 1);
 532                   brcbgyend = int_ceildivpow2 (brprcyend, 1);
 533                   cbgwidthexpn = pdx - 1;
 534                   cbgheightexpn = pdy - 1;
 535                 }
 536
 537               cblkwidthexpn = int_min (tccp->cblkw, cbgwidthexpn);
 538               cblkheightexpn = int_min (tccp->cblkh, cbgheightexpn);
 539
 540               for (bandno = 0; bandno < res->numbands; bandno++)
 541                 {
 542                   int x0b, y0b;
 543                   int gain, numbps;
 544                   j2k_stepsize_t *ss;
 545                   band = &res->bands[bandno];
 546                   band->bandno = resno == 0 ? 0 : bandno + 1;
 547                   x0b = (band->bandno == 1) || (band->bandno == 3) ? 1 : 0;
 548                   y0b = (band->bandno == 2) || (band->bandno == 3) ? 1 : 0;
 549
 550                   if (band->bandno == 0)
 551                     {
 552                       /* band border */
 553                       band->x0 = int_ceildivpow2 (tilec->x0, levelno);
 554                       band->y0 = int_ceildivpow2 (tilec->y0, levelno);
 555                       band->x1 = int_ceildivpow2 (tilec->x1, levelno);
 556                       band->y1 = int_ceildivpow2 (tilec->y1, levelno);
 557                     }
 558                   else
 559                     {
 560                       band->x0 =
 561                         int_ceildivpow2 (tilec->x0 -
 562                                          (1 << levelno) * x0b, levelno + 1);
 563                       band->y0 =
 564                         int_ceildivpow2 (tilec->y0 -
 565                                          (1 << levelno) * y0b, levelno + 1);
 566                       band->x1 =
 567                         int_ceildivpow2 (tilec->x1 -
 568                                          (1 << levelno) * x0b, levelno + 1);
 569                       band->y1 =
 570                         int_ceildivpow2 (tilec->y1 -
 571                                          (1 << levelno) * y0b, levelno + 1);
 572                     }
 573
 574                   ss = &tccp->stepsizes[resno ==
 575                                         0 ? 0 : 3 * (resno - 1) + bandno + 1];
 576                   gain =
 577                     tccp->qmfbid ==
 578                     0 ? dwt_getgain_real (band->bandno) : dwt_getgain (band->
 579                                                                        bandno);
 580                   numbps = img->comps[compno].prec + gain;
 581                   band->stepsize =
 582                     (int) floor ((1.0 + ss->mant / 2048.0) *
 583                                  pow (2.0, numbps - ss->expn) * 8192.0);
 584                   band->numbps = ss->expn + tccp->numgbits - 1; /* WHY -1 ? */
 585
 586                   for (precno = 0; precno < res->pw * res->ph; precno++)
 587                     {
 588                       int tlcblkxstart, tlcblkystart, brcblkxend, brcblkyend;
 589                       int cbgxstart =
 590                         tlcbgxstart +
 591                         (precno % res->pw) * (1 << cbgwidthexpn);
 592                       int cbgystart =
 593                         tlcbgystart +
 594                         (precno / res->pw) * (1 << cbgheightexpn);
 595                       int cbgxend = cbgxstart + (1 << cbgwidthexpn);
 596                       int cbgyend = cbgystart + (1 << cbgheightexpn);
 597
 598                       prc = &band->precincts[precno];
 599                       /* precinct size (global) */
 600                       prc->x0 = int_max (cbgxstart, band->x0);
 601                       prc->y0 = int_max (cbgystart, band->y0);
 602                       prc->x1 = int_min (cbgxend, band->x1);
 603                       prc->y1 = int_min (cbgyend, band->y1);
 604
 605                       tlcblkxstart =
 606                         int_floordivpow2 (prc->x0,
 607                                           cblkwidthexpn) << cblkwidthexpn;
 608                       tlcblkystart =
 609                         int_floordivpow2 (prc->y0,
 610                                           cblkheightexpn) << cblkheightexpn;
 611                       brcblkxend =
 612                         int_ceildivpow2 (prc->x1,
 613                                          cblkwidthexpn) << cblkwidthexpn;
 614                       brcblkyend =
 615                         int_ceildivpow2 (prc->y1,
 616                                          cblkheightexpn) << cblkheightexpn;
 617                       prc->cw = (brcblkxend - tlcblkxstart) >> cblkwidthexpn;
 618                       prc->ch = (brcblkyend - tlcblkystart) >> cblkheightexpn;
 619
 620                       free (prc->cblks);
 621                       prc->cblks =
 622                         (tcd_cblk_t *) malloc (prc->cw * prc->ch *
 623                                                sizeof (tcd_cblk_t));
 624
 625                       if (prc->incltree != NULL)
 626                         tgt_destroy (prc->incltree);
 627                       if (prc->imsbtree != NULL)
 628                         tgt_destroy (prc->imsbtree);
 629
 630                       prc->incltree = tgt_create (prc->cw, prc->ch);
 631                       prc->imsbtree = tgt_create (prc->cw, prc->ch);
 632
 633                       for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++)
 634                         {
 635                           int cblkxstart =
 636                             tlcblkxstart +
 637                             (cblkno % prc->cw) * (1 << cblkwidthexpn);
 638                           int cblkystart =
 639                             tlcblkystart +
 640                             (cblkno / prc->cw) * (1 << cblkheightexpn);
 641                           int cblkxend = cblkxstart + (1 << cblkwidthexpn);
 642                           int cblkyend = cblkystart + (1 << cblkheightexpn);
 643                           cblk = &prc->cblks[cblkno];
 644
 645                           /* code-block size (global) */
 646                           cblk->x0 = int_max (cblkxstart, prc->x0);
 647                           cblk->y0 = int_max (cblkystart, prc->y0);
 648                           cblk->x1 = int_min (cblkxend, prc->x1);
 649                           cblk->y1 = int_min (cblkyend, prc->y1);
 650
 651                         }
 652                     }
 653                 }
 654             }
 655         }
 656     }
 657   /* tcd_dump(&tcd_image,0); */
 658 }
 659
 660 void
 661 tcd_init (j2k_image_t * img, j2k_cp_t * cp)
 662 {
 663   int tileno, compno, resno, bandno, precno, cblkno, i, j;
 664   unsigned int x0 = 0, y0 = 0, x1 = 0, y1 = 0, w, h, p, q;
 665   tcd_img = img;
 666   tcd_cp = cp;
 667   tcd_image.tw = cp->tw;
 668   tcd_image.th = cp->th;
 669   tcd_image.tiles =
 670     (tcd_tile_t *) malloc (cp->tw * cp->th * sizeof (tcd_tile_t));
 671
 672   /*for (tileno = 0; tileno < cp->tw * cp->th; tileno++) {
 673      j2k_tcp_t *tcp = &cp->tcps[tileno];
 674      tcd_tile_t *tile = &tcd_image.tiles[tileno]; */
 675
 676   for (i = 0; i < cp->tileno_size; i++)
 677     {
 678       j2k_tcp_t *tcp = &cp->tcps[cp->tileno[i]];
 679       tcd_tile_t *tile = &tcd_image.tiles[cp->tileno[i]];
 680       tileno = cp->tileno[i];
 681
 682
 683       //              int previous_x0, previous_x1, previous_y0, previous_y1;
 684       /* cfr p59 ISO/IEC FDIS15444-1 : 2000 (18 august 2000) */
 685       p = tileno % cp->tw;      /* si numerotation matricielle .. */
 686       q = tileno / cp->tw;      /* .. coordonnees de la tile (q,p) q pour ligne et p pour colonne */
 687
 688       /* 4 borders of the tile rescale on the image if necessary */
 689       tile->x0 = int_max (cp->tx0 + p * cp->tdx, img->x0);
 690       tile->y0 = int_max (cp->ty0 + q * cp->tdy, img->y0);
 691       tile->x1 = int_min (cp->tx0 + (p + 1) * cp->tdx, img->x1);
 692       tile->y1 = int_min (cp->ty0 + (q + 1) * cp->tdy, img->y1);
 693
 694       tile->numcomps = img->numcomps;
 695       tile->comps =
 696         (tcd_tilecomp_t *) malloc (img->numcomps * sizeof (tcd_tilecomp_t));
 697       for (compno = 0; compno < tile->numcomps; compno++)
 698         {
 699           j2k_tccp_t *tccp = &tcp->tccps[compno];
 700           tcd_tilecomp_t *tilec = &tile->comps[compno];
 701           /* border of each tile component (global) */
 702           tilec->x0 = int_ceildiv (tile->x0, img->comps[compno].dx);
 703           tilec->y0 = int_ceildiv (tile->y0, img->comps[compno].dy);
 704           tilec->x1 = int_ceildiv (tile->x1, img->comps[compno].dx);
 705           tilec->y1 = int_ceildiv (tile->y1, img->comps[compno].dy);
 706
 707           tilec->data =
 708             (int *) malloc ((tilec->x1 - tilec->x0) *
 709                             (tilec->y1 - tilec->y0) * sizeof (int));
 710           tilec->numresolutions = tccp->numresolutions;
 711           tilec->resolutions =
 712             (tcd_resolution_t *) malloc (tilec->numresolutions *
 713                                          sizeof (tcd_resolution_t));
 714           for (resno = 0; resno < tilec->numresolutions; resno++)
 715             {
 716               int pdx, pdy;
 717               int levelno = tilec->numresolutions - 1 - resno;
 718               int tlprcxstart, tlprcystart, brprcxend, brprcyend;
 719               int tlcbgxstart, tlcbgystart, brcbgxend, brcbgyend;
 720               int cbgwidthexpn, cbgheightexpn;
 721               int cblkwidthexpn, cblkheightexpn;
 722               tcd_resolution_t *res = &tilec->resolutions[resno];
 723
 724               /* border for each resolution level (global) */
 725               res->x0 = int_ceildivpow2 (tilec->x0, levelno);
 726               res->y0 = int_ceildivpow2 (tilec->y0, levelno);
 727               res->x1 = int_ceildivpow2 (tilec->x1, levelno);
 728               res->y1 = int_ceildivpow2 (tilec->y1, levelno);
 729
 730               res->numbands = resno == 0 ? 1 : 3;
 731               /* p. 35, table A-23, ISO/IEC FDIS154444-1 : 2000 (18 august 2000) */
 732               if (tccp->csty & J2K_CCP_CSTY_PRT)
 733                 {
 734                   pdx = tccp->prcw[resno];
 735                   pdy = tccp->prch[resno];
 736                 }
 737               else
 738                 {
 739                   pdx = 15;
 740                   pdy = 15;
 741                 }
 742               /* p. 64, B.6, ISO/IEC FDIS15444-1 : 2000 (18 august 2000)  */
 743               tlprcxstart = int_floordivpow2 (res->x0, pdx) << pdx;
 744               tlprcystart = int_floordivpow2 (res->y0, pdy) << pdy;
 745               brprcxend = int_ceildivpow2 (res->x1, pdx) << pdx;
 746               brprcyend = int_ceildivpow2 (res->y1, pdy) << pdy;
 747               res->pw = (res->x0 == res->x1) ? 0 : ((brprcxend - tlprcxstart) >> pdx);  // Mod Antonin : sizebug1
 748               res->ph = (res->y0 == res->y1) ? 0 : ((brprcyend - tlprcystart) >> pdy);  // Mod Antonin : sizebug1
 749
 750               if (resno == 0)
 751                 {
 752                   tlcbgxstart = tlprcxstart;
 753                   tlcbgystart = tlprcystart;
 754                   brcbgxend = brprcxend;
 755                   brcbgyend = brprcyend;
 756                   cbgwidthexpn = pdx;
 757                   cbgheightexpn = pdy;
 758                 }
 759               else
 760                 {
 761                   tlcbgxstart = int_ceildivpow2 (tlprcxstart, 1);
 762                   tlcbgystart = int_ceildivpow2 (tlprcystart, 1);
 763                   brcbgxend = int_ceildivpow2 (brprcxend, 1);
 764                   brcbgyend = int_ceildivpow2 (brprcyend, 1);
 765                   cbgwidthexpn = pdx - 1;
 766                   cbgheightexpn = pdy - 1;
 767                 }
 768
 769               cblkwidthexpn = int_min (tccp->cblkw, cbgwidthexpn);
 770               cblkheightexpn = int_min (tccp->cblkh, cbgheightexpn);
 771
 772               for (bandno = 0; bandno < res->numbands; bandno++)
 773                 {
 774                   int x0b, y0b;
 775                   int gain, numbps;
 776                   j2k_stepsize_t *ss;
 777                   tcd_band_t *band = &res->bands[bandno];
 778                   band->bandno = resno == 0 ? 0 : bandno + 1;
 779                   x0b = (band->bandno == 1) || (band->bandno == 3) ? 1 : 0;
 780                   y0b = (band->bandno == 2) || (band->bandno == 3) ? 1 : 0;
 781
 782                   if (band->bandno == 0)
 783                     {
 784                       /* band border (global) */
 785                       band->x0 = int_ceildivpow2 (tilec->x0, levelno);
 786                       band->y0 = int_ceildivpow2 (tilec->y0, levelno);
 787                       band->x1 = int_ceildivpow2 (tilec->x1, levelno);
 788                       band->y1 = int_ceildivpow2 (tilec->y1, levelno);
 789                     }
 790                   else
 791                     {
 792                       /* band border (global) */
 793                       band->x0 =
 794                         int_ceildivpow2 (tilec->x0 -
 795                                          (1 << levelno) * x0b, levelno + 1);
 796                       band->y0 =
 797                         int_ceildivpow2 (tilec->y0 -
 798                                          (1 << levelno) * y0b, levelno + 1);
 799                       band->x1 =
 800                         int_ceildivpow2 (tilec->x1 -
 801                                          (1 << levelno) * x0b, levelno + 1);
 802                       band->y1 =
 803                         int_ceildivpow2 (tilec->y1 -
 804                                          (1 << levelno) * y0b, levelno + 1);
 805                     }
 806
 807                   ss = &tccp->stepsizes[resno ==
 808                                         0 ? 0 : 3 * (resno - 1) + bandno + 1];
 809                   gain =
 810                     tccp->qmfbid ==
 811                     0 ? dwt_getgain_real (band->bandno) : dwt_getgain (band->
 812                                                                        bandno);
 813                   numbps = img->comps[compno].prec + gain;
 814                   band->stepsize =
 815                     (int) floor ((1.0 + ss->mant / 2048.0) *
 816                                  pow (2.0, numbps - ss->expn) * 8192.0);
 817                   band->numbps = ss->expn + tccp->numgbits - 1; /* WHY -1 ? */
 818
 819                   band->precincts =
 820                     (tcd_precinct_t *) malloc (res->pw * res->ph *
 821                                                sizeof (tcd_precinct_t));
 822
 823                   for (precno = 0; precno < res->pw * res->ph; precno++)
 824                     {
 825                       int tlcblkxstart, tlcblkystart, brcblkxend, brcblkyend;
 826                       int cbgxstart =
 827                         tlcbgxstart +
 828                         (precno % res->pw) * (1 << cbgwidthexpn);
 829                       int cbgystart =
 830                         tlcbgystart +
 831                         (precno / res->pw) * (1 << cbgheightexpn);
 832                       int cbgxend = cbgxstart + (1 << cbgwidthexpn);
 833                       int cbgyend = cbgystart + (1 << cbgheightexpn);
 834                       tcd_precinct_t *prc = &band->precincts[precno];
 835                       /* precinct size (global) */
 836                       prc->x0 = int_max (cbgxstart, band->x0);
 837                       prc->y0 = int_max (cbgystart, band->y0);
 838                       prc->x1 = int_min (cbgxend, band->x1);
 839                       prc->y1 = int_min (cbgyend, band->y1);
 840
 841                       tlcblkxstart =
 842                         int_floordivpow2 (prc->x0,
 843                                           cblkwidthexpn) << cblkwidthexpn;
 844                       tlcblkystart =
 845                         int_floordivpow2 (prc->y0,
 846                                           cblkheightexpn) << cblkheightexpn;
 847                       brcblkxend =
 848                         int_ceildivpow2 (prc->x1,
 849                                          cblkwidthexpn) << cblkwidthexpn;
 850                       brcblkyend =
 851                         int_ceildivpow2 (prc->y1,
 852                                          cblkheightexpn) << cblkheightexpn;
 853                       prc->cw = (brcblkxend - tlcblkxstart) >> cblkwidthexpn;
 854                       prc->ch = (brcblkyend - tlcblkystart) >> cblkheightexpn;
 855
 856                       prc->cblks =
 857                         (tcd_cblk_t *) malloc (prc->cw * prc->ch *
 858                                                sizeof (tcd_cblk_t));
 859
 860                       prc->incltree = tgt_create (prc->cw, prc->ch);
 861                       prc->imsbtree = tgt_create (prc->cw, prc->ch);
 862
 863                       for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++)
 864                         {
 865                           int cblkxstart =
 866                             tlcblkxstart +
 867                             (cblkno % prc->cw) * (1 << cblkwidthexpn);
 868                           int cblkystart =
 869                             tlcblkystart +
 870                             (cblkno / prc->cw) * (1 << cblkheightexpn);
 871                           int cblkxend = cblkxstart + (1 << cblkwidthexpn);
 872                           int cblkyend = cblkystart + (1 << cblkheightexpn);
 873                           tcd_cblk_t *cblk = &prc->cblks[cblkno];
 874                           /* code-block size (global) */
 875                           cblk->x0 = int_max (cblkxstart, prc->x0);
 876                           cblk->y0 = int_max (cblkystart, prc->y0);
 877                           cblk->x1 = int_min (cblkxend, prc->x1);
 878                           cblk->y1 = int_min (cblkyend, prc->y1);
 879
 880
 881
 882                           cblk->lastbp = 0;     // Add Antonin : quantizbug1
 883                         }
 884                     }
 885                 }
 886             }
 887         }
 888     }
 889   //tcd_dump(&tcd_image,0);
 890
 891
 892   /* Allocate place to store the data decoded = final image */
 893   /* Place limited by the tile really present in the codestream */
 894
 895
 896   for (i = 0; i < img->numcomps; i++)
 897     {
 898       for (j = 0; j < cp->tileno_size; j++)
 899         {
 900           tileno = cp->tileno[j];
 901           x0 = j == 0 ? tcd_image.tiles[tileno].comps[i].x0 : int_min (x0,
 902                                                                        tcd_image.
 903                                                                        tiles
 904                                                                        [tileno].
 905                                                                        comps
 906                                                                        [i].
 907                                                                        x0);
 908           y0 =
 909             j == 0 ? tcd_image.tiles[tileno].comps[i].y0 : int_min (y0,
 910                                                                     tcd_image.
 911                                                                     tiles
 912                                                                     [tileno].
 913                                                                     comps[i].
 914                                                                     y0);
 915           x1 =
 916             j == 0 ? tcd_image.tiles[tileno].comps[i].x1 : int_max (x1,
 917                                                                     tcd_image.
 918                                                                     tiles
 919                                                                     [tileno].
 920                                                                     comps[i].
 921                                                                     x1);
 922           y1 =
 923             j == 0 ? tcd_image.tiles[tileno].comps[i].y1 : int_max (y1,
 924                                                                     tcd_image.
 925                                                                     tiles
 926                                                                     [tileno].
 927                                                                     comps[i].
 928                                                                     y1);
 929         }
 930       //w = int_ceildiv(x1 - x0, img->comps[i].dx);
 931       //h = int_ceildiv(y1 - y0, img->comps[i].dy);
 932
 933       w = x1 - x0;
 934
 935       h = y1 - y0;
 936       img->comps[i].data = (int *) calloc (w * h, sizeof (int));
 937       img->comps[i].w = w;
 938       img->comps[i].h = h;
 939       img->comps[i].x0 = x0;
 940       img->comps[i].y0 = y0;
 941     }
 942 }
 943
 944 void
 945 tcd_makelayer_fixed (int layno, int final)
 946 {
 947   int compno, resno, bandno, precno, cblkno;
 948   int value;                    //, matrice[tcd_tcp->numlayers][tcd_tile->comps[0].numresolutions][3];
 949   int matrice[10][10][3];
 950   int i, j, k;
 951
 952   /*matrice=(int*)malloc(tcd_tcp->numlayers*tcd_tile->comps[0].numresolutions*3*sizeof(int)); */
 953
 954   for (compno = 0; compno < tcd_tile->numcomps; compno++)
 955     {
 956       tcd_tilecomp_t *tilec = &tcd_tile->comps[compno];
 957       for (i = 0; i < tcd_tcp->numlayers; i++)
 958         {
 959           for (j = 0; j < tilec->numresolutions; j++)
 960             {
 961               for (k = 0; k < 3; k++)
 962                 {
 963                   matrice[i][j][k] =
 964                     (int) (tcd_cp->
 965                            matrice[i * tilec->numresolutions * 3 +
 966                                    j * 3 +
 967                                    k] *
 968                            (float) (tcd_img->comps[compno].prec / 16.0));
 969         }}}
 970
 971       for (resno = 0; resno < tilec->numresolutions; resno++)
 972         {
 973           tcd_resolution_t *res = &tilec->resolutions[resno];
 974           for (bandno = 0; bandno < res->numbands; bandno++)
 975             {
 976               tcd_band_t *band = &res->bands[bandno];
 977               for (precno = 0; precno < res->pw * res->ph; precno++)
 978                 {
 979                   tcd_precinct_t *prc = &band->precincts[precno];
 980                   for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++)
 981                     {
 982                       tcd_cblk_t *cblk = &prc->cblks[cblkno];
 983                       tcd_layer_t *layer = &cblk->layers[layno];
 984                       int n;
 985                       int imsb = tcd_img->comps[compno].prec - cblk->numbps;    /* number of bit-plan equal to zero */
 986                       /* Correction of the matrix of coefficient to include the IMSB information */
 987
 988                       if (layno == 0)
 989                         {
 990                           value = matrice[layno][resno][bandno];
 991                           if (imsb >= value)
 992                             value = 0;
 993                           else
 994                             value -= imsb;
 995                         }
 996                       else
 997                         {
 998                           value =
 999                             matrice[layno][resno][bandno] -
1000                             matrice[layno - 1][resno][bandno];
1001                           if (imsb >= matrice[layno - 1][resno][bandno])
1002                             {
1003                               value -=
1004                                 (imsb - matrice[layno - 1][resno][bandno]);
1005                               if (value < 0)
1006                                 value = 0;
1007                             }
1008                         }
1009
1010                       if (layno == 0)
1011                         cblk->numpassesinlayers = 0;
1012
1013                       n = cblk->numpassesinlayers;
1014                       if (cblk->numpassesinlayers == 0)
1015                         {
1016                           if (value != 0)
1017                             n = 3 * value - 2 + cblk->numpassesinlayers;
1018                           else
1019                             n = cblk->numpassesinlayers;
1020                         }
1021                       else
1022                         n = 3 * value + cblk->numpassesinlayers;
1023
1024                       layer->numpasses = n - cblk->numpassesinlayers;
1025
1026                       if (!layer->numpasses)
1027                         continue;
1028
1029                       if (cblk->numpassesinlayers == 0)
1030                         {
1031                           layer->len = cblk->passes[n - 1].rate;
1032                           layer->data = cblk->data;
1033                         }
1034                       else
1035                         {
1036                           layer->len =
1037                             cblk->passes[n - 1].rate -
1038                             cblk->passes[cblk->numpassesinlayers - 1].rate;
1039                           layer->data =
1040                             cblk->data +
1041                             cblk->passes[cblk->numpassesinlayers - 1].rate;
1042                         }
1043                       if (final)
1044                         cblk->numpassesinlayers = n;
1045                     }
1046                 }
1047             }
1048         }
1049     }
1050 }
1051
1052 void
1053 tcd_rateallocate_fixed ()
1054 {
1055   int layno;
1056
1057   for (layno = 0; layno < tcd_tcp->numlayers; layno++)
1058     {
1059       tcd_makelayer_fixed (layno, 1);
1060     }
1061 }
1062
1063 void
1064 tcd_makelayer (int layno, double thresh, int final)
1065 {
1066   int compno, resno, bandno, precno, cblkno, passno;
1067
1068   tcd_tile->distolayer[layno] = 0;      //add fixed_quality
1069
1070   for (compno = 0; compno < tcd_tile->numcomps; compno++)
1071     {
1072       tcd_tilecomp_t *tilec = &tcd_tile->comps[compno];
1073       for (resno = 0; resno < tilec->numresolutions; resno++)
1074         {
1075           tcd_resolution_t *res = &tilec->resolutions[resno];
1076           for (bandno = 0; bandno < res->numbands; bandno++)
1077             {
1078               tcd_band_t *band = &res->bands[bandno];
1079               for (precno = 0; precno < res->pw * res->ph; precno++)
1080                 {
1081                   tcd_precinct_t *prc = &band->precincts[precno];
1082                   for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++)
1083                     {
1084                       tcd_cblk_t *cblk = &prc->cblks[cblkno];
1085                       tcd_layer_t *layer = &cblk->layers[layno];
1086                       int n;
1087
1088                       if (layno == 0)
1089                         {
1090                           cblk->numpassesinlayers = 0;
1091                         }
1092                       n = cblk->numpassesinlayers;
1093                       for (passno = cblk->numpassesinlayers;
1094                            passno < cblk->totalpasses; passno++)
1095                         {
1096                           int dr;
1097                           double dd;
1098                           tcd_pass_t *pass = &cblk->passes[passno];
1099                           if (n == 0)
1100                             {
1101                               dr = pass->rate;
1102                               dd = pass->distortiondec;
1103                             }
1104                           else
1105                             {
1106                               dr = pass->rate - cblk->passes[n - 1].rate;
1107                               dd = pass->distortiondec - cblk->passes[n -
1108                                                                       1].
1109                                 distortiondec;
1110                             }
1111                           if (!dr)
1112                             {
1113                               if (dd)
1114                                 n = passno + 1;
1115                               continue;
1116                             }
1117
1118                           if (dd / dr >= thresh)
1119                             n = passno + 1;
1120                         }
1121                       layer->numpasses = n - cblk->numpassesinlayers;
1122
1123                       if (!layer->numpasses)
1124                         {
1125                           layer->disto = 0;
1126                           continue;
1127                         }
1128
1129                       if (cblk->numpassesinlayers == 0)
1130                         {
1131                           layer->len = cblk->passes[n - 1].rate;
1132                           layer->data = cblk->data;
1133                           layer->disto = cblk->passes[n - 1].distortiondec;
1134                         }
1135                       else
1136                         {
1137                           layer->len = cblk->passes[n - 1].rate -
1138                             cblk->passes[cblk->numpassesinlayers - 1].rate;
1139                           layer->data =
1140                             cblk->data +
1141                             cblk->passes[cblk->numpassesinlayers - 1].rate;
1142                           layer->disto =
1143                             cblk->passes[n - 1].distortiondec -
1144                             cblk->passes[cblk->numpassesinlayers -
1145                                          1].distortiondec;
1146                         }
1147
1148                       tcd_tile->distolayer[layno] += layer->disto;      //add fixed_quality
1149
1150                       if (final)
1151                         cblk->numpassesinlayers = n;
1152                     }
1153                 }
1154             }
1155         }
1156     }
1157 }
1158
1159 void
1160 tcd_rateallocate (unsigned char *dest, int len, info_image * info_IM)
1161 {
1162   int compno, resno, bandno, precno, cblkno, passno, layno;
1163   double min, max;
1164   double cumdisto[100];         //add fixed_quality
1165   const double K = 1;           // 1.1; //add fixed_quality
1166
1167   double maxSE = 0;
1168   min = DBL_MAX;
1169   max = 0;
1170
1171   tcd_tile->nbpix = 0;          //add fixed_quality
1172
1173   for (compno = 0; compno < tcd_tile->numcomps; compno++)
1174     {
1175       tcd_tilecomp_t *tilec = &tcd_tile->comps[compno];
1176
1177       tilec->nbpix = 0;
1178       for (resno = 0; resno < tilec->numresolutions; resno++)
1179         {
1180           tcd_resolution_t *res = &tilec->resolutions[resno];
1181           for (bandno = 0; bandno < res->numbands; bandno++)
1182             {
1183               tcd_band_t *band = &res->bands[bandno];
1184               for (precno = 0; precno < res->pw * res->ph; precno++)
1185                 {
1186                   tcd_precinct_t *prc = &band->precincts[precno];
1187                   for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++)
1188                     {
1189                       tcd_cblk_t *cblk = &prc->cblks[cblkno];
1190                       for (passno = 0; passno < cblk->totalpasses; passno++)
1191                         {
1192                           tcd_pass_t *pass = &cblk->passes[passno];
1193                           int dr;
1194                           double dd, rdslope;
1195                           if (passno == 0)
1196                             {
1197                               dr = pass->rate;
1198                               dd = pass->distortiondec;
1199                             }
1200                           else
1201                             {
1202                               dr = pass->rate - cblk->passes[passno - 1].rate;
1203                               dd = pass->distortiondec -
1204                                 cblk->passes[passno - 1].distortiondec;
1205                             }
1206                           if (dr == 0)
1207                             {
1208                               continue;
1209                             }
1210
1211                           rdslope = dd / dr;
1212
1213                           if (rdslope < min)
1214                             {
1215                               min = rdslope;
1216                             }
1217                           if (rdslope > max)
1218                             {
1219                               max = rdslope;
1220                             }
1221                         }       /* passno */
1222
1223                       tcd_tile->nbpix += ((cblk->x1 - cblk->x0) * (cblk->y1 - cblk->y0));       //add fixed_quality
1224
1225                       tilec->nbpix += ((cblk->x1 - cblk->x0) * (cblk->y1 - cblk->y0));  //add fixed_quality
1226
1227                     }           /* cbklno */
1228                 }               /* precno */
1229             }                   /* bandno */
1230         }                       /* resno */
1231
1232       maxSE +=
1233         (double) (((1 << tcd_img->comps[compno].prec) -
1234                    1) * ((1 << tcd_img->comps[compno].prec) -
1235                          1)) * (tilec->nbpix);
1236     }                           /* compno */
1237
1238   /* add antonin index */
1239   if (info_IM->index_on)
1240     {
1241       info_tile *info_TL = &info_IM->tile[tcd_tileno];
1242       info_TL->nbpix = tcd_tile->nbpix;
1243       info_TL->distotile = tcd_tile->distotile;
1244       info_TL->thresh =
1245         (double *) malloc (tcd_tcp->numlayers * sizeof (double));
1246     }
1247   /* dda */
1248
1249   for (layno = 0; layno < tcd_tcp->numlayers; layno++)
1250     {
1251       volatile double lo = min;
1252       volatile double hi = max;
1253       volatile int success = 0;
1254       volatile int maxlen = tcd_tcp->rates[layno] ? int_min (tcd_tcp->rates[layno], len) : len; //Mod antonin losslessbug
1255       volatile double goodthresh;
1256       volatile int i;
1257       double distotarget;       //add fixed_quality
1258
1259       distotarget = tcd_tile->distotile - ((K * maxSE) / pow (10, tcd_tcp->distoratio[layno] / 10));    // add fixed_quality
1260
1261       if (tcd_tcp->rates[layno])
1262         {
1263           for (i = 0; i < 32; i++)
1264             {
1265               volatile double thresh = (lo + hi) / 2;
1266               int l = 0;
1267               double distoachieved = 0; // add fixed_quality
1268
1269               tcd_makelayer (layno, thresh, 0);
1270
1271               if (tcd_cp->fixed_quality)
1272                 {               // add fixed_quality
1273                   distoachieved =
1274                     layno ==
1275                     0 ? tcd_tile->distolayer[0] : cumdisto[layno - 1] +
1276                     tcd_tile->distolayer[layno];
1277                   if (distoachieved < distotarget)
1278                     {
1279                       hi = thresh;
1280                       continue;
1281                     }
1282                   lo = thresh;
1283                 }
1284               else
1285                 {
1286                   l =
1287                     t2_encode_packets (tcd_img, tcd_cp, tcd_tileno, tcd_tile,
1288                                        layno + 1, dest, maxlen, info_IM);
1289                   /* fprintf(stderr, "rate alloc: len=%d, max=%d\n", l, maxlen); */
1290                   if (l == -999)
1291                     {
1292                       lo = thresh;
1293                       continue;
1294                     }
1295                   hi = thresh;
1296                 }
1297
1298               success = 1;
1299               goodthresh = thresh;
1300             }
1301         }
1302       else
1303         {
1304           success = 1;
1305           goodthresh = min;
1306         }
1307
1308       if (!success)
1309         {
1310           longjmp (j2k_error, 1);
1311         }
1312
1313       if (info_IM->index_on)
1314         {                       /* Threshold for Marcela Index */
1315           info_IM->tile[tcd_tileno].thresh[layno] = goodthresh;
1316         }
1317       tcd_makelayer (layno, goodthresh, 1);
1318
1319       cumdisto[layno] = layno == 0 ? tcd_tile->distolayer[0] : cumdisto[layno - 1] + tcd_tile->distolayer[layno];       // add fixed_quality
1320     }
1321 }
1322
1323 int
1324 tcd_encode_tile_pxm (int tileno, unsigned char *dest, int len,
1325                      info_image * info_IM)
1326 {
1327   int compno;
1328   int l, i;
1329   clock_t time7;
1330   tcd_tile_t *tile;
1331   j2k_tcp_t *tcp = &tcd_cp->tcps[0];
1332   j2k_tccp_t *tccp = &tcp->tccps[0];
1333
1334   tcd_tileno = tileno;
1335   tcd_tile = tcd_image.tiles;
1336   tcd_tcp = &tcd_cp->tcps[tileno];
1337   tile = tcd_tile;
1338   /* INDEX >> "Precinct_nb_X et Precinct_nb_Y" */
1339   if (info_IM->index_on)
1340     {
1341       tcd_tilecomp_t *tilec_idx = &tile->comps[0];      //Based on Component 0
1342
1343       for (i = 0; i < tilec_idx->numresolutions; i++)
1344         {
1345
1346           tcd_resolution_t *res_idx = &tilec_idx->resolutions[i];
1347
1348           info_IM->tile[tileno].pw[i] = res_idx->pw;
1349           info_IM->tile[tileno].ph[i] = res_idx->ph;
1350
1351           info_IM->tile[tileno].pdx[i] = tccp->prcw[i];
1352           info_IM->tile[tileno].pdy[i] = tccp->prch[i];
1353
1354         }
1355     }
1356   /* << INDEX */
1357
1358 /*---------------TILE-------------------*/
1359
1360   time7 = clock ();
1361
1362   for (compno = 0; compno < tile->numcomps; compno++)
1363     {
1364       FILE *src;
1365       char tmp[256];
1366       int k;
1367       unsigned char elmt;
1368       int i, j;
1369       int tw, w;
1370       tcd_tilecomp_t *tilec = &tile->comps[compno];
1371       int adjust =
1372         tcd_img->comps[compno].sgnd ? 0 : 1 << (tcd_img->comps[compno].
1373                                                 prec - 1);
1374       int offset_x, offset_y;
1375
1376       offset_x = int_ceildiv (tcd_img->x0, tcd_img->comps[compno].dx);
1377       offset_y = int_ceildiv (tcd_img->y0, tcd_img->comps[compno].dy);
1378       tw = tilec->x1 - tilec->x0;
1379       w = int_ceildiv (tcd_img->x1 - tcd_img->x0, tcd_img->comps[compno].dx);
1380       sprintf (tmp, "Compo%d", compno); /* component file */
1381       src = fopen (tmp, "rb");
1382       if (!src)
1383         {
1384           fprintf (stderr, "failed to open %s for reading\n", tmp);
1385           return 1;
1386         }
1387
1388       /* read the Compo file to extract data of the tile */
1389       k = 0;
1390       fseek (src, (tilec->x0 - offset_x) + (tilec->y0 - offset_y) * w,
1391              SEEK_SET);
1392       k = (tilec->x0 - offset_x) + (tilec->y0 - offset_y) * w;
1393       for (j = tilec->y0; j < tilec->y1; j++)
1394         {
1395           for (i = tilec->x0; i < tilec->x1; i++)
1396             {
1397               if (tcd_tcp->tccps[compno].qmfbid == 1)
1398                 {
1399                   elmt = fgetc (src);
1400                   tilec->data[i - tilec->x0 + (j - tilec->y0) * tw] =
1401                     elmt - adjust;
1402                   k++;
1403                 }
1404               else if (tcd_tcp->tccps[compno].qmfbid == 0)
1405                 {
1406                   elmt = fgetc (src);
1407                   tilec->data[i - tilec->x0 + (j - tilec->y0) * tw] =
1408                     (elmt - adjust) << 13;
1409                   k++;
1410                 }
1411             }
1412           fseek (src, (tilec->x0 - offset_x) + (j + 1 - offset_y) * w - k,
1413                  SEEK_CUR);
1414           k = tilec->x0 - offset_x + (j + 1 - offset_y) * w;
1415
1416         }
1417       fclose (src);
1418     }
1419
1420 /*----------------MCT-------------------*/
1421
1422   if (tcd_tcp->mct)
1423     {
1424       if (tcd_tcp->tccps[0].qmfbid == 0)
1425         {
1426           mct_encode_real (tile->comps[0].data, tile->comps[1].data,
1427                            tile->comps[2].data,
1428                            (tile->comps[0].x1 -
1429                             tile->comps[0].x0) * (tile->comps[0].y1 -
1430                                                   tile->comps[0].y0));
1431         }
1432       else
1433         {
1434           mct_encode (tile->comps[0].data, tile->comps[1].data,
1435                       tile->comps[2].data,
1436                       (tile->comps[0].x1 -
1437                        tile->comps[0].x0) * (tile->comps[0].y1 -
1438                                              tile->comps[0].y0));
1439         }
1440     }
1441 /*----------------DWT---------------------*/
1442
1443   /* time3=clock(); */
1444   for (compno = 0; compno < tile->numcomps; compno++)
1445     {
1446       tcd_tilecomp_t *tilec = &tile->comps[compno];
1447       if (tcd_tcp->tccps[compno].qmfbid == 1)
1448         {
1449           dwt_encode (tilec->data, tilec->x1 - tilec->x0,
1450                       tilec->y1 - tilec->y0, tilec,
1451                       tilec->numresolutions - 1);
1452         }
1453       else if (tcd_tcp->tccps[compno].qmfbid == 0)
1454         {
1455           dwt_encode_real (tilec->data, tilec->x1 - tilec->x0,
1456                            tilec->y1 - tilec->y0, tilec,
1457                            tilec->numresolutions - 1);
1458         }
1459     }
1460 /*------------------TIER1-----------------*/
1461
1462   t1_init_luts ();
1463   t1_encode_cblks (tile, tcd_tcp);
1464
1465 /*-----------RATE-ALLOCATE------------------*/
1466   info_IM->index_write = 0;     /* INDEX     */
1467
1468   if (tcd_cp->disto_alloc || tcd_cp->fixed_quality)     // mod fixed_quality
1469     /* Normal Rate/distortion allocation */
1470     tcd_rateallocate (dest, len, info_IM);
1471   else
1472     /* Fixed layer allocation */
1473     tcd_rateallocate_fixed ();
1474
1475 /*--------------TIER2------------------*/
1476   info_IM->index_write = 1;     /* INDEX     */
1477   l = t2_encode_packets (tcd_img, tcd_cp, tileno, tile,
1478                          tcd_tcp->numlayers, dest, len, info_IM);
1479 /*---------------CLEAN-------------------*/
1480
1481   time7 = clock () - time7;
1482   printf ("total:     %ld.%.3ld s\n", time7 / CLOCKS_PER_SEC,
1483           (time7 % CLOCKS_PER_SEC) * 1000 / CLOCKS_PER_SEC);
1484
1485   /* cleaning memory */
1486   for (compno = 0; compno < tile->numcomps; compno++)
1487     {
1488       tilec = &tile->comps[compno];
1489       free (tilec->data);
1490     }
1491
1492   return l;
1493 }
1494
1495 int
1496 tcd_encode_tile_pgx (int tileno, unsigned char *dest, int len,
1497                      info_image * info_IM)
1498 {
1499   int compno;
1500   int l, i;
1501   clock_t time;
1502   tcd_tile_t *tile;
1503   j2k_tcp_t *tcp = &tcd_cp->tcps[0];
1504   j2k_tccp_t *tccp = &tcp->tccps[0];
1505
1506   tcd_tileno = tileno;
1507   tcd_tile = tcd_image.tiles;
1508   tcd_tcp = &tcd_cp->tcps[tileno];
1509   tile = tcd_tile;
1510   /* INDEX >> "Precinct_nb_X et Precinct_nb_Y" */
1511
1512   if (info_IM->index_on)
1513     {
1514
1515       tcd_tilecomp_t *tilec_idx = &tile->comps[0];      //Based on Component 0
1516
1517       for (i = 0; i < tilec_idx->numresolutions; i++)
1518         {
1519
1520           tcd_resolution_t *res_idx = &tilec_idx->resolutions[i];
1521
1522
1523
1524           info_IM->tile[tileno].pw[i] = res_idx->pw;
1525
1526           info_IM->tile[tileno].ph[i] = res_idx->ph;
1527
1528
1529
1530           info_IM->tile[tileno].pdx[i] = tccp->prcw[i];
1531
1532           info_IM->tile[tileno].pdy[i] = tccp->prch[i];
1533
1534         }
1535
1536     }
1537
1538   /* << INDEX */
1539 /*---------------TILE-------------------*/
1540   time = clock ();
1541
1542   for (compno = 0; compno < tile->numcomps; compno++)
1543     {
1544       FILE *src;
1545       char tmp[256];
1546       int k;
1547       int elmt;
1548       int i, j;
1549       int tw, w;
1550       tcd_tilecomp_t *tilec = &tile->comps[compno];
1551       int adjust =
1552         tcd_img->comps[compno].sgnd ? 0 : 1 << (tcd_img->comps[compno].
1553                                                 prec - 1);
1554       int offset_x, offset_y;
1555
1556       offset_x = int_ceildiv (tcd_img->x0, tcd_img->comps[compno].dx);
1557       offset_y = int_ceildiv (tcd_img->y0, tcd_img->comps[compno].dy);
1558       tw = tilec->x1 - tilec->x0;
1559       w = int_ceildiv (tcd_img->x1 - tcd_img->x0, tcd_img->comps[compno].dx);
1560       sprintf (tmp, "bandtile%d", tileno / tcd_cp->tw + 1);     /* bandtile file opening */
1561       src = fopen (tmp, "rb");
1562       if (!src)
1563         {
1564           fprintf (stderr, "failed to open %s for reading\n", tmp);
1565           return 1;
1566         }
1567       /* Extract data from bandtile file limited to the current tile */
1568       k = 0;
1569       while (k < tilec->x0 - offset_x)
1570         {
1571           k++;
1572           fscanf (src, "%d", &elmt);
1573         }
1574
1575       for (j = 0; j < tilec->y1 - tilec->y0; j++)
1576         {
1577           for (i = tilec->x0; i < tilec->x1; i++)
1578             {
1579               if (tcd_tcp->tccps[compno].qmfbid == 1)
1580                 {
1581                   fscanf (src, "%d", &elmt);
1582                   tilec->data[i - tilec->x0 + (j) * tw] = elmt - adjust;
1583                   k++;
1584                 }
1585               else if (tcd_tcp->tccps[compno].qmfbid == 0)
1586                 {
1587                   fscanf (src, "%d", &elmt);
1588                   tilec->data[i - tilec->x0 + (j) * tw] =
1589                     (elmt - adjust) << 13;
1590                   k++;
1591                 }
1592             }
1593           while (k < tilec->x0 - offset_x + (j + 1) * w)
1594             {
1595               k++;
1596               fscanf (src, "%d", &elmt);
1597             }
1598         }
1599       fclose (src);
1600     }
1601
1602 /*----------------MCT-------------------*/
1603
1604   if (tcd_tcp->mct)
1605     {
1606       if (tcd_tcp->tccps[0].qmfbid == 0)
1607         {
1608           mct_encode_real (tile->comps[0].data, tile->comps[1].data,
1609                            tile->comps[2].data,
1610                            (tile->comps[0].x1 -
1611                             tile->comps[0].x0) * (tile->comps[0].y1 -
1612                                                   tile->comps[0].y0));
1613         }
1614       else
1615         {
1616           mct_encode (tile->comps[0].data, tile->comps[1].data,
1617                       tile->comps[2].data,
1618                       (tile->comps[0].x1 -
1619                        tile->comps[0].x0) * (tile->comps[0].y1 -
1620                                              tile->comps[0].y0));
1621         }
1622     }
1623
1624 /*----------------DWT---------------------*/
1625
1626   for (compno = 0; compno < tile->numcomps; compno++)
1627     {
1628       tcd_tilecomp_t *tilec = &tile->comps[compno];
1629       if (tcd_tcp->tccps[compno].qmfbid == 1)
1630         {
1631           dwt_encode (tilec->data, tilec->x1 - tilec->x0,
1632                       tilec->y1 - tilec->y0, tilec,
1633                       tilec->numresolutions - 1);
1634         }
1635       else if (tcd_tcp->tccps[compno].qmfbid == 0)
1636         {
1637           dwt_encode_real (tilec->data, tilec->x1 - tilec->x0,
1638                            tilec->y1 - tilec->y0, tilec,
1639                            tilec->numresolutions - 1);
1640         }
1641     }
1642
1643 /*------------------TIER1-----------------*/
1644
1645   t1_init_luts ();
1646   t1_encode_cblks (tile, tcd_tcp);
1647
1648 /*-----------RATE-ALLOCATE------------------*/
1649
1650   info_IM->index_write = 0;     /* INDEX */
1651
1652   if (tcd_cp->disto_alloc || tcd_cp->fixed_quality)     // mod fixed_quality
1653
1654     /* Normal Rate/distortion allocation */
1655
1656     tcd_rateallocate (dest, len, info_IM);
1657
1658   else
1659     /* Fixed layer allocation */
1660
1661     tcd_rateallocate_fixed ();
1662
1663 /*--------------TIER2------------------*/
1664   info_IM->index_write = 1;     /* INDEX */
1665
1666   l = t2_encode_packets (tcd_img, tcd_cp, tileno, tile,
1667                          tcd_tcp->numlayers, dest, len, info_IM);
1668
1669  /*---------------CLEAN-------------------*/
1670   time = clock () - time;
1671   printf ("total:     %ld.%.3ld s\n", time / CLOCKS_PER_SEC,
1672           (time % CLOCKS_PER_SEC) * 1000 / CLOCKS_PER_SEC);
1673
1674   for (compno = 0; compno < tile->numcomps; compno++)
1675     {
1676       tilec = &tile->comps[compno];
1677       free (tilec->data);
1678     }
1679
1680   return l;
1681 }
1682
1683
1684 int
1685 tcd_decode_tile (unsigned char *src, int len, int tileno)
1686 {
1687   int l;
1688   int compno;
1689   int eof = 0;
1690   clock_t time;
1691   tcd_tile_t *tile;
1692
1693   tcd_tileno = tileno;
1694   tcd_tile = &tcd_image.tiles[tileno];
1695   tcd_tcp = &tcd_cp->tcps[tileno];
1696   tile = tcd_tile;
1697
1698   time = clock ();
1699
1700   fprintf (stderr, "tile decoding time %d/%d: ", tileno + 1,
1701            tcd_cp->tw * tcd_cp->th);
1702
1703         /*--------------TIER2------------------*/
1704
1705   l = t2_decode_packets (src, len, tcd_img, tcd_cp, tileno, tile);
1706
1707   if (l == -999)
1708     {
1709       eof = 1;
1710       fprintf (stderr, "tcd_decode: incomplete bistream\n");
1711     }
1712
1713         /*------------------TIER1-----------------*/
1714   t1_init_luts ();
1715   t1_decode_cblks (tile, tcd_tcp);
1716
1717         /*----------------DWT---------------------*/
1718
1719   for (compno = 0; compno < tile->numcomps; compno++)
1720     {
1721       tcd_tilecomp_t *tilec = &tile->comps[compno];
1722       if (tcd_cp->reduce_on == 1)
1723         {
1724           tcd_img->comps[compno].resno_decoded =
1725             tile->comps[compno].numresolutions - tcd_cp->reduce_value - 1;
1726         }
1727
1728
1729       if (tcd_tcp->tccps[compno].qmfbid == 1)
1730         {
1731           dwt_decode (tilec->data, tilec->x1 - tilec->x0,
1732                       tilec->y1 - tilec->y0, tilec,
1733                       tilec->numresolutions - 1,
1734                       tilec->numresolutions - 1 -
1735                       tcd_img->comps[compno].resno_decoded);
1736         }
1737       else
1738         {
1739           dwt_decode_real (tilec->data, tilec->x1 - tilec->x0,
1740                            tilec->y1 - tilec->y0, tilec,
1741                            tilec->numresolutions - 1,
1742                            tilec->numresolutions - 1 -
1743                            tcd_img->comps[compno].resno_decoded);
1744         }
1745
1746       if (tile->comps[compno].numresolutions > 0)
1747         tcd_img->comps[compno].factor =
1748           tile->comps[compno].numresolutions -
1749           (tcd_img->comps[compno].resno_decoded + 1);
1750     }
1751
1752         /*----------------MCT-------------------*/
1753
1754   if (tcd_tcp->mct)
1755     {
1756       if (tcd_tcp->tccps[0].qmfbid == 1)
1757         {
1758           mct_decode (tile->comps[0].data, tile->comps[1].data,
1759                       tile->comps[2].data,
1760                       (tile->comps[0].x1 -
1761                        tile->comps[0].x0) * (tile->comps[0].y1 -
1762                                              tile->comps[0].y0));
1763         }
1764       else
1765         {
1766           mct_decode_real (tile->comps[0].data, tile->comps[1].data,
1767                            tile->comps[2].data,
1768                            (tile->comps[0].x1 -
1769                             tile->comps[0].x0) * (tile->comps[0].y1 -
1770                                                   tile->comps[0].y0));
1771         }
1772     }
1773
1774         /*---------------TILE-------------------*/
1775
1776   for (compno = 0; compno < tile->numcomps; compno++)
1777     {
1778       tcd_tilecomp_t *tilec = &tile->comps[compno];
1779       tcd_resolution_t *res =
1780         &tilec->resolutions[tcd_img->comps[compno].resno_decoded];
1781       int adjust =
1782         tcd_img->comps[compno].sgnd ? 0 : 1 << (tcd_img->comps[compno].
1783                                                 prec - 1);
1784       int min =
1785         tcd_img->comps[compno].
1786         sgnd ? -(1 << (tcd_img->comps[compno].prec - 1)) : 0;
1787       int max =
1788         tcd_img->comps[compno].
1789         sgnd ? (1 << (tcd_img->comps[compno].prec - 1)) -
1790         1 : (1 << tcd_img->comps[compno].prec) - 1;
1791
1792       int tw = tilec->x1 - tilec->x0;
1793       int w = tcd_img->comps[compno].w;
1794
1795       int i, j;
1796       int offset_x = int_ceildivpow2 (tcd_img->comps[compno].x0,
1797                                       tcd_img->comps[compno].factor);
1798       int offset_y = int_ceildivpow2 (tcd_img->comps[compno].y0,
1799                                       tcd_img->comps[compno].factor);
1800
1801       for (j = res->y0; j < res->y1; j++)
1802         {
1803           for (i = res->x0; i < res->x1; i++)
1804             {
1805
1806               int v;
1807
1808               double tmp =
1809                 (double) tilec->data[i - res->x0 + (j - res->y0) * tw];
1810               if (tcd_tcp->tccps[compno].qmfbid == 1)
1811                 {
1812                   v = (int) tmp;
1813                 }
1814               else
1815                 {
1816
1817                   //v = (int) tmp >> 13;
1818
1819                   //Mod antonin : multbug1
1820                   v =
1821                     (int) ((fabs (tmp / 8192.0) >=
1822                             floor (fabs (tmp / 8192.0)) +
1823                             0.5) ? fabs (tmp / 8192.0) +
1824                            1.0 : fabs (tmp / 8192.0));
1825
1826                   v = (tmp < 0) ? -v : v;
1827
1828                   //doM
1829                 }
1830               v += adjust;
1831
1832               tcd_img->comps[compno].data[(i - offset_x) +
1833                                           (j - offset_y) * w] =
1834                 int_clamp (v, min, max);
1835             }
1836         }
1837     }
1838
1839   time = clock () - time;
1840   fprintf (stderr, "total:     %ld.%.3ld s\n", time / CLOCKS_PER_SEC,
1841            (time % CLOCKS_PER_SEC) * 1000 / CLOCKS_PER_SEC);
1842
1843   for (compno = 0; compno < tile->numcomps; compno++)
1844     {
1845       free (tcd_image.tiles[tileno].comps[compno].data);
1846     }
1847
1848   if (eof)
1849     {
1850       longjmp (j2k_error, 1);
1851     }
1852
1853   return l;
1854 }