src/lib/openjp3d/tcd.c

   1 /*
   2  * The copyright in this software is being made available under the 2-clauses
   3  * BSD License, included below. This software may be subject to other third
   4  * party and contributor rights, including patent rights, and no such rights
   5  * are granted under this license.
   6  *
   7  * Copyright (c) 2001-2003, David Janssens
   8  * Copyright (c) 2002-2003, Yannick Verschueren
   9  * Copyright (c) 2003-2005, Francois Devaux and Antonin Descampe
  10  * Copyright (c) 2005, Herve Drolon, FreeImage Team
  11  * Copyright (c) 2002-2005, Communications and remote sensing Laboratory, Universite catholique de Louvain, Belgium
  12  * Copyright (c) 2006, Mónica Díez, LPI-UVA, Spain
  13  * All rights reserved.
  14  *
  15  * Redistribution and use in source and binary forms, with or without
  16  * modification, are permitted provided that the following conditions
  17  * are met:
  18  * 1. Redistributions of source code must retain the above copyright
  19  *    notice, this list of conditions and the following disclaimer.
  20  * 2. Redistributions in binary form must reproduce the above copyright
  21  *    notice, this list of conditions and the following disclaimer in the
  22  *    documentation and/or other materials provided with the distribution.
  23  *
  24  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS `AS IS'
  25  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  26  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  27  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  28  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  29  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  30  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  31  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  32  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  33  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  34  * POSSIBILITY OF SUCH DAMAGE.
  35  */
  36
  37 #include "opj_includes.h"
  38
  39 void tcd_dump(FILE *fd, opj_tcd_t *tcd, opj_tcd_volume_t * vol)
  40 {
  41     int tileno, compno, resno, bandno, precno, cblkno;
  42
  43     fprintf(fd, "volume {\n");
  44     fprintf(fd, "  tw=%d, th=%d, tl=%d, x0=%d x1=%d y0=%d y1=%d z0=%d z1=%d\n",
  45             vol->tw, vol->th, vol->tl, tcd->volume->x0, tcd->volume->x1, tcd->volume->y0,
  46             tcd->volume->y1, tcd->volume->z0, tcd->volume->z1);
  47
  48     for (tileno = 0; tileno < vol->th * vol->tw * vol->tl; tileno++) {
  49         opj_tcd_tile_t *tile = &tcd->tcd_volume->tiles[tileno];
  50         fprintf(fd, "  tile {\n");
  51         fprintf(fd, "    x0=%d, y0=%d, z0=%d, x1=%d, y1=%d, z1=%d, numcomps=%d\n",
  52                 tile->x0, tile->y0, tile->z0, tile->x1, tile->y1, tile->z1, tile->numcomps);
  53         for (compno = 0; compno < tile->numcomps; compno++) {
  54             opj_tcd_tilecomp_t *tilec = &tile->comps[compno];
  55             fprintf(fd, "    tilecomp %d {\n", compno);
  56             fprintf(fd,
  57                     "     x0=%d, y0=%d, z0=%d, x1=%d, y1=%d, z1=%d, numresx=%d, numresy=%d, numresz=%d\n",
  58                     tilec->x0, tilec->y0, tilec->z0, tilec->x1, tilec->y1, tilec->z1,
  59                     tilec->numresolution[0], tilec->numresolution[1], tilec->numresolution[2]);
  60             for (resno = 0; resno < tilec->numresolution[0]; resno++) {
  61                 opj_tcd_resolution_t *res = &tilec->resolutions[resno];
  62                 fprintf(fd, "     res %d{\n", resno);
  63                 fprintf(fd,
  64                         "      x0=%d, y0=%d, z0=%d, x1=%d, y1=%d, z1=%d, pw=%d, ph=%d, pl=%d, numbands=%d\n",
  65                         res->x0, res->y0, res->z0, res->x1, res->y1, res->z1, res->prctno[0],
  66                         res->prctno[1], res->prctno[2], res->numbands);
  67                 for (bandno = 0; bandno < res->numbands; bandno++) {
  68                     opj_tcd_band_t *band = &res->bands[bandno];
  69                     fprintf(fd, "       band %d{\n", bandno);
  70                     fprintf(fd,
  71                             "            x0=%d, y0=%d, z0=%d, x1=%d, y1=%d, z1=%d, stepsize=%f, numbps=%d\n",
  72                             band->x0, band->y0, band->z0, band->x1, band->y1, band->z1, band->stepsize,
  73                             band->numbps);
  74                     for (precno = 0; precno < (res->prctno[0] * res->prctno[1] * res->prctno[2]);
  75                             precno++) {
  76                         opj_tcd_precinct_t *prec = &band->precincts[precno];
  77                         fprintf(fd, "             prec %d{\n", precno);
  78                         fprintf(fd,
  79                                 "                  x0=%d, y0=%d, z0=%d, x1=%d, y1=%d, z1=%d, cw=%d, ch=%d, cl=%d,\n",
  80                                 prec->x0, prec->y0, prec->z0, prec->x1, prec->y1, prec->z1, prec->cblkno[0],
  81                                 prec->cblkno[1], prec->cblkno[2]);
  82                         for (cblkno = 0; cblkno < (prec->cblkno[0] * prec->cblkno[1] * prec->cblkno[2]);
  83                                 cblkno++) {
  84                             opj_tcd_cblk_t *cblk = &prec->cblks[cblkno];
  85                             fprintf(fd, "                   cblk %d{\n", cblkno);
  86                             fprintf(fd, "                    x0=%d, y0=%d, z0=%d, x1=%d, y1=%d, z1=%d\n", cblk->x0,
  87                                     cblk->y0, cblk->z0, cblk->x1, cblk->y1, cblk->z1);
  88                             fprintf(fd, "            }\n");
  89                         }
  90                         fprintf(fd, "          }\n");
  91                     }
  92                     fprintf(fd, "        }\n");
  93                 }
  94                 fprintf(fd, "      }\n");
  95             }
  96             fprintf(fd, "    }\n");
  97         }
  98         fprintf(fd, "  }\n");
  99     }
 100     fprintf(fd, "}\n");
 101 }
 102
 103 static void tilec_dump(FILE *fd, opj_tcd_tilecomp_t *tilec)
 104 {
 105
 106     int i = 0, k;
 107     int datalen;
 108     int *a;
 109
 110     fprintf(fd, "    tilecomp{\n");
 111     fprintf(fd,
 112             "     x0=%d, y0=%d, z0=%d, x1=%d, y1=%d, z1=%d, numresx=%d, numresy=%d, numresz=%d\n",
 113             tilec->x0, tilec->y0, tilec->z0, tilec->x1, tilec->y1, tilec->z1,
 114             tilec->numresolution[0], tilec->numresolution[1], tilec->numresolution[2]);
 115     fprintf(fd, "     data {\n");
 116     datalen = (tilec->z1 - tilec->z0) * (tilec->y1 - tilec->y0) *
 117               (tilec->x1 - tilec->x0);
 118     a = tilec->data;
 119     for (k = 0; k < datalen; k++) {
 120         if (!(k % tilec->x1)) {
 121             fprintf(fd, "\n");
 122         }
 123         if (!(k % (tilec->y1 * tilec->x1))) {
 124             fprintf(fd, "Slice %d\n", i++);
 125         }
 126         fprintf(fd, " %d", a[k]);
 127
 128
 129     }
 130     fprintf(fd, "     }\n");
 131     /*i=0;
 132     fprintf(fd, "Slice %d\n");
 133     if (tilec->prediction->prederr) {
 134         fprintf(fd, "     prederror {\n");
 135         a = tilec->prediction->prederr;
 136         for (k = 0; k < datalen; k++) {
 137             fprintf(fd," %d",*(a++));
 138             if (!(k % (tilec->y1 - tilec->y0) * (tilec->x1 - tilec->x0))){
 139                 fprintf(fd, "\n");fprintf(fd, "Slice %d\n",i++);
 140             }
 141             if (!(k % (tilec->x1 - tilec->x0))){
 142                 fprintf(fd, "\n");
 143             }
 144         }
 145     }
 146     fprintf(fd, "     }\n");*/
 147     fprintf(fd, "}\n");
 148 }
 149
 150 /* ----------------------------------------------------------------------- */
 151
 152 /**
 153 Create a new TCD handle
 154 */
 155 opj_tcd_t* tcd_create(opj_common_ptr cinfo)
 156 {
 157     /* create the tcd structure */
 158     opj_tcd_t *tcd = (opj_tcd_t*)opj_malloc(sizeof(opj_tcd_t));
 159     if (!tcd) {
 160         return NULL;
 161     }
 162     tcd->cinfo = cinfo;
 163     tcd->tcd_volume = (opj_tcd_volume_t*)opj_malloc(sizeof(opj_tcd_volume_t));
 164     if (!tcd->tcd_volume) {
 165         opj_free(tcd);
 166         return NULL;
 167     }
 168
 169     return tcd;
 170 }
 171
 172 /**
 173 Destroy a previously created TCD handle
 174 */
 175 void tcd_destroy(opj_tcd_t *tcd)
 176 {
 177     if (tcd) {
 178         opj_free(tcd->tcd_volume);
 179         opj_free(tcd);
 180     }
 181 }
 182
 183 /* ----------------------------------------------------------------------- */
 184 void tcd_malloc_encode(opj_tcd_t *tcd, opj_volume_t * volume, opj_cp_t * cp,
 185                        int curtileno)
 186 {
 187     int compno, resno, bandno, precno, cblkno, i, j;/*, k;*/
 188
 189     opj_tcd_tile_t      *tile = NULL;       /* pointer to tcd->tile */
 190     opj_tcd_tilecomp_t  *tilec = NULL;      /* pointer to tcd->tilec */
 191     opj_tcd_resolution_t    *res = NULL;        /* pointer to tcd->res */
 192     opj_tcd_band_t      *band = NULL;       /* pointer to tcd->band */
 193     opj_tcd_precinct_t  *prc = NULL;        /* pointer to tcd->prc */
 194     opj_tcd_cblk_t      *cblk = NULL;       /* pointer to tcd->cblk */
 195     opj_tcp_t       *tcp = &cp->tcps[curtileno];
 196     int p, q, r;
 197
 198     tcd->volume = volume;
 199     tcd->cp = cp;
 200     tcd->tcd_volume->tw = cp->tw;
 201     tcd->tcd_volume->th = cp->th;
 202     tcd->tcd_volume->tl = cp->tl;
 203     tcd->tcd_volume->tiles = (opj_tcd_tile_t *) opj_malloc(sizeof(opj_tcd_tile_t));
 204     tcd->tile = tcd->tcd_volume->tiles;
 205     tile = tcd->tile;
 206
 207
 208     /* p61 ISO/IEC IS15444-1 : 2002 */
 209     /* curtileno --> raster scanned index of tiles */
 210     /* p,q,r --> matricial index of tiles */
 211     p = curtileno % cp->tw;
 212     q = curtileno / cp->tw;
 213     r = curtileno / (cp->tw * cp->th); /* extension to 3-D */
 214
 215     /* 4 borders of the tile rescale on the volume if necessary (B.3)*/
 216     tile->x0 = int_max(cp->tx0 + p * cp->tdx, volume->x0);
 217     tile->y0 = int_max(cp->ty0 + q * cp->tdy, volume->y0);
 218     tile->z0 = int_max(cp->tz0 + r * cp->tdz, volume->z0);
 219     tile->x1 = int_min(cp->tx0 + (p + 1) * cp->tdx, volume->x1);
 220     tile->y1 = int_min(cp->ty0 + (q + 1) * cp->tdy, volume->y1);
 221     tile->z1 = int_min(cp->tz0 + (r + 1) * cp->tdz, volume->z1);
 222     tile->numcomps = volume->numcomps;
 223
 224     /* Modification of the RATE >> */
 225     for (j = 0; j < tcp->numlayers; j++) {
 226         if (tcp->rates[j] <= 1) {
 227             tcp->rates[j] = 0;
 228         } else {
 229             float num = (float)(tile->numcomps * (tile->x1 - tile->x0) *
 230                                 (tile->y1 - tile->y0) * (tile->z1 - tile->z0) * volume->comps[0].prec);
 231             float den = (float)(8 * volume->comps[0].dx * volume->comps[0].dy *
 232                                 volume->comps[0].dz);
 233             den = tcp->rates[j] * den;
 234             tcp->rates[j] = (num + den - 1) / den;
 235         }
 236         /*tcp->rates[j] = tcp->rates[j] ? int_ceildiv(
 237             tile->numcomps * (tile->x1 - tile->x0) * (tile->y1 - tile->y0) * (tile->z1 - tile->z0) * volume->comps[0].prec,
 238             (tcp->rates[j] * 8 * volume->comps[0].dx * volume->comps[0].dy * volume->comps[0].dz)) : 0;*/
 239         if (tcp->rates[j]) {
 240             if (j && tcp->rates[j] < tcp->rates[j - 1] + 10) {
 241                 tcp->rates[j] = tcp->rates[j - 1] + 20;
 242             } else if (!j && tcp->rates[j] < 30) {
 243                 tcp->rates[j] = 30;
 244             }
 245         }
 246     }
 247     /* << Modification of the RATE */
 248
 249     tile->comps = (opj_tcd_tilecomp_t *) opj_malloc(volume->numcomps * sizeof(
 250                       opj_tcd_tilecomp_t));
 251     for (compno = 0; compno < tile->numcomps; compno++) {
 252         opj_tccp_t *tccp = &tcp->tccps[compno];
 253         int res_max;
 254         int prevnumbands = 0;
 255
 256         /* opj_tcd_tilecomp_t *tilec=&tile->comps[compno]; */
 257         tcd->tilec = &tile->comps[compno];
 258         tilec = tcd->tilec;
 259
 260         /* border of each tile component (global) (B.3) */
 261         tilec->x0 = int_ceildiv(tile->x0, volume->comps[compno].dx);
 262         tilec->y0 = int_ceildiv(tile->y0, volume->comps[compno].dy);
 263         tilec->z0 = int_ceildiv(tile->z0, volume->comps[compno].dz);
 264         tilec->x1 = int_ceildiv(tile->x1, volume->comps[compno].dx);
 265         tilec->y1 = int_ceildiv(tile->y1, volume->comps[compno].dy);
 266         tilec->z1 = int_ceildiv(tile->z1, volume->comps[compno].dz);
 267
 268         tilec->data = (int *) opj_malloc((tilec->x1 - tilec->x0) *
 269                                          (tilec->y1 - tilec->y0) * (tilec->z1 - tilec->z0) * sizeof(int));
 270
 271         res_max = 0;
 272         for (i = 0; i < 3; i++) {
 273             tilec->numresolution[i] = tccp->numresolution[i];
 274             /*Greater of 3 resolutions contains all information*/
 275             res_max = (tilec->numresolution[i] > res_max) ? tilec->numresolution[i] :
 276                       res_max;
 277         }
 278
 279
 280         tilec->resolutions = (opj_tcd_resolution_t *) opj_malloc(res_max * sizeof(
 281                                  opj_tcd_resolution_t));
 282         for (resno = 0; resno < res_max; resno++) {
 283
 284             int pdx, pdy, pdz;
 285             int tlprcxstart, tlprcystart, tlprczstart;
 286             int brprcxend, brprcyend, brprczend;
 287             int tlcbgxstart, tlcbgystart, tlcbgzstart;
 288             int brcbgxend, brcbgyend, brcbgzend;
 289             int cbgwidthexpn, cbgheightexpn, cbglengthexpn;
 290             int cblkwidthexpn, cblkheightexpn, cblklengthexpn;
 291
 292             int diff = tccp->numresolution[0] - tccp->numresolution[2];
 293             int levelnox = tilec->numresolution[0] - 1 - resno;
 294             int levelnoy = tilec->numresolution[1] - 1 - resno;
 295             int levelnoz = tilec->numresolution[2] - 1 - ((resno <= diff) ? 0 :
 296                            (resno - diff));
 297             if (levelnoz < 0) {
 298                 levelnoz = 0;
 299             }
 300
 301             /* opj_tcd_resolution_t *res=&tilec->resolutions[resno]; */
 302             tcd->res = &tilec->resolutions[resno];
 303             res = tcd->res;
 304
 305             /* border for each resolution level (global) (B.14)*/
 306             res->x0 = int_ceildivpow2(tilec->x0, levelnox);
 307             res->y0 = int_ceildivpow2(tilec->y0, levelnoy);
 308             res->z0 = int_ceildivpow2(tilec->z0, levelnoz);
 309             res->x1 = int_ceildivpow2(tilec->x1, levelnox);
 310             res->y1 = int_ceildivpow2(tilec->y1, levelnoy);
 311             res->z1 = int_ceildivpow2(tilec->z1, levelnoz);
 312             /*if (res->z1 < 0)fprintf(stdout,"Res: %d       %d/%d --> %d\n",resno,tilec->z1, levelnoz, int_ceildivpow2(tilec->z1, levelnoz));*/
 313
 314             res->numbands = (resno == 0) ? 1 : (resno <= diff) ? 3 : 7; /* --> 3D */
 315
 316             /* p. 30, table A-13, ISO/IEC IS154444-1 : 2002 */
 317             if (tccp->csty & J3D_CCP_CSTY_PRT) {
 318                 pdx = tccp->prctsiz[0][resno];
 319                 pdy = tccp->prctsiz[1][resno];
 320                 pdz = tccp->prctsiz[2][resno];
 321             } else {
 322                 pdx = 15;
 323                 pdy = 15;
 324                 pdz = 15;
 325             }
 326
 327             /* p. 66, B.16, ISO/IEC IS15444-1 : 2002  */
 328             tlprcxstart = int_floordivpow2(res->x0, pdx) << pdx;
 329             tlprcystart = int_floordivpow2(res->y0, pdy) << pdy;
 330             tlprczstart = int_floordivpow2(res->z0, pdz) << pdz;
 331             brprcxend = int_ceildivpow2(res->x1, pdx) << pdx;
 332             brprcyend = int_ceildivpow2(res->y1, pdy) << pdy;
 333             brprczend = int_ceildivpow2(res->z1, pdz) << pdz;
 334
 335             res->prctno[0] = (brprcxend - tlprcxstart) >> pdx;
 336             res->prctno[1] = (brprcyend - tlprcystart) >> pdy;
 337             res->prctno[2] = (brprczend - tlprczstart) >> pdz;
 338             if (res->prctno[2] == 0) {
 339                 res->prctno[2] = 1;
 340             }
 341
 342             /* p. 67, B.17 & B.18, ISO/IEC IS15444-1 : 2002  */
 343             if (resno == 0) {
 344                 tlcbgxstart = tlprcxstart;
 345                 tlcbgystart = tlprcystart;
 346                 tlcbgzstart = tlprczstart;
 347                 brcbgxend = brprcxend;
 348                 brcbgyend = brprcyend;
 349                 brcbgzend = brprczend;
 350                 cbgwidthexpn = pdx;
 351                 cbgheightexpn = pdy;
 352                 cbglengthexpn = pdz;
 353             } else {
 354                 tlcbgxstart = int_ceildivpow2(tlprcxstart, 1);
 355                 tlcbgystart = int_ceildivpow2(tlprcystart, 1);
 356                 tlcbgzstart = int_ceildivpow2(tlprczstart, 1);
 357                 brcbgxend = int_ceildivpow2(brprcxend, 1);
 358                 brcbgyend = int_ceildivpow2(brprcyend, 1);
 359                 brcbgzend = int_ceildivpow2(brprczend, 1);
 360                 cbgwidthexpn = pdx - 1;
 361                 cbgheightexpn = pdy - 1;
 362                 cbglengthexpn = pdz - 1;
 363             }
 364
 365             cblkwidthexpn = int_min(tccp->cblk[0], cbgwidthexpn); /*6*/
 366             cblkheightexpn = int_min(tccp->cblk[1], cbgheightexpn); /*6*/
 367             cblklengthexpn = int_min(tccp->cblk[2], cbglengthexpn); /*6*/
 368
 369             res->bands = (opj_tcd_band_t *) opj_malloc(res->numbands * sizeof(
 370                              opj_tcd_band_t));
 371             for (bandno = 0; bandno < res->numbands; bandno++) {
 372                 int x0b, y0b, z0b, i;
 373                 int gain, numbps;
 374                 opj_stepsize_t *ss = NULL;
 375
 376                 tcd->band = &res->bands[bandno];
 377                 band = tcd->band;
 378
 379                 band->bandno = (resno == 0) ? 0 : bandno + 1;
 380                 /* Bandno:  0 - LLL     2 - LHL
 381                             1 - HLL     3 - HHL
 382                             4 - LLH     6 - LHH
 383                             5 - HLH     7 - HHH     */
 384                 x0b = (band->bandno == 1) || (band->bandno == 3) || (band->bandno == 5) ||
 385                       (band->bandno == 7) ? 1 : 0;
 386                 y0b = (band->bandno == 2) || (band->bandno == 3) || (band->bandno == 6) ||
 387                       (band->bandno == 7) ? 1 : 0;
 388                 z0b = (band->bandno == 4) || (band->bandno == 5) || (band->bandno == 6) ||
 389                       (band->bandno == 7) ? 1 : 0;
 390
 391                 /* p. 65, B.15, ISO/IEC IS15444-1 : 2002  */
 392                 if (band->bandno == 0) {
 393                     /* band border (global) */
 394                     band->x0 = int_ceildivpow2(tilec->x0, levelnox);
 395                     band->y0 = int_ceildivpow2(tilec->y0, levelnoy);
 396                     band->z0 = int_ceildivpow2(tilec->z0, levelnoz);
 397                     band->x1 = int_ceildivpow2(tilec->x1, levelnox);
 398                     band->y1 = int_ceildivpow2(tilec->y1, levelnoy);
 399                     band->z1 = int_ceildivpow2(tilec->z1, levelnoz);
 400                 } else {
 401                     band->x0 = int_ceildivpow2(tilec->x0 - (1 << levelnox) * x0b, levelnox + 1);
 402                     band->y0 = int_ceildivpow2(tilec->y0 - (1 << levelnoy) * y0b, levelnoy + 1);
 403                     band->z0 = int_ceildivpow2(tilec->z0 - (1 << levelnoz) * z0b,
 404                                                (resno <= diff) ? levelnoz : levelnoz + 1);
 405                     band->x1 = int_ceildivpow2(tilec->x1 - (1 << levelnox) * x0b, levelnox + 1);
 406                     band->y1 = int_ceildivpow2(tilec->y1 - (1 << levelnoy) * y0b, levelnoy + 1);
 407                     band->z1 = int_ceildivpow2(tilec->z1 - (1 << levelnoz) * z0b,
 408                                                (resno <= diff) ? levelnoz : levelnoz + 1);
 409                 }
 410
 411                 ss = &tccp->stepsizes[(resno == 0) ? 0 : (prevnumbands + bandno + 1)];
 412                 if (bandno == (res->numbands - 1)) {
 413                     prevnumbands += (resno == 0) ? 0 : res->numbands;
 414                 }
 415                 gain = dwt_getgain(band->bandno, tccp->reversible);
 416                 numbps = volume->comps[compno].prec + gain;
 417                 band->stepsize = (float)((1.0 + ss->mant / 2048.0) * pow(2.0,
 418                                          numbps - ss->expn));
 419                 band->numbps = ss->expn + tccp->numgbits - 1;   /* WHY -1 ? */
 420
 421                 band->precincts = (opj_tcd_precinct_t *) opj_malloc((res->prctno[0] *
 422                                   res->prctno[1] * res->prctno[2]) * sizeof(opj_tcd_precinct_t));
 423
 424                 for (i = 0; i < (res->prctno[0] * res->prctno[1] * res->prctno[2]); i++) {
 425                     band->precincts[i].imsbtree = NULL;
 426                     band->precincts[i].incltree = NULL;
 427                 }
 428
 429                 for (precno = 0; precno < (res->prctno[0] * res->prctno[1] * res->prctno[2]);
 430                         precno++) {
 431                     int tlcblkxstart, tlcblkystart, tlcblkzstart, brcblkxend, brcblkyend,
 432                         brcblkzend;
 433                     int cbgxstart, cbgystart, cbgzstart, cbgxend, cbgyend, cbgzend;
 434
 435                     cbgxstart = tlcbgxstart + (precno % res->prctno[0]) * (1 << cbgwidthexpn);
 436                     cbgystart = tlcbgystart + ((precno % (res->prctno[0] * res->prctno[1])) /
 437                                                res->prctno[0]) * (1 << cbgheightexpn);
 438                     cbgzstart = tlcbgzstart + (precno / (res->prctno[0] * res->prctno[1])) *
 439                                 (1 << cbglengthexpn);
 440                     cbgxend = cbgxstart + (1 << cbgwidthexpn);
 441                     cbgyend = cbgystart + (1 << cbgheightexpn);
 442                     cbgzend = cbgzstart + (1 << cbglengthexpn);
 443
 444                     tcd->prc = &band->precincts[precno];
 445                     prc = tcd->prc;
 446
 447                     /* precinct size (global) */
 448                     prc->x0 = int_max(cbgxstart, band->x0);
 449                     prc->y0 = int_max(cbgystart, band->y0);
 450                     prc->z0 = int_max(cbgzstart, band->z0);
 451                     prc->x1 = int_min(cbgxend, band->x1);
 452                     prc->y1 = int_min(cbgyend, band->y1);
 453                     prc->z1 = int_min(cbgzend, band->z1);
 454
 455                     tlcblkxstart = int_floordivpow2(prc->x0, cblkwidthexpn) << cblkwidthexpn;
 456                     tlcblkystart = int_floordivpow2(prc->y0, cblkheightexpn) << cblkheightexpn;
 457                     tlcblkzstart = int_floordivpow2(prc->z0, cblklengthexpn) << cblklengthexpn;
 458                     brcblkxend = int_ceildivpow2(prc->x1, cblkwidthexpn) << cblkwidthexpn;
 459                     brcblkyend = int_ceildivpow2(prc->y1, cblkheightexpn) << cblkheightexpn;
 460                     brcblkzend = int_ceildivpow2(prc->z1, cblklengthexpn) << cblklengthexpn;
 461                     prc->cblkno[0] = (brcblkxend - tlcblkxstart) >> cblkwidthexpn;
 462                     prc->cblkno[1] = (brcblkyend - tlcblkystart) >> cblkheightexpn;
 463                     prc->cblkno[2] = (brcblkzend - tlcblkzstart) >> cblklengthexpn;
 464                     prc->cblkno[2] = (prc->cblkno[2] == 0) ? 1 : prc->cblkno[2];
 465
 466                     prc->cblks = (opj_tcd_cblk_t *) opj_malloc((prc->cblkno[0] * prc->cblkno[1] *
 467                                  prc->cblkno[2]) * sizeof(opj_tcd_cblk_t));
 468                     prc->incltree = tgt_create(prc->cblkno[0], prc->cblkno[1], prc->cblkno[2]);
 469                     prc->imsbtree = tgt_create(prc->cblkno[0], prc->cblkno[1], prc->cblkno[2]);
 470                     /*tgt_tree_dump(stdout,prc->incltree);*/
 471                     for (cblkno = 0; cblkno < (prc->cblkno[0] * prc->cblkno[1] * prc->cblkno[2]);
 472                             cblkno++) {
 473                         int cblkxstart = tlcblkxstart + (cblkno % prc->cblkno[0]) *
 474                                          (1 << cblkwidthexpn);
 475                         int cblkystart = tlcblkystart + ((cblkno % (prc->cblkno[0] * prc->cblkno[1])) /
 476                                                          prc->cblkno[0]) * (1 << cblkheightexpn);
 477                         int cblkzstart = tlcblkzstart + (cblkno / (prc->cblkno[0] * prc->cblkno[1])) *
 478                                          (1 << cblklengthexpn);
 479                         int cblkxend = cblkxstart + (1 << cblkwidthexpn);
 480                         int cblkyend = cblkystart + (1 << cblkheightexpn);
 481                         int cblkzend = cblkzstart + (1 << cblklengthexpn);
 482                         int prec = ((tilec->bpp > 16) ? 3 : ((tilec->bpp > 8) ? 2 : 1));
 483
 484                         tcd->cblk = &prc->cblks[cblkno];
 485                         cblk = tcd->cblk;
 486
 487                         /* code-block size (global) */
 488                         cblk->x0 = int_max(cblkxstart, prc->x0);
 489                         cblk->y0 = int_max(cblkystart, prc->y0);
 490                         cblk->z0 = int_max(cblkzstart, prc->z0);
 491                         cblk->x1 = int_min(cblkxend, prc->x1);
 492                         cblk->y1 = int_min(cblkyend, prc->y1);
 493                         cblk->z1 = int_min(cblkzend, prc->z1);
 494                     }
 495                 }
 496             }
 497         }
 498     }
 499     /*tcd_dump(stdout, tcd, tcd->tcd_volume);*/
 500
 501 }
 502 void tcd_init_encode(opj_tcd_t *tcd, opj_volume_t * volume, opj_cp_t * cp,
 503                      int curtileno)
 504 {
 505     int compno, resno, bandno, precno, cblkno;
 506     int j, p, q, r;
 507
 508     opj_tcd_tile_t      *tile = NULL;       /* pointer to tcd->tile */
 509     opj_tcd_tilecomp_t  *tilec = NULL;      /* pointer to tcd->tilec */
 510     opj_tcd_resolution_t    *res = NULL;    /* pointer to tcd->res */
 511     opj_tcd_band_t      *band = NULL;       /* pointer to tcd->band */
 512     opj_tcd_precinct_t  *prc = NULL;        /* pointer to tcd->prc */
 513     opj_tcd_cblk_t      *cblk = NULL;       /* pointer to tcd->cblk */
 514     opj_tcp_t *tcp = &cp->tcps[curtileno];
 515
 516     tcd->tile = tcd->tcd_volume->tiles;
 517     tile = tcd->tile;
 518
 519     /* p61 ISO/IEC IS15444-1 : 2002 */
 520     /* curtileno --> raster scanned index of tiles */
 521     /* p,q,r --> matricial index of tiles */
 522     p = curtileno % cp->tw;
 523     q = curtileno / cp->tw;
 524     r = curtileno / (cp->tw * cp->th); /* extension to 3-D */
 525
 526     /* 4 borders of the tile rescale on the volume if necessary (B.3)*/
 527     tile->x0 = int_max(cp->tx0 + p * cp->tdx, volume->x0);
 528     tile->y0 = int_max(cp->ty0 + q * cp->tdy, volume->y0);
 529     tile->z0 = int_max(cp->tz0 + r * cp->tdz, volume->z0);
 530     tile->x1 = int_min(cp->tx0 + (p + 1) * cp->tdx, volume->x1);
 531     tile->y1 = int_min(cp->ty0 + (q + 1) * cp->tdy, volume->y1);
 532     tile->z1 = int_min(cp->tz0 + (r + 1) * cp->tdz, volume->z1);
 533     tile->numcomps = volume->numcomps;
 534
 535     /* Modification of the RATE >> */
 536     for (j = 0; j < tcp->numlayers; j++) {
 537         if (tcp->rates[j] <= 1) {
 538             tcp->rates[j] = 0;
 539         } else {
 540             float num = (float)(tile->numcomps * (tile->x1 - tile->x0) *
 541                                 (tile->y1 - tile->y0) * (tile->z1 - tile->z0) * volume->comps[0].prec);
 542             float den = (float)(8 * volume->comps[0].dx * volume->comps[0].dy *
 543                                 volume->comps[0].dz);
 544             den = tcp->rates[j] * den;
 545             tcp->rates[j] = (num + den - 1) / den;
 546         }
 547         /*tcp->rates[j] = tcp->rates[j] ? int_ceildiv(
 548             tile->numcomps * (tile->x1 - tile->x0) * (tile->y1 - tile->y0) * (tile->z1 - tile->z0) * volume->comps[0].prec,
 549             (tcp->rates[j] * 8 * volume->comps[0].dx * volume->comps[0].dy * volume->comps[0].dz)) : 0;*/
 550         if (tcp->rates[j]) {
 551             if (j && tcp->rates[j] < tcp->rates[j - 1] + 10) {
 552                 tcp->rates[j] = tcp->rates[j - 1] + 20;
 553             } else if (!j && tcp->rates[j] < 30) {
 554                 tcp->rates[j] = 30;
 555             }
 556         }
 557     }
 558     /* << Modification of the RATE */
 559
 560     for (compno = 0; compno < tile->numcomps; compno++) {
 561         opj_tccp_t *tccp = &tcp->tccps[compno];
 562         int res_max, i;
 563         int prevnumbands = 0;
 564
 565         /* opj_tcd_tilecomp_t *tilec=&tile->comps[compno]; */
 566         tcd->tilec = &tile->comps[compno];
 567         tilec = tcd->tilec;
 568
 569         /* border of each tile component (global) (B.3) */
 570         tilec->x0 = int_ceildiv(tile->x0, volume->comps[compno].dx);
 571         tilec->y0 = int_ceildiv(tile->y0, volume->comps[compno].dy);
 572         tilec->z0 = int_ceildiv(tile->z0, volume->comps[compno].dz);
 573         tilec->x1 = int_ceildiv(tile->x1, volume->comps[compno].dx);
 574         tilec->y1 = int_ceildiv(tile->y1, volume->comps[compno].dy);
 575         tilec->z1 = int_ceildiv(tile->z1, volume->comps[compno].dz);
 576
 577         tilec->data = (int *) opj_malloc((tilec->x1 - tilec->x0) *
 578                                          (tilec->y1 - tilec->y0) * (tilec->z1 - tilec->z0) * sizeof(int));
 579
 580         res_max = 0;
 581         for (i = 0; i < 3; i++) {
 582             tilec->numresolution[i] = tccp->numresolution[i];
 583             /*Greater of 3 resolutions contains all information*/
 584             res_max = (tilec->numresolution[i] > res_max) ? tilec->numresolution[i] :
 585                       res_max;
 586         }
 587
 588         tilec->resolutions = (opj_tcd_resolution_t *) opj_malloc(res_max * sizeof(
 589                                  opj_tcd_resolution_t));
 590         for (resno = 0; resno < res_max; resno++) {
 591             int pdx, pdy, pdz;
 592             int tlprcxstart, tlprcystart, tlprczstart, brprcxend, brprcyend, brprczend;
 593             int tlcbgxstart, tlcbgystart, tlcbgzstart, brcbgxend, brcbgyend, brcbgzend;
 594             int cbgwidthexpn, cbgheightexpn, cbglengthexpn;
 595             int cblkwidthexpn, cblkheightexpn, cblklengthexpn;
 596
 597             int levelnox = tilec->numresolution[0] - 1 - resno;
 598             int levelnoy = tilec->numresolution[1] - 1 - resno;
 599             int diff = tccp->numresolution[0] - tccp->numresolution[2];
 600             int levelnoz = tilec->numresolution[2] - 1 - ((resno <= diff) ? 0 :
 601                            (resno - diff));
 602             if (levelnoz < 0) {
 603                 levelnoz = 0;
 604             }
 605
 606             tcd->res = &tilec->resolutions[resno];
 607             res = tcd->res;
 608
 609             /* border for each resolution level (global) (B.14)*/
 610             res->x0 = int_ceildivpow2(tilec->x0, levelnox);
 611             res->y0 = int_ceildivpow2(tilec->y0, levelnoy);
 612             res->z0 = int_ceildivpow2(tilec->z0, levelnoz);
 613             res->x1 = int_ceildivpow2(tilec->x1, levelnox);
 614             res->y1 = int_ceildivpow2(tilec->y1, levelnoy);
 615             res->z1 = int_ceildivpow2(tilec->z1, levelnoz);
 616
 617             /* res->numbands = resno == 0 ? 1 : 3; *//* --> 2D */
 618
 619             res->numbands = (resno == 0) ? 1 : (resno <= diff) ? 3 : 7; /* --> 3D */
 620
 621             /* p. 30, table A-13, ISO/IEC IS154444-1 : 2002 */
 622             if (tccp->csty & J3D_CCP_CSTY_PRT) {
 623                 pdx = tccp->prctsiz[0][resno];
 624                 pdy = tccp->prctsiz[1][resno];
 625                 pdz = tccp->prctsiz[2][resno];
 626             } else {
 627                 pdx = 15;
 628                 pdy = 15;
 629                 pdz = 15;
 630             }
 631             /* p. 66, B.16, ISO/IEC IS15444-1 : 2002  */
 632             tlprcxstart = int_floordivpow2(res->x0, pdx) << pdx;
 633             tlprcystart = int_floordivpow2(res->y0, pdy) << pdy;
 634             tlprczstart = int_floordivpow2(res->z0, pdz) << pdz;
 635             brprcxend = int_ceildivpow2(res->x1, pdx) << pdx;
 636             brprcyend = int_ceildivpow2(res->y1, pdy) << pdy;
 637             brprczend = int_ceildivpow2(res->z1, pdz) << pdz;
 638
 639             res->prctno[0] = (brprcxend - tlprcxstart) >> pdx;
 640             res->prctno[1] = (brprcyend - tlprcystart) >> pdy;
 641             res->prctno[2] = (brprczend - tlprczstart) >> pdz;
 642             if (res->prctno[2] == 0) {
 643                 res->prctno[2] = 1;
 644             }
 645
 646             /* p. 67, B.17 & B.18, ISO/IEC IS15444-1 : 2002  */
 647             if (resno == 0) {
 648                 tlcbgxstart = tlprcxstart;
 649                 tlcbgystart = tlprcystart;
 650                 tlcbgzstart = tlprczstart;
 651                 brcbgxend = brprcxend;
 652                 brcbgyend = brprcyend;
 653                 brcbgzend = brprczend;
 654                 cbgwidthexpn = pdx;
 655                 cbgheightexpn = pdy;
 656                 cbglengthexpn = pdz;
 657             } else {
 658                 tlcbgxstart = int_ceildivpow2(tlprcxstart, 1);
 659                 tlcbgystart = int_ceildivpow2(tlprcystart, 1);
 660                 tlcbgzstart = int_ceildivpow2(tlprczstart, 1);
 661                 brcbgxend = int_ceildivpow2(brprcxend, 1);
 662                 brcbgyend = int_ceildivpow2(brprcyend, 1);
 663                 brcbgzend = int_ceildivpow2(brprczend, 1);
 664                 cbgwidthexpn = pdx - 1;
 665                 cbgheightexpn = pdy - 1;
 666                 cbglengthexpn = pdz - 1;
 667             }
 668
 669             cblkwidthexpn = int_min(tccp->cblk[0], cbgwidthexpn);
 670             cblkheightexpn = int_min(tccp->cblk[1], cbgheightexpn);
 671             cblklengthexpn = int_min(tccp->cblk[2], cbglengthexpn);
 672
 673             res->bands = (opj_tcd_band_t *) opj_malloc(res->numbands * sizeof(
 674                              opj_tcd_band_t));
 675             for (bandno = 0; bandno < res->numbands; bandno++) {
 676                 int x0b, y0b, z0b;
 677                 int gain, numbps;
 678                 opj_stepsize_t *ss = NULL;
 679
 680                 tcd->band = &res->bands[bandno];
 681                 band = tcd->band;
 682
 683                 band->bandno = resno == 0 ? 0 : bandno + 1;
 684                 /* Bandno:  0 - LLL     2 - LHL
 685                             1 - HLL     3 - HHL
 686                             4 - LLH     6 - LHH
 687                             5 - HLH     7 - HHH     */
 688                 x0b = (band->bandno == 1) || (band->bandno == 3) || (band->bandno == 5) ||
 689                       (band->bandno == 7) ? 1 : 0;
 690                 y0b = (band->bandno == 2) || (band->bandno == 3) || (band->bandno == 6) ||
 691                       (band->bandno == 7) ? 1 : 0;
 692                 z0b = (band->bandno == 4) || (band->bandno == 5) || (band->bandno == 6) ||
 693                       (band->bandno == 7) ? 1 : 0;
 694
 695                 /* p. 65, B.15, ISO/IEC IS15444-1 : 2002  */
 696                 if (band->bandno == 0) {
 697                     /* band border (global) */
 698                     band->x0 = int_ceildivpow2(tilec->x0, levelnox);
 699                     band->y0 = int_ceildivpow2(tilec->y0, levelnoy);
 700                     band->z0 = int_ceildivpow2(tilec->z0, levelnoz);
 701                     band->x1 = int_ceildivpow2(tilec->x1, levelnox);
 702                     band->y1 = int_ceildivpow2(tilec->y1, levelnoy);
 703                     band->z1 = int_ceildivpow2(tilec->z1, levelnoz);
 704                 } else {
 705                     band->x0 = int_ceildivpow2(tilec->x0 - (1 << levelnox) * x0b, levelnox + 1);
 706                     band->y0 = int_ceildivpow2(tilec->y0 - (1 << levelnoy) * y0b, levelnoy + 1);
 707                     band->z0 = int_ceildivpow2(tilec->z0 - (1 << levelnoz) * z0b,
 708                                                (resno <= diff) ? levelnoz : levelnoz + 1);
 709                     band->x1 = int_ceildivpow2(tilec->x1 - (1 << levelnox) * x0b, levelnox + 1);
 710                     band->y1 = int_ceildivpow2(tilec->y1 - (1 << levelnoy) * y0b, levelnoy + 1);
 711                     band->z1 = int_ceildivpow2(tilec->z1 - (1 << levelnoz) * z0b,
 712                                                (resno <= diff) ? levelnoz : levelnoz + 1);
 713                 }
 714
 715                 ss = &tccp->stepsizes[(resno == 0) ? 0 : (prevnumbands + bandno + 1)];
 716                 if (bandno == (res->numbands - 1)) {
 717                     prevnumbands += (resno == 0) ? 0 : res->numbands;
 718                 }
 719                 gain = dwt_getgain(band->bandno, tccp->reversible);
 720                 numbps = volume->comps[compno].prec + gain;
 721
 722                 band->stepsize = (float)((1.0 + ss->mant / 2048.0) * pow(2.0,
 723                                          numbps - ss->expn));
 724                 band->numbps = ss->expn + tccp->numgbits - 1;   /* WHY -1 ? */
 725
 726                 for (precno = 0; precno < res->prctno[0] * res->prctno[1] * res->prctno[2];
 727                         precno++) {
 728                     int tlcblkxstart, tlcblkystart, tlcblkzstart, brcblkxend, brcblkyend,
 729                         brcblkzend;
 730
 731                     int cbgxstart = tlcbgxstart + (precno % res->prctno[0]) * (1 << cbgwidthexpn);
 732                     int cbgystart = tlcbgystart + ((precno / (res->prctno[0] * res->prctno[1])) /
 733                                                    res->prctno[0]) * (1 << cbgheightexpn);
 734                     int cbgzstart = tlcbgzstart + (precno / (res->prctno[0] * res->prctno[1])) *
 735                                     (1 << cbglengthexpn);
 736                     int cbgxend = cbgxstart + (1 << cbgwidthexpn);
 737                     int cbgyend = cbgystart + (1 << cbgheightexpn);
 738                     int cbgzend = cbgzstart + (1 << cbglengthexpn);
 739
 740                     /* opj_tcd_precinct_t *prc=&band->precincts[precno]; */
 741                     tcd->prc = &band->precincts[precno];
 742                     prc = tcd->prc;
 743
 744                     /* precinct size (global) */
 745                     prc->x0 = int_max(cbgxstart, band->x0);
 746                     prc->y0 = int_max(cbgystart, band->y0);
 747                     prc->z0 = int_max(cbgzstart, band->z0);
 748                     prc->x1 = int_min(cbgxend, band->x1);
 749                     prc->y1 = int_min(cbgyend, band->y1);
 750                     prc->z1 = int_min(cbgzend, band->z1);
 751
 752                     tlcblkxstart = int_floordivpow2(prc->x0, cblkwidthexpn) << cblkwidthexpn;
 753                     tlcblkystart = int_floordivpow2(prc->y0, cblkheightexpn) << cblkheightexpn;
 754                     tlcblkzstart = int_floordivpow2(prc->z0, cblklengthexpn) << cblklengthexpn;
 755                     brcblkxend = int_ceildivpow2(prc->x1, cblkwidthexpn) << cblkwidthexpn;
 756                     brcblkyend = int_ceildivpow2(prc->y1, cblkheightexpn) << cblkheightexpn;
 757                     brcblkzend = int_ceildivpow2(prc->z1, cblklengthexpn) << cblklengthexpn;
 758                     prc->cblkno[0] = (brcblkxend - tlcblkxstart) >> cblkwidthexpn;
 759                     prc->cblkno[1] = (brcblkyend - tlcblkystart) >> cblkheightexpn;
 760                     prc->cblkno[2] = (brcblkzend - tlcblkzstart) >> cblklengthexpn;
 761                     prc->cblkno[2] = (prc->cblkno[2] == 0) ? 1 : prc->cblkno[2];
 762
 763                     opj_free(prc->cblks);
 764                     prc->cblks = (opj_tcd_cblk_t *) opj_malloc((prc->cblkno[0] * prc->cblkno[1] *
 765                                  prc->cblkno[2]) * sizeof(opj_tcd_cblk_t));
 766                     prc->incltree = tgt_create(prc->cblkno[0], prc->cblkno[1], prc->cblkno[2]);
 767                     prc->imsbtree = tgt_create(prc->cblkno[0], prc->cblkno[1], prc->cblkno[2]);
 768
 769                     for (cblkno = 0; cblkno < (prc->cblkno[0] * prc->cblkno[1] * prc->cblkno[2]);
 770                             cblkno++) {
 771                         int cblkxstart = tlcblkxstart + (cblkno % prc->cblkno[0]) *
 772                                          (1 << cblkwidthexpn);
 773                         int cblkystart = tlcblkystart + ((cblkno % (prc->cblkno[0] * prc->cblkno[1])) /
 774                                                          prc->cblkno[0]) * (1 << cblkheightexpn);
 775                         int cblkzstart = tlcblkzstart + (cblkno / (prc->cblkno[0] * prc->cblkno[1])) *
 776                                          (1 << cblklengthexpn);
 777                         int cblkxend = cblkxstart + (1 << cblkwidthexpn);
 778                         int cblkyend = cblkystart + (1 << cblkheightexpn);
 779                         int cblkzend = cblkzstart + (1 << cblklengthexpn);
 780                         int prec = ((tilec->bpp > 16) ? 3 : ((tilec->bpp > 8) ? 2 : 1));
 781
 782                         tcd->cblk = &prc->cblks[cblkno];
 783                         cblk = tcd->cblk;
 784
 785                         /* code-block size (global) */
 786                         cblk->x0 = int_max(cblkxstart, prc->x0);
 787                         cblk->y0 = int_max(cblkystart, prc->y0);
 788                         cblk->z0 = int_max(cblkzstart, prc->z0);
 789                         cblk->x1 = int_min(cblkxend, prc->x1);
 790                         cblk->y1 = int_min(cblkyend, prc->y1);
 791                         cblk->z1 = int_min(cblkzend, prc->z1);
 792                     }
 793                 } /* precno */
 794             } /* bandno */
 795         } /* resno */
 796     } /* compno */
 797     /*tcd_dump(stdout, tcd, tcd->tcd_volume);*/
 798 }
 799
 800
 801 void tcd_free_encode(opj_tcd_t *tcd)
 802 {
 803     int tileno, compno, resno, bandno, precno;
 804
 805     opj_tcd_tile_t *tile = NULL;        /* pointer to tcd->tile     */
 806     /*  opj_tcd_slice_t *slice = NULL; */       /* pointer to tcd->slice */
 807     opj_tcd_tilecomp_t *tilec = NULL;   /* pointer to tcd->tilec    */
 808     opj_tcd_resolution_t *res = NULL;   /* pointer to tcd->res      */
 809     opj_tcd_band_t *band = NULL;        /* pointer to tcd->band     */
 810     opj_tcd_precinct_t *prc = NULL;     /* pointer to tcd->prc      */
 811
 812     for (tileno = 0; tileno < 1; tileno++) {
 813         tcd->tile = tcd->tcd_volume->tiles;
 814         tile = tcd->tile;
 815
 816         for (compno = 0; compno < tile->numcomps; compno++) {
 817             tcd->tilec = &tile->comps[compno];
 818             tilec = tcd->tilec;
 819
 820             for (resno = 0; resno < tilec->numresolution[0]; resno++) {
 821                 tcd->res = &tilec->resolutions[resno];
 822                 res = tcd->res;
 823
 824                 for (bandno = 0; bandno < res->numbands; bandno++) {
 825                     tcd->band = &res->bands[bandno];
 826                     band = tcd->band;
 827
 828                     for (precno = 0; precno < res->prctno[0] * res->prctno[1] * res->prctno[2];
 829                             precno++) {
 830                         tcd->prc = &band->precincts[precno];
 831                         prc = tcd->prc;
 832
 833                         if (prc->incltree != NULL) {
 834                             tgt_destroy(prc->incltree);
 835                             prc->incltree = NULL;
 836                         }
 837                         if (prc->imsbtree != NULL) {
 838                             tgt_destroy(prc->imsbtree);
 839                             prc->imsbtree = NULL;
 840                         }
 841                         opj_free(prc->cblks);
 842                         prc->cblks = NULL;
 843                     } /* for (precno */
 844                     opj_free(band->precincts);
 845                     band->precincts = NULL;
 846                 } /* for (bandno */
 847             } /* for (resno */
 848             opj_free(tilec->resolutions);
 849             tilec->resolutions = NULL;
 850         } /* for (compno */
 851         opj_free(tile->comps);
 852         tile->comps = NULL;
 853     } /* for (tileno */
 854     opj_free(tcd->tcd_volume->tiles);
 855     tcd->tcd_volume->tiles = NULL;
 856 }
 857
 858 /* ----------------------------------------------------------------------- */
 859 void tcd_malloc_decode(opj_tcd_t *tcd, opj_volume_t * volume, opj_cp_t * cp)
 860 {
 861     int tileno, compno, resno, bandno, precno, cblkno, res_max,
 862         i, j, p, q, r;
 863     unsigned int x0 = 0, y0 = 0, z0 = 0,
 864                  x1 = 0, y1 = 0, z1 = 0,
 865                  w, h, l;
 866
 867     tcd->volume = volume;
 868     tcd->cp = cp;
 869     tcd->tcd_volume->tw = cp->tw;
 870     tcd->tcd_volume->th = cp->th;
 871     tcd->tcd_volume->tl = cp->tl;
 872     tcd->tcd_volume->tiles = (opj_tcd_tile_t *) opj_malloc(cp->tw * cp->th * cp->tl
 873                              * sizeof(opj_tcd_tile_t));
 874
 875     for (i = 0; i < cp->tileno_size; i++) {
 876         opj_tcp_t *tcp = &(cp->tcps[cp->tileno[i]]);
 877         opj_tcd_tile_t *tile = &(tcd->tcd_volume->tiles[cp->tileno[i]]);
 878
 879         /* p61 ISO/IEC IS15444-1 : 2002 */
 880         /* curtileno --> raster scanned index of tiles */
 881         /* p,q,r --> matricial index of tiles */
 882         tileno = cp->tileno[i];
 883         p = tileno % cp->tw;
 884         q = tileno / cp->tw;
 885         r = tileno / (cp->tw * cp->th); /* extension to 3-D */
 886
 887         /* 4 borders of the tile rescale on the volume if necessary (B.3)*/
 888         tile->x0 = int_max(cp->tx0 + p * cp->tdx, volume->x0);
 889         tile->y0 = int_max(cp->ty0 + q * cp->tdy, volume->y0);
 890         tile->z0 = int_max(cp->tz0 + r * cp->tdz, volume->z0);
 891         tile->x1 = int_min(cp->tx0 + (p + 1) * cp->tdx, volume->x1);
 892         tile->y1 = int_min(cp->ty0 + (q + 1) * cp->tdy, volume->y1);
 893         tile->z1 = int_min(cp->tz0 + (r + 1) * cp->tdz, volume->z1);
 894         tile->numcomps = volume->numcomps;
 895
 896         tile->comps = (opj_tcd_tilecomp_t *) opj_malloc(volume->numcomps * sizeof(
 897                           opj_tcd_tilecomp_t));
 898         for (compno = 0; compno < tile->numcomps; compno++) {
 899             opj_tccp_t *tccp = &tcp->tccps[compno];
 900             opj_tcd_tilecomp_t *tilec = &tile->comps[compno];
 901             int prevnumbands = 0;
 902
 903             /* border of each tile component (global) */
 904             tilec->x0 = int_ceildiv(tile->x0, volume->comps[compno].dx);
 905             tilec->y0 = int_ceildiv(tile->y0, volume->comps[compno].dy);
 906             tilec->z0 = int_ceildiv(tile->z0, volume->comps[compno].dz);
 907             tilec->x1 = int_ceildiv(tile->x1, volume->comps[compno].dx);
 908             tilec->y1 = int_ceildiv(tile->y1, volume->comps[compno].dy);
 909             tilec->z1 = int_ceildiv(tile->z1, volume->comps[compno].dz);
 910
 911             tilec->data = (int *) opj_malloc((tilec->x1 - tilec->x0) *
 912                                              (tilec->y1 - tilec->y0) * (tilec->z1 - tilec->z0) * sizeof(int));
 913
 914             res_max = 0;
 915             for (i = 0; i < 3; i++) {
 916                 tilec->numresolution[i] = tccp->numresolution[i];
 917                 /*Greater of 3 resolutions contains all information*/
 918                 res_max = (tilec->numresolution[i] > res_max) ? tilec->numresolution[i] :
 919                           res_max;
 920             }
 921
 922             tilec->resolutions = (opj_tcd_resolution_t *) opj_malloc(res_max * sizeof(
 923                                      opj_tcd_resolution_t));
 924
 925             for (resno = 0; resno < res_max; resno++) {
 926                 opj_tcd_resolution_t *res = &tilec->resolutions[resno];
 927                 int pdx, pdy, pdz;
 928                 int tlprcxstart, tlprcystart, tlprczstart, brprcxend, brprcyend, brprczend;
 929                 int tlcbgxstart, tlcbgystart, tlcbgzstart, brcbgxend, brcbgyend, brcbgzend;
 930                 int cbgwidthexpn, cbgheightexpn, cbglengthexpn;
 931                 int cblkwidthexpn, cblkheightexpn, cblklengthexpn;
 932                 int levelnox = tilec->numresolution[0] - 1 - resno;
 933                 int levelnoy = tilec->numresolution[1] - 1 - resno;
 934                 int diff = tccp->numresolution[0] - tccp->numresolution[2];
 935                 int levelnoz = tilec->numresolution[2] - 1 - ((resno <= diff) ? 0 :
 936                                (resno - diff));
 937                 if (levelnoz < 0) {
 938                     levelnoz = 0;
 939                 }
 940
 941                 /* border for each resolution level (global) */
 942                 res->x0 = int_ceildivpow2(tilec->x0, levelnox);
 943                 res->y0 = int_ceildivpow2(tilec->y0, levelnoy);
 944                 res->z0 = int_ceildivpow2(tilec->z0, levelnoz);
 945                 res->x1 = int_ceildivpow2(tilec->x1, levelnox);
 946                 res->y1 = int_ceildivpow2(tilec->y1, levelnoy);
 947                 res->z1 = int_ceildivpow2(tilec->z1, levelnoz);
 948                 res->numbands = (resno == 0) ? 1 : (resno <= diff) ? 3 : 7; /* --> 3D */
 949
 950                 /* p. 30, table A-13, ISO/IEC IS154444-1 : 2002 */
 951                 if (tccp->csty & J3D_CCP_CSTY_PRT) {
 952                     pdx = tccp->prctsiz[0][resno];
 953                     pdy = tccp->prctsiz[1][resno];
 954                     pdz = tccp->prctsiz[2][resno];
 955                 } else {
 956                     pdx = 15;
 957                     pdy = 15;
 958                     pdz = 15;
 959                 }
 960
 961                 /* p. 66, B.16, ISO/IEC IS15444-1 : 2002  */
 962                 tlprcxstart = int_floordivpow2(res->x0, pdx) << pdx;
 963                 tlprcystart = int_floordivpow2(res->y0, pdy) << pdy;
 964                 tlprczstart = int_floordivpow2(res->z0, pdz) << pdz;
 965                 brprcxend = int_ceildivpow2(res->x1, pdx) << pdx;
 966                 brprcyend = int_ceildivpow2(res->y1, pdy) << pdy;
 967                 brprczend = int_ceildivpow2(res->z1, pdz) << pdz;
 968
 969                 res->prctno[0] = (brprcxend - tlprcxstart) >> pdx;
 970                 res->prctno[1] = (brprcyend - tlprcystart) >> pdy;
 971                 res->prctno[2] = (brprczend - tlprczstart) >> pdz;
 972
 973                 /* p. 67, B.17 & B.18, ISO/IEC IS15444-1 : 2002  */
 974                 if (resno == 0) {
 975                     tlcbgxstart = tlprcxstart;/*0*/
 976                     tlcbgystart = tlprcystart;
 977                     tlcbgzstart = tlprczstart;
 978                     brcbgxend = brprcxend;/*1*/
 979                     brcbgyend = brprcyend;
 980                     brcbgzend = brprczend;
 981                     cbgwidthexpn = pdx; /*15*/
 982                     cbgheightexpn = pdy;
 983                     cbglengthexpn = pdz;
 984                 } else {
 985                     tlcbgxstart = int_ceildivpow2(tlprcxstart, 1);
 986                     tlcbgystart = int_ceildivpow2(tlprcystart, 1);
 987                     tlcbgzstart = int_ceildivpow2(tlprczstart, 1);
 988                     brcbgxend = int_ceildivpow2(brprcxend, 1);
 989                     brcbgyend = int_ceildivpow2(brprcyend, 1);
 990                     brcbgzend = int_ceildivpow2(brprczend, 1);
 991                     cbgwidthexpn = pdx - 1;
 992                     cbgheightexpn = pdy - 1;
 993                     cbglengthexpn = pdz - 1;
 994                 }
 995
 996                 cblkwidthexpn = int_min(tccp->cblk[0], cbgwidthexpn); /*6*/
 997                 cblkheightexpn = int_min(tccp->cblk[1], cbgheightexpn); /*6*/
 998                 cblklengthexpn = int_min(tccp->cblk[2], cbglengthexpn); /*6*/
 999
1000                 res->bands = (opj_tcd_band_t *) opj_malloc(res->numbands * sizeof(
1001                                  opj_tcd_band_t));
1002                 for (bandno = 0; bandno < res->numbands; bandno++) {
1003                     int x0b, y0b, z0b;
1004                     int gain, numbps;
1005                     opj_stepsize_t *ss = NULL;
1006
1007                     opj_tcd_band_t *band = &res->bands[bandno];
1008                     band->bandno = resno == 0 ? 0 : bandno + 1;
1009                     /* Bandno:  0 - LLL     2 - LHL
1010                                 1 - HLL     3 - HHL
1011                                 4 - LLH     6 - LHH
1012                                 5 - HLH     7 - HHH     */
1013                     x0b = (band->bandno == 1) || (band->bandno == 3) || (band->bandno == 5) ||
1014                           (band->bandno == 7) ? 1 : 0;
1015                     y0b = (band->bandno == 2) || (band->bandno == 3) || (band->bandno == 6) ||
1016                           (band->bandno == 7) ? 1 : 0;
1017                     z0b = (band->bandno == 4) || (band->bandno == 5) || (band->bandno == 6) ||
1018                           (band->bandno == 7) ? 1 : 0;
1019
1020                     /* p. 65, B.15, ISO/IEC IS15444-1 : 2002  */
1021                     if (band->bandno == 0) {
1022                         /* band border (global) */
1023                         band->x0 = int_ceildivpow2(tilec->x0, levelnox);
1024                         band->y0 = int_ceildivpow2(tilec->y0, levelnoy);
1025                         band->z0 = int_ceildivpow2(tilec->z0, levelnoz);
1026                         band->x1 = int_ceildivpow2(tilec->x1, levelnox);
1027                         band->y1 = int_ceildivpow2(tilec->y1, levelnoy);
1028                         band->z1 = int_ceildivpow2(tilec->z1, levelnoz);
1029                     } else {
1030                         band->x0 = int_ceildivpow2(tilec->x0 - (1 << levelnox) * x0b, levelnox + 1);
1031                         band->y0 = int_ceildivpow2(tilec->y0 - (1 << levelnoy) * y0b, levelnoy + 1);
1032                         band->z0 = int_ceildivpow2(tilec->z0 - (1 << levelnoz) * z0b,
1033                                                    (resno <= diff) ? levelnoz : levelnoz + 1);
1034                         band->x1 = int_ceildivpow2(tilec->x1 - (1 << levelnox) * x0b, levelnox + 1);
1035                         band->y1 = int_ceildivpow2(tilec->y1 - (1 << levelnoy) * y0b, levelnoy + 1);
1036                         band->z1 = int_ceildivpow2(tilec->z1 - (1 << levelnoz) * z0b,
1037                                                    (resno <= diff) ? levelnoz : levelnoz + 1);
1038                     }
1039
1040                     ss = &tccp->stepsizes[(resno == 0) ? 0 : (prevnumbands + bandno + 1)];
1041                     if (bandno == (res->numbands - 1)) {
1042                         prevnumbands += (resno == 0) ? 0 : res->numbands;
1043                     }
1044                     gain = dwt_getgain(band->bandno, tccp->reversible);
1045                     numbps = volume->comps[compno].prec + gain;
1046
1047                     band->stepsize = (float)((1.0 + ss->mant / 2048.0) * pow(2.0,
1048                                              numbps - ss->expn));
1049                     band->numbps = ss->expn + tccp->numgbits - 1;   /* WHY -1 ? */
1050
1051                     band->precincts = (opj_tcd_precinct_t *) opj_malloc(res->prctno[0] *
1052                                       res->prctno[1] * res->prctno[2] * sizeof(opj_tcd_precinct_t));
1053
1054                     for (precno = 0; precno < res->prctno[0] * res->prctno[1] * res->prctno[2];
1055                             precno++) {
1056                         int tlcblkxstart, tlcblkystart, tlcblkzstart, brcblkxend, brcblkyend,
1057                             brcblkzend;
1058
1059                         int cbgxstart = tlcbgxstart + (precno % res->prctno[0]) * (1 << cbgwidthexpn);
1060                         int cbgystart = tlcbgystart + (precno / res->prctno[0]) * (1 << cbgheightexpn);
1061                         int cbgzstart = tlcbgzstart + (precno / (res->prctno[0] * res->prctno[1])) *
1062                                         (1 << cbglengthexpn);
1063                         int cbgxend = cbgxstart + (1 << cbgwidthexpn);
1064                         int cbgyend = cbgystart + (1 << cbgheightexpn);
1065                         int cbgzend = cbgzstart + (1 << cbglengthexpn);
1066
1067                         opj_tcd_precinct_t *prc = &band->precincts[precno];
1068                         /* precinct size (global) */
1069                         prc->x0 = int_max(cbgxstart, band->x0);
1070                         prc->y0 = int_max(cbgystart, band->y0);
1071                         prc->z0 = int_max(cbgzstart, band->z0);
1072                         prc->x1 = int_min(cbgxend, band->x1);
1073                         prc->y1 = int_min(cbgyend, band->y1);
1074                         prc->z1 = int_min(cbgzend, band->z1);
1075
1076                         tlcblkxstart = int_floordivpow2(prc->x0, cblkwidthexpn) << cblkwidthexpn;
1077                         tlcblkystart = int_floordivpow2(prc->y0, cblkheightexpn) << cblkheightexpn;
1078                         tlcblkzstart = int_floordivpow2(prc->z0, cblklengthexpn) << cblklengthexpn;
1079                         brcblkxend = int_ceildivpow2(prc->x1, cblkwidthexpn) << cblkwidthexpn;
1080                         brcblkyend = int_ceildivpow2(prc->y1, cblkheightexpn) << cblkheightexpn;
1081                         brcblkzend = int_ceildivpow2(prc->z1, cblklengthexpn) << cblklengthexpn;
1082                         prc->cblkno[0] = (brcblkxend - tlcblkxstart) >> cblkwidthexpn;
1083                         prc->cblkno[1] = (brcblkyend - tlcblkystart) >> cblkheightexpn;
1084                         prc->cblkno[2] = (brcblkzend - tlcblkzstart) >> cblklengthexpn;
1085                         prc->cblkno[2] = (prc->cblkno[2] == 0) ? 1 : prc->cblkno[2];
1086
1087                         prc->cblks = (opj_tcd_cblk_t *) opj_malloc((prc->cblkno[0] * prc->cblkno[1] *
1088                                      prc->cblkno[2]) * sizeof(opj_tcd_cblk_t));
1089                         prc->incltree = tgt_create(prc->cblkno[0], prc->cblkno[1], prc->cblkno[2]);
1090                         prc->imsbtree = tgt_create(prc->cblkno[0], prc->cblkno[1], prc->cblkno[2]);
1091
1092                         for (cblkno = 0; cblkno < prc->cblkno[0] * prc->cblkno[1] * prc->cblkno[2];
1093                                 cblkno++) {
1094                             int cblkxstart = tlcblkxstart + (cblkno % prc->cblkno[0]) *
1095                                              (1 << cblkwidthexpn);
1096                             int cblkystart = tlcblkystart + ((cblkno % (prc->cblkno[0] * prc->cblkno[1])) /
1097                                                              prc->cblkno[0]) * (1 << cblkheightexpn);
1098                             int cblkzstart = tlcblkzstart + (cblkno / (prc->cblkno[0] * prc->cblkno[1])) *
1099                                              (1 << cblklengthexpn);
1100                             int cblkxend = cblkxstart + (1 << cblkwidthexpn);
1101                             int cblkyend = cblkystart + (1 << cblkheightexpn);
1102                             int cblkzend = cblkzstart + (1 << cblklengthexpn);
1103                             int prec = ((tilec->bpp > 16) ? 3 : ((tilec->bpp > 8) ? 2 : 1));
1104                             /* code-block size (global) */
1105                             opj_tcd_cblk_t *cblk = &prc->cblks[cblkno];
1106
1107                             /* code-block size (global) */
1108                             cblk->x0 = int_max(cblkxstart, prc->x0);
1109                             cblk->y0 = int_max(cblkystart, prc->y0);
1110                             cblk->z0 = int_max(cblkzstart, prc->z0);
1111                             cblk->x1 = int_min(cblkxend, prc->x1);
1112                             cblk->y1 = int_min(cblkyend, prc->y1);
1113                             cblk->z1 = int_min(cblkzend, prc->z1);
1114                         }
1115                     } /* precno */
1116                 } /* bandno */
1117             } /* resno */
1118         } /* compno */
1119     } /* i = 0..cp->tileno_size */
1120
1121     /*tcd_dump(stdout, tcd, tcd->tcd_volume);*/
1122
1123     /*
1124     Allocate place to store the decoded data = final volume
1125     Place limited by the tile really present in the codestream
1126     */
1127
1128     for (i = 0; i < volume->numcomps; i++) {
1129         for (j = 0; j < cp->tileno_size; j++) {
1130             tileno = cp->tileno[j];
1131             x0 = (j == 0) ? tcd->tcd_volume->tiles[tileno].comps[i].x0 : int_min(x0,
1132                     (unsigned int) tcd->tcd_volume->tiles[tileno].comps[i].x0);
1133             y0 = (j == 0) ? tcd->tcd_volume->tiles[tileno].comps[i].y0 : int_min(y0,
1134                     (unsigned int) tcd->tcd_volume->tiles[tileno].comps[i].y0);
1135             z0 = (j == 0) ? tcd->tcd_volume->tiles[tileno].comps[i].z0 : int_min(z0,
1136                     (unsigned int) tcd->tcd_volume->tiles[tileno].comps[i].z0);
1137             x1 = (j == 0) ? tcd->tcd_volume->tiles[tileno].comps[i].x1 : int_max(x1,
1138                     (unsigned int) tcd->tcd_volume->tiles[tileno].comps[i].x1);
1139             y1 = (j == 0) ? tcd->tcd_volume->tiles[tileno].comps[i].y1 : int_max(y1,
1140                     (unsigned int) tcd->tcd_volume->tiles[tileno].comps[i].y1);
1141             z1 = (j == 0) ? tcd->tcd_volume->tiles[tileno].comps[i].z1 : int_max(z1,
1142                     (unsigned int) tcd->tcd_volume->tiles[tileno].comps[i].z1);
1143         }
1144
1145         w = x1 - x0;
1146         h = y1 - y0;
1147         l = z1 - z0;
1148
1149         volume->comps[i].data = (int *) opj_malloc(w * h * l * sizeof(int));
1150         volume->comps[i].w = w;
1151         volume->comps[i].h = h;
1152         volume->comps[i].l = l;
1153         volume->comps[i].x0 = x0;
1154         volume->comps[i].y0 = y0;
1155         volume->comps[i].z0 = z0;
1156         volume->comps[i].bigendian = cp->bigendian;
1157     }
1158 }
1159
1160 void tcd_free_decode(opj_tcd_t *tcd)
1161 {
1162     int tileno, compno, resno, bandno, precno;
1163
1164     opj_tcd_volume_t *tcd_volume = tcd->tcd_volume;
1165
1166     for (tileno = 0; tileno < tcd_volume->tw * tcd_volume->th * tcd_volume->tl;
1167             tileno++) {
1168         opj_tcd_tile_t *tile = &tcd_volume->tiles[tileno];
1169         for (compno = 0; compno < tile->numcomps; compno++) {
1170             opj_tcd_tilecomp_t *tilec = &tile->comps[compno];
1171             for (resno = 0; resno < tilec->numresolution[0]; resno++) {
1172                 opj_tcd_resolution_t *res = &tilec->resolutions[resno];
1173                 for (bandno = 0; bandno < res->numbands; bandno++) {
1174                     opj_tcd_band_t *band = &res->bands[bandno];
1175                     for (precno = 0; precno < res->prctno[1] * res->prctno[0] * res->prctno[2];
1176                             precno++) {
1177                         opj_tcd_precinct_t *prec = &band->precincts[precno];
1178                         if (prec->cblks != NULL) {
1179                             opj_free(prec->cblks);
1180                         }
1181                         if (prec->imsbtree != NULL) {
1182                             tgt_destroy(prec->imsbtree);
1183                         }
1184                         if (prec->incltree != NULL) {
1185                             tgt_destroy(prec->incltree);
1186                         }
1187                         /*for (treeno = 0; treeno < prec->numtrees; treeno++){
1188                             if (prec->imsbtree[treeno] != NULL) tgt_destroy(prec->imsbtree[treeno]);
1189                             if (prec->incltree[treeno] != NULL) tgt_destroy(prec->incltree[treeno]);
1190                         }*/
1191                     }
1192                     if (band->precincts != NULL) {
1193                         opj_free(band->precincts);
1194                     }
1195                 }
1196             }
1197             if (tilec->resolutions != NULL) {
1198                 opj_free(tilec->resolutions);
1199             }
1200         }
1201         if (tile->comps != NULL) {
1202             opj_free(tile->comps);
1203         }
1204     }
1205
1206     if (tcd_volume->tiles != NULL) {
1207         opj_free(tcd_volume->tiles);
1208     }
1209 }
1210
1211
1212
1213 /* ----------------------------------------------------------------------- */
1214 void tcd_makelayer_fixed(opj_tcd_t *tcd, int layno, int final)
1215 {
1216     int compno, resno, bandno, precno, cblkno;
1217     int value;          /*, matrice[tcd_tcp->numlayers][tcd_tile->comps[0].numresolution[0]][3]; */
1218     int matrice[10][10][3];
1219     int i, j, k;
1220
1221     opj_cp_t *cp = tcd->cp;
1222     opj_tcd_tile_t *tcd_tile = tcd->tcd_tile;
1223     opj_tcp_t *tcd_tcp = tcd->tcp;
1224
1225     /*matrice=(int*)opj_malloc(tcd_tcp->numlayers*tcd_tile->comps[0].numresolution[0]*3*sizeof(int)); */
1226
1227     for (compno = 0; compno < tcd_tile->numcomps; compno++) {
1228         opj_tcd_tilecomp_t *tilec = &tcd_tile->comps[compno];
1229         for (i = 0; i < tcd_tcp->numlayers; i++) {
1230             for (j = 0; j < tilec->numresolution[0]; j++) {
1231                 for (k = 0; k < 3; k++) {
1232                     matrice[i][j][k] =
1233                         (int)(cp->matrice[i * tilec->numresolution[0] * 3 + j * 3 + k]
1234                               * (float)(tcd->volume->comps[compno].prec / 16.0));
1235                 }
1236             }
1237         }
1238
1239         for (resno = 0; resno < tilec->numresolution[0]; resno++) {
1240             opj_tcd_resolution_t *res = &tilec->resolutions[resno];
1241             for (bandno = 0; bandno < res->numbands; bandno++) {
1242                 opj_tcd_band_t *band = &res->bands[bandno];
1243                 for (precno = 0; precno < res->prctno[0] * res->prctno[1] * res->prctno[2];
1244                         precno++) {
1245                     opj_tcd_precinct_t *prc = &band->precincts[precno];
1246                     for (cblkno = 0; cblkno < prc->cblkno[0] * prc->cblkno[1] * prc->cblkno[2];
1247                             cblkno++) {
1248                         opj_tcd_cblk_t *cblk = &prc->cblks[cblkno];
1249                         opj_tcd_layer_t *layer = &cblk->layers[layno];
1250                         int n;
1251                         int imsb = tcd->volume->comps[compno].prec -
1252                                    cblk->numbps;  /* number of bit-plan equal to zero */
1253                         /* Correction of the matrix of coefficient to include the IMSB information */
1254                         if (layno == 0) {
1255                             value = matrice[layno][resno][bandno];
1256                             if (imsb >= value) {
1257                                 value = 0;
1258                             } else {
1259                                 value -= imsb;
1260                             }
1261                         } else {
1262                             value = matrice[layno][resno][bandno] - matrice[layno - 1][resno][bandno];
1263                             if (imsb >= matrice[layno - 1][resno][bandno]) {
1264                                 value -= (imsb - matrice[layno - 1][resno][bandno]);
1265                                 if (value < 0) {
1266                                     value = 0;
1267                                 }
1268                             }
1269                         }
1270
1271                         if (layno == 0) {
1272                             cblk->numpassesinlayers = 0;
1273                         }
1274
1275                         n = cblk->numpassesinlayers;
1276                         if (cblk->numpassesinlayers == 0) {
1277                             if (value != 0) {
1278                                 n = 3 * value - 2 + cblk->numpassesinlayers;
1279                             } else {
1280                                 n = cblk->numpassesinlayers;
1281                             }
1282                         } else {
1283                             n = 3 * value + cblk->numpassesinlayers;
1284                         }
1285
1286                         layer->numpasses = n - cblk->numpassesinlayers;
1287
1288                         if (!layer->numpasses) {
1289                             continue;
1290                         }
1291
1292                         if (cblk->numpassesinlayers == 0) {
1293                             layer->len = cblk->passes[n - 1].rate;
1294                             layer->data = cblk->data;
1295                         } else {
1296                             layer->len = cblk->passes[n - 1].rate - cblk->passes[cblk->numpassesinlayers -
1297                                          1].rate;
1298                             layer->data = cblk->data + cblk->passes[cblk->numpassesinlayers - 1].rate;
1299                         }
1300                         if (final) {
1301                             cblk->numpassesinlayers = n;
1302                         }
1303                     }
1304                 }
1305             }
1306         }
1307     }
1308 }
1309
1310 void tcd_rateallocate_fixed(opj_tcd_t *tcd)
1311 {
1312     int layno;
1313     for (layno = 0; layno < tcd->tcp->numlayers; layno++) {
1314         tcd_makelayer_fixed(tcd, layno, 1);
1315     }
1316 }
1317
1318 void tcd_makelayer(opj_tcd_t *tcd, int layno, double thresh, int final)
1319 {
1320     int compno, resno, bandno, precno, cblkno, passno;
1321
1322     opj_tcd_tile_t *tcd_tile = tcd->tcd_tile;
1323
1324     tcd_tile->distolayer[layno] = 0;    /* fixed_quality */
1325
1326     for (compno = 0; compno < tcd_tile->numcomps; compno++) {
1327         opj_tcd_tilecomp_t *tilec = &tcd_tile->comps[compno];
1328         for (resno = 0; resno < tilec->numresolution[0]; resno++) {
1329             opj_tcd_resolution_t *res = &tilec->resolutions[resno];
1330             for (bandno = 0; bandno < res->numbands; bandno++) {
1331                 opj_tcd_band_t *band = &res->bands[bandno];
1332                 for (precno = 0; precno < res->prctno[0] * res->prctno[1] * res->prctno[2];
1333                         precno++) {
1334                     opj_tcd_precinct_t *prc = &band->precincts[precno];
1335                     for (cblkno = 0; cblkno < prc->cblkno[0] * prc->cblkno[1] * prc->cblkno[2];
1336                             cblkno++) {
1337                         opj_tcd_cblk_t *cblk = &prc->cblks[cblkno];
1338                         opj_tcd_layer_t *layer = &cblk->layers[layno];
1339
1340                         int n;
1341                         if (layno == 0) {
1342                             cblk->numpassesinlayers = 0;
1343                         }
1344                         n = cblk->numpassesinlayers;
1345                         for (passno = cblk->numpassesinlayers; passno < cblk->totalpasses; passno++) {
1346                             int dr;
1347                             double dd;
1348                             opj_tcd_pass_t *pass = &cblk->passes[passno];
1349                             if (n == 0) {
1350                                 dr = pass->rate;
1351                                 dd = pass->distortiondec;
1352                             } else {
1353                                 dr = pass->rate - cblk->passes[n - 1].rate;
1354                                 dd = pass->distortiondec - cblk->passes[n - 1].distortiondec;
1355                             }
1356                             if (!dr) {
1357                                 if (dd) {
1358                                     n = passno + 1;
1359                                 }
1360                                 continue;
1361                             }
1362                             if (dd / dr >= thresh) {
1363                                 n = passno + 1;
1364                             }
1365                         }
1366                         layer->numpasses = n - cblk->numpassesinlayers;
1367
1368                         if (!layer->numpasses) {
1369                             layer->disto = 0;
1370                             continue;
1371                         }
1372                         if (cblk->numpassesinlayers == 0) {
1373                             layer->len = cblk->passes[n - 1].rate;
1374                             layer->data = cblk->data;
1375                             layer->disto = cblk->passes[n - 1].distortiondec;
1376                         } else {
1377                             layer->len = cblk->passes[n - 1].rate - cblk->passes[cblk->numpassesinlayers -
1378                                          1].rate;
1379                             layer->data = cblk->data + cblk->passes[cblk->numpassesinlayers - 1].rate;
1380                             layer->disto = cblk->passes[n - 1].distortiondec -
1381                                            cblk->passes[cblk->numpassesinlayers - 1].distortiondec;
1382                         }
1383
1384                         tcd_tile->distolayer[layno] += layer->disto;    /* fixed_quality */
1385
1386                         if (final) {
1387                             cblk->numpassesinlayers = n;
1388                         }
1389
1390                         /*  fprintf(stdout,"MakeLayer : %d %f %d %d \n",layer->len, layer->disto, layer->numpasses, n);*/
1391                     }
1392                 }
1393             }
1394         }
1395     }
1396 }
1397
1398 bool tcd_rateallocate(opj_tcd_t *tcd, unsigned char *dest, int len,
1399                       opj_volume_info_t * volume_info)
1400 {
1401     int compno, resno, bandno, precno, cblkno, passno, layno;
1402     double min, max;
1403     double cumdisto[100];   /* fixed_quality */
1404     const double K = 1;     /* 1.1; // fixed_quality */
1405     double maxSE = 0;
1406
1407     opj_cp_t *cp = tcd->cp;
1408     opj_tcd_tile_t *tcd_tile = tcd->tcd_tile;
1409     opj_tcp_t *tcd_tcp = tcd->tcp;
1410
1411     min = DBL_MAX;
1412     max = 0;
1413
1414     tcd_tile->nbpix = 0;        /* fixed_quality */
1415
1416     for (compno = 0; compno < tcd_tile->numcomps; compno++) {
1417         opj_tcd_tilecomp_t *tilec = &tcd_tile->comps[compno];
1418         tilec->nbpix = 0;
1419         for (resno = 0; resno < tilec->numresolution[0]; resno++) {
1420             opj_tcd_resolution_t *res = &tilec->resolutions[resno];
1421             for (bandno = 0; bandno < res->numbands; bandno++) {
1422                 opj_tcd_band_t *band = &res->bands[bandno];
1423                 for (precno = 0; precno < res->prctno[0] * res->prctno[1] * res->prctno[2];
1424                         precno++) {
1425                     opj_tcd_precinct_t *prc = &band->precincts[precno];
1426                     for (cblkno = 0; cblkno < prc->cblkno[0] * prc->cblkno[1] * prc->cblkno[2];
1427                             cblkno++) {
1428                         opj_tcd_cblk_t *cblk = &prc->cblks[cblkno];
1429                         for (passno = 0; passno < cblk->totalpasses; passno++) {
1430                             opj_tcd_pass_t *pass = &cblk->passes[passno];
1431                             int dr;
1432                             double dd, rdslope;
1433                             if (passno == 0) {
1434                                 dr = pass->rate;
1435                                 dd = pass->distortiondec;
1436                             } else {
1437                                 dr = pass->rate - cblk->passes[passno - 1].rate;
1438                                 dd = pass->distortiondec - cblk->passes[passno - 1].distortiondec;
1439                             }
1440                             if (dr == 0) {
1441                                 continue;
1442                             }
1443                             rdslope = dd / dr;
1444                             if (rdslope < min) {
1445                                 min = rdslope;
1446                             }
1447                             if (rdslope > max) {
1448                                 max = rdslope;
1449                             }
1450
1451                         } /* passno */
1452
1453                         /* fixed_quality */
1454                         tcd_tile->nbpix += ((cblk->x1 - cblk->x0) * (cblk->y1 - cblk->y0) *
1455                                             (cblk->z1 - cblk->z0));
1456                         tilec->nbpix += ((cblk->x1 - cblk->x0) * (cblk->y1 - cblk->y0) *
1457                                          (cblk->z1 - cblk->z0));
1458                     } /* cbklno */
1459                 } /* precno */
1460             } /* bandno */
1461         } /* resno */
1462
1463         maxSE += (((double)(1 << tcd->volume->comps[compno].prec) - 1.0)
1464                   * ((double)(1 << tcd->volume->comps[compno].prec) - 1.0))
1465                  * ((double)(tilec->nbpix));
1466     } /* compno */
1467
1468     /* add antonin index */
1469     if (volume_info && volume_info->index_on) {
1470         opj_tile_info_t *info_TL = &volume_info->tile[tcd->tcd_tileno];
1471         info_TL->nbpix = tcd_tile->nbpix;
1472         info_TL->distotile = tcd_tile->distotile;
1473         info_TL->thresh = (double *) opj_malloc(tcd_tcp->numlayers * sizeof(double));
1474     }
1475     /* dda */
1476
1477     for (layno = 0; layno < tcd_tcp->numlayers; layno++) {
1478         double lo = min;
1479         double hi = max;
1480         int success = 0;
1481         int maxlen = tcd_tcp->rates[layno] ? int_min(((int) tcd_tcp->rates[layno]),
1482                      len) : len;
1483         double goodthresh;
1484         double distotarget;     /* fixed_quality */
1485         int i = 0;
1486
1487         /* fixed_quality */
1488         distotarget = tcd_tile->distotile - ((K * maxSE) / pow((float)10,
1489                                              tcd_tcp->distoratio[layno] / 10));
1490
1491         if ((tcd_tcp->rates[layno]) || (cp->disto_alloc == 0)) {
1492             opj_t2_t *t2 = t2_create(tcd->cinfo, tcd->volume, cp);
1493             int oldl = 0, oldoldl = 0;
1494             for (i = 0; i < 128; i++) {
1495                 double thresh = (lo + hi) / 2;
1496                 int l = 0;
1497                 double distoachieved = 0;   /* fixed_quality -q */
1498
1499                 tcd_makelayer(tcd, layno, thresh, 0);
1500
1501                 if (cp->fixed_quality) {    /* fixed_quality -q */
1502                     distoachieved = (layno == 0) ? tcd_tile->distolayer[0] : cumdisto[layno - 1] +
1503                                     tcd_tile->distolayer[layno];
1504                     if (distoachieved < distotarget) {
1505                         hi = thresh;
1506                         continue;
1507                     }
1508                     lo = thresh;
1509                 } else {        /* disto_alloc -r, fixed_alloc -f */
1510                     l = t2_encode_packets(t2, tcd->tcd_tileno, tcd_tile, layno + 1, dest, maxlen,
1511                                           volume_info);
1512                     /*fprintf(stdout, "layno %d i %d len=%d max=%d \n",layno,i,l,maxlen);*/
1513                     if (l == -999) {
1514                         lo = thresh;
1515                         continue;
1516                     } else if (l == oldl && oldl == oldoldl && tcd_tile->distolayer[layno] > 0.0 &&
1517                                i > 32) {
1518                         break;
1519                     }
1520                     hi = thresh;
1521                     oldoldl = oldl;
1522                     oldl = l;
1523                 }
1524                 success = 1;
1525                 goodthresh = thresh;
1526             }
1527             t2_destroy(t2);
1528         } else {
1529             success = 1;
1530             goodthresh = min;
1531         }
1532         if (!success) {
1533             return false;
1534         }
1535
1536         if (volume_info && volume_info->index_on) { /* Threshold for Marcela Index */
1537             volume_info->tile[tcd->tcd_tileno].thresh[layno] = goodthresh;
1538         }
1539         tcd_makelayer(tcd, layno, goodthresh, 1);
1540
1541         /* fixed_quality */
1542         cumdisto[layno] = (layno == 0) ? tcd_tile->distolayer[0] : cumdisto[layno - 1] +
1543                           tcd_tile->distolayer[layno];
1544     }
1545
1546     return true;
1547 }
1548
1549 /* ----------------------------------------------------------------------- */
1550 int tcd_encode_tile(opj_tcd_t *tcd, int tileno, unsigned char *dest, int len,
1551                     opj_volume_info_t * volume_info)
1552 {
1553     int compno;
1554     int l = 0, i, npck = 0;
1555     double encoding_time;
1556
1557     opj_tcd_tile_t  *tile = NULL;
1558     opj_tcp_t       *tcd_tcp = NULL;
1559     opj_cp_t        *cp = NULL;
1560
1561     opj_tcp_t       *tcp = &tcd->cp->tcps[0];
1562     opj_tccp_t      *tccp = &tcp->tccps[0];
1563     opj_volume_t    *volume = tcd->volume;
1564     opj_t2_t        *t2 = NULL;     /* T2 component */
1565
1566     tcd->tcd_tileno = tileno;           /* current encoded/decoded tile */
1567
1568     tcd->tcd_tile = tcd->tcd_volume->tiles; /* tile information */
1569     tile = tcd->tcd_tile;
1570
1571     tcd->tcp = &tcd->cp->tcps[tileno];  /* coding/decoding params of tileno */
1572     tcd_tcp = tcd->tcp;
1573
1574     cp = tcd->cp;       /* coding parameters */
1575
1576     /* INDEX >> */
1577     if (volume_info && volume_info->index_on) {
1578         opj_tcd_tilecomp_t *tilec_idx = &tile->comps[0];    /* based on component 0 */
1579         for (i = 0; i < tilec_idx->numresolution[0]; i++) {
1580             opj_tcd_resolution_t *res_idx = &tilec_idx->resolutions[i];
1581
1582             volume_info->tile[tileno].prctno[0][i] = res_idx->prctno[0];
1583             volume_info->tile[tileno].prctno[1][i] = res_idx->prctno[1];
1584             volume_info->tile[tileno].prctno[2][i] = res_idx->prctno[2];
1585
1586             npck += res_idx->prctno[0] * res_idx->prctno[1] * res_idx->prctno[2];
1587
1588             volume_info->tile[tileno].prctsiz[0][i] = tccp->prctsiz[0][i];
1589             volume_info->tile[tileno].prctsiz[1][i] = tccp->prctsiz[1][i];
1590             volume_info->tile[tileno].prctsiz[2][i] = tccp->prctsiz[2][i];
1591         }
1592         volume_info->tile[tileno].packet = (opj_packet_info_t *) opj_malloc(
1593                                                volume_info->comp * volume_info->layer * npck * sizeof(opj_packet_info_t));
1594     }
1595     /* << INDEX */
1596
1597     /*---------------TILE-------------------*/
1598     encoding_time = opj_clock();    /* time needed to encode a tile */
1599
1600     for (compno = 0; compno < tile->numcomps; compno++) {
1601         int x, y, z;
1602         opj_tcd_tilecomp_t *tilec = &tile->comps[compno];
1603
1604         int adjust;
1605         int offset_x = int_ceildiv(volume->x0,
1606                                    volume->comps[compno].dx); /*ceil(x0 / subsampling_dx)*/
1607         int offset_y = int_ceildiv(volume->y0, volume->comps[compno].dy);
1608         int offset_z = int_ceildiv(volume->z0, volume->comps[compno].dz);
1609
1610         int tw = tilec->x1 - tilec->x0;
1611         int w = int_ceildiv(volume->x1 - volume->x0, volume->comps[compno].dx);
1612         int th = tilec->y1 - tilec->y0;
1613         int h = int_ceildiv(volume->y1 - volume->y0, volume->comps[compno].dy);
1614         int tl = tilec->z1 - tilec->z0;
1615         int l = int_ceildiv(volume->z1 - volume->z0, volume->comps[compno].dz);
1616
1617
1618
1619         /* extract tile data from volume.comps[0].data to tile.comps[0].data */
1620         /*fprintf(stdout,"[INFO] Extract tile data\n");*/
1621         if (tcd->cp->transform_format == TRF_3D_RLS ||
1622                 tcd->cp->transform_format == TRF_3D_LSE) {
1623             adjust = 0;
1624         } else {
1625             adjust = volume->comps[compno].sgnd ? 0 : 1 << (volume->comps[compno].prec -
1626                      1); /*sign=='+' --> 2^(prec-1)*/
1627             if (volume->comps[compno].dcoffset != 0) {
1628                 adjust += volume->comps[compno].dcoffset;
1629                 fprintf(stdout, "[INFO] DC Offset applied: DCO = %d -> adjust = %d\n",
1630                         volume->comps[compno].dcoffset, adjust);
1631             }
1632         }
1633
1634         if (tcd_tcp->tccps[compno].reversible ==
1635                 1) { /*IF perfect reconstruction (DWT.5-3)*/
1636             for (z = tilec->z0; z < tilec->z1; z++) {
1637                 for (y = tilec->y0; y < tilec->y1; y++) {
1638                     /* start of the src tile scanline */
1639                     int *data = &volume->comps[compno].data[(tilec->x0 - offset_x) +
1640                                                             (y - offset_y) * w + (z - offset_z) * w * h];
1641                     /* start of the dst tile scanline */
1642                     int *tile_data = &tilec->data[(y - tilec->y0) * tw + (z - tilec->z0) * tw * th];
1643                     for (x = tilec->x0; x < tilec->x1; x++) {
1644                         *tile_data++ = *data++ - adjust;
1645                     }
1646                 }
1647             }
1648         } else if (tcd_tcp->tccps[compno].reversible == 0) { /*IF not (DWT.9-7)*/
1649             for (z = tilec->z0; z < tilec->z1; z++) {
1650                 for (y = tilec->y0; y < tilec->y1; y++) {
1651                     /* start of the src tile scanline */
1652                     int *data = &volume->comps[compno].data[(tilec->x0 - offset_x) +
1653                                                             (y - offset_y) * w + (z - offset_z) * w * h];
1654                     /* start of the dst tile scanline */
1655                     int *tile_data = &tilec->data[(y - tilec->y0) * tw + (z - tilec->z0) * tw * th];
1656                     for (x = tilec->x0; x < tilec->x1; x++) {
1657                         *tile_data++ = (*data++ - adjust) << 13;
1658                     }
1659                 }
1660             }
1661         }
1662
1663     }
1664
1665     /*----------------MCT-------------------*/
1666     if (tcd_tcp->mct) {
1667         int samples = (tile->comps[0].x1 - tile->comps[0].x0) *
1668                       (tile->comps[0].y1 - tile->comps[0].y0) * (tile->comps[0].z1 -
1669                               tile->comps[0].z0);
1670         fprintf(stdout, "[INFO] Tcd_encode_tile: mct\n");
1671         if (tcd_tcp->tccps[0].reversible == 0) {
1672             mct_encode_real(tile->comps[0].data, tile->comps[1].data, tile->comps[2].data,
1673                             samples);
1674         } else {
1675             mct_encode(tile->comps[0].data, tile->comps[1].data, tile->comps[2].data,
1676                        samples);
1677         }
1678     }
1679     /*----------------TRANSFORM---------------------------------*/
1680     fprintf(stdout, "[INFO] Tcd_encode_tile: Transform\n");
1681     for (compno = 0; compno < tile->numcomps; compno++) {
1682         opj_tcd_tilecomp_t *tilec = &tile->comps[compno];
1683         dwt_encode(tilec, tcd_tcp->tccps[compno].dwtid);
1684     }
1685
1686     /*-------------------ENTROPY CODING-----------------------------*/
1687     fprintf(stdout, "[INFO] Tcd_encode_tile: Entropy coding\n");
1688     if ((cp->encoding_format == ENCOD_2EB) || (cp->encoding_format == ENCOD_3EB)) {
1689         if (cp->encoding_format == ENCOD_2EB) {
1690             opj_t1_t *t1 = NULL;
1691             t1 = t1_create(tcd->cinfo);
1692             t1_encode_cblks(t1, tile, tcd_tcp);
1693             t1_destroy(t1);
1694         } else if (cp->encoding_format == ENCOD_3EB) {
1695             opj_t1_3d_t *t1 = NULL;
1696             t1 = t1_3d_create(tcd->cinfo);
1697             t1_3d_encode_cblks(t1, tile, tcd_tcp);
1698             t1_3d_destroy(t1);
1699         }
1700         /*-----------RATE-ALLOCATE------------------*/
1701         /* INDEX */
1702         if (volume_info) {
1703             volume_info->index_write = 0;
1704         }
1705         if (cp->disto_alloc || cp->fixed_quality) {
1706             fprintf(stdout, "[INFO] Tcd_encode_tile: Rate-allocate\n");
1707             tcd_rateallocate(tcd, dest, len,
1708                              volume_info);          /* Normal Rate/distortion allocation */
1709         } else {/* fixed_alloc */
1710             fprintf(stdout, "[INFO] Tcd_encode_tile: Rate-allocate fixed\n");
1711             tcd_rateallocate_fixed(
1712                 tcd);                            /* Fixed layer allocation */
1713         }
1714
1715         /*--------------TIER2------------------*/
1716         /* INDEX */
1717         if (volume_info) {
1718             volume_info->index_write = 1;
1719         }
1720         fprintf(stdout, "[INFO] Tcd_encode_tile: Tier - 2\n");
1721         t2 = t2_create(tcd->cinfo, volume, cp);
1722         l = t2_encode_packets(t2, tileno, tile, tcd_tcp->numlayers, dest, len,
1723                               volume_info);
1724         t2_destroy(t2);
1725     } else if ((cp->encoding_format == ENCOD_2GR) ||
1726                (cp->encoding_format == ENCOD_3GR)) {
1727         /*if(volume_info) {
1728             volume_info->index_write = 1;
1729         }
1730         gr = golomb_create(tcd->cinfo, volume, cp);
1731         l = golomb_encode(gr, tileno, tile, dest, len, volume_info);
1732         golomb_destroy(gr);*/
1733     }
1734
1735
1736     /*---------------CLEAN-------------------*/
1737     fprintf(stdout, "[INFO] Tcd_encode_tile: %d bytes coded\n", l);
1738     encoding_time = opj_clock() - encoding_time;
1739     opj_event_msg(tcd->cinfo, EVT_INFO, "- tile encoded in %f s\n", encoding_time);
1740
1741     /* cleaning memory */
1742     for (compno = 0; compno < tile->numcomps; compno++) {
1743         tcd->tilec = &tile->comps[compno];
1744         opj_free(tcd->tilec->data);
1745     }
1746
1747     if (l == -999) {
1748         fprintf(stdout, "[ERROR] Unable to perform T2 tier. Return -999.\n");
1749         return 0;
1750     }
1751
1752     return l;
1753 }
1754
1755
1756 bool tcd_decode_tile(opj_tcd_t *tcd, unsigned char *src, int len, int tileno)
1757 {
1758     int l, i;
1759     int compno, eof = 0;
1760     double tile_time, t1_time, dwt_time;
1761
1762     opj_tcd_tile_t *tile = NULL;
1763     opj_t2_t *t2 = NULL;        /* T2 component */
1764
1765     tcd->tcd_tileno = tileno;
1766     tcd->tcd_tile = &(tcd->tcd_volume->tiles[tileno]);
1767     tcd->tcp = &(tcd->cp->tcps[tileno]);
1768     tile = tcd->tcd_tile;
1769
1770     tile_time = opj_clock();    /* time needed to decode a tile */
1771     opj_event_msg(tcd->cinfo, EVT_INFO, "tile %d / %d\n", tileno + 1,
1772                   tcd->cp->tw * tcd->cp->th * tcd->cp->tl);
1773
1774     if ((tcd->cp->encoding_format == ENCOD_2EB) ||
1775             (tcd->cp->encoding_format == ENCOD_3EB)) {
1776         /*--------------TIER2------------------*/
1777         t2 = t2_create(tcd->cinfo, tcd->volume, tcd->cp);
1778         l = t2_decode_packets(t2, src, len, tileno, tile);
1779         t2_destroy(t2);
1780         opj_event_msg(tcd->cinfo, EVT_INFO, "Tcd_decode_tile: %d bytes decoded\n", l);
1781
1782         if (l == -999) {
1783             eof = 1;
1784             opj_event_msg(tcd->cinfo, EVT_ERROR, "Tcd_decode_tile: incomplete bitstream\n");
1785         }
1786
1787         /*------------------TIER1-----------------*/
1788         opj_event_msg(tcd->cinfo, EVT_INFO, "Tcd_decode_tile: Entropy decoding %d \n",
1789                       tcd->cp->encoding_format);
1790         t1_time = opj_clock();  /* time needed to decode a tile */
1791         if (tcd->cp->encoding_format == ENCOD_2EB) {
1792             opj_t1_t *t1 = NULL;        /* T1 component */
1793             t1 = t1_create(tcd->cinfo);
1794             t1_decode_cblks(t1, tile, tcd->tcp);
1795             t1_destroy(t1);
1796         } else if (tcd->cp->encoding_format == ENCOD_3EB) {
1797             opj_t1_3d_t *t1 = NULL;     /* T1 component */
1798             t1 = t1_3d_create(tcd->cinfo);
1799             t1_3d_decode_cblks(t1, tile, tcd->tcp);
1800             t1_3d_destroy(t1);
1801         }
1802
1803         t1_time = opj_clock() - t1_time;
1804 #ifdef VERBOSE
1805         opj_event_msg(tcd->cinfo, EVT_INFO, "- tier-1 took %f s\n", t1_time);
1806 #endif
1807     } else if ((tcd->cp->encoding_format == ENCOD_2GR) ||
1808                (tcd->cp->encoding_format == ENCOD_3GR)) {
1809         opj_event_msg(tcd->cinfo, EVT_INFO,
1810                       "Tcd_decode_tile: Entropy decoding -- Does nothing :-D\n");
1811         /*
1812         gr = golomb_create(tcd->cinfo, tcd->volume, tcd->cp);
1813         l = golomb_decode(gr, tileno, tile, src, len);
1814         golomb_destroy(gr);
1815         if (l == -999) {
1816             eof = 1;
1817             opj_event_msg(tcd->cinfo, EVT_ERROR, "Tcd_decode_tile: incomplete bitstream\n");
1818         }
1819         */
1820     }
1821
1822     /*----------------DWT---------------------*/
1823     fprintf(stdout, "[INFO] Tcd_decode_tile: Inverse DWT\n");
1824     dwt_time = opj_clock(); /* time needed to decode a tile */
1825     for (compno = 0; compno < tile->numcomps; compno++) {
1826         opj_tcd_tilecomp_t *tilec = &tile->comps[compno];
1827         int stops[3], dwtid[3];
1828
1829         for (i = 0; i < 3; i++) {
1830             if (tcd->cp->reduce[i] != 0) {
1831                 tcd->volume->comps[compno].resno_decoded[i] =
1832                     tile->comps[compno].numresolution[i] - tcd->cp->reduce[i] - 1;
1833             }
1834             stops[i] = tilec->numresolution[i] - 1 -
1835                        tcd->volume->comps[compno].resno_decoded[i];
1836             if (stops[i] < 0) {
1837                 stops[i] = 0;
1838             }
1839             dwtid[i] = tcd->cp->tcps->tccps[compno].dwtid[i];
1840         }
1841
1842         dwt_decode(tilec, stops, dwtid);
1843
1844         for (i = 0; i < 3; i++) {
1845             if (tile->comps[compno].numresolution[i] > 0) {
1846                 tcd->volume->comps[compno].factor[i] = tile->comps[compno].numresolution[i] -
1847                                                        (tcd->volume->comps[compno].resno_decoded[i] + 1);
1848                 if ((tcd->volume->comps[compno].factor[i]) < 0) {
1849                     tcd->volume->comps[compno].factor[i] = 0;
1850                 }
1851             }
1852         }
1853     }
1854     dwt_time = opj_clock() - dwt_time;
1855 #ifdef VERBOSE
1856     opj_event_msg(tcd->cinfo, EVT_INFO, "- dwt took %f s\n", dwt_time);
1857 #endif
1858
1859     /*----------------MCT-------------------*/
1860
1861     if (tcd->tcp->mct) {
1862         if (tcd->tcp->tccps[0].reversible == 1) {
1863             mct_decode(tile->comps[0].data, tile->comps[1].data, tile->comps[2].data,
1864                        (tile->comps[0].x1 - tile->comps[0].x0) * (tile->comps[0].y1 -
1865                                tile->comps[0].y0) * (tile->comps[0].z1 - tile->comps[0].z0));
1866         } else {
1867             mct_decode_real(tile->comps[0].data, tile->comps[1].data, tile->comps[2].data,
1868                             (tile->comps[0].x1 - tile->comps[0].x0) * (tile->comps[0].y1 -
1869                                     tile->comps[0].y0) * (tile->comps[0].z1 - tile->comps[0].z0));
1870         }
1871     }
1872
1873     /*---------------TILE-------------------*/
1874
1875     for (compno = 0; compno < tile->numcomps; compno++) {
1876         opj_tcd_tilecomp_t *tilec = &tile->comps[compno];
1877         opj_tcd_resolution_t *res =
1878             &tilec->resolutions[tcd->volume->comps[compno].resno_decoded[0]];
1879         int adjust;
1880         int minval = tcd->volume->comps[compno].sgnd ? -(1 <<
1881                      (tcd->volume->comps[compno].prec - 1)) : 0;
1882         int maxval = tcd->volume->comps[compno].sgnd ? (1 <<
1883                      (tcd->volume->comps[compno].prec - 1)) - 1 : (1 <<
1884                              tcd->volume->comps[compno].prec) - 1;
1885
1886         int tw = tilec->x1 - tilec->x0;
1887         int w = tcd->volume->comps[compno].w;
1888         int th = tilec->y1 - tilec->y0;
1889         int h = tcd->volume->comps[compno].h;
1890
1891         int i, j, k;
1892         int offset_x = int_ceildivpow2(tcd->volume->comps[compno].x0,
1893                                        tcd->volume->comps[compno].factor[0]);
1894         int offset_y = int_ceildivpow2(tcd->volume->comps[compno].y0,
1895                                        tcd->volume->comps[compno].factor[1]);
1896         int offset_z = int_ceildivpow2(tcd->volume->comps[compno].z0,
1897                                        tcd->volume->comps[compno].factor[2]);
1898
1899         if (tcd->cp->transform_format == TRF_3D_RLS ||
1900                 tcd->cp->transform_format == TRF_3D_LSE) {
1901             adjust = 0;
1902         } else {
1903             adjust = tcd->volume->comps[compno].sgnd ? 0 : 1 <<
1904                      (tcd->volume->comps[compno].prec - 1); /*sign=='+' --> 2^(prec-1)*/
1905             if (tcd->volume->comps[compno].dcoffset != 0) {
1906                 adjust += tcd->volume->comps[compno].dcoffset;
1907                 fprintf(stdout, "[INFO] DC Offset applied: DCO = %d -> adjust = %d\n",
1908                         tcd->volume->comps[compno].dcoffset, adjust);
1909             }
1910         }
1911
1912         for (k = res->z0; k < res->z1; k++) {
1913             for (j = res->y0; j < res->y1; j++) {
1914                 for (i = res->x0; i < res->x1; i++) {
1915                     int v;
1916                     float tmp = (float)((tilec->data[i - res->x0 + (j - res->y0) * tw +
1917                                                        (k - res->z0) * tw * th]) / 8192.0);
1918
1919                     if (tcd->tcp->tccps[compno].reversible == 1) {
1920                         v = tilec->data[i - res->x0 + (j - res->y0) * tw + (k - res->z0) * tw * th];
1921                     } else {
1922                         int tmp2 = ((int)(floor(fabs(tmp)))) + ((int) floor(fabs(tmp * 2)) % 2);
1923                         v = ((tmp < 0) ? -tmp2 : tmp2);
1924                     }
1925                     v += adjust;
1926
1927                     tcd->volume->comps[compno].data[(i - offset_x) + (j - offset_y) * w +
1928                                                     (k - offset_z) * w * h] = int_clamp(v, minval, maxval);
1929                 }
1930             }
1931         }
1932     }
1933
1934     tile_time = opj_clock() - tile_time;    /* time needed to decode a tile */
1935     opj_event_msg(tcd->cinfo, EVT_INFO, "- tile decoded in %f s\n", tile_time);
1936
1937     for (compno = 0; compno < tile->numcomps; compno++) {
1938         opj_free(tcd->tcd_volume->tiles[tileno].comps[compno].data);
1939         tcd->tcd_volume->tiles[tileno].comps[compno].data = NULL;
1940     }
1941
1942     if (eof) {
1943         return false;
1944     }
1945
1946     return true;
1947 }
1948