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     for (resno = 0; resno < tilec->numresolutions; resno++) {
 859       tcd_resolution_t *res = &tilec->resolutions[resno];
 860       for (bandno = 0; bandno < res->numbands; bandno++) {
 861    tcd_band_t *band = &res->bands[bandno];
 862    for (precno = 0; precno < res->pw * res->ph; precno++) {
 863      tcd_precinct_t *prc = &band->precincts[precno];
 864      for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++) {
 865        tcd_cblk_t *cblk = &prc->cblks[cblkno];
 866        tcd_layer_t *layer = &cblk->layers[layno];
 867        int n;
 868        int imsb = tcd_img->comps[compno].prec - cblk->numbps;   /* number of bit-plan equal to zero */
 869        /* Correction of the matrix of coefficient to include the IMSB information */
 870
 871        if (layno == 0) {
 872          value = matrice[layno][resno][bandno];
 873          if (imsb >= value)
 874       value = 0;
 875          else
 876       value -= imsb;
 877        } else {
 878          value =
 879       matrice[layno][resno][bandno] -
 880       matrice[layno - 1][resno][bandno];
 881          if (imsb >= matrice[layno - 1][resno][bandno]) {
 882       value -= (imsb - matrice[layno - 1][resno][bandno]);
 883       if (value < 0)
 884         value = 0;
 885          }
 886        }
 887
 888        if (layno == 0)
 889          cblk->numpassesinlayers = 0;
 890
 891        n = cblk->numpassesinlayers;
 892        if (cblk->numpassesinlayers == 0) {
 893          if (value != 0)
 894       n = 3 * value - 2 + cblk->numpassesinlayers;
 895          else
 896       n = cblk->numpassesinlayers;
 897        } else
 898          n = 3 * value + cblk->numpassesinlayers;
 899
 900        layer->numpasses = n - cblk->numpassesinlayers;
 901
 902        if (!layer->numpasses)
 903          continue;
 904
 905        if (cblk->numpassesinlayers == 0) {
 906          layer->len = cblk->passes[n - 1].rate;
 907          layer->data = cblk->data;
 908        } else {
 909          layer->len =
 910       cblk->passes[n - 1].rate -
 911       cblk->passes[cblk->numpassesinlayers - 1].rate;
 912          layer->data =
 913       cblk->data +
 914       cblk->passes[cblk->numpassesinlayers - 1].rate;
 915        }
 916        if (final)
 917          cblk->numpassesinlayers = n;
 918      }
 919    }
 920       }
 921     }
 922   }
 923 }
 924
 925 void tcd_rateallocate_fixed()
 926 {
 927   int layno;
 928
 929   for (layno = 0; layno < tcd_tcp->numlayers; layno++) {
 930     tcd_makelayer_fixed(layno, 1);
 931   }
 932 }
 933
 934 void tcd_makelayer(int layno, double thresh, int final)
 935 {
 936   int compno, resno, bandno, precno, cblkno, passno;
 937
 938   tcd_tile->distolayer[layno] = 0;   /*add fixed_quality*/
 939
 940   for (compno = 0; compno < tcd_tile->numcomps; compno++) {
 941     tcd_tilecomp_t *tilec = &tcd_tile->comps[compno];
 942     for (resno = 0; resno < tilec->numresolutions; resno++) {
 943       tcd_resolution_t *res = &tilec->resolutions[resno];
 944       for (bandno = 0; bandno < res->numbands; bandno++) {
 945    tcd_band_t *band = &res->bands[bandno];
 946    for (precno = 0; precno < res->pw * res->ph; precno++) {
 947      tcd_precinct_t *prc = &band->precincts[precno];
 948      for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++) {
 949        tcd_cblk_t *cblk = &prc->cblks[cblkno];
 950        tcd_layer_t *layer = &cblk->layers[layno];
 951        int n;
 952
 953        if (layno == 0) {
 954          cblk->numpassesinlayers = 0;
 955        }
 956        n = cblk->numpassesinlayers;
 957        for (passno = cblk->numpassesinlayers;
 958        passno < cblk->totalpasses; passno++) {
 959          int dr;
 960          double dd;
 961          tcd_pass_t *pass = &cblk->passes[passno];
 962          if (n == 0) {
 963       dr = pass->rate;
 964       dd = pass->distortiondec;
 965          } else {
 966       dr = pass->rate - cblk->passes[n - 1].rate;
 967       dd = pass->distortiondec - cblk->passes[n -
 968                      1].distortiondec;
 969          }
 970          if (!dr) {
 971       if (dd)
 972         n = passno + 1;
 973       continue;
 974          }
 975
 976          if (dd / dr >= thresh)
 977       n = passno + 1;
 978        }
 979        layer->numpasses = n - cblk->numpassesinlayers;
 980
 981        if (!layer->numpasses) {
 982          layer->disto = 0;
 983          continue;
 984        }
 985
 986        if (cblk->numpassesinlayers == 0) {
 987          layer->len = cblk->passes[n - 1].rate;
 988          layer->data = cblk->data;
 989          layer->disto = cblk->passes[n - 1].distortiondec;
 990        } else {
 991          layer->len = cblk->passes[n - 1].rate -
 992       cblk->passes[cblk->numpassesinlayers - 1].rate;
 993          layer->data =
 994       cblk->data +
 995       cblk->passes[cblk->numpassesinlayers - 1].rate;
 996          layer->disto =
 997       cblk->passes[n - 1].distortiondec -
 998       cblk->passes[cblk->numpassesinlayers - 1].distortiondec;
 999        }
1000
1001        tcd_tile->distolayer[layno] += layer->disto;   /*add fixed_quality*/
1002
1003        if (final)
1004          cblk->numpassesinlayers = n;
1005      }
1006    }
1007       }
1008     }
1009   }
1010 }
1011
1012 void tcd_rateallocate(unsigned char *dest, int len, info_image * info_IM)
1013 {
1014   int compno, resno, bandno, precno, cblkno, passno, layno;
1015   double min, max;
1016   double cumdisto[100];      /*add fixed_quality*/
1017   const double K = 1;      /* 1.1; //add fixed_quality*/
1018
1019   double maxSE = 0;
1020   min = DBL_MAX;
1021   max = 0;
1022
1023   tcd_tile->nbpix = 0;      /*add fixed_quality*/
1024
1025   for (compno = 0; compno < tcd_tile->numcomps; compno++) {
1026     tcd_tilecomp_t *tilec = &tcd_tile->comps[compno];
1027
1028     tilec->nbpix = 0;
1029     for (resno = 0; resno < tilec->numresolutions; resno++) {
1030       tcd_resolution_t *res = &tilec->resolutions[resno];
1031       for (bandno = 0; bandno < res->numbands; bandno++) {
1032    tcd_band_t *band = &res->bands[bandno];
1033    for (precno = 0; precno < res->pw * res->ph; precno++) {
1034      tcd_precinct_t *prc = &band->precincts[precno];
1035      for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++) {
1036        tcd_cblk_t *cblk = &prc->cblks[cblkno];
1037        for (passno = 0; passno < cblk->totalpasses; passno++) {
1038          tcd_pass_t *pass = &cblk->passes[passno];
1039          int dr;
1040          double dd, rdslope;
1041          if (passno == 0) {
1042       dr = pass->rate;
1043       dd = pass->distortiondec;
1044          } else {
1045       dr = pass->rate - cblk->passes[passno - 1].rate;
1046       dd = pass->distortiondec -
1047         cblk->passes[passno - 1].distortiondec;
1048          }
1049          if (dr == 0) {
1050       continue;
1051          }
1052
1053          rdslope = dd / dr;
1054
1055          if (rdslope < min) {
1056       min = rdslope;
1057          }
1058          if (rdslope > max) {
1059       max = rdslope;
1060          }
1061        }         /* passno */
1062
1063        tcd_tile->nbpix += ((cblk->x1 - cblk->x0) * (cblk->y1 - cblk->y0));   /*add fixed_quality*/
1064
1065        tilec->nbpix += ((cblk->x1 - cblk->x0) * (cblk->y1 - cblk->y0));   /*add fixed_quality*/
1066
1067      }         /* cbklno */
1068    }         /* precno */
1069       }            /* bandno */
1070     }            /* resno */
1071     maxSE += (((double)(1 << tcd_img->comps[compno].prec) - 1.0) * ((double)(1 << tcd_img->comps[compno].prec) -1.0)) * ((double)(tilec->nbpix));
1072   }            /* compno */
1073
1074   /* add antonin index */
1075   if (info_IM->index_on) {
1076     info_tile *info_TL = &info_IM->tile[tcd_tileno];
1077     info_TL->nbpix = tcd_tile->nbpix;
1078     info_TL->distotile = tcd_tile->distotile;
1079     info_TL->thresh =
1080       (double *) malloc(tcd_tcp->numlayers * sizeof(double));
1081   }
1082   /* dda */
1083
1084   for (layno = 0; layno < tcd_tcp->numlayers; layno++) {
1085     volatile double lo = min;
1086     volatile double hi = max;
1087     volatile int success = 0;
1088     volatile int maxlen = tcd_tcp->rates[layno] ? int_min(tcd_tcp->rates[layno], len) : len;   /*Mod antonin losslessbug*/
1089     volatile double goodthresh;
1090     volatile int i;
1091     double distotarget;      /*add fixed_quality*/
1092
1093     distotarget = tcd_tile->distotile - ((K * maxSE) / pow(10, tcd_tcp->distoratio[layno] / 10));   /* add fixed_quality*/
1094
1095     if ((tcd_tcp->rates[layno]) || (tcd_cp->disto_alloc==0)) {
1096       for (i = 0; i < 32; i++) {
1097    volatile double thresh = (lo + hi) / 2;
1098    int l = 0;
1099    double distoachieved = 0;   /* add fixed_quality*/
1100
1101    tcd_makelayer(layno, thresh, 0);
1102
1103    if (tcd_cp->fixed_quality) {   /* add fixed_quality*/
1104      distoachieved =
1105        layno ==
1106        0 ? tcd_tile->distolayer[0] : cumdisto[layno - 1] +
1107        tcd_tile->distolayer[layno];
1108      if (distoachieved < distotarget) {
1109        hi = thresh;
1110        continue;
1111      }
1112      lo = thresh;
1113    } else {
1114      l =
1115        t2_encode_packets(tcd_img, tcd_cp, tcd_tileno, tcd_tile,
1116                layno + 1, dest, maxlen, info_IM);
1117      /* fprintf(stderr, "rate alloc: len=%d, max=%d\n", l, maxlen); */
1118      if (l == -999) {
1119        lo = thresh;
1120        continue;
1121      }
1122      hi = thresh;
1123    }
1124
1125    success = 1;
1126    goodthresh = thresh;
1127       }
1128     } else {
1129       success = 1;
1130       goodthresh = min;
1131     }
1132
1133     if (!success) {
1134       longjmp(j2k_error, 1);
1135     }
1136
1137     if (info_IM->index_on) {   /* Threshold for Marcela Index */
1138       info_IM->tile[tcd_tileno].thresh[layno] = goodthresh;
1139     }
1140     tcd_makelayer(layno, goodthresh, 1);
1141
1142     cumdisto[layno] = layno == 0 ? tcd_tile->distolayer[0] : cumdisto[layno - 1] + tcd_tile->distolayer[layno];   /* add fixed_quality*/
1143   }
1144 }
1145
1146 int
1147 tcd_encode_tile_pxm(int tileno, unsigned char *dest, int len,
1148           info_image * info_IM)
1149 {
1150   int compno;
1151   int l, i, npck=0;
1152   clock_t time7;
1153   tcd_tile_t *tile;
1154   j2k_tcp_t *tcp = &tcd_cp->tcps[0];
1155   j2k_tccp_t *tccp = &tcp->tccps[0];
1156
1157   tcd_tileno = tileno;
1158   tcd_tile = tcd_image.tiles;
1159   tcd_tcp = &tcd_cp->tcps[tileno];
1160   tile = tcd_tile;
1161   /* INDEX >> "Precinct_nb_X et Precinct_nb_Y" */
1162   if (info_IM->index_on) {
1163     tcd_tilecomp_t *tilec_idx = &tile->comps[0];   /*Based on Component 0*/
1164
1165     for (i = 0; i < tilec_idx->numresolutions; i++) {
1166
1167       tcd_resolution_t *res_idx = &tilec_idx->resolutions[i];
1168
1169       info_IM->tile[tileno].pw[i] = res_idx->pw;
1170       info_IM->tile[tileno].ph[i] = res_idx->ph;
1171
1172       npck+=res_idx->pw * res_idx->ph;
1173
1174       info_IM->tile[tileno].pdx[i] = tccp->prcw[i];
1175       info_IM->tile[tileno].pdy[i] = tccp->prch[i];
1176
1177     }
1178     info_IM->tile[tileno].packet = (info_packet *) calloc(info_IM->Comp * info_IM->Layer * npck, sizeof(info_packet));
1179   }
1180   /* << INDEX */
1181
1182 /*---------------TILE-------------------*/
1183
1184   time7 = clock();
1185
1186   for (compno = 0; compno < tile->numcomps; compno++) {
1187     FILE *src;
1188     char tmp[256];
1189     int k;
1190     unsigned char elmt;
1191     int i, j;
1192     int tw, w;
1193     tcd_tilecomp_t *tilec = &tile->comps[compno];
1194     int adjust =
1195       tcd_img->comps[compno].sgnd ? 0 : 1 << (tcd_img->comps[compno].
1196                      prec - 1);
1197     int offset_x, offset_y;
1198
1199     offset_x = int_ceildiv(tcd_img->x0, tcd_img->comps[compno].dx);
1200     offset_y = int_ceildiv(tcd_img->y0, tcd_img->comps[compno].dy);
1201     tw = tilec->x1 - tilec->x0;
1202     w = int_ceildiv(tcd_img->x1 - tcd_img->x0, tcd_img->comps[compno].dx);
1203     sprintf(tmp, "Compo%d", compno);   /* component file */
1204     src = fopen(tmp, "rb");
1205     if (!src) {
1206       fprintf(stderr, "failed to open %s for reading\n", tmp);
1207       return 1;
1208     }
1209
1210     /* read the Compo file to extract data of the tile */
1211     k = 0;
1212     fseek(src, (tilec->x0 - offset_x) + (tilec->y0 - offset_y) * w,
1213      SEEK_SET);
1214     k = (tilec->x0 - offset_x) + (tilec->y0 - offset_y) * w;
1215     for (j = tilec->y0; j < tilec->y1; j++) {
1216       for (i = tilec->x0; i < tilec->x1; i++) {
1217    if (tcd_tcp->tccps[compno].qmfbid == 1) {
1218      elmt = fgetc(src);
1219      tilec->data[i - tilec->x0 + (j - tilec->y0) * tw] =
1220        elmt - adjust;
1221      k++;
1222    } else if (tcd_tcp->tccps[compno].qmfbid == 0) {
1223      elmt = fgetc(src);
1224      tilec->data[i - tilec->x0 + (j - tilec->y0) * tw] =
1225        (elmt - adjust) << 13;
1226      k++;
1227    }
1228       }
1229       fseek(src, (tilec->x0 - offset_x) + (j + 1 - offset_y) * w - k,
1230        SEEK_CUR);
1231       k = tilec->x0 - offset_x + (j + 1 - offset_y) * w;
1232
1233     }
1234     fclose(src);
1235   }
1236
1237 /*----------------MCT-------------------*/
1238
1239   if (tcd_tcp->mct) {
1240     if (tcd_tcp->tccps[0].qmfbid == 0) {
1241       mct_encode_real(tile->comps[0].data, tile->comps[1].data,
1242             tile->comps[2].data,
1243             (tile->comps[0].x1 -
1244              tile->comps[0].x0) * (tile->comps[0].y1 -
1245                     tile->comps[0].y0));
1246     } else {
1247       mct_encode(tile->comps[0].data, tile->comps[1].data,
1248        tile->comps[2].data,
1249        (tile->comps[0].x1 -
1250         tile->comps[0].x0) * (tile->comps[0].y1 -
1251                tile->comps[0].y0));
1252     }
1253   }
1254 /*----------------DWT---------------------*/
1255
1256 /* mod Ive*/
1257 for (compno = 0; compno < tile->numcomps; compno++) {
1258   tcd_tilecomp_t *tilec = &tile->comps[compno];
1259   if (tcd_tcp->tccps[compno].qmfbid == 1) {
1260     dwt_encode(tilec);
1261   } else if (tcd_tcp->tccps[compno].qmfbid == 0) {
1262     dwt_encode_real(tilec);
1263   }
1264 }
1265 /* /mod Ive*/
1266 /*------------------TIER1-----------------*/
1267
1268   t1_init_luts();
1269   t1_encode_cblks(tile, tcd_tcp);
1270
1271 /*-----------RATE-ALLOCATE------------------*/
1272   info_IM->index_write = 0;   /* INDEX     */
1273
1274   if (tcd_cp->disto_alloc || tcd_cp->fixed_quality)   /* mod fixed_quality*/
1275     /* Normal Rate/distortion allocation */
1276     tcd_rateallocate(dest, len, info_IM);
1277   else
1278     /* Fixed layer allocation */
1279     tcd_rateallocate_fixed();
1280
1281 /*--------------TIER2------------------*/
1282   info_IM->index_write = 1;   /* INDEX     */
1283   l = t2_encode_packets(tcd_img, tcd_cp, tileno, tile,
1284          tcd_tcp->numlayers, dest, len, info_IM);
1285 /*---------------CLEAN-------------------*/
1286
1287   time7 = clock() - time7;
1288   fprintf(stdout,"total:     %ld.%.3ld s\n", time7 / CLOCKS_PER_SEC,
1289     (time7 % (int)CLOCKS_PER_SEC) * 1000 / CLOCKS_PER_SEC);
1290
1291   /* cleaning memory */
1292   for (compno = 0; compno < tile->numcomps; compno++) {
1293     tilec = &tile->comps[compno];
1294     free(tilec->data);
1295   }
1296
1297   return l;
1298 }
1299
1300 int
1301 tcd_encode_tile_pgx(int tileno, unsigned char *dest, int len,
1302           info_image * info_IM)
1303 {
1304   int compno;
1305   int l, i, npck=0;
1306   clock_t time;
1307   tcd_tile_t *tile;
1308   j2k_tcp_t *tcp = &tcd_cp->tcps[0];
1309   j2k_tccp_t *tccp = &tcp->tccps[0];
1310
1311   tcd_tileno = tileno;
1312   tcd_tile = tcd_image.tiles;
1313   tcd_tcp = &tcd_cp->tcps[tileno];
1314   tile = tcd_tile;
1315   /* INDEX >> "Precinct_nb_X et Precinct_nb_Y" */
1316   if (info_IM->index_on) {
1317     tcd_tilecomp_t *tilec_idx = &tile->comps[0];   /*Based on Component 0*/
1318
1319     for (i = 0; i < tilec_idx->numresolutions; i++) {
1320
1321       tcd_resolution_t *res_idx = &tilec_idx->resolutions[i];
1322
1323       info_IM->tile[tileno].pw[i] = res_idx->pw;
1324       info_IM->tile[tileno].ph[i] = res_idx->ph;
1325
1326       npck+=res_idx->pw * res_idx->ph;
1327
1328       info_IM->tile[tileno].pdx[i] = tccp->prcw[i];
1329       info_IM->tile[tileno].pdy[i] = tccp->prch[i];
1330
1331     }
1332     info_IM->tile[tileno].packet = (info_packet *) calloc(info_IM->Comp * info_IM->Layer * npck, sizeof(info_packet));
1333   }
1334   /* << INDEX */
1335 /*---------------TILE-------------------*/
1336   time = clock();
1337
1338   for (compno = 0; compno < tile->numcomps; compno++) {
1339     FILE *src;
1340     char tmp[256];
1341     int k;
1342     int elmt;
1343     int i, j;
1344     int tw, w;
1345     tcd_tilecomp_t *tilec = &tile->comps[compno];
1346     int adjust =
1347       tcd_img->comps[compno].sgnd ? 0 : 1 << (tcd_img->comps[compno].
1348                      prec - 1);
1349     int offset_x, offset_y;
1350
1351     offset_x = int_ceildiv(tcd_img->x0, tcd_img->comps[compno].dx);
1352     offset_y = int_ceildiv(tcd_img->y0, tcd_img->comps[compno].dy);
1353     tw = tilec->x1 - tilec->x0;
1354     w = int_ceildiv(tcd_img->x1 - tcd_img->x0, tcd_img->comps[compno].dx);
1355     sprintf(tmp, "bandtile%d", tileno / tcd_cp->tw + 1);   /* bandtile file opening */
1356     src = fopen(tmp, "rb");
1357     if (!src) {
1358       fprintf(stderr, "failed to open %s for reading\n", tmp);
1359       return 1;
1360     }
1361     /* Extract data from bandtile file limited to the current tile */
1362     k = 0;
1363     while (k < tilec->x0 - offset_x) {
1364       k++;
1365       fscanf(src, "%d", &elmt);
1366     }
1367
1368     for (j = 0; j < tilec->y1 - tilec->y0; j++) {
1369       for (i = tilec->x0; i < tilec->x1; i++) {
1370    if (tcd_tcp->tccps[compno].qmfbid == 1) {
1371      fscanf(src, "%d", &elmt);
1372      tilec->data[i - tilec->x0 + (j) * tw] = elmt - adjust;
1373      k++;
1374    } else if (tcd_tcp->tccps[compno].qmfbid == 0) {
1375      fscanf(src, "%d", &elmt);
1376      tilec->data[i - tilec->x0 + (j) * tw] = (elmt - adjust) << 13;
1377      k++;
1378    }
1379       }
1380       while (k < tilec->x0 - offset_x + (j + 1) * w) {
1381    k++;
1382    fscanf(src, "%d", &elmt);
1383       }
1384     }
1385     fclose(src);
1386   }
1387
1388 /*----------------MCT-------------------*/
1389
1390   if (tcd_tcp->mct) {
1391     if (tcd_tcp->tccps[0].qmfbid == 0) {
1392       mct_encode_real(tile->comps[0].data, tile->comps[1].data,
1393             tile->comps[2].data,
1394             (tile->comps[0].x1 -
1395              tile->comps[0].x0) * (tile->comps[0].y1 -
1396                     tile->comps[0].y0));
1397     } else {
1398       mct_encode(tile->comps[0].data, tile->comps[1].data,
1399        tile->comps[2].data,
1400        (tile->comps[0].x1 -
1401         tile->comps[0].x0) * (tile->comps[0].y1 -
1402                tile->comps[0].y0));
1403     }
1404   }
1405
1406 /*----------------DWT---------------------*/
1407
1408 /* mod Ive*/
1409 for (compno = 0; compno < tile->numcomps; compno++) {
1410   tcd_tilecomp_t *tilec = &tile->comps[compno];
1411   if (tcd_tcp->tccps[compno].qmfbid == 1) {
1412     dwt_encode(tilec);
1413   } else if (tcd_tcp->tccps[compno].qmfbid == 0) {
1414     dwt_encode_real(tilec);
1415   }
1416 }
1417 /* /mod Ive*/
1418
1419 /*------------------TIER1-----------------*/
1420
1421   t1_init_luts();
1422   t1_encode_cblks(tile, tcd_tcp);
1423
1424 /*-----------RATE-ALLOCATE------------------*/
1425
1426   info_IM->index_write = 0;   /* INDEX */
1427
1428   if (tcd_cp->disto_alloc || tcd_cp->fixed_quality)   /* mod fixed_quality*/
1429
1430     /* Normal Rate/distortion allocation */
1431
1432     tcd_rateallocate(dest, len, info_IM);
1433
1434   else
1435     /* Fixed layer allocation */
1436
1437     tcd_rateallocate_fixed();
1438
1439 /*--------------TIER2------------------*/
1440   info_IM->index_write = 1;   /* INDEX */
1441
1442   l = t2_encode_packets(tcd_img, tcd_cp, tileno, tile,
1443          tcd_tcp->numlayers, dest, len, info_IM);
1444
1445  /*---------------CLEAN-------------------*/
1446   time = clock() - time;
1447   fprintf(stdout,"total:     %ld.%.3ld s\n", time / CLOCKS_PER_SEC,
1448     (time % (int)CLOCKS_PER_SEC) * 1000 / CLOCKS_PER_SEC);
1449
1450   for (compno = 0; compno < tile->numcomps; compno++) {
1451     tilec = &tile->comps[compno];
1452     free(tilec->data);
1453   }
1454
1455   return l;
1456 }
1457
1458
1459 int tcd_decode_tile(unsigned char *src, int len, int tileno)
1460 {
1461   int l;
1462   int compno;
1463   int eof = 0;
1464   clock_t time;
1465   tcd_tile_t *tile;
1466
1467   tcd_tileno = tileno;
1468   tcd_tile = &tcd_image.tiles[tileno];
1469   tcd_tcp = &tcd_cp->tcps[tileno];
1470   tile = tcd_tile;
1471
1472   time = clock();
1473
1474   fprintf(stdout, "Tile %d of %d decoded in ", tileno + 1,
1475      tcd_cp->tw * tcd_cp->th);
1476
1477    /*--------------TIER2------------------*/
1478
1479   l = t2_decode_packets(src, len, tcd_img, tcd_cp, tileno, tile);
1480
1481   if (l == -999) {
1482     eof = 1;
1483     fprintf(stderr, "tcd_decode: incomplete bistream\n");
1484   }
1485
1486    /*------------------TIER1-----------------*/
1487   t1_init_luts();
1488   t1_decode_cblks(tile, tcd_tcp);
1489
1490    /*----------------DWT---------------------*/
1491
1492   for (compno = 0; compno < tile->numcomps; compno++) {
1493     tcd_tilecomp_t *tilec = &tile->comps[compno];
1494     if (tcd_cp->reduce != 0) {
1495       tcd_img->comps[compno].resno_decoded =
1496    tile->comps[compno].numresolutions - tcd_cp->reduce - 1;
1497     }
1498
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
1590
1591   for (compno = 0; compno < tile->numcomps; compno++) {
1592     free(tcd_image.tiles[tileno].comps[compno].data);
1593   }
1594
1595   if (eof) {
1596     longjmp(j2k_error, 1);
1597   }
1598
1599   return l;
1600 }
1601
1602
1603
1604 void tcd_dec_release()
1605
1606 {
1607
1608   int tileno,compno,resno,bandno,precno;
1609
1610   for (tileno=0;tileno<tcd_image.tw*tcd_image.th;tileno++) {
1611
1612     tcd_tile_t tile=tcd_image.tiles[tileno];
1613
1614     for (compno=0;compno<tile.numcomps;compno++) {
1615
1616       tcd_tilecomp_t tilec=tile.comps[compno];
1617
1618       for (resno=0;resno<tilec.numresolutions;resno++) {
1619
1620    tcd_resolution_t res=tilec.resolutions[resno];
1621
1622    for (bandno=0;bandno<res.numbands;bandno++) {
1623
1624      tcd_band_t band=res.bands[bandno];
1625
1626      for (precno=0;precno<res.ph*res.pw;precno++) {
1627
1628        tcd_precinct_t prec=band.precincts[precno];
1629
1630        if (prec.cblks!=NULL) free(prec.cblks);
1631
1632        if (prec.imsbtree!=NULL) free(prec.imsbtree);
1633
1634        if (prec.incltree!=NULL) free(prec.incltree);
1635
1636      }
1637
1638      if (band.precincts!=NULL) free(band.precincts);
1639
1640    }
1641
1642       }
1643
1644       if (tilec.resolutions!=NULL) free(tilec.resolutions);
1645
1646     }
1647
1648     if (tile.comps!=NULL) free(tile.comps);
1649
1650   }
1651
1652   if (tcd_image.tiles!=NULL) free(tcd_image.tiles);
1653
1654 }