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
1080     maxSE +=
1081       (double) (((1 << tcd_img->comps[compno].prec) -
1082                  1) * ((1 << tcd_img->comps[compno].prec) -
1083                        1)) * (tilec->nbpix);
1084   }                             /* compno */
1085
1086   /* add antonin index */
1087   if (info_IM->index_on) {
1088     info_tile *info_TL = &info_IM->tile[tcd_tileno];
1089     info_TL->nbpix = tcd_tile->nbpix;
1090     info_TL->distotile = tcd_tile->distotile;
1091     info_TL->thresh =
1092       (double *) malloc(tcd_tcp->numlayers * sizeof(double));
1093   }
1094   /* dda */
1095
1096   for (layno = 0; layno < tcd_tcp->numlayers; layno++) {
1097     volatile double lo = min;
1098     volatile double hi = max;
1099     volatile int success = 0;
1100     volatile int maxlen = tcd_tcp->rates[layno] ? int_min(tcd_tcp->rates[layno], len) : len;    //Mod antonin losslessbug
1101     volatile double goodthresh;
1102     volatile int i;
1103     double distotarget;         //add fixed_quality
1104
1105     distotarget = tcd_tile->distotile - ((K * maxSE) / pow(10, tcd_tcp->distoratio[layno] / 10));       // add fixed_quality
1106
1107     if (tcd_tcp->rates[layno]) {
1108       for (i = 0; i < 32; i++) {
1109         volatile double thresh = (lo + hi) / 2;
1110         int l = 0;
1111         double distoachieved = 0;       // add fixed_quality
1112
1113         tcd_makelayer(layno, thresh, 0);
1114
1115         if (tcd_cp->fixed_quality) {    // add fixed_quality
1116           distoachieved =
1117             layno ==
1118             0 ? tcd_tile->distolayer[0] : cumdisto[layno - 1] +
1119             tcd_tile->distolayer[layno];
1120           if (distoachieved < distotarget) {
1121             hi = thresh;
1122             continue;
1123           }
1124           lo = thresh;
1125         } else {
1126           l =
1127             t2_encode_packets(tcd_img, tcd_cp, tcd_tileno, tcd_tile,
1128                               layno + 1, dest, maxlen, info_IM);
1129           /* fprintf(stderr, "rate alloc: len=%d, max=%d\n", l, maxlen); */
1130           if (l == -999) {
1131             lo = thresh;
1132             continue;
1133           }
1134           hi = thresh;
1135         }
1136
1137         success = 1;
1138         goodthresh = thresh;
1139       }
1140     } else {
1141       success = 1;
1142       goodthresh = min;
1143     }
1144
1145     if (!success) {
1146       longjmp(j2k_error, 1);
1147     }
1148
1149     if (info_IM->index_on) {    /* Threshold for Marcela Index */
1150       info_IM->tile[tcd_tileno].thresh[layno] = goodthresh;
1151     }
1152     tcd_makelayer(layno, goodthresh, 1);
1153
1154     cumdisto[layno] = layno == 0 ? tcd_tile->distolayer[0] : cumdisto[layno - 1] + tcd_tile->distolayer[layno]; // add fixed_quality
1155   }
1156 }
1157
1158 int
1159 tcd_encode_tile_pxm(int tileno, unsigned char *dest, int len,
1160                     info_image * info_IM)
1161 {
1162   int compno;
1163   int l, i;
1164   clock_t time7;
1165   tcd_tile_t *tile;
1166   j2k_tcp_t *tcp = &tcd_cp->tcps[0];
1167   j2k_tccp_t *tccp = &tcp->tccps[0];
1168
1169   tcd_tileno = tileno;
1170   tcd_tile = tcd_image.tiles;
1171   tcd_tcp = &tcd_cp->tcps[tileno];
1172   tile = tcd_tile;
1173   /* INDEX >> "Precinct_nb_X et Precinct_nb_Y" */
1174   if (info_IM->index_on) {
1175     tcd_tilecomp_t *tilec_idx = &tile->comps[0];        //Based on Component 0
1176
1177     for (i = 0; i < tilec_idx->numresolutions; i++) {
1178
1179       tcd_resolution_t *res_idx = &tilec_idx->resolutions[i];
1180
1181       info_IM->tile[tileno].pw[i] = res_idx->pw;
1182       info_IM->tile[tileno].ph[i] = res_idx->ph;
1183
1184       info_IM->tile[tileno].pdx[i] = tccp->prcw[i];
1185       info_IM->tile[tileno].pdy[i] = tccp->prch[i];
1186
1187     }
1188   }
1189   /* << INDEX */
1190
1191 /*---------------TILE-------------------*/
1192
1193   time7 = clock();
1194
1195   for (compno = 0; compno < tile->numcomps; compno++) {
1196     FILE *src;
1197     char tmp[256];
1198     int k;
1199     unsigned char elmt;
1200     int i, j;
1201     int tw, w;
1202     tcd_tilecomp_t *tilec = &tile->comps[compno];
1203     int adjust =
1204       tcd_img->comps[compno].sgnd ? 0 : 1 << (tcd_img->comps[compno].
1205                                               prec - 1);
1206     int offset_x, offset_y;
1207
1208     offset_x = int_ceildiv(tcd_img->x0, tcd_img->comps[compno].dx);
1209     offset_y = int_ceildiv(tcd_img->y0, tcd_img->comps[compno].dy);
1210     tw = tilec->x1 - tilec->x0;
1211     w = int_ceildiv(tcd_img->x1 - tcd_img->x0, tcd_img->comps[compno].dx);
1212     sprintf(tmp, "Compo%d", compno);    /* component file */
1213     src = fopen(tmp, "rb");
1214     if (!src) {
1215       fprintf(stderr, "failed to open %s for reading\n", tmp);
1216       return 1;
1217     }
1218
1219     /* read the Compo file to extract data of the tile */
1220     k = 0;
1221     fseek(src, (tilec->x0 - offset_x) + (tilec->y0 - offset_y) * w,
1222           SEEK_SET);
1223     k = (tilec->x0 - offset_x) + (tilec->y0 - offset_y) * w;
1224     for (j = tilec->y0; j < tilec->y1; j++) {
1225       for (i = tilec->x0; i < tilec->x1; i++) {
1226         if (tcd_tcp->tccps[compno].qmfbid == 1) {
1227           elmt = fgetc(src);
1228           tilec->data[i - tilec->x0 + (j - tilec->y0) * tw] =
1229             elmt - adjust;
1230           k++;
1231         } else if (tcd_tcp->tccps[compno].qmfbid == 0) {
1232           elmt = fgetc(src);
1233           tilec->data[i - tilec->x0 + (j - tilec->y0) * tw] =
1234             (elmt - adjust) << 13;
1235           k++;
1236         }
1237       }
1238       fseek(src, (tilec->x0 - offset_x) + (j + 1 - offset_y) * w - k,
1239             SEEK_CUR);
1240       k = tilec->x0 - offset_x + (j + 1 - offset_y) * w;
1241
1242     }
1243     fclose(src);
1244   }
1245
1246 /*----------------MCT-------------------*/
1247
1248   if (tcd_tcp->mct) {
1249     if (tcd_tcp->tccps[0].qmfbid == 0) {
1250       mct_encode_real(tile->comps[0].data, tile->comps[1].data,
1251                       tile->comps[2].data,
1252                       (tile->comps[0].x1 -
1253                        tile->comps[0].x0) * (tile->comps[0].y1 -
1254                                              tile->comps[0].y0));
1255     } else {
1256       mct_encode(tile->comps[0].data, tile->comps[1].data,
1257                  tile->comps[2].data,
1258                  (tile->comps[0].x1 -
1259                   tile->comps[0].x0) * (tile->comps[0].y1 -
1260                                         tile->comps[0].y0));
1261     }
1262   }
1263 /*----------------DWT---------------------*/
1264
1265   /* time3=clock(); */
1266   for (compno = 0; compno < tile->numcomps; compno++) {
1267     tcd_tilecomp_t *tilec = &tile->comps[compno];
1268     if (tcd_tcp->tccps[compno].qmfbid == 1) {
1269       dwt_encode(tilec->data, tilec->x1 - tilec->x0,
1270                  tilec->y1 - tilec->y0, tilec, tilec->numresolutions - 1);
1271     } else if (tcd_tcp->tccps[compno].qmfbid == 0) {
1272       dwt_encode_real(tilec->data, tilec->x1 - tilec->x0,
1273                       tilec->y1 - tilec->y0, tilec,
1274                       tilec->numresolutions - 1);
1275     }
1276   }
1277 /*------------------TIER1-----------------*/
1278
1279   t1_init_luts();
1280   t1_encode_cblks(tile, tcd_tcp);
1281
1282 /*-----------RATE-ALLOCATE------------------*/
1283   info_IM->index_write = 0;     /* INDEX     */
1284
1285   if (tcd_cp->disto_alloc || tcd_cp->fixed_quality)     // mod fixed_quality
1286     /* Normal Rate/distortion allocation */
1287     tcd_rateallocate(dest, len, info_IM);
1288   else
1289     /* Fixed layer allocation */
1290     tcd_rateallocate_fixed();
1291
1292 /*--------------TIER2------------------*/
1293   info_IM->index_write = 1;     /* INDEX     */
1294   l = t2_encode_packets(tcd_img, tcd_cp, tileno, tile,
1295                         tcd_tcp->numlayers, dest, len, info_IM);
1296 /*---------------CLEAN-------------------*/
1297
1298   time7 = clock() - time7;
1299   fprintf(stdout,"total:     %ld.%.3ld s\n", time7 / CLOCKS_PER_SEC,
1300          (time7 % CLOCKS_PER_SEC) * 1000 / CLOCKS_PER_SEC);
1301
1302   /* cleaning memory */
1303   for (compno = 0; compno < tile->numcomps; compno++) {
1304     tilec = &tile->comps[compno];
1305     free(tilec->data);
1306   }
1307
1308   return l;
1309 }
1310
1311 int
1312 tcd_encode_tile_pgx(int tileno, unsigned char *dest, int len,
1313                     info_image * info_IM)
1314 {
1315   int compno;
1316   int l, i;
1317   clock_t time;
1318   tcd_tile_t *tile;
1319   j2k_tcp_t *tcp = &tcd_cp->tcps[0];
1320   j2k_tccp_t *tccp = &tcp->tccps[0];
1321
1322   tcd_tileno = tileno;
1323   tcd_tile = tcd_image.tiles;
1324   tcd_tcp = &tcd_cp->tcps[tileno];
1325   tile = tcd_tile;
1326   /* INDEX >> "Precinct_nb_X et Precinct_nb_Y" */
1327
1328   if (info_IM->index_on) {
1329
1330     tcd_tilecomp_t *tilec_idx = &tile->comps[0];        //Based on Component 0
1331
1332     for (i = 0; i < tilec_idx->numresolutions; i++) {
1333
1334       tcd_resolution_t *res_idx = &tilec_idx->resolutions[i];
1335
1336
1337
1338       info_IM->tile[tileno].pw[i] = res_idx->pw;
1339
1340       info_IM->tile[tileno].ph[i] = res_idx->ph;
1341
1342
1343
1344       info_IM->tile[tileno].pdx[i] = tccp->prcw[i];
1345
1346       info_IM->tile[tileno].pdy[i] = tccp->prch[i];
1347
1348     }
1349
1350   }
1351
1352   /* << INDEX */
1353 /*---------------TILE-------------------*/
1354   time = clock();
1355
1356   for (compno = 0; compno < tile->numcomps; compno++) {
1357     FILE *src;
1358     char tmp[256];
1359     int k;
1360     int elmt;
1361     int i, j;
1362     int tw, w;
1363     tcd_tilecomp_t *tilec = &tile->comps[compno];
1364     int adjust =
1365       tcd_img->comps[compno].sgnd ? 0 : 1 << (tcd_img->comps[compno].
1366                                               prec - 1);
1367     int offset_x, offset_y;
1368
1369     offset_x = int_ceildiv(tcd_img->x0, tcd_img->comps[compno].dx);
1370     offset_y = int_ceildiv(tcd_img->y0, tcd_img->comps[compno].dy);
1371     tw = tilec->x1 - tilec->x0;
1372     w = int_ceildiv(tcd_img->x1 - tcd_img->x0, tcd_img->comps[compno].dx);
1373     sprintf(tmp, "bandtile%d", tileno / tcd_cp->tw + 1);        /* bandtile file opening */
1374     src = fopen(tmp, "rb");
1375     if (!src) {
1376       fprintf(stderr, "failed to open %s for reading\n", tmp);
1377       return 1;
1378     }
1379     /* Extract data from bandtile file limited to the current tile */
1380     k = 0;
1381     while (k < tilec->x0 - offset_x) {
1382       k++;
1383       fscanf(src, "%d", &elmt);
1384     }
1385
1386     for (j = 0; j < tilec->y1 - tilec->y0; j++) {
1387       for (i = tilec->x0; i < tilec->x1; i++) {
1388         if (tcd_tcp->tccps[compno].qmfbid == 1) {
1389           fscanf(src, "%d", &elmt);
1390           tilec->data[i - tilec->x0 + (j) * tw] = elmt - adjust;
1391           k++;
1392         } else if (tcd_tcp->tccps[compno].qmfbid == 0) {
1393           fscanf(src, "%d", &elmt);
1394           tilec->data[i - tilec->x0 + (j) * tw] = (elmt - adjust) << 13;
1395           k++;
1396         }
1397       }
1398       while (k < tilec->x0 - offset_x + (j + 1) * w) {
1399         k++;
1400         fscanf(src, "%d", &elmt);
1401       }
1402     }
1403     fclose(src);
1404   }
1405
1406 /*----------------MCT-------------------*/
1407
1408   if (tcd_tcp->mct) {
1409     if (tcd_tcp->tccps[0].qmfbid == 0) {
1410       mct_encode_real(tile->comps[0].data, tile->comps[1].data,
1411                       tile->comps[2].data,
1412                       (tile->comps[0].x1 -
1413                        tile->comps[0].x0) * (tile->comps[0].y1 -
1414                                              tile->comps[0].y0));
1415     } else {
1416       mct_encode(tile->comps[0].data, tile->comps[1].data,
1417                  tile->comps[2].data,
1418                  (tile->comps[0].x1 -
1419                   tile->comps[0].x0) * (tile->comps[0].y1 -
1420                                         tile->comps[0].y0));
1421     }
1422   }
1423
1424 /*----------------DWT---------------------*/
1425
1426   for (compno = 0; compno < tile->numcomps; compno++) {
1427     tcd_tilecomp_t *tilec = &tile->comps[compno];
1428     if (tcd_tcp->tccps[compno].qmfbid == 1) {
1429       dwt_encode(tilec->data, tilec->x1 - tilec->x0,
1430                  tilec->y1 - tilec->y0, tilec, tilec->numresolutions - 1);
1431     } else if (tcd_tcp->tccps[compno].qmfbid == 0) {
1432       dwt_encode_real(tilec->data, tilec->x1 - tilec->x0,
1433                       tilec->y1 - tilec->y0, tilec,
1434                       tilec->numresolutions - 1);
1435     }
1436   }
1437
1438 /*------------------TIER1-----------------*/
1439
1440   t1_init_luts();
1441   t1_encode_cblks(tile, tcd_tcp);
1442
1443 /*-----------RATE-ALLOCATE------------------*/
1444
1445   info_IM->index_write = 0;     /* INDEX */
1446
1447   if (tcd_cp->disto_alloc || tcd_cp->fixed_quality)     // mod fixed_quality
1448
1449     /* Normal Rate/distortion allocation */
1450
1451     tcd_rateallocate(dest, len, info_IM);
1452
1453   else
1454     /* Fixed layer allocation */
1455
1456     tcd_rateallocate_fixed();
1457
1458 /*--------------TIER2------------------*/
1459   info_IM->index_write = 1;     /* INDEX */
1460
1461   l = t2_encode_packets(tcd_img, tcd_cp, tileno, tile,
1462                         tcd_tcp->numlayers, dest, len, info_IM);
1463
1464  /*---------------CLEAN-------------------*/
1465   time = clock() - time;
1466   fprintf(stdout,"total:     %ld.%.3ld s\n", time / CLOCKS_PER_SEC,
1467          (time % CLOCKS_PER_SEC) * 1000 / CLOCKS_PER_SEC);
1468
1469   for (compno = 0; compno < tile->numcomps; compno++) {
1470     tilec = &tile->comps[compno];
1471     free(tilec->data);
1472   }
1473
1474   return l;
1475 }
1476
1477
1478 int tcd_decode_tile(unsigned char *src, int len, int tileno)
1479 {
1480   int l;
1481   int compno;
1482   int eof = 0;
1483   clock_t time;
1484   tcd_tile_t *tile;
1485
1486   tcd_tileno = tileno;
1487   tcd_tile = &tcd_image.tiles[tileno];
1488   tcd_tcp = &tcd_cp->tcps[tileno];
1489   tile = tcd_tile;
1490
1491   time = clock();
1492
1493   fprintf(stdout, "Tile %d of %d decoded in ", tileno + 1,
1494           tcd_cp->tw * tcd_cp->th);
1495
1496         /*--------------TIER2------------------*/
1497
1498   l = t2_decode_packets(src, len, tcd_img, tcd_cp, tileno, tile);
1499
1500   if (l == -999) {
1501     eof = 1;
1502     fprintf(stderr, "tcd_decode: incomplete bistream\n");
1503   }
1504
1505         /*------------------TIER1-----------------*/
1506   t1_init_luts();
1507   t1_decode_cblks(tile, tcd_tcp);
1508
1509         /*----------------DWT---------------------*/
1510
1511   for (compno = 0; compno < tile->numcomps; compno++) {
1512     tcd_tilecomp_t *tilec = &tile->comps[compno];
1513     if (tcd_cp->reduce != 0) {
1514       tcd_img->comps[compno].resno_decoded =
1515         tile->comps[compno].numresolutions - tcd_cp->reduce - 1;
1516     }
1517
1518
1519     if (tcd_tcp->tccps[compno].qmfbid == 1) {
1520       dwt_decode(tilec->data, tilec->x1 - tilec->x0,
1521                  tilec->y1 - tilec->y0, tilec,
1522                  tilec->numresolutions - 1,
1523                  tilec->numresolutions - 1 -
1524                  tcd_img->comps[compno].resno_decoded);
1525     } else {
1526       dwt_decode_real(tilec->data, tilec->x1 - tilec->x0,
1527                       tilec->y1 - tilec->y0, tilec,
1528                       tilec->numresolutions - 1,
1529                       tilec->numresolutions - 1 -
1530                       tcd_img->comps[compno].resno_decoded);
1531     }
1532
1533     if (tile->comps[compno].numresolutions > 0)
1534       tcd_img->comps[compno].factor =
1535         tile->comps[compno].numresolutions -
1536         (tcd_img->comps[compno].resno_decoded + 1);
1537   }
1538
1539         /*----------------MCT-------------------*/
1540
1541   if (tcd_tcp->mct) {
1542     if (tcd_tcp->tccps[0].qmfbid == 1) {
1543       mct_decode(tile->comps[0].data, tile->comps[1].data,
1544                  tile->comps[2].data,
1545                  (tile->comps[0].x1 -
1546                   tile->comps[0].x0) * (tile->comps[0].y1 -
1547                                         tile->comps[0].y0));
1548     } else {
1549       mct_decode_real(tile->comps[0].data, tile->comps[1].data,
1550                       tile->comps[2].data,
1551                       (tile->comps[0].x1 -
1552                        tile->comps[0].x0) * (tile->comps[0].y1 -
1553                                              tile->comps[0].y0));
1554     }
1555   }
1556
1557         /*---------------TILE-------------------*/
1558
1559   for (compno = 0; compno < tile->numcomps; compno++) {
1560     tcd_tilecomp_t *tilec = &tile->comps[compno];
1561     tcd_resolution_t *res =
1562       &tilec->resolutions[tcd_img->comps[compno].resno_decoded];
1563     int adjust =
1564       tcd_img->comps[compno].sgnd ? 0 : 1 << (tcd_img->comps[compno].
1565                                               prec - 1);
1566     int min =
1567       tcd_img->comps[compno].
1568       sgnd ? -(1 << (tcd_img->comps[compno].prec - 1)) : 0;
1569     int max =
1570       tcd_img->comps[compno].
1571       sgnd ? (1 << (tcd_img->comps[compno].prec - 1)) -
1572       1 : (1 << tcd_img->comps[compno].prec) - 1;
1573
1574     int tw = tilec->x1 - tilec->x0;
1575     int w = tcd_img->comps[compno].w;
1576
1577     int i, j;
1578     int offset_x = int_ceildivpow2(tcd_img->comps[compno].x0,
1579                                    tcd_img->comps[compno].factor);
1580     int offset_y = int_ceildivpow2(tcd_img->comps[compno].y0,
1581                                    tcd_img->comps[compno].factor);
1582
1583     for (j = res->y0; j < res->y1; j++) {
1584       for (i = res->x0; i < res->x1; i++) {
1585
1586         int v;
1587
1588         double tmp =
1589           (double) tilec->data[i - res->x0 + (j - res->y0) * tw];
1590         if (tcd_tcp->tccps[compno].qmfbid == 1) {
1591           v = (int) tmp;
1592         } else {
1593
1594           //v = (int) tmp >> 13;
1595
1596           //Mod antonin : multbug1
1597           v =
1598             (int) ((fabs(tmp / 8192.0) >=
1599                     floor(fabs(tmp / 8192.0)) +
1600                     0.5) ? fabs(tmp / 8192.0) + 1.0 : fabs(tmp / 8192.0));
1601
1602           v = (tmp < 0) ? -v : v;
1603
1604           //doM
1605         }
1606         v += adjust;
1607
1608         tcd_img->comps[compno].data[(i - offset_x) +
1609                                     (j - offset_y) * w] =
1610           int_clamp(v, min, max);
1611       }
1612     }
1613   }
1614
1615   time = clock() - time;
1616   fprintf(stdout, "%ld.%.3ld s\n", time / CLOCKS_PER_SEC,
1617           (time % CLOCKS_PER_SEC) * 1000 / CLOCKS_PER_SEC);
1618
1619
1620
1621   for (compno = 0; compno < tile->numcomps; compno++) {
1622     free(tcd_image.tiles[tileno].comps[compno].data);
1623   }
1624
1625   if (eof) {
1626     longjmp(j2k_error, 1);
1627   }
1628
1629   return l;
1630 }
1631
1632
1633
1634 void tcd_dec_release()
1635
1636 {
1637
1638   int tileno,compno,resno,bandno,precno;
1639
1640   for (tileno=0;tileno<tcd_image.tw*tcd_image.th;tileno++) {
1641
1642     tcd_tile_t tile=tcd_image.tiles[tileno];
1643
1644     for (compno=0;compno<tile.numcomps;compno++) {
1645
1646       tcd_tilecomp_t tilec=tile.comps[compno];
1647
1648       for (resno=0;resno<tilec.numresolutions;resno++) {
1649
1650         tcd_resolution_t res=tilec.resolutions[resno];
1651
1652         for (bandno=0;bandno<res.numbands;bandno++) {
1653
1654           tcd_band_t band=res.bands[bandno];
1655
1656           for (precno=0;precno<res.ph*res.pw;precno++) {
1657
1658             tcd_precinct_t prec=band.precincts[precno];
1659
1660             if (prec.cblks!=NULL) free(prec.cblks);
1661
1662             if (prec.imsbtree!=NULL) free(prec.imsbtree);
1663
1664             if (prec.incltree!=NULL) free(prec.incltree);
1665
1666           }
1667
1668           if (band.precincts!=NULL) free(band.precincts);
1669
1670         }
1671
1672       }
1673
1674       if (tilec.resolutions!=NULL) free(tilec.resolutions);
1675
1676     }
1677
1678     if (tile.comps!=NULL) free(tile.comps);
1679
1680   }
1681
1682   if (tcd_image.tiles!=NULL) free(tcd_image.tiles);
1683
1684 }