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