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=%d, 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 =
 274             (int) floor((1.0 + ss->mant / 2048.0) *
 275                         pow(2.0, numbps - ss->expn) * 8192.0);
 276           band->numbps = ss->expn + tccp->numgbits - 1; /* WHY -1 ? */
 277
 278           band->precincts =
 279             (tcd_precinct_t *) malloc(3 * res->pw * res->ph *
 280                                       sizeof(tcd_precinct_t));
 281
 282           for (i = 0; i < res->pw * res->ph * 3; i++) {
 283             band->precincts[i].imsbtree = NULL;
 284             band->precincts[i].incltree = NULL;
 285           }
 286
 287           for (precno = 0; precno < res->pw * res->ph; precno++) {
 288             int tlcblkxstart, tlcblkystart, brcblkxend, brcblkyend;
 289             int cbgxstart =
 290               tlcbgxstart + (precno % res->pw) * (1 << cbgwidthexpn);
 291             int cbgystart =
 292               tlcbgystart + (precno / res->pw) * (1 << cbgheightexpn);
 293             int cbgxend = cbgxstart + (1 << cbgwidthexpn);
 294             int cbgyend = cbgystart + (1 << cbgheightexpn);
 295             /* tcd_precinct_t *prc=&band->precincts[precno]; */
 296             prc = &band->precincts[precno];
 297             /* precinct size (global) */
 298             prc->x0 = int_max(cbgxstart, band->x0);
 299             prc->y0 = int_max(cbgystart, band->y0);
 300             prc->x1 = int_min(cbgxend, band->x1);
 301             prc->y1 = int_min(cbgyend, band->y1);
 302
 303             tlcblkxstart =
 304               int_floordivpow2(prc->x0, cblkwidthexpn) << cblkwidthexpn;
 305             tlcblkystart =
 306               int_floordivpow2(prc->y0, cblkheightexpn) << cblkheightexpn;
 307             brcblkxend =
 308               int_ceildivpow2(prc->x1, cblkwidthexpn) << cblkwidthexpn;
 309             brcblkyend =
 310               int_ceildivpow2(prc->y1, cblkheightexpn) << cblkheightexpn;
 311             prc->cw = (brcblkxend - tlcblkxstart) >> cblkwidthexpn;
 312             prc->ch = (brcblkyend - tlcblkystart) >> cblkheightexpn;
 313
 314             prc->cblks =
 315               (tcd_cblk_t *) malloc((prc->cw * prc->ch) *
 316                                     sizeof(tcd_cblk_t));
 317             prc->incltree = tgt_create(prc->cw, prc->ch);
 318             prc->imsbtree = tgt_create(prc->cw, prc->ch);
 319
 320             for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++) {
 321               int cblkxstart =
 322                 tlcblkxstart + (cblkno % prc->cw) * (1 << cblkwidthexpn);
 323               int cblkystart =
 324                 tlcblkystart + (cblkno / prc->cw) * (1 << cblkheightexpn);
 325               int cblkxend = cblkxstart + (1 << cblkwidthexpn);
 326               int cblkyend = cblkystart + (1 << cblkheightexpn);
 327
 328               cblk = &prc->cblks[cblkno];
 329               /* code-block size (global) */
 330               cblk->x0 = int_max(cblkxstart, prc->x0);
 331               cblk->y0 = int_max(cblkystart, prc->y0);
 332               cblk->x1 = int_min(cblkxend, prc->x1);
 333               cblk->y1 = int_min(cblkyend, prc->y1);
 334             }
 335           }
 336         }
 337       }
 338     }
 339   }
 340   /* tcd_dump(&tcd_image,curtileno); */
 341 }
 342
 343 void tcd_free_encode(j2k_image_t * img, j2k_cp_t * cp, int curtileno)
 344 {
 345   int tileno, compno, resno, bandno, precno;
 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 =
 518             (int) floor((1.0 + ss->mant / 2048.0) *
 519                         pow(2.0, numbps - ss->expn) * 8192.0);
 520           band->numbps = ss->expn + tccp->numgbits - 1; /* WHY -1 ? */
 521
 522           for (precno = 0; precno < res->pw * res->ph; precno++) {
 523             int tlcblkxstart, tlcblkystart, brcblkxend, brcblkyend;
 524             int cbgxstart =
 525               tlcbgxstart + (precno % res->pw) * (1 << cbgwidthexpn);
 526             int cbgystart =
 527               tlcbgystart + (precno / res->pw) * (1 << cbgheightexpn);
 528             int cbgxend = cbgxstart + (1 << cbgwidthexpn);
 529             int cbgyend = cbgystart + (1 << cbgheightexpn);
 530
 531             prc = &band->precincts[precno];
 532             /* precinct size (global) */
 533             prc->x0 = int_max(cbgxstart, band->x0);
 534             prc->y0 = int_max(cbgystart, band->y0);
 535             prc->x1 = int_min(cbgxend, band->x1);
 536             prc->y1 = int_min(cbgyend, band->y1);
 537
 538             tlcblkxstart =
 539               int_floordivpow2(prc->x0, cblkwidthexpn) << cblkwidthexpn;
 540             tlcblkystart =
 541               int_floordivpow2(prc->y0, cblkheightexpn) << cblkheightexpn;
 542             brcblkxend =
 543               int_ceildivpow2(prc->x1, cblkwidthexpn) << cblkwidthexpn;
 544             brcblkyend =
 545               int_ceildivpow2(prc->y1, cblkheightexpn) << cblkheightexpn;
 546             prc->cw = (brcblkxend - tlcblkxstart) >> cblkwidthexpn;
 547             prc->ch = (brcblkyend - tlcblkystart) >> cblkheightexpn;
 548
 549             free(prc->cblks);
 550             prc->cblks =
 551               (tcd_cblk_t *) malloc(prc->cw * prc->ch *
 552                                     sizeof(tcd_cblk_t));
 553
 554             if (prc->incltree != NULL)
 555               tgt_destroy(prc->incltree);
 556             if (prc->imsbtree != NULL)
 557               tgt_destroy(prc->imsbtree);
 558
 559             prc->incltree = tgt_create(prc->cw, prc->ch);
 560             prc->imsbtree = tgt_create(prc->cw, prc->ch);
 561
 562             for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++) {
 563               int cblkxstart =
 564                 tlcblkxstart + (cblkno % prc->cw) * (1 << cblkwidthexpn);
 565               int cblkystart =
 566                 tlcblkystart + (cblkno / prc->cw) * (1 << cblkheightexpn);
 567               int cblkxend = cblkxstart + (1 << cblkwidthexpn);
 568               int cblkyend = cblkystart + (1 << cblkheightexpn);
 569               cblk = &prc->cblks[cblkno];
 570
 571               /* code-block size (global) */
 572               cblk->x0 = int_max(cblkxstart, prc->x0);
 573               cblk->y0 = int_max(cblkystart, prc->y0);
 574               cblk->x1 = int_min(cblkxend, prc->x1);
 575               cblk->y1 = int_min(cblkyend, prc->y1);
 576
 577             }
 578           }
 579         }
 580       }
 581     }
 582   }
 583   /* tcd_dump(&tcd_image,0); */
 584 }
 585
 586 void tcd_init(j2k_image_t * img, j2k_cp_t * cp)
 587 {
 588   int tileno, compno, resno, bandno, precno, cblkno, i, j, p, q;
 589   unsigned int x0 = 0, y0 = 0, x1 = 0, y1 = 0, w, h;
 590   tcd_img = img;
 591   tcd_cp = cp;
 592   tcd_image.tw = cp->tw;
 593   tcd_image.th = cp->th;
 594   tcd_image.tiles =
 595     (tcd_tile_t *) malloc(cp->tw * cp->th * sizeof(tcd_tile_t));
 596
 597   /*for (tileno = 0; tileno < cp->tw * cp->th; tileno++) {
 598      j2k_tcp_t *tcp = &cp->tcps[tileno];
 599      tcd_tile_t *tile = &tcd_image.tiles[tileno]; */
 600
 601   for (i = 0; i < cp->tileno_size; i++) {
 602     j2k_tcp_t *tcp = &cp->tcps[cp->tileno[i]];
 603     tcd_tile_t *tile = &tcd_image.tiles[cp->tileno[i]];
 604     tileno = cp->tileno[i];
 605
 606
 607     //              int previous_x0, previous_x1, previous_y0, previous_y1;
 608     /* cfr p59 ISO/IEC FDIS15444-1 : 2000 (18 august 2000) */
 609     p = tileno % cp->tw;        /* si numerotation matricielle .. */
 610     q = tileno / cp->tw;        /* .. coordonnees de la tile (q,p) q pour ligne et p pour colonne */
 611
 612     /* 4 borders of the tile rescale on the image if necessary */
 613     tile->x0 = int_max(cp->tx0 + p * cp->tdx, img->x0);
 614     tile->y0 = int_max(cp->ty0 + q * cp->tdy, img->y0);
 615     tile->x1 = int_min(cp->tx0 + (p + 1) * cp->tdx, img->x1);
 616     tile->y1 = int_min(cp->ty0 + (q + 1) * cp->tdy, img->y1);
 617
 618     tile->numcomps = img->numcomps;
 619     tile->comps =
 620       (tcd_tilecomp_t *) malloc(img->numcomps * sizeof(tcd_tilecomp_t));
 621     for (compno = 0; compno < tile->numcomps; compno++) {
 622       j2k_tccp_t *tccp = &tcp->tccps[compno];
 623       tcd_tilecomp_t *tilec = &tile->comps[compno];
 624       /* border of each tile component (global) */
 625       tilec->x0 = int_ceildiv(tile->x0, img->comps[compno].dx);
 626       tilec->y0 = int_ceildiv(tile->y0, img->comps[compno].dy);
 627       tilec->x1 = int_ceildiv(tile->x1, img->comps[compno].dx);
 628       tilec->y1 = int_ceildiv(tile->y1, img->comps[compno].dy);
 629
 630       tilec->data =
 631         (int *) malloc((tilec->x1 - tilec->x0) *
 632                        (tilec->y1 - tilec->y0) * sizeof(int));
 633       tilec->numresolutions = tccp->numresolutions;
 634       tilec->resolutions =
 635         (tcd_resolution_t *) malloc(tilec->numresolutions *
 636                                     sizeof(tcd_resolution_t));
 637       for (resno = 0; resno < tilec->numresolutions; resno++) {
 638         int pdx, pdy;
 639         int levelno = tilec->numresolutions - 1 - resno;
 640         int tlprcxstart, tlprcystart, brprcxend, brprcyend;
 641         int tlcbgxstart, tlcbgystart, brcbgxend, brcbgyend;
 642         int cbgwidthexpn, cbgheightexpn;
 643         int cblkwidthexpn, cblkheightexpn;
 644         tcd_resolution_t *res = &tilec->resolutions[resno];
 645
 646         /* border for each resolution level (global) */
 647         res->x0 = int_ceildivpow2(tilec->x0, levelno);
 648         res->y0 = int_ceildivpow2(tilec->y0, levelno);
 649         res->x1 = int_ceildivpow2(tilec->x1, levelno);
 650         res->y1 = int_ceildivpow2(tilec->y1, levelno);
 651
 652         res->numbands = resno == 0 ? 1 : 3;
 653         /* p. 35, table A-23, ISO/IEC FDIS154444-1 : 2000 (18 august 2000) */
 654         if (tccp->csty & J2K_CCP_CSTY_PRT) {
 655           pdx = tccp->prcw[resno];
 656           pdy = tccp->prch[resno];
 657         } else {
 658           pdx = 15;
 659           pdy = 15;
 660         }
 661         /* p. 64, B.6, ISO/IEC FDIS15444-1 : 2000 (18 august 2000)  */
 662         tlprcxstart = int_floordivpow2(res->x0, pdx) << pdx;
 663         tlprcystart = int_floordivpow2(res->y0, pdy) << pdy;
 664         brprcxend = int_ceildivpow2(res->x1, pdx) << pdx;
 665         brprcyend = int_ceildivpow2(res->y1, pdy) << pdy;
 666         res->pw = (res->x0 == res->x1) ? 0 : ((brprcxend - tlprcxstart) >> pdx);        // Mod Antonin : sizebug1
 667         res->ph = (res->y0 == res->y1) ? 0 : ((brprcyend - tlprcystart) >> pdy);        // Mod Antonin : sizebug1
 668
 669         if (resno == 0) {
 670           tlcbgxstart = tlprcxstart;
 671           tlcbgystart = tlprcystart;
 672           brcbgxend = brprcxend;
 673           brcbgyend = brprcyend;
 674           cbgwidthexpn = pdx;
 675           cbgheightexpn = pdy;
 676         } else {
 677           tlcbgxstart = int_ceildivpow2(tlprcxstart, 1);
 678           tlcbgystart = int_ceildivpow2(tlprcystart, 1);
 679           brcbgxend = int_ceildivpow2(brprcxend, 1);
 680           brcbgyend = int_ceildivpow2(brprcyend, 1);
 681           cbgwidthexpn = pdx - 1;
 682           cbgheightexpn = pdy - 1;
 683         }
 684
 685         cblkwidthexpn = int_min(tccp->cblkw, cbgwidthexpn);
 686         cblkheightexpn = int_min(tccp->cblkh, cbgheightexpn);
 687
 688         for (bandno = 0; bandno < res->numbands; bandno++) {
 689           int x0b, y0b;
 690           int gain, numbps;
 691           j2k_stepsize_t *ss;
 692           tcd_band_t *band = &res->bands[bandno];
 693           band->bandno = resno == 0 ? 0 : bandno + 1;
 694           x0b = (band->bandno == 1) || (band->bandno == 3) ? 1 : 0;
 695           y0b = (band->bandno == 2) || (band->bandno == 3) ? 1 : 0;
 696
 697           if (band->bandno == 0) {
 698             /* band border (global) */
 699             band->x0 = int_ceildivpow2(tilec->x0, levelno);
 700             band->y0 = int_ceildivpow2(tilec->y0, levelno);
 701             band->x1 = int_ceildivpow2(tilec->x1, levelno);
 702             band->y1 = int_ceildivpow2(tilec->y1, levelno);
 703           } else {
 704             /* band border (global) */
 705             band->x0 =
 706               int_ceildivpow2(tilec->x0 -
 707                               (1 << levelno) * x0b, levelno + 1);
 708             band->y0 =
 709               int_ceildivpow2(tilec->y0 -
 710                               (1 << levelno) * y0b, levelno + 1);
 711             band->x1 =
 712               int_ceildivpow2(tilec->x1 -
 713                               (1 << levelno) * x0b, levelno + 1);
 714             band->y1 =
 715               int_ceildivpow2(tilec->y1 -
 716                               (1 << levelno) * y0b, levelno + 1);
 717           }
 718
 719           ss = &tccp->stepsizes[resno ==
 720                                 0 ? 0 : 3 * (resno - 1) + bandno + 1];
 721           gain =
 722             tccp->qmfbid ==
 723             0 ? dwt_getgain_real(band->bandno) : dwt_getgain(band->bandno);
 724           numbps = img->comps[compno].prec + gain;
 725           band->stepsize =
 726             (int) floor((1.0 + ss->mant / 2048.0) *
 727                         pow(2.0, numbps - ss->expn) * 8192.0);
 728           band->numbps = ss->expn + tccp->numgbits - 1; /* WHY -1 ? */
 729
 730           band->precincts =
 731             (tcd_precinct_t *) malloc(res->pw * res->ph *
 732                                       sizeof(tcd_precinct_t));
 733
 734           for (precno = 0; precno < res->pw * res->ph; precno++) {
 735             int tlcblkxstart, tlcblkystart, brcblkxend, brcblkyend;
 736             int cbgxstart =
 737               tlcbgxstart + (precno % res->pw) * (1 << cbgwidthexpn);
 738             int cbgystart =
 739               tlcbgystart + (precno / res->pw) * (1 << cbgheightexpn);
 740             int cbgxend = cbgxstart + (1 << cbgwidthexpn);
 741             int cbgyend = cbgystart + (1 << cbgheightexpn);
 742             tcd_precinct_t *prc = &band->precincts[precno];
 743             /* precinct size (global) */
 744             prc->x0 = int_max(cbgxstart, band->x0);
 745             prc->y0 = int_max(cbgystart, band->y0);
 746             prc->x1 = int_min(cbgxend, band->x1);
 747             prc->y1 = int_min(cbgyend, band->y1);
 748
 749             tlcblkxstart =
 750               int_floordivpow2(prc->x0, cblkwidthexpn) << cblkwidthexpn;
 751             tlcblkystart =
 752               int_floordivpow2(prc->y0, cblkheightexpn) << cblkheightexpn;
 753             brcblkxend =
 754               int_ceildivpow2(prc->x1, cblkwidthexpn) << cblkwidthexpn;
 755             brcblkyend =
 756               int_ceildivpow2(prc->y1, cblkheightexpn) << cblkheightexpn;
 757             prc->cw = (brcblkxend - tlcblkxstart) >> cblkwidthexpn;
 758             prc->ch = (brcblkyend - tlcblkystart) >> cblkheightexpn;
 759
 760             prc->cblks =
 761               (tcd_cblk_t *) malloc(prc->cw * prc->ch *
 762                                     sizeof(tcd_cblk_t));
 763
 764             prc->incltree = tgt_create(prc->cw, prc->ch);
 765             prc->imsbtree = tgt_create(prc->cw, prc->ch);
 766
 767             for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++) {
 768               int cblkxstart =
 769                 tlcblkxstart + (cblkno % prc->cw) * (1 << cblkwidthexpn);
 770               int cblkystart =
 771                 tlcblkystart + (cblkno / prc->cw) * (1 << cblkheightexpn);
 772               int cblkxend = cblkxstart + (1 << cblkwidthexpn);
 773               int cblkyend = cblkystart + (1 << cblkheightexpn);
 774               tcd_cblk_t *cblk = &prc->cblks[cblkno];
 775               /* code-block size (global) */
 776               cblk->x0 = int_max(cblkxstart, prc->x0);
 777               cblk->y0 = int_max(cblkystart, prc->y0);
 778               cblk->x1 = int_min(cblkxend, prc->x1);
 779               cblk->y1 = int_min(cblkyend, prc->y1);
 780
 781
 782
 783               cblk->lastbp = 0; // Add Antonin : quantizbug1
 784             }
 785           }
 786         }
 787       }
 788     }
 789   }
 790   //tcd_dump(&tcd_image,0);
 791
 792
 793   /* Allocate place to store the data decoded = final image */
 794   /* Place limited by the tile really present in the codestream */
 795
 796
 797   for (i = 0; i < img->numcomps; i++) {
 798     for (j = 0; j < cp->tileno_size; j++) {
 799       tileno = cp->tileno[j];
 800       x0 = j == 0 ? tcd_image.tiles[tileno].comps[i].x0 : int_min(x0,
 801                                                                   (unsigned int)
 802                                                                   tcd_image.
 803                                                                   tiles
 804                                                                   [tileno].
 805                                                                   comps
 806                                                                   [i].x0);
 807       y0 =
 808         j == 0 ? tcd_image.tiles[tileno].comps[i].y0 : int_min(y0,
 809                                                                (unsigned int)
 810                                                                tcd_image.
 811                                                                tiles
 812                                                                [tileno].
 813                                                                comps[i].
 814                                                                y0);
 815       x1 =
 816         j == 0 ? tcd_image.tiles[tileno].comps[i].x1 : int_max(x1,
 817                                                                (unsigned int)
 818                                                                tcd_image.
 819                                                                tiles
 820                                                                [tileno].
 821                                                                comps[i].
 822                                                                x1);
 823       y1 =
 824         j == 0 ? tcd_image.tiles[tileno].comps[i].y1 : int_max(y1,
 825                                                                (unsigned int)
 826                                                                tcd_image.
 827                                                                tiles
 828                                                                [tileno].
 829                                                                comps[i].
 830                                                                y1);
 831     }
 832
 833     w = x1 - x0;
 834
 835     h = y1 - y0;
 836     img->comps[i].data = (int *) calloc(w * h, sizeof(int));
 837     img->comps[i].w = w;
 838     img->comps[i].h = h;
 839     img->comps[i].x0 = x0;
 840     img->comps[i].y0 = y0;
 841   }
 842 }
 843
 844 void tcd_makelayer_fixed(int layno, int final)
 845 {
 846   int compno, resno, bandno, precno, cblkno;
 847   int value;                    //, matrice[tcd_tcp->numlayers][tcd_tile->comps[0].numresolutions][3];
 848   int matrice[10][10][3];
 849   int i, j, k;
 850
 851   /*matrice=(int*)malloc(tcd_tcp->numlayers*tcd_tile->comps[0].numresolutions*3*sizeof(int)); */
 852
 853   for (compno = 0; compno < tcd_tile->numcomps; compno++) {
 854     tcd_tilecomp_t *tilec = &tcd_tile->comps[compno];
 855     for (i = 0; i < tcd_tcp->numlayers; i++) {
 856       for (j = 0; j < tilec->numresolutions; j++) {
 857         for (k = 0; k < 3; k++) {
 858           matrice[i][j][k] =
 859             (int) (tcd_cp->
 860                    matrice[i * tilec->numresolutions * 3 +
 861                            j * 3 +
 862                            k] *
 863                    (float) (tcd_img->comps[compno].prec / 16.0));
 864     }}}
 865
 866     for (resno = 0; resno < tilec->numresolutions; resno++) {
 867       tcd_resolution_t *res = &tilec->resolutions[resno];
 868       for (bandno = 0; bandno < res->numbands; bandno++) {
 869         tcd_band_t *band = &res->bands[bandno];
 870         for (precno = 0; precno < res->pw * res->ph; precno++) {
 871           tcd_precinct_t *prc = &band->precincts[precno];
 872           for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++) {
 873             tcd_cblk_t *cblk = &prc->cblks[cblkno];
 874             tcd_layer_t *layer = &cblk->layers[layno];
 875             int n;
 876             int imsb = tcd_img->comps[compno].prec - cblk->numbps;      /* number of bit-plan equal to zero */
 877             /* Correction of the matrix of coefficient to include the IMSB information */
 878
 879             if (layno == 0) {
 880               value = matrice[layno][resno][bandno];
 881               if (imsb >= value)
 882                 value = 0;
 883               else
 884                 value -= imsb;
 885             } else {
 886               value =
 887                 matrice[layno][resno][bandno] -
 888                 matrice[layno - 1][resno][bandno];
 889               if (imsb >= matrice[layno - 1][resno][bandno]) {
 890                 value -= (imsb - matrice[layno - 1][resno][bandno]);
 891                 if (value < 0)
 892                   value = 0;
 893               }
 894             }
 895
 896             if (layno == 0)
 897               cblk->numpassesinlayers = 0;
 898
 899             n = cblk->numpassesinlayers;
 900             if (cblk->numpassesinlayers == 0) {
 901               if (value != 0)
 902                 n = 3 * value - 2 + cblk->numpassesinlayers;
 903               else
 904                 n = cblk->numpassesinlayers;
 905             } else
 906               n = 3 * value + cblk->numpassesinlayers;
 907
 908             layer->numpasses = n - cblk->numpassesinlayers;
 909
 910             if (!layer->numpasses)
 911               continue;
 912
 913             if (cblk->numpassesinlayers == 0) {
 914               layer->len = cblk->passes[n - 1].rate;
 915               layer->data = cblk->data;
 916             } else {
 917               layer->len =
 918                 cblk->passes[n - 1].rate -
 919                 cblk->passes[cblk->numpassesinlayers - 1].rate;
 920               layer->data =
 921                 cblk->data +
 922                 cblk->passes[cblk->numpassesinlayers - 1].rate;
 923             }
 924             if (final)
 925               cblk->numpassesinlayers = n;
 926           }
 927         }
 928       }
 929     }
 930   }
 931 }
 932
 933 void tcd_rateallocate_fixed()
 934 {
 935   int layno;
 936
 937   for (layno = 0; layno < tcd_tcp->numlayers; layno++) {
 938     tcd_makelayer_fixed(layno, 1);
 939   }
 940 }
 941
 942 void tcd_makelayer(int layno, double thresh, int final)
 943 {
 944   int compno, resno, bandno, precno, cblkno, passno;
 945
 946   tcd_tile->distolayer[layno] = 0;      //add fixed_quality
 947
 948   for (compno = 0; compno < tcd_tile->numcomps; compno++) {
 949     tcd_tilecomp_t *tilec = &tcd_tile->comps[compno];
 950     for (resno = 0; resno < tilec->numresolutions; resno++) {
 951       tcd_resolution_t *res = &tilec->resolutions[resno];
 952       for (bandno = 0; bandno < res->numbands; bandno++) {
 953         tcd_band_t *band = &res->bands[bandno];
 954         for (precno = 0; precno < res->pw * res->ph; precno++) {
 955           tcd_precinct_t *prc = &band->precincts[precno];
 956           for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++) {
 957             tcd_cblk_t *cblk = &prc->cblks[cblkno];
 958             tcd_layer_t *layer = &cblk->layers[layno];
 959             int n;
 960
 961             if (layno == 0) {
 962               cblk->numpassesinlayers = 0;
 963             }
 964             n = cblk->numpassesinlayers;
 965             for (passno = cblk->numpassesinlayers;
 966                  passno < cblk->totalpasses; passno++) {
 967               int dr;
 968               double dd;
 969               tcd_pass_t *pass = &cblk->passes[passno];
 970               if (n == 0) {
 971                 dr = pass->rate;
 972                 dd = pass->distortiondec;
 973               } else {
 974                 dr = pass->rate - cblk->passes[n - 1].rate;
 975                 dd = pass->distortiondec - cblk->passes[n -
 976                                                         1].distortiondec;
 977               }
 978               if (!dr) {
 979                 if (dd)
 980                   n = passno + 1;
 981                 continue;
 982               }
 983
 984               if (dd / dr >= thresh)
 985                 n = passno + 1;
 986             }
 987             layer->numpasses = n - cblk->numpassesinlayers;
 988
 989             if (!layer->numpasses) {
 990               layer->disto = 0;
 991               continue;
 992             }
 993
 994             if (cblk->numpassesinlayers == 0) {
 995               layer->len = cblk->passes[n - 1].rate;
 996               layer->data = cblk->data;
 997               layer->disto = cblk->passes[n - 1].distortiondec;
 998             } else {
 999               layer->len = cblk->passes[n - 1].rate -
1000                 cblk->passes[cblk->numpassesinlayers - 1].rate;
1001               layer->data =
1002                 cblk->data +
1003                 cblk->passes[cblk->numpassesinlayers - 1].rate;
1004               layer->disto =
1005                 cblk->passes[n - 1].distortiondec -
1006                 cblk->passes[cblk->numpassesinlayers - 1].distortiondec;
1007             }
1008
1009             tcd_tile->distolayer[layno] += layer->disto;        //add fixed_quality
1010
1011             if (final)
1012               cblk->numpassesinlayers = n;
1013           }
1014         }
1015       }
1016     }
1017   }
1018 }
1019
1020 void tcd_rateallocate(unsigned char *dest, int len, info_image * info_IM)
1021 {
1022   int compno, resno, bandno, precno, cblkno, passno, layno;
1023   double min, max;
1024   double cumdisto[100];         //add fixed_quality
1025   const double K = 1;           // 1.1; //add fixed_quality
1026
1027   double maxSE = 0;
1028   min = DBL_MAX;
1029   max = 0;
1030
1031   tcd_tile->nbpix = 0;          //add fixed_quality
1032
1033   for (compno = 0; compno < tcd_tile->numcomps; compno++) {
1034     tcd_tilecomp_t *tilec = &tcd_tile->comps[compno];
1035
1036     tilec->nbpix = 0;
1037     for (resno = 0; resno < tilec->numresolutions; resno++) {
1038       tcd_resolution_t *res = &tilec->resolutions[resno];
1039       for (bandno = 0; bandno < res->numbands; bandno++) {
1040         tcd_band_t *band = &res->bands[bandno];
1041         for (precno = 0; precno < res->pw * res->ph; precno++) {
1042           tcd_precinct_t *prc = &band->precincts[precno];
1043           for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++) {
1044             tcd_cblk_t *cblk = &prc->cblks[cblkno];
1045             for (passno = 0; passno < cblk->totalpasses; passno++) {
1046               tcd_pass_t *pass = &cblk->passes[passno];
1047               int dr;
1048               double dd, rdslope;
1049               if (passno == 0) {
1050                 dr = pass->rate;
1051                 dd = pass->distortiondec;
1052               } else {
1053                 dr = pass->rate - cblk->passes[passno - 1].rate;
1054                 dd = pass->distortiondec -
1055                   cblk->passes[passno - 1].distortiondec;
1056               }
1057               if (dr == 0) {
1058                 continue;
1059               }
1060
1061               rdslope = dd / dr;
1062
1063               if (rdslope < min) {
1064                 min = rdslope;
1065               }
1066               if (rdslope > max) {
1067                 max = rdslope;
1068               }
1069             }                   /* passno */
1070
1071             tcd_tile->nbpix += ((cblk->x1 - cblk->x0) * (cblk->y1 - cblk->y0)); //add fixed_quality
1072
1073             tilec->nbpix += ((cblk->x1 - cblk->x0) * (cblk->y1 - cblk->y0));    //add fixed_quality
1074
1075           }                     /* cbklno */
1076         }                       /* precno */
1077       }                         /* bandno */
1078     }                           /* resno */
1079     maxSE += (((double)(1 << tcd_img->comps[compno].prec) - 1.0) * ((double)(1 << tcd_img->comps[compno].prec) -1.0)) * ((double)(tilec->nbpix));
1080   }                             /* compno */
1081
1082   /* add antonin index */
1083   if (info_IM->index_on) {
1084     info_tile *info_TL = &info_IM->tile[tcd_tileno];
1085     info_TL->nbpix = tcd_tile->nbpix;
1086     info_TL->distotile = tcd_tile->distotile;
1087     info_TL->thresh =
1088       (double *) malloc(tcd_tcp->numlayers * sizeof(double));
1089   }
1090   /* dda */
1091
1092   for (layno = 0; layno < tcd_tcp->numlayers; layno++) {
1093     volatile double lo = min;
1094     volatile double hi = max;
1095     volatile int success = 0;
1096     volatile int maxlen = tcd_tcp->rates[layno] ? int_min(tcd_tcp->rates[layno], len) : len;    //Mod antonin losslessbug
1097     volatile double goodthresh;
1098     volatile int i;
1099     double distotarget;         //add fixed_quality
1100
1101     distotarget = tcd_tile->distotile - ((K * maxSE) / pow(10, tcd_tcp->distoratio[layno] / 10));       // add fixed_quality
1102
1103     if ((tcd_tcp->rates[layno]) || (tcd_cp->disto_alloc==0)) {
1104       for (i = 0; i < 32; i++) {
1105         volatile double thresh = (lo + hi) / 2;
1106         int l = 0;
1107         double distoachieved = 0;       // add fixed_quality
1108
1109         tcd_makelayer(layno, thresh, 0);
1110
1111         if (tcd_cp->fixed_quality) {    // add fixed_quality
1112           distoachieved =
1113             layno ==
1114             0 ? tcd_tile->distolayer[0] : cumdisto[layno - 1] +
1115             tcd_tile->distolayer[layno];
1116           if (distoachieved < distotarget) {
1117             hi = thresh;
1118             continue;
1119           }
1120           lo = thresh;
1121         } else {
1122           l =
1123             t2_encode_packets(tcd_img, tcd_cp, tcd_tileno, tcd_tile,
1124                               layno + 1, dest, maxlen, info_IM);
1125           /* fprintf(stderr, "rate alloc: len=%d, max=%d\n", l, maxlen); */
1126           if (l == -999) {
1127             lo = thresh;
1128             continue;
1129           }
1130           hi = thresh;
1131         }
1132
1133         success = 1;
1134         goodthresh = thresh;
1135       }
1136     } else {
1137       success = 1;
1138       goodthresh = min;
1139     }
1140
1141     if (!success) {
1142       longjmp(j2k_error, 1);
1143     }
1144
1145     if (info_IM->index_on) {    /* Threshold for Marcela Index */
1146       info_IM->tile[tcd_tileno].thresh[layno] = goodthresh;
1147     }
1148     tcd_makelayer(layno, goodthresh, 1);
1149
1150     cumdisto[layno] = layno == 0 ? tcd_tile->distolayer[0] : cumdisto[layno - 1] + tcd_tile->distolayer[layno]; // add fixed_quality
1151   }
1152 }
1153
1154 int
1155 tcd_encode_tile_pxm(int tileno, unsigned char *dest, int len,
1156                     info_image * info_IM)
1157 {
1158   int compno;
1159   int l, i, npck=0;
1160   clock_t time7;
1161   tcd_tile_t *tile;
1162   j2k_tcp_t *tcp = &tcd_cp->tcps[0];
1163   j2k_tccp_t *tccp = &tcp->tccps[0];
1164
1165   tcd_tileno = tileno;
1166   tcd_tile = tcd_image.tiles;
1167   tcd_tcp = &tcd_cp->tcps[tileno];
1168   tile = tcd_tile;
1169   /* INDEX >> "Precinct_nb_X et Precinct_nb_Y" */
1170   if (info_IM->index_on) {
1171     tcd_tilecomp_t *tilec_idx = &tile->comps[0];        //Based on Component 0
1172
1173     for (i = 0; i < tilec_idx->numresolutions; i++) {
1174
1175       tcd_resolution_t *res_idx = &tilec_idx->resolutions[i];
1176
1177       info_IM->tile[tileno].pw[i] = res_idx->pw;
1178       info_IM->tile[tileno].ph[i] = res_idx->ph;
1179
1180       npck+=res_idx->pw * res_idx->ph;
1181
1182       info_IM->tile[tileno].pdx[i] = tccp->prcw[i];
1183       info_IM->tile[tileno].pdy[i] = tccp->prch[i];
1184
1185     }
1186     info_IM->tile[tileno].packet = (info_packet *) calloc(info_IM->Comp * info_IM->Layer * npck, sizeof(info_packet));
1187   }
1188   /* << INDEX */
1189
1190 /*---------------TILE-------------------*/
1191
1192   time7 = clock();
1193
1194   for (compno = 0; compno < tile->numcomps; compno++) {
1195     FILE *src;
1196     char tmp[256];
1197     int k;
1198     unsigned char elmt;
1199     int i, j;
1200     int tw, w;
1201     tcd_tilecomp_t *tilec = &tile->comps[compno];
1202     int adjust =
1203       tcd_img->comps[compno].sgnd ? 0 : 1 << (tcd_img->comps[compno].
1204                                               prec - 1);
1205     int offset_x, offset_y;
1206
1207     offset_x = int_ceildiv(tcd_img->x0, tcd_img->comps[compno].dx);
1208     offset_y = int_ceildiv(tcd_img->y0, tcd_img->comps[compno].dy);
1209     tw = tilec->x1 - tilec->x0;
1210     w = int_ceildiv(tcd_img->x1 - tcd_img->x0, tcd_img->comps[compno].dx);
1211     sprintf(tmp, "Compo%d", compno);    /* component file */
1212     src = fopen(tmp, "rb");
1213     if (!src) {
1214       fprintf(stderr, "failed to open %s for reading\n", tmp);
1215       return 1;
1216     }
1217
1218     /* read the Compo file to extract data of the tile */
1219     k = 0;
1220     fseek(src, (tilec->x0 - offset_x) + (tilec->y0 - offset_y) * w,
1221           SEEK_SET);
1222     k = (tilec->x0 - offset_x) + (tilec->y0 - offset_y) * w;
1223     for (j = tilec->y0; j < tilec->y1; j++) {
1224       for (i = tilec->x0; i < tilec->x1; i++) {
1225         if (tcd_tcp->tccps[compno].qmfbid == 1) {
1226           elmt = fgetc(src);
1227           tilec->data[i - tilec->x0 + (j - tilec->y0) * tw] =
1228             elmt - adjust;
1229           k++;
1230         } else if (tcd_tcp->tccps[compno].qmfbid == 0) {
1231           elmt = fgetc(src);
1232           tilec->data[i - tilec->x0 + (j - tilec->y0) * tw] =
1233             (elmt - adjust) << 13;
1234           k++;
1235         }
1236       }
1237       fseek(src, (tilec->x0 - offset_x) + (j + 1 - offset_y) * w - k,
1238             SEEK_CUR);
1239       k = tilec->x0 - offset_x + (j + 1 - offset_y) * w;
1240
1241     }
1242     fclose(src);
1243   }
1244
1245 /*----------------MCT-------------------*/
1246
1247   if (tcd_tcp->mct) {
1248     if (tcd_tcp->tccps[0].qmfbid == 0) {
1249       mct_encode_real(tile->comps[0].data, tile->comps[1].data,
1250                       tile->comps[2].data,
1251                       (tile->comps[0].x1 -
1252                        tile->comps[0].x0) * (tile->comps[0].y1 -
1253                                              tile->comps[0].y0));
1254     } else {
1255       mct_encode(tile->comps[0].data, tile->comps[1].data,
1256                  tile->comps[2].data,
1257                  (tile->comps[0].x1 -
1258                   tile->comps[0].x0) * (tile->comps[0].y1 -
1259                                         tile->comps[0].y0));
1260     }
1261   }
1262 /*----------------DWT---------------------*/
1263
1264   /* time3=clock(); */
1265   for (compno = 0; compno < tile->numcomps; compno++) {
1266     tcd_tilecomp_t *tilec = &tile->comps[compno];
1267     if (tcd_tcp->tccps[compno].qmfbid == 1) {
1268       dwt_encode(tilec->data, tilec->x1 - tilec->x0,
1269                  tilec->y1 - tilec->y0, tilec, tilec->numresolutions - 1);
1270     } else if (tcd_tcp->tccps[compno].qmfbid == 0) {
1271       dwt_encode_real(tilec->data, tilec->x1 - tilec->x0,
1272                       tilec->y1 - tilec->y0, tilec,
1273                       tilec->numresolutions - 1);
1274     }
1275   }
1276 /*------------------TIER1-----------------*/
1277
1278   t1_init_luts();
1279   t1_encode_cblks(tile, tcd_tcp);
1280
1281 /*-----------RATE-ALLOCATE------------------*/
1282   info_IM->index_write = 0;     /* INDEX     */
1283
1284   if (tcd_cp->disto_alloc || tcd_cp->fixed_quality)     // mod fixed_quality
1285     /* Normal Rate/distortion allocation */
1286     tcd_rateallocate(dest, len, info_IM);
1287   else
1288     /* Fixed layer allocation */
1289     tcd_rateallocate_fixed();
1290
1291 /*--------------TIER2------------------*/
1292   info_IM->index_write = 1;     /* INDEX     */
1293   l = t2_encode_packets(tcd_img, tcd_cp, tileno, tile,
1294                         tcd_tcp->numlayers, dest, len, info_IM);
1295 /*---------------CLEAN-------------------*/
1296
1297   time7 = clock() - time7;
1298   fprintf(stdout,"total:     %ld.%.3ld s\n", time7 / CLOCKS_PER_SEC,
1299          (time7 % CLOCKS_PER_SEC) * 1000 / CLOCKS_PER_SEC);
1300
1301   /* cleaning memory */
1302   for (compno = 0; compno < tile->numcomps; compno++) {
1303     tilec = &tile->comps[compno];
1304     free(tilec->data);
1305   }
1306
1307   return l;
1308 }
1309
1310 int
1311 tcd_encode_tile_pgx(int tileno, unsigned char *dest, int len,
1312                     info_image * info_IM)
1313 {
1314   int compno;
1315   int l, i;
1316   clock_t time;
1317   tcd_tile_t *tile;
1318   j2k_tcp_t *tcp = &tcd_cp->tcps[0];
1319   j2k_tccp_t *tccp = &tcp->tccps[0];
1320
1321   tcd_tileno = tileno;
1322   tcd_tile = tcd_image.tiles;
1323   tcd_tcp = &tcd_cp->tcps[tileno];
1324   tile = tcd_tile;
1325   /* INDEX >> "Precinct_nb_X et Precinct_nb_Y" */
1326
1327   if (info_IM->index_on) {
1328
1329     tcd_tilecomp_t *tilec_idx = &tile->comps[0];        //Based on Component 0
1330
1331     for (i = 0; i < tilec_idx->numresolutions; i++) {
1332
1333       tcd_resolution_t *res_idx = &tilec_idx->resolutions[i];
1334
1335
1336
1337       info_IM->tile[tileno].pw[i] = res_idx->pw;
1338
1339       info_IM->tile[tileno].ph[i] = res_idx->ph;
1340
1341
1342
1343       info_IM->tile[tileno].pdx[i] = tccp->prcw[i];
1344
1345       info_IM->tile[tileno].pdy[i] = tccp->prch[i];
1346
1347     }
1348
1349   }
1350
1351   /* << INDEX */
1352 /*---------------TILE-------------------*/
1353   time = clock();
1354
1355   for (compno = 0; compno < tile->numcomps; compno++) {
1356     FILE *src;
1357     char tmp[256];
1358     int k;
1359     int elmt;
1360     int i, j;
1361     int tw, w;
1362     tcd_tilecomp_t *tilec = &tile->comps[compno];
1363     int adjust =
1364       tcd_img->comps[compno].sgnd ? 0 : 1 << (tcd_img->comps[compno].
1365                                               prec - 1);
1366     int offset_x, offset_y;
1367
1368     offset_x = int_ceildiv(tcd_img->x0, tcd_img->comps[compno].dx);
1369     offset_y = int_ceildiv(tcd_img->y0, tcd_img->comps[compno].dy);
1370     tw = tilec->x1 - tilec->x0;
1371     w = int_ceildiv(tcd_img->x1 - tcd_img->x0, tcd_img->comps[compno].dx);
1372     sprintf(tmp, "bandtile%d", tileno / tcd_cp->tw + 1);        /* bandtile file opening */
1373     src = fopen(tmp, "rb");
1374     if (!src) {
1375       fprintf(stderr, "failed to open %s for reading\n", tmp);
1376       return 1;
1377     }
1378     /* Extract data from bandtile file limited to the current tile */
1379     k = 0;
1380     while (k < tilec->x0 - offset_x) {
1381       k++;
1382       fscanf(src, "%d", &elmt);
1383     }
1384
1385     for (j = 0; j < tilec->y1 - tilec->y0; j++) {
1386       for (i = tilec->x0; i < tilec->x1; i++) {
1387         if (tcd_tcp->tccps[compno].qmfbid == 1) {
1388           fscanf(src, "%d", &elmt);
1389           tilec->data[i - tilec->x0 + (j) * tw] = elmt - adjust;
1390           k++;
1391         } else if (tcd_tcp->tccps[compno].qmfbid == 0) {
1392           fscanf(src, "%d", &elmt);
1393           tilec->data[i - tilec->x0 + (j) * tw] = (elmt - adjust) << 13;
1394           k++;
1395         }
1396       }
1397       while (k < tilec->x0 - offset_x + (j + 1) * w) {
1398         k++;
1399         fscanf(src, "%d", &elmt);
1400       }
1401     }
1402     fclose(src);
1403   }
1404
1405 /*----------------MCT-------------------*/
1406
1407   if (tcd_tcp->mct) {
1408     if (tcd_tcp->tccps[0].qmfbid == 0) {
1409       mct_encode_real(tile->comps[0].data, tile->comps[1].data,
1410                       tile->comps[2].data,
1411                       (tile->comps[0].x1 -
1412                        tile->comps[0].x0) * (tile->comps[0].y1 -
1413                                              tile->comps[0].y0));
1414     } else {
1415       mct_encode(tile->comps[0].data, tile->comps[1].data,
1416                  tile->comps[2].data,
1417                  (tile->comps[0].x1 -
1418                   tile->comps[0].x0) * (tile->comps[0].y1 -
1419                                         tile->comps[0].y0));
1420     }
1421   }
1422
1423 /*----------------DWT---------------------*/
1424
1425   for (compno = 0; compno < tile->numcomps; compno++) {
1426     tcd_tilecomp_t *tilec = &tile->comps[compno];
1427     if (tcd_tcp->tccps[compno].qmfbid == 1) {
1428       dwt_encode(tilec->data, tilec->x1 - tilec->x0,
1429                  tilec->y1 - tilec->y0, tilec, tilec->numresolutions - 1);
1430     } else if (tcd_tcp->tccps[compno].qmfbid == 0) {
1431       dwt_encode_real(tilec->data, tilec->x1 - tilec->x0,
1432                       tilec->y1 - tilec->y0, tilec,
1433                       tilec->numresolutions - 1);
1434     }
1435   }
1436
1437 /*------------------TIER1-----------------*/
1438
1439   t1_init_luts();
1440   t1_encode_cblks(tile, tcd_tcp);
1441
1442 /*-----------RATE-ALLOCATE------------------*/
1443
1444   info_IM->index_write = 0;     /* INDEX */
1445
1446   if (tcd_cp->disto_alloc || tcd_cp->fixed_quality)     // mod fixed_quality
1447
1448     /* Normal Rate/distortion allocation */
1449
1450     tcd_rateallocate(dest, len, info_IM);
1451
1452   else
1453     /* Fixed layer allocation */
1454
1455     tcd_rateallocate_fixed();
1456
1457 /*--------------TIER2------------------*/
1458   info_IM->index_write = 1;     /* INDEX */
1459
1460   l = t2_encode_packets(tcd_img, tcd_cp, tileno, tile,
1461                         tcd_tcp->numlayers, dest, len, info_IM);
1462
1463  /*---------------CLEAN-------------------*/
1464   time = clock() - time;
1465   fprintf(stdout,"total:     %ld.%.3ld s\n", time / CLOCKS_PER_SEC,
1466          (time % CLOCKS_PER_SEC) * 1000 / CLOCKS_PER_SEC);
1467
1468   for (compno = 0; compno < tile->numcomps; compno++) {
1469     tilec = &tile->comps[compno];
1470     free(tilec->data);
1471   }
1472
1473   return l;
1474 }
1475
1476
1477 int tcd_decode_tile(unsigned char *src, int len, int tileno)
1478 {
1479   int l;
1480   int compno;
1481   int eof = 0;
1482   clock_t time;
1483   tcd_tile_t *tile;
1484
1485   tcd_tileno = tileno;
1486   tcd_tile = &tcd_image.tiles[tileno];
1487   tcd_tcp = &tcd_cp->tcps[tileno];
1488   tile = tcd_tile;
1489
1490   time = clock();
1491
1492   fprintf(stdout, "Tile %d of %d decoded in ", tileno + 1,
1493           tcd_cp->tw * tcd_cp->th);
1494
1495         /*--------------TIER2------------------*/
1496
1497   l = t2_decode_packets(src, len, tcd_img, tcd_cp, tileno, tile);
1498
1499   if (l == -999) {
1500     eof = 1;
1501     fprintf(stderr, "tcd_decode: incomplete bistream\n");
1502   }
1503
1504         /*------------------TIER1-----------------*/
1505   t1_init_luts();
1506   t1_decode_cblks(tile, tcd_tcp);
1507
1508         /*----------------DWT---------------------*/
1509
1510   for (compno = 0; compno < tile->numcomps; compno++) {
1511     tcd_tilecomp_t *tilec = &tile->comps[compno];
1512     if (tcd_cp->reduce != 0) {
1513       tcd_img->comps[compno].resno_decoded =
1514         tile->comps[compno].numresolutions - tcd_cp->reduce - 1;
1515     }
1516
1517
1518     if (tcd_tcp->tccps[compno].qmfbid == 1) {
1519       dwt_decode(tilec->data, tilec->x1 - tilec->x0,
1520                  tilec->y1 - tilec->y0, tilec,
1521                  tilec->numresolutions - 1,
1522                  tilec->numresolutions - 1 -
1523                  tcd_img->comps[compno].resno_decoded);
1524     } else {
1525       dwt_decode_real(tilec->data, tilec->x1 - tilec->x0,
1526                       tilec->y1 - tilec->y0, tilec,
1527                       tilec->numresolutions - 1,
1528                       tilec->numresolutions - 1 -
1529                       tcd_img->comps[compno].resno_decoded);
1530     }
1531
1532     if (tile->comps[compno].numresolutions > 0)
1533       tcd_img->comps[compno].factor =
1534         tile->comps[compno].numresolutions -
1535         (tcd_img->comps[compno].resno_decoded + 1);
1536   }
1537
1538         /*----------------MCT-------------------*/
1539
1540   if (tcd_tcp->mct) {
1541     if (tcd_tcp->tccps[0].qmfbid == 1) {
1542       mct_decode(tile->comps[0].data, tile->comps[1].data,
1543                  tile->comps[2].data,
1544                  (tile->comps[0].x1 -
1545                   tile->comps[0].x0) * (tile->comps[0].y1 -
1546                                         tile->comps[0].y0));
1547     } else {
1548       mct_decode_real(tile->comps[0].data, tile->comps[1].data,
1549                       tile->comps[2].data,
1550                       (tile->comps[0].x1 -
1551                        tile->comps[0].x0) * (tile->comps[0].y1 -
1552                                              tile->comps[0].y0));
1553     }
1554   }
1555
1556         /*---------------TILE-------------------*/
1557
1558   for (compno = 0; compno < tile->numcomps; compno++) {
1559     tcd_tilecomp_t *tilec = &tile->comps[compno];
1560     tcd_resolution_t *res =
1561       &tilec->resolutions[tcd_img->comps[compno].resno_decoded];
1562     int adjust =
1563       tcd_img->comps[compno].sgnd ? 0 : 1 << (tcd_img->comps[compno].
1564                                               prec - 1);
1565     int min =
1566       tcd_img->comps[compno].
1567       sgnd ? -(1 << (tcd_img->comps[compno].prec - 1)) : 0;
1568     int max =
1569       tcd_img->comps[compno].
1570       sgnd ? (1 << (tcd_img->comps[compno].prec - 1)) -
1571       1 : (1 << tcd_img->comps[compno].prec) - 1;
1572
1573     int tw = tilec->x1 - tilec->x0;
1574     int w = tcd_img->comps[compno].w;
1575
1576     int i, j;
1577     int offset_x = int_ceildivpow2(tcd_img->comps[compno].x0,
1578                                    tcd_img->comps[compno].factor);
1579     int offset_y = int_ceildivpow2(tcd_img->comps[compno].y0,
1580                                    tcd_img->comps[compno].factor);
1581
1582     for (j = res->y0; j < res->y1; j++) {
1583       for (i = res->x0; i < res->x1; i++) {
1584
1585         int v;
1586
1587         double tmp =
1588           (double) tilec->data[i - res->x0 + (j - res->y0) * tw];
1589         if (tcd_tcp->tccps[compno].qmfbid == 1) {
1590           v = (int) tmp;
1591         } else {
1592
1593           //v = (int) tmp >> 13;
1594
1595           //Mod antonin : multbug1
1596           v =
1597             (int) ((fabs(tmp / 8192.0) >=
1598                     floor(fabs(tmp / 8192.0)) +
1599                     0.5) ? fabs(tmp / 8192.0) + 1.0 : fabs(tmp / 8192.0));
1600
1601           v = (tmp < 0) ? -v : v;
1602
1603           //doM
1604         }
1605         v += adjust;
1606
1607         tcd_img->comps[compno].data[(i - offset_x) +
1608                                     (j - offset_y) * w] =
1609           int_clamp(v, min, max);
1610       }
1611     }
1612   }
1613
1614   time = clock() - time;
1615   fprintf(stdout, "%ld.%.3ld s\n", time / CLOCKS_PER_SEC,
1616           (time % CLOCKS_PER_SEC) * 1000 / CLOCKS_PER_SEC);
1617
1618
1619
1620   for (compno = 0; compno < tile->numcomps; compno++) {
1621     free(tcd_image.tiles[tileno].comps[compno].data);
1622   }
1623
1624   if (eof) {
1625     longjmp(j2k_error, 1);
1626   }
1627
1628   return l;
1629 }
1630
1631
1632
1633 void tcd_dec_release()
1634
1635 {
1636
1637   int tileno,compno,resno,bandno,precno;
1638
1639   for (tileno=0;tileno<tcd_image.tw*tcd_image.th;tileno++) {
1640
1641     tcd_tile_t tile=tcd_image.tiles[tileno];
1642
1643     for (compno=0;compno<tile.numcomps;compno++) {
1644
1645       tcd_tilecomp_t tilec=tile.comps[compno];
1646
1647       for (resno=0;resno<tilec.numresolutions;resno++) {
1648
1649         tcd_resolution_t res=tilec.resolutions[resno];
1650
1651         for (bandno=0;bandno<res.numbands;bandno++) {
1652
1653           tcd_band_t band=res.bands[bandno];
1654
1655           for (precno=0;precno<res.ph*res.pw;precno++) {
1656
1657             tcd_precinct_t prec=band.precincts[precno];
1658
1659             if (prec.cblks!=NULL) free(prec.cblks);
1660
1661             if (prec.imsbtree!=NULL) free(prec.imsbtree);
1662
1663             if (prec.incltree!=NULL) free(prec.incltree);
1664
1665           }
1666
1667           if (band.precincts!=NULL) free(band.precincts);
1668
1669         }
1670
1671       }
1672
1673       if (tilec.resolutions!=NULL) free(tilec.resolutions);
1674
1675     }
1676
1677     if (tile.comps!=NULL) free(tile.comps);
1678
1679   }
1680
1681   if (tcd_image.tiles!=NULL) free(tcd_image.tiles);
1682
1683 }