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