src/lib/openjp3d/t1.c

   1 /*
   2  * Copyright (c) 2001-2003, David Janssens
   3  * Copyright (c) 2002-2003, Yannick Verschueren
   4  * Copyright (c) 2003-2005, Francois Devaux and Antonin Descampe
   5  * Copyright (c) 2005, Herve Drolon, FreeImage Team
   6  * Copyright (c) 2002-2005, Communications and remote sensing Laboratory, Universite catholique de Louvain, Belgium
   7  * All rights reserved.
   8  *
   9  * Redistribution and use in source and binary forms, with or without
  10  * modification, are permitted provided that the following conditions
  11  * are met:
  12  * 1. Redistributions of source code must retain the above copyright
  13  *    notice, this list of conditions and the following disclaimer.
  14  * 2. Redistributions in binary form must reproduce the above copyright
  15  *    notice, this list of conditions and the following disclaimer in the
  16  *    documentation and/or other materials provided with the distribution.
  17  *
  18  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS `AS IS'
  19  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  20  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  21  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  22  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  23  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  24  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  25  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  26  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  27  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  28  * POSSIBILITY OF SUCH DAMAGE.
  29  */
  30
  31 #include "opj_includes.h"
  32
  33 /** @defgroup T1 T1 - Implementation of the tier-1 coding */
  34 /*@{*/
  35
  36 /** @name Local static functions */
  37 /*@{*/
  38
  39 static int t1_getctxno_zc(opj_t1_t *t1, int f, int orient);
  40 static int t1_getctxno_sc(opj_t1_t *t1, int f);
  41 static int t1_getctxno_mag(opj_t1_t *t1, int f);
  42 static int t1_getspb(opj_t1_t *t1, int f);
  43 static int t1_getnmsedec_sig(opj_t1_t *t1, int x, int bitpos);
  44 static int t1_getnmsedec_ref(opj_t1_t *t1, int x, int bitpos);
  45 static void t1_updateflags(int *fp, int s);
  46 /**
  47 Encode significant pass
  48 */
  49 static void t1_enc_sigpass_step(opj_t1_t *t1, int *fp, int *dp, int orient, int bpno, int one, int *nmsedec, char type, int vsc);
  50 /**
  51 Decode significant pass
  52 */
  53 static void t1_dec_sigpass_step(opj_t1_t *t1, int *fp, int *dp, int orient, int oneplushalf, char type, int vsc);
  54 /**
  55 Encode significant pass
  56 */
  57 static void t1_enc_sigpass(opj_t1_t *t1, int w, int h, int l, int bpno, int orient, int *nmsedec, char type, int cblksty);
  58 /**
  59 Decode significant pass
  60 */
  61 static void t1_dec_sigpass(opj_t1_t *t1, int w, int h, int l, int bpno, int orient, char type, int cblksty);
  62 /**
  63 Encode refinement pass
  64 */
  65 static void t1_enc_refpass_step(opj_t1_t *t1, int *fp, int *dp, int bpno, int one, int *nmsedec, char type, int vsc);
  66 /**
  67 Decode refinement pass
  68 */
  69 static void t1_dec_refpass_step(opj_t1_t *t1, int *fp, int *dp, int poshalf, int neghalf, char type, int vsc);
  70 /**
  71 Encode refinement pass
  72 */
  73 static void t1_enc_refpass(opj_t1_t *t1, int w, int h, int l, int bpno, int *nmsedec, char type, int cblksty);
  74 /**
  75 Decode refinement pass
  76 */
  77 static void t1_dec_refpass(opj_t1_t *t1, int w, int h, int l, int bpno, char type, int cblksty);
  78 /**
  79 Encode clean-up pass
  80 */
  81 static void t1_enc_clnpass_step(opj_t1_t *t1, int *fp, int *dp, int orient, int bpno, int one, int *nmsedec, int partial, int vsc);
  82 /**
  83 Decode clean-up pass
  84 */
  85 static void t1_dec_clnpass_step(opj_t1_t *t1, int *fp, int *dp, int orient, int oneplushalf, int partial, int vsc);
  86 /**
  87 Encode clean-up pass
  88 */
  89 static void t1_enc_clnpass(opj_t1_t *t1, int w, int h, int l, int bpno, int orient, int *nmsedec, int cblksty);
  90 /**
  91 Decode clean-up pass
  92 */
  93 static void t1_dec_clnpass(opj_t1_t *t1, int w, int h, int l, int bpno, int orient, int cblksty);
  94 /**
  95 Encode 1 code-block
  96 @param t1 T1 handle
  97 @param cblk Code-block coding parameters
  98 @param orient
  99 @param compno Component number
 100 @param level
 101 @param dwtid
 102 @param stepsize
 103 @param cblksty Code-block style
 104 @param numcomps
 105 @param tile
 106 */
 107 static void t1_encode_cblk(opj_t1_t *t1, opj_tcd_cblk_t * cblk, int orient, int compno, int level[3], int dwtid[3], double stepsize, int cblksty, int numcomps, opj_tcd_tile_t * tile);
 108 /**
 109 Decode 1 code-block
 110 @param t1 T1 handle
 111 @param cblk Code-block coding parameters
 112 @param orient
 113 @param roishift Region of interest shifting value
 114 @param cblksty Code-block style
 115 */
 116 static void t1_decode_cblk(opj_t1_t *t1, opj_tcd_cblk_t * cblk, int orient, int roishift, int cblksty);
 117
 118 static int t1_init_ctxno_zc(int f, int orient);
 119 static int t1_init_ctxno_sc(int f);
 120 static int t1_init_ctxno_mag(int f);
 121 static int t1_init_spb(int f);
 122 /**
 123 Initialize the look-up tables of the Tier-1 coder/decoder
 124 @param t1 T1 handle
 125 */
 126 static void t1_init_luts(opj_t1_t *t1);
 127
 128 /*@}*/
 129
 130 /*@}*/
 131
 132 /* ----------------------------------------------------------------------- */
 133
 134 static int t1_getctxno_zc(opj_t1_t *t1, int f, int orient) {
 135         return t1->lut_ctxno_zc[(orient << 8) | (f & T1_SIG_OTH)];
 136 }
 137
 138 static int t1_getctxno_sc(opj_t1_t *t1, int f) {
 139         return t1->lut_ctxno_sc[(f & (T1_SIG_PRIM | T1_SGN)) >> 4];
 140 }
 141
 142 static int t1_getctxno_mag(opj_t1_t *t1, int f) {
 143         return t1->lut_ctxno_mag[(f & T1_SIG_OTH) | (((f & T1_REFINE) != 0) << 11)];
 144 }
 145
 146 static int t1_getspb(opj_t1_t *t1, int f) {
 147         return t1->lut_spb[(f & (T1_SIG_PRIM | T1_SGN)) >> 4];
 148 }
 149
 150 static int t1_getnmsedec_sig(opj_t1_t *t1, int x, int bitpos) {
 151         if (bitpos > T1_NMSEDEC_FRACBITS) {
 152                 return t1->lut_nmsedec_sig[(x >> (bitpos - T1_NMSEDEC_FRACBITS)) & ((1 << T1_NMSEDEC_BITS) - 1)];
 153         }
 154
 155         return t1->lut_nmsedec_sig0[x & ((1 << T1_NMSEDEC_BITS) - 1)];
 156 }
 157
 158 static int t1_getnmsedec_ref(opj_t1_t *t1, int x, int bitpos) {
 159         if (bitpos > T1_NMSEDEC_FRACBITS) {
 160                 return t1->lut_nmsedec_ref[(x >> (bitpos - T1_NMSEDEC_FRACBITS)) & ((1 << T1_NMSEDEC_BITS) - 1)];
 161         }
 162
 163     return t1->lut_nmsedec_ref0[x & ((1 << T1_NMSEDEC_BITS) - 1)];
 164 }
 165
 166 static void t1_updateflags(int *fp, int s) {
 167         int *np = fp - (T1_MAXCBLKW + 2);
 168         int *sp = fp + (T1_MAXCBLKW + 2);
 169         np[-1] |= T1_SIG_SE;
 170         np[1] |= T1_SIG_SW;
 171         sp[-1] |= T1_SIG_NE;
 172         sp[1] |= T1_SIG_NW;
 173         *np |= T1_SIG_S;
 174         *sp |= T1_SIG_N;
 175         fp[-1] |= T1_SIG_E;
 176         fp[1] |= T1_SIG_W;
 177         if (s) {
 178                 *np |= T1_SGN_S;
 179                 *sp |= T1_SGN_N;
 180                 fp[-1] |= T1_SGN_E;
 181                 fp[1] |= T1_SGN_W;
 182         }
 183 }
 184
 185 static void t1_enc_sigpass_step(opj_t1_t *t1, int *fp, int *dp, int orient, int bpno, int one, int *nmsedec, char type, int vsc) {
 186         int v, flag;
 187
 188         opj_mqc_t *mqc = t1->mqc;       /* MQC component */
 189
 190         flag = vsc ? ((*fp) & (~(T1_SIG_S | T1_SIG_SE | T1_SIG_SW | T1_SGN_S))) : (*fp);
 191         if ((flag & T1_SIG_OTH) && !(flag & (T1_SIG | T1_VISIT))) {
 192                 v = int_abs(*dp) & one ? 1 : 0;
 193                 if (type == T1_TYPE_RAW) {      /* BYPASS/LAZY MODE */
 194                         mqc_setcurctx(mqc, t1_getctxno_zc(t1, flag, orient));   /* ESSAI */
 195                         mqc_bypass_enc(mqc, v);
 196                 } else {
 197                         mqc_setcurctx(mqc, t1_getctxno_zc(t1, flag, orient));
 198                         mqc_encode(mqc, v);
 199                 }
 200                 if (v) {
 201                         v = *dp < 0 ? 1 : 0;
 202                         *nmsedec +=     t1_getnmsedec_sig(t1, int_abs(*dp), bpno + T1_NMSEDEC_FRACBITS);
 203                         if (type == T1_TYPE_RAW) {      /* BYPASS/LAZY MODE */
 204                                 mqc_setcurctx(mqc, t1_getctxno_sc(t1, flag));   /* ESSAI */
 205                                 mqc_bypass_enc(mqc, v);
 206                         } else {
 207                                 mqc_setcurctx(mqc, t1_getctxno_sc(t1, flag));
 208                                 mqc_encode(mqc, v ^ t1_getspb(t1, flag));
 209                         }
 210                         t1_updateflags(fp, v);
 211                         *fp |= T1_SIG;
 212                 }
 213                 *fp |= T1_VISIT;
 214         }
 215 }
 216
 217 static void t1_dec_sigpass_step(opj_t1_t *t1, int *fp, int *dp, int orient, int oneplushalf, char type, int vsc) {
 218         int v, flag;
 219
 220         opj_raw_t *raw = t1->raw;       /* RAW component */
 221         opj_mqc_t *mqc = t1->mqc;       /* MQC component */
 222
 223         flag = vsc ? ((*fp) & (~(T1_SIG_S | T1_SIG_SE | T1_SIG_SW | T1_SGN_S))) : (*fp);
 224         if ((flag & T1_SIG_OTH) && !(flag & (T1_SIG | T1_VISIT))) {
 225                 if (type == T1_TYPE_RAW) {
 226                         if (raw_decode(raw)) {
 227                                 v = raw_decode(raw);    /* ESSAI */
 228                                 *dp = v ? -oneplushalf : oneplushalf;
 229                                 t1_updateflags(fp, v);
 230                                 *fp |= T1_SIG;
 231                         }
 232                 } else {
 233                         mqc_setcurctx(mqc, t1_getctxno_zc(t1, flag, orient));
 234                         if (mqc_decode(mqc)) {
 235                                 mqc_setcurctx(mqc, t1_getctxno_sc(t1, flag));
 236                                 v = mqc_decode(mqc) ^ t1_getspb(t1, flag);
 237                                 *dp = v ? -oneplushalf : oneplushalf;
 238                                 t1_updateflags(fp, v);
 239                                 *fp |= T1_SIG;
 240                         }
 241                 }
 242                 *fp |= T1_VISIT;
 243         }
 244 }                               /* VSC and  BYPASS by Antonin */
 245
 246 static void t1_enc_sigpass(opj_t1_t *t1, int w, int h, int l, int bpno, int orient, int *nmsedec, char type, int cblksty) {
 247         int i, j, k, m, one, vsc;
 248         *nmsedec = 0;
 249         one = 1 << (bpno + T1_NMSEDEC_FRACBITS);
 250         for (m = 0; m < l; m++) {
 251         for (k = 0; k < h; k += 4) {
 252                 for (i = 0; i < w; i++) {
 253                         for (j = k; j < k + 4 && j < h; j++) {
 254                                 vsc = ((cblksty & J3D_CCP_CBLKSTY_VSC) && (j == k + 3 || j == h - 1)) ? 1 : 0;
 255                                 t1_enc_sigpass_step(t1, &t1->flags[1 + m][1 + j][1 + i], &t1->data[m][j][i], orient, bpno, one, nmsedec, type, vsc);
 256                         }
 257                 }
 258         }
 259         }
 260 }
 261
 262 static void t1_dec_sigpass(opj_t1_t *t1, int w, int h, int l, int bpno, int orient, char type, int cblksty) {
 263         int i, j, k, m, one, half, oneplushalf, vsc;
 264         one = 1 << bpno;
 265         half = one >> 1;
 266         oneplushalf = one | half;
 267         for (m = 0; m < l; m++) {
 268         for (k = 0; k < h; k += 4) {
 269                 for (i = 0; i < w; i++) {
 270                         for (j = k; j < k + 4 && j < h; j++) {
 271                                 vsc = ((cblksty & J3D_CCP_CBLKSTY_VSC) && (j == k + 3 || j == h - 1)) ? 1 : 0;
 272                                 t1_dec_sigpass_step(t1, &t1->flags[1 + m][1 + j][1 + i], &t1->data[m][j][i], orient, oneplushalf, type, vsc);
 273                         }
 274                 }
 275         }
 276         }
 277 }                               /* VSC and  BYPASS by Antonin */
 278
 279 static void t1_enc_refpass_step(opj_t1_t *t1, int *fp, int *dp, int bpno, int one, int *nmsedec, char type, int vsc) {
 280         int v, flag;
 281
 282         opj_mqc_t *mqc = t1->mqc;       /* MQC component */
 283
 284         flag = vsc ? ((*fp) & (~(T1_SIG_S | T1_SIG_SE | T1_SIG_SW | T1_SGN_S))) : (*fp);
 285         if ((flag & (T1_SIG | T1_VISIT)) == T1_SIG) {
 286                 *nmsedec += t1_getnmsedec_ref(t1, int_abs(*dp), bpno + T1_NMSEDEC_FRACBITS);
 287                 v = int_abs(*dp) & one ? 1 : 0;
 288                 if (type == T1_TYPE_RAW) {      /* BYPASS/LAZY MODE */
 289                         mqc_setcurctx(mqc, t1_getctxno_mag(t1, flag));  /* ESSAI */
 290                         mqc_bypass_enc(mqc, v);
 291                 } else {
 292                         mqc_setcurctx(mqc, t1_getctxno_mag(t1, flag));
 293                         mqc_encode(mqc, v);
 294                 }
 295                 *fp |= T1_REFINE;
 296         }
 297 }
 298
 299 static void t1_dec_refpass_step(opj_t1_t *t1, int *fp, int *dp, int poshalf, int neghalf, char type, int vsc) {
 300         int v, t, flag;
 301
 302         opj_mqc_t *mqc = t1->mqc;       /* MQC component */
 303         opj_raw_t *raw = t1->raw;       /* RAW component */
 304
 305         flag = vsc ? ((*fp) & (~(T1_SIG_S | T1_SIG_SE | T1_SIG_SW | T1_SGN_S))) : (*fp);
 306         if ((flag & (T1_SIG | T1_VISIT)) == T1_SIG) {
 307                 if (type == T1_TYPE_RAW) {
 308                         mqc_setcurctx(mqc, t1_getctxno_mag(t1, flag));  /* ESSAI */
 309                         v = raw_decode(raw);
 310                 } else {
 311                         mqc_setcurctx(mqc, t1_getctxno_mag(t1, flag));
 312                         v = mqc_decode(mqc);
 313                 }
 314                 t = v ? poshalf : neghalf;
 315                 *dp += *dp < 0 ? -t : t;
 316                 *fp |= T1_REFINE;
 317         }
 318 }                               /* VSC and  BYPASS by Antonin  */
 319
 320 static void t1_enc_refpass(opj_t1_t *t1, int w, int h, int l, int bpno, int *nmsedec, char type, int cblksty) {
 321         int i, j, k, m, one, vsc;
 322         *nmsedec = 0;
 323         one = 1 << (bpno + T1_NMSEDEC_FRACBITS);
 324         for (m = 0; m < l; m++) {
 325                 for (k = 0; k < h; k += 4) {
 326                 for (i = 0; i < w; i++) {
 327                         for (j = k; j < k + 4 && j < h; j++) {
 328                                 vsc = ((cblksty & J3D_CCP_CBLKSTY_VSC) && (j == k + 3 || j == h - 1)) ? 1 : 0;
 329                                 t1_enc_refpass_step(t1, &t1->flags[1 + m][1 + j][1 + i], &t1->data[m][j][i], bpno, one, nmsedec, type, vsc);
 330                         }
 331                 }
 332         }
 333         }
 334 }
 335
 336 static void t1_dec_refpass(opj_t1_t *t1, int w, int h, int l, int bpno, char type, int cblksty) {
 337         int i, j, k, m, one, poshalf, neghalf;
 338         int vsc;
 339         one = 1 << bpno;
 340         poshalf = one >> 1;
 341         neghalf = bpno > 0 ? -poshalf : -1;
 342         for (m = 0; m < l; m++) {
 343                 for (k = 0; k < h; k += 4) {
 344                 for (i = 0; i < w; i++) {
 345                         for (j = k; j < k + 4 && j < h; j++) {
 346                                 vsc = ((cblksty & J3D_CCP_CBLKSTY_VSC) && (j == k + 3 || j == h - 1)) ? 1 : 0;
 347                                 t1_dec_refpass_step(t1, &t1->flags[1 + m][1 + j][1 + i], &t1->data[m][j][i], poshalf, neghalf, type, vsc);
 348                         }
 349                 }
 350         }
 351         }
 352 }                               /* VSC and  BYPASS by Antonin */
 353
 354 static void t1_enc_clnpass_step(opj_t1_t *t1, int *fp, int *dp, int orient, int bpno, int one, int *nmsedec, int partial, int vsc) {
 355         int v, flag;
 356
 357         opj_mqc_t *mqc = t1->mqc;       /* MQC component */
 358
 359         flag = vsc ? ((*fp) & (~(T1_SIG_S | T1_SIG_SE | T1_SIG_SW | T1_SGN_S))) : (*fp);
 360         if (partial) {
 361                 goto LABEL_PARTIAL;
 362         }
 363         if (!(*fp & (T1_SIG | T1_VISIT))) {
 364                 mqc_setcurctx(mqc, t1_getctxno_zc(t1, flag, orient));
 365                 v = int_abs(*dp) & one ? 1 : 0;
 366                 mqc_encode(mqc, v);
 367                 if (v) {
 368 LABEL_PARTIAL:
 369                         *nmsedec += t1_getnmsedec_sig(t1, int_abs(*dp), bpno + T1_NMSEDEC_FRACBITS);
 370                         mqc_setcurctx(mqc, t1_getctxno_sc(t1, flag));
 371                         v = *dp < 0 ? 1 : 0;
 372                         mqc_encode(mqc, v ^ t1_getspb(t1, flag));
 373                         t1_updateflags(fp, v);
 374                         *fp |= T1_SIG;
 375                 }
 376         }
 377         *fp &= ~T1_VISIT;
 378 }
 379
 380 static void t1_dec_clnpass_step(opj_t1_t *t1, int *fp, int *dp, int orient, int oneplushalf, int partial, int vsc) {
 381         int v, flag;
 382
 383         opj_mqc_t *mqc = t1->mqc;       /* MQC component */
 384
 385         flag = vsc ? ((*fp) & (~(T1_SIG_S | T1_SIG_SE | T1_SIG_SW | T1_SGN_S))) : (*fp);
 386         if (partial) {
 387                 goto LABEL_PARTIAL;
 388         }
 389         if (!(flag & (T1_SIG | T1_VISIT))) {
 390                 mqc_setcurctx(mqc, t1_getctxno_zc(t1, flag, orient));
 391                 if (mqc_decode(mqc)) {
 392 LABEL_PARTIAL:
 393                         mqc_setcurctx(mqc, t1_getctxno_sc(t1, flag));
 394                         v = mqc_decode(mqc) ^ t1_getspb(t1, flag);
 395                         *dp = v ? -oneplushalf : oneplushalf;
 396                         t1_updateflags(fp, v);
 397                         *fp |= T1_SIG;
 398                 }
 399         }
 400         *fp &= ~T1_VISIT;
 401 }                               /* VSC and  BYPASS by Antonin */
 402
 403 static void t1_enc_clnpass(opj_t1_t *t1, int w, int h, int l, int bpno, int orient, int *nmsedec, int cblksty) {
 404         int i, j, k, m, one, agg, runlen, vsc;
 405
 406         opj_mqc_t *mqc = t1->mqc;       /* MQC component */
 407
 408         *nmsedec = 0;
 409         one = 1 << (bpno + T1_NMSEDEC_FRACBITS);
 410         for (m = 0; m < l; m++) {
 411                 for (k = 0; k < h; k += 4) {
 412                         for (i = 0; i < w; i++) {
 413                                 if (k + 3 < h) {
 414                                         if (cblksty & J3D_CCP_CBLKSTY_VSC) {
 415                                                 agg = !(t1->flags[1 + m][1 + k][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH)
 416                                                         || t1->flags[1 + m][1 + k + 1][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH)
 417                                                         || t1->flags[1 + m][1 + k + 2][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH)
 418                                                         || (t1->flags[1 + m][1 + k + 3][1 + i]
 419                                                         & (~(T1_SIG_S | T1_SIG_SE | T1_SIG_SW | T1_SGN_S))) & (T1_SIG | T1_VISIT | T1_SIG_OTH));
 420                                         } else {
 421                                                 agg = !(t1->flags[1 + m][1 + k][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH)
 422                                                         || t1->flags[1 + m][1 + k + 1][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH)
 423                                                         || t1->flags[1 + m][1 + k + 2][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH)
 424                                                         || t1->flags[1 + m][1 + k + 3][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH));
 425                                         }
 426                                 } else {
 427                                         agg = 0;
 428                                 }
 429                                 if (agg) {
 430                                         for (runlen = 0; runlen < 4; runlen++) {
 431                                                 if (int_abs(t1->data[m][k + runlen][i]) & one)
 432                                                         break;
 433                                         }
 434                                         mqc_setcurctx(mqc, T1_CTXNO_AGG);
 435                                         mqc_encode(mqc, runlen != 4);
 436                                         if (runlen == 4) {
 437                                                 continue;
 438                                         }
 439                                         mqc_setcurctx(mqc, T1_CTXNO_UNI);
 440                                         mqc_encode(mqc, runlen >> 1);
 441                                         mqc_encode(mqc, runlen & 1);
 442                                 } else {
 443                                         runlen = 0;
 444                                 }
 445                                 for (j = k + runlen; j < k + 4 && j < h; j++) {
 446                                         vsc = ((cblksty & J3D_CCP_CBLKSTY_VSC) && (j == k + 3 || j == h - 1)) ? 1 : 0;
 447                                         t1_enc_clnpass_step(t1, &(t1->flags[1 + m][1 + j][1 + i]), &(t1->data[m][j][i]), orient, bpno, one, nmsedec, agg && (j == k + runlen), vsc);
 448                                 }
 449                         }
 450         }
 451         }
 452 }
 453
 454 static void t1_dec_clnpass(opj_t1_t *t1, int w, int h, int l, int bpno, int orient, int cblksty) {
 455         int i, j, k, m, one, half, oneplushalf, agg, runlen, vsc;
 456         int segsym = cblksty & J3D_CCP_CBLKSTY_SEGSYM;
 457
 458         opj_mqc_t *mqc = t1->mqc;       /* MQC component */
 459
 460         one = 1 << bpno;
 461         half = one >> 1;
 462         oneplushalf = one | half;
 463         for (m = 0; m < l; m++) {
 464                 for (k = 0; k < h; k += 4) {
 465                 for (i = 0; i < w; i++) {
 466                         if (k + 3 < h) {
 467                                 if (cblksty & J3D_CCP_CBLKSTY_VSC) {
 468                                         agg = !(t1->flags[1 + m][1 + k][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH)
 469                                                 || t1->flags[1 + m][1 + k + 1][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH)
 470                                                 || t1->flags[1 + m][1 + k + 2][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH)
 471                                                 || (t1->flags[1 + m][1 + k + 3][1 + i]
 472                                                 & (~(T1_SIG_S | T1_SIG_SE | T1_SIG_SW | T1_SGN_S))) & (T1_SIG | T1_VISIT | T1_SIG_OTH));
 473                                 } else {
 474                                         agg = !(t1->flags[1 + m][1 + k][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH)
 475                                                 || t1->flags[1 + m][1 + k + 1][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH)
 476                                                 || t1->flags[1 + m][1 + k + 2][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH)
 477                                                 || t1->flags[1 + m][1 + k + 3][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH));
 478                                 }
 479                         } else {
 480                                 agg = 0;
 481                         }
 482                         if (agg) {
 483                                 mqc_setcurctx(mqc, T1_CTXNO_AGG);
 484                                 if (!mqc_decode(mqc)) {
 485                                         continue;
 486                                 }
 487                                 mqc_setcurctx(mqc, T1_CTXNO_UNI);
 488                                 runlen = mqc_decode(mqc);
 489                                 runlen = (runlen << 1) | mqc_decode(mqc);
 490                         } else {
 491                                 runlen = 0;
 492                         }
 493                         for (j = k + runlen; j < k + 4 && j < h; j++) {
 494                                 vsc = ((cblksty & J3D_CCP_CBLKSTY_VSC) && (j == k + 3 || j == h - 1)) ? 1 : 0;
 495                                 t1_dec_clnpass_step(t1, &t1->flags[1 + m][1 + j][1 + i], &t1->data[m][j][i], orient, oneplushalf, agg && (j == k + runlen), vsc);
 496                         }
 497                 }
 498         }
 499         }
 500         if (segsym) {
 501                 int v = 0;
 502                 mqc_setcurctx(mqc, T1_CTXNO_UNI);
 503                 v = mqc_decode(mqc);
 504                 v = (v << 1) | mqc_decode(mqc);
 505                 v = (v << 1) | mqc_decode(mqc);
 506                 v = (v << 1) | mqc_decode(mqc);
 507                 /*
 508                 if (v!=0xa) {
 509                         opj_event_msg(t1->cinfo, EVT_WARNING, "Bad segmentation symbol %x\n", v);
 510                 }
 511                 */
 512         }
 513 }                               /* VSC and  BYPASS by Antonin */
 514
 515
 516 static void t1_encode_cblk(opj_t1_t *t1, opj_tcd_cblk_t * cblk, int orient, int compno, int level[3], int dwtid[3], double stepsize, int cblksty, int numcomps, opj_tcd_tile_t * tile) {
 517         int i, j, k;
 518         int w, h, l;
 519         int passno;
 520         int bpno, passtype;
 521         int max;
 522         int nmsedec = 0;
 523         double cumwmsedec = 0;
 524         char type = T1_TYPE_MQ;
 525
 526         opj_mqc_t *mqc = t1->mqc;       /* MQC component */
 527
 528         w = cblk->x1 - cblk->x0;
 529         h = cblk->y1 - cblk->y0;
 530         l = cblk->z1 - cblk->z0;
 531
 532         max = 0;
 533         for (k = 0; k < l; k++) {
 534                 for (j = 0; j < h; j++) {
 535                         for (i = 0; i < w; i++) {
 536                                 max = int_max(max, int_abs(t1->data[k][j][i]));
 537                         }
 538                 }
 539         }
 540         for (k = 0; k <= l; k++) {
 541                 for (j = 0; j <= h; j++) {
 542                         for (i = 0; i <= w; i++) {
 543                                 t1->flags[k][j][i] = 0;
 544                         }
 545                 }
 546         }
 547
 548         cblk->numbps = max ? (int_floorlog2(max) + 1) - T1_NMSEDEC_FRACBITS : 0;
 549
 550         bpno = cblk->numbps - 1;
 551         passtype = 2;
 552
 553         mqc_reset_enc(mqc);
 554         mqc_init_enc(mqc, cblk->data);
 555
 556         for (passno = 0; bpno >= 0; passno++) {
 557                 opj_tcd_pass_t *pass = &cblk->passes[passno];
 558                 int correction = 3;
 559                 double tmpwmsedec;
 560                 type = ((bpno < (cblk->numbps - 4)) && (passtype < 2) && (cblksty & J3D_CCP_CBLKSTY_LAZY)) ? T1_TYPE_RAW : T1_TYPE_MQ;
 561                 /*fprintf(stdout,"passno %d passtype %d w %d h %d l %d bpno %d orient %d type %d cblksty %d\n",passno,passtype,w,h,l,bpno,orient,type,cblksty);*/
 562
 563                 switch (passtype) {
 564                         case 0:
 565                                 t1_enc_sigpass(t1, w, h, l, bpno, orient, &nmsedec, type, cblksty);
 566                                 break;
 567                         case 1:
 568                                 t1_enc_refpass(t1, w, h, l, bpno, &nmsedec, type, cblksty);
 569                                 break;
 570                         case 2:
 571                                 /*fprintf(stdout,"w %d h %d l %d bpno %d orient %d \n",w,h,l,bpno,orient);*/
 572                                 t1_enc_clnpass(t1, w, h, l, bpno, orient, &nmsedec, cblksty);
 573                                 /* code switch SEGMARK (i.e. SEGSYM) */
 574                                 if (cblksty & J3D_CCP_CBLKSTY_SEGSYM)
 575                                         mqc_segmark_enc(mqc);
 576                                 break;
 577                 }
 578
 579                 /* fixed_quality */
 580                 tmpwmsedec = t1_getwmsedec(nmsedec, compno, level, orient, bpno, stepsize, numcomps, dwtid);
 581                 cumwmsedec += tmpwmsedec;
 582                 tile->distotile += tmpwmsedec;
 583
 584                 /* Code switch "RESTART" (i.e. TERMALL) */
 585                 if ((cblksty & J3D_CCP_CBLKSTY_TERMALL) && !((passtype == 2) && (bpno - 1 < 0))) {
 586                         if (type == T1_TYPE_RAW) {
 587                                 mqc_flush(mqc);
 588                                 correction = 1;
 589                                 /* correction = mqc_bypass_flush_enc(); */
 590                         } else {                        /* correction = mqc_restart_enc(); */
 591                                 mqc_flush(mqc);
 592                                 correction = 1;
 593                         }
 594                         pass->term = 1;
 595                 } else {
 596                         if (((bpno < (cblk->numbps - 4) && (passtype > 0))
 597                                 || ((bpno == (cblk->numbps - 4)) && (passtype == 2))) && (cblksty & J3D_CCP_CBLKSTY_LAZY)) {
 598                                 if (type == T1_TYPE_RAW) {
 599                                         mqc_flush(mqc);
 600                                         correction = 1;
 601                                         /* correction = mqc_bypass_flush_enc(); */
 602                                 } else {                /* correction = mqc_restart_enc(); */
 603                                         mqc_flush(mqc);
 604                                         correction = 1;
 605                                 }
 606                                 pass->term = 1;
 607                         } else {
 608                                 pass->term = 0;
 609                         }
 610                 }
 611
 612                 if (++passtype == 3) {
 613                         passtype = 0;
 614                         bpno--;
 615                 }
 616
 617                 if (pass->term && bpno > 0) {
 618                         type = ((bpno < (cblk->numbps - 4)) && (passtype < 2) && (cblksty & J3D_CCP_CBLKSTY_LAZY)) ? T1_TYPE_RAW : T1_TYPE_MQ;
 619                         if (type == T1_TYPE_RAW)
 620                                 mqc_bypass_init_enc(mqc);
 621                         else
 622                                 mqc_restart_init_enc(mqc);
 623                 }
 624
 625                 pass->distortiondec = cumwmsedec;
 626                 pass->rate = mqc_numbytes(mqc) + correction;    /* FIXME */
 627                 pass->len = pass->rate - (passno == 0 ? 0 : cblk->passes[passno - 1].rate);
 628
 629                 /* Code-switch "RESET" */
 630                 if (cblksty & J3D_CCP_CBLKSTY_RESET)
 631                         mqc_reset_enc(mqc);
 632         }
 633
 634         /* Code switch "ERTERM" (i.e. PTERM) */
 635         if (cblksty & J3D_CCP_CBLKSTY_PTERM)
 636                 mqc_erterm_enc(mqc);
 637         else /* Default coding */ if (!(cblksty & J3D_CCP_CBLKSTY_LAZY))
 638                 mqc_flush(mqc);
 639
 640         cblk->totalpasses = passno;
 641 }
 642
 643 static void t1_decode_cblk(opj_t1_t *t1, opj_tcd_cblk_t * cblk, int orient, int roishift, int cblksty) {
 644         int i, j, k, w, h, l;
 645         int bpno, passtype;
 646         int segno, passno;
 647         char type = T1_TYPE_MQ; /* BYPASS mode */
 648
 649         opj_raw_t *raw = t1->raw;       /* RAW component */
 650         opj_mqc_t *mqc = t1->mqc;       /* MQC component */
 651
 652         w = cblk->x1 - cblk->x0;
 653         h = cblk->y1 - cblk->y0;
 654         l = cblk->z1 - cblk->z0;
 655
 656     for (k = 0; k < l; k++) {
 657                 for (j = 0; j < h; j++) {
 658                         for (i = 0; i < w; i++) {
 659                                 t1->data[k][j][i] = 0;
 660                         }
 661                 }
 662         }
 663
 664         for (k = 0; k <= l; k++) {
 665                 for (j = 0; j <= h; j++) {
 666                         for (i = 0; i <= w; i++) {
 667                                 t1->flags[k][j][i] = 0;
 668                         }
 669                 }
 670         }
 671
 672         bpno = roishift + cblk->numbps - 1;
 673         passtype = 2;
 674
 675         mqc_reset_enc(mqc);
 676
 677         for (segno = 0; segno < cblk->numsegs; segno++) {
 678                 opj_tcd_seg_t *seg = &cblk->segs[segno];
 679
 680                 /* BYPASS mode */
 681                 type = ((bpno <= (cblk->numbps - 1) - 4) && (passtype < 2) && (cblksty & J3D_CCP_CBLKSTY_LAZY)) ? T1_TYPE_RAW : T1_TYPE_MQ;
 682                 if (type == T1_TYPE_RAW) {
 683                         raw_init_dec(raw, seg->data, seg->len);
 684                 } else {
 685                         mqc_init_dec(mqc, seg->data, seg->len);
 686                 }
 687
 688                 for (passno = 0; passno < seg->numpasses; passno++) {
 689                         switch (passtype) {
 690                                 case 0:
 691                                         t1_dec_sigpass(t1, w, h, l, bpno+1, orient, type, cblksty);
 692                                         break;
 693                                 case 1:
 694                                         t1_dec_refpass(t1, w, h, l, bpno+1, type, cblksty);
 695                                         break;
 696                                 case 2:
 697                                         t1_dec_clnpass(t1, w, h, l, bpno+1, orient, cblksty);
 698                                         break;
 699                         }
 700
 701                         if ((cblksty & J3D_CCP_CBLKSTY_RESET) && type == T1_TYPE_MQ) {
 702                                 mqc_reset_enc(mqc);
 703                         }
 704                         if (++passtype == 3) {
 705                                 passtype = 0;
 706                                 bpno--;
 707                         }
 708                 }
 709         }
 710 }
 711
 712 static int t1_init_ctxno_zc(int f, int orient) {
 713         int h, v, d, n, t, hv;
 714         n = 0;
 715         h = ((f & T1_SIG_W) != 0) + ((f & T1_SIG_E) != 0);
 716         v = ((f & T1_SIG_N) != 0) + ((f & T1_SIG_S) != 0);
 717         d = ((f & T1_SIG_NW) != 0) + ((f & T1_SIG_NE) != 0) + ((f & T1_SIG_SE) != 0) + ((f & T1_SIG_SW) != 0);
 718
 719         switch (orient) {
 720                 case 2:
 721                         t = h;
 722                         h = v;
 723                         v = t;
 724                 case 0:
 725                 case 1:
 726                         if (!h) {
 727                                 if (!v) {
 728                                         if (!d)
 729                                                 n = 0;
 730                                         else if (d == 1)
 731                                                 n = 1;
 732                                         else
 733                                                 n = 2;
 734                                 } else if (v == 1) {
 735                                         n = 3;
 736                                 } else {
 737                                         n = 4;
 738                                 }
 739                         } else if (h == 1) {
 740                                 if (!v) {
 741                                         if (!d)
 742                                                 n = 5;
 743                                         else
 744                                                 n = 6;
 745                                 } else {
 746                                         n = 7;
 747                                 }
 748                         } else
 749                                 n = 8;
 750                         break;
 751                 case 3:
 752                         hv = h + v;
 753                         if (!d) {
 754                                 if (!hv) {
 755                                         n = 0;
 756                                 } else if (hv == 1) {
 757                                         n = 1;
 758                                 } else {
 759                                         n = 2;
 760                                 }
 761                         } else if (d == 1) {
 762                                 if (!hv) {
 763                                         n = 3;
 764                                 } else if (hv == 1) {
 765                                         n = 4;
 766                                 } else {
 767                                         n = 5;
 768                                 }
 769                         } else if (d == 2) {
 770                                 if (!hv) {
 771                                         n = 6;
 772                                 } else {
 773                                         n = 7;
 774                                 }
 775                         } else {
 776                                 n = 8;
 777                         }
 778                         break;
 779         }
 780
 781         return (T1_CTXNO_ZC + n);
 782 }
 783
 784 static int t1_init_ctxno_sc(int f) {
 785         int hc, vc, n;
 786         n = 0;
 787
 788         hc = int_min(((f & (T1_SIG_E | T1_SGN_E)) ==
 789                 T1_SIG_E) + ((f & (T1_SIG_W | T1_SGN_W)) == T1_SIG_W),
 790                1) - int_min(((f & (T1_SIG_E | T1_SGN_E)) ==
 791                    (T1_SIG_E | T1_SGN_E)) +
 792                    ((f & (T1_SIG_W | T1_SGN_W)) ==
 793                    (T1_SIG_W | T1_SGN_W)), 1);
 794
 795         vc = int_min(((f & (T1_SIG_N | T1_SGN_N)) ==
 796                 T1_SIG_N) + ((f & (T1_SIG_S | T1_SGN_S)) == T1_SIG_S),
 797                1) - int_min(((f & (T1_SIG_N | T1_SGN_N)) ==
 798                    (T1_SIG_N | T1_SGN_N)) +
 799                    ((f & (T1_SIG_S | T1_SGN_S)) ==
 800                    (T1_SIG_S | T1_SGN_S)), 1);
 801
 802         if (hc < 0) {
 803                 hc = -hc;
 804                 vc = -vc;
 805         }
 806         if (!hc) {
 807                 if (vc == -1)
 808                         n = 1;
 809                 else if (!vc)
 810                         n = 0;
 811                 else
 812                         n = 1;
 813         } else if (hc == 1) {
 814                 if (vc == -1)
 815                         n = 2;
 816                 else if (!vc)
 817                         n = 3;
 818                 else
 819                         n = 4;
 820         }
 821
 822         return (T1_CTXNO_SC + n);
 823 }
 824
 825 static int t1_init_ctxno_mag(int f) {
 826         int n;
 827         if (!(f & T1_REFINE))
 828                 n = (f & (T1_SIG_OTH)) ? 1 : 0;
 829         else
 830                 n = 2;
 831
 832         return (T1_CTXNO_MAG + n);
 833 }
 834
 835 static int t1_init_spb(int f) {
 836         int hc, vc, n;
 837
 838         hc = int_min(((f & (T1_SIG_E | T1_SGN_E)) ==
 839                 T1_SIG_E) + ((f & (T1_SIG_W | T1_SGN_W)) == T1_SIG_W),
 840                1) - int_min(((f & (T1_SIG_E | T1_SGN_E)) ==
 841                    (T1_SIG_E | T1_SGN_E)) +
 842                    ((f & (T1_SIG_W | T1_SGN_W)) ==
 843                    (T1_SIG_W | T1_SGN_W)), 1);
 844
 845         vc = int_min(((f & (T1_SIG_N | T1_SGN_N)) ==
 846                 T1_SIG_N) + ((f & (T1_SIG_S | T1_SGN_S)) == T1_SIG_S),
 847                1) - int_min(((f & (T1_SIG_N | T1_SGN_N)) ==
 848                    (T1_SIG_N | T1_SGN_N)) +
 849                    ((f & (T1_SIG_S | T1_SGN_S)) ==
 850                    (T1_SIG_S | T1_SGN_S)), 1);
 851
 852         if (!hc && !vc)
 853                 n = 0;
 854         else
 855                 n = (!(hc > 0 || (!hc && vc > 0)));
 856
 857         return n;
 858 }
 859
 860 static void t1_init_luts(opj_t1_t *t1) {
 861         int i, j;
 862         double u, v, t;
 863         for (j = 0; j < 4; j++) {
 864                 for (i = 0; i < 256; ++i) {
 865                         t1->lut_ctxno_zc[(j << 8) | i] = t1_init_ctxno_zc(i, j);
 866                 }
 867         }
 868         for (i = 0; i < 256; i++) {
 869                 t1->lut_ctxno_sc[i] = t1_init_ctxno_sc(i << 4);
 870         }
 871         for (j = 0; j < 2; j++) {
 872                 for (i = 0; i < 2048; ++i) {
 873                         t1->lut_ctxno_mag[(j << 11) + i] = t1_init_ctxno_mag((j ? T1_REFINE : 0) | i);
 874                 }
 875         }
 876         for (i = 0; i < 256; ++i) {
 877                 t1->lut_spb[i] = t1_init_spb(i << 4);
 878         }
 879         /* FIXME FIXME FIXME */
 880         /* fprintf(stdout,"nmsedec luts:\n"); */
 881         for (i = 0; i < (1 << T1_NMSEDEC_BITS); i++) {
 882                 t = i / pow(2, T1_NMSEDEC_FRACBITS);
 883                 u = t;
 884                 v = t - 1.5;
 885                 t1->lut_nmsedec_sig[i] =
 886                         int_max(0,
 887                         (int) (floor((u * u - v * v) * pow(2, T1_NMSEDEC_FRACBITS) + 0.5) / pow(2, T1_NMSEDEC_FRACBITS) * 8192.0));
 888                 t1->lut_nmsedec_sig0[i] =
 889                         int_max(0,
 890                         (int) (floor((u * u) * pow(2, T1_NMSEDEC_FRACBITS) + 0.5) / pow(2, T1_NMSEDEC_FRACBITS) * 8192.0));
 891                 u = t - 1.0;
 892                 if (i & (1 << (T1_NMSEDEC_BITS - 1))) {
 893                         v = t - 1.5;
 894                 } else {
 895                         v = t - 0.5;
 896                 }
 897                 t1->lut_nmsedec_ref[i] =
 898                         int_max(0,
 899                         (int) (floor((u * u - v * v) * pow(2, T1_NMSEDEC_FRACBITS) + 0.5) / pow(2, T1_NMSEDEC_FRACBITS) * 8192.0));
 900                 t1->lut_nmsedec_ref0[i] =
 901                         int_max(0,
 902                         (int) (floor((u * u) * pow(2, T1_NMSEDEC_FRACBITS) + 0.5) / pow(2, T1_NMSEDEC_FRACBITS) * 8192.0));
 903         }
 904 }
 905
 906 /* ----------------------------------------------------------------------- */
 907
 908 opj_t1_t* t1_create(opj_common_ptr cinfo) {
 909         opj_t1_t *t1 = (opj_t1_t*)opj_malloc(sizeof(opj_t1_t));
 910         if(t1) {
 911                 t1->cinfo = cinfo;
 912                 /* create MQC and RAW handles */
 913                 t1->mqc = mqc_create();
 914                 t1->raw = raw_create();
 915                 /* initialize the look-up tables of the Tier-1 coder/decoder */
 916                 t1_init_luts(t1);
 917         }
 918         return t1;
 919 }
 920
 921 void t1_destroy(opj_t1_t *t1) {
 922         if(t1) {
 923                 /* destroy MQC and RAW handles */
 924                 mqc_destroy(t1->mqc);
 925                 raw_destroy(t1->raw);
 926                 /*opj_free(t1->data);*/
 927                 /*opj_free(t1->flags);*/
 928                 opj_free(t1);
 929         }
 930 }
 931
 932 void t1_encode_cblks(opj_t1_t *t1, opj_tcd_tile_t *tile, opj_tcp_t *tcp) {
 933         int compno, resno, bandno, precno, cblkno;
 934         int x, y, z, i, j, k, orient;
 935         int n=0;
 936         int level[3];
 937         FILE *fid = NULL;
 938 /*      char filename[10];*/
 939         tile->distotile = 0;            /* fixed_quality */
 940
 941         for (compno = 0; compno < tile->numcomps; compno++) {
 942                 opj_tcd_tilecomp_t *tilec = &tile->comps[compno];
 943
 944                 for (resno = 0; resno < tilec->numresolution[0]; resno++) {
 945                         opj_tcd_resolution_t *res = &tilec->resolutions[resno];
 946
 947                         /* Weighted first order entropy
 948                         sprintf(filename,"res%d.txt",resno);
 949                         if ((fid = fopen(filename,"w")) == 0){
 950                                 fprintf(stdout,"Error while opening %s\n", filename);
 951                                 exit(1);
 952                         }
 953                         */
 954                         for (bandno = 0; bandno < res->numbands; bandno++) {
 955                                 opj_tcd_band_t *band = &res->bands[bandno];
 956                                 for (precno = 0; precno < res->prctno[0] * res->prctno[1] * res->prctno[2]; precno++) {
 957                                         opj_tcd_precinct_t *prc = &band->precincts[precno];
 958
 959                                         for (cblkno = 0; cblkno < prc->cblkno[0] * prc->cblkno[1] * prc->cblkno[2]; cblkno++) {
 960                                                 opj_tcd_cblk_t *cblk = &prc->cblks[cblkno];
 961
 962                                                 /*fprintf(stdout,"Precno %d Cblkno %d \n",precno,cblkno);*/
 963                                                 if (band->bandno == 0) {
 964                                                         x = cblk->x0 - band->x0;
 965                                                         y = cblk->y0 - band->y0;
 966                                                         z = cblk->z0 - band->z0;
 967                                                 } else if (band->bandno == 1) {
 968                                                         opj_tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
 969                                                         x = pres->x1 - pres->x0 + cblk->x0 - band->x0;
 970                                                         y = cblk->y0 - band->y0;
 971                                                         z = cblk->z0 - band->z0;
 972                                                 } else if (band->bandno == 2) {
 973                                                         opj_tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
 974                                                         x = cblk->x0 - band->x0;
 975                                                         y = pres->y1 - pres->y0 + cblk->y0 - band->y0;
 976                                                         z = cblk->z0 - band->z0;
 977                                                 } else if (band->bandno == 3) {
 978                                                         opj_tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
 979                                                         x = pres->x1 - pres->x0 + cblk->x0 - band->x0;
 980                                                         y = pres->y1 - pres->y0 + cblk->y0 - band->y0;
 981                                                         z = cblk->z0 - band->z0;
 982                                                 } else if (band->bandno == 4) {
 983                                                         opj_tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
 984                                                         x = cblk->x0 - band->x0;
 985                                                         y = cblk->y0 - band->y0;
 986                                                         z = pres->z1 - pres->z0 + cblk->z0 - band->z0;
 987                                                 } else if (band->bandno == 5) {
 988                                                         opj_tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
 989                                                         x = pres->x1 - pres->x0 + cblk->x0 - band->x0;
 990                                                         y = cblk->y0 - band->y0;
 991                                                         z = pres->z1 - pres->z0 + cblk->z0 - band->z0;
 992                                                 } else if (band->bandno == 6) {
 993                                                         opj_tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
 994                                                         x = cblk->x0 - band->x0;
 995                                                         y = pres->y1 - pres->y0 + cblk->y0 - band->y0;
 996                                                         z = pres->z1 - pres->z0 + cblk->z0 - band->z0;
 997                                                 } else if (band->bandno == 7) {
 998                                                         opj_tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
 999                                                         x = pres->x1 - pres->x0 + cblk->x0 - band->x0;
1000                                                         y = pres->y1 - pres->y0 + cblk->y0 - band->y0;
1001                                                         z = pres->z1 - pres->z0 + cblk->z0 - band->z0;
1002                                                 }
1003
1004                                                 if (tcp->tccps[compno].reversible == 1) {
1005                                                         for (k = 0; k < cblk->z1 - cblk->z0; k++) {
1006                                                                 for (j = 0; j < cblk->y1 - cblk->y0; j++) {
1007                                     for (i = 0; i < cblk->x1 - cblk->x0; i++) {
1008                                         t1->data[k][j][i] =
1009                                                                                 tilec->data[(x + i) + (y + j) * (tilec->x1 - tilec->x0) + (z + k) * (tilec->x1 - tilec->x0) * (tilec->y1 - tilec->y0)] << T1_NMSEDEC_FRACBITS;
1010 /*fprintf(fid," %d",t1->data[k][j][i]);*/
1011                                                                         }
1012                                                                 }
1013                                                         }
1014                                                 } else if (tcp->tccps[compno].reversible == 0) {
1015                                                         for (k = 0; k < cblk->z1 - cblk->z0; k++) {
1016                                                                 for (j = 0; j < cblk->y1 - cblk->y0; j++) {
1017                                     for (i = 0; i < cblk->x1 - cblk->x0; i++) {
1018                                         t1->data[k][j][i] = fix_mul(
1019                                                                                 tilec->data[(x + i) + (y + j) * (tilec->x1 - tilec->x0) + (z + k) * (tilec->x1 - tilec->x0) * (tilec->y1 - tilec->y0)],
1020                                                                                 8192 * 8192 / ((int) floor(band->stepsize * 8192))) >> (13 - T1_NMSEDEC_FRACBITS);
1021                                                                         }
1022                                                                 }
1023                                                         }
1024                                                 }
1025
1026                                                 orient = band->bandno;  /* FIXME */
1027                                                 if (orient == 2) {
1028                                                         orient = 1;
1029                                                 } else if (orient == 1) {
1030                                                         orient = 2;
1031                                                 }
1032                                                 for (i = 0; i < 3; i++)
1033                                                         level[i] = tilec->numresolution[i] - 1 - resno;
1034                                                 /*fprintf(stdout,"t1_encode_cblk(t1, cblk, %d, %d, %d %d %d, %d, %f, %d, %d, tile);\n", orient, compno, level[0], level[1], level[2], tcp->tccps[compno].reversible, band->stepsize, tcp->tccps[compno].cblksty, tile->numcomps);*/
1035                                                 t1_encode_cblk(t1, cblk, orient, compno, level, tcp->tccps[compno].dwtid, band->stepsize, tcp->tccps[compno].cblksty, tile->numcomps, tile);
1036
1037                                         } /* cblkno */
1038                                 } /* precno */
1039 /*fprintf(fid,"\n");*/
1040                         } /* bandno */
1041 /*fclose(fid);*/
1042                 } /* resno  */
1043         } /* compno  */
1044 }
1045
1046 void t1_decode_cblks(opj_t1_t *t1, opj_tcd_tile_t *tile, opj_tcp_t *tcp) {
1047         int compno, resno, bandno, precno, cblkno;
1048
1049         for (compno = 0; compno < tile->numcomps; compno++) {
1050                 opj_tcd_tilecomp_t *tilec = &tile->comps[compno];
1051
1052                 for (resno = 0; resno < tilec->numresolution[0]; resno++) {
1053                         opj_tcd_resolution_t *res = &tilec->resolutions[resno];
1054
1055                         for (bandno = 0; bandno < res->numbands; bandno++) {
1056                                 opj_tcd_band_t *band = &res->bands[bandno];
1057
1058                                 for (precno = 0; precno < res->prctno[0] * res->prctno[1] * res->prctno[2]; precno++) {
1059                                         opj_tcd_precinct_t *prc = &band->precincts[precno];
1060
1061                                         for (cblkno = 0; cblkno < prc->cblkno[0] * prc->cblkno[1] * prc->cblkno[2]; cblkno++) {
1062                                                 int x, y, k, i, j, z, orient;
1063                                                 opj_tcd_cblk_t *cblk = &prc->cblks[cblkno];
1064
1065                                                 orient = band->bandno;  /* FIXME */
1066                                                 if (orient == 2) {
1067                                                         orient = 1;
1068                                                 } else if (orient == 1) {
1069                                                         orient = 2;
1070                                                 }
1071
1072                                                 t1_decode_cblk(t1, cblk, orient, tcp->tccps[compno].roishift, tcp->tccps[compno].cblksty);
1073
1074                                                 if (band->bandno == 0) {
1075                                                         x = cblk->x0 - band->x0;
1076                                                         y = cblk->y0 - band->y0;
1077                                                         z = cblk->z0 - band->z0;
1078                                                 } else if (band->bandno == 1) {
1079                                                         opj_tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
1080                                                         x = pres->x1 - pres->x0 + cblk->x0 - band->x0;
1081                                                         y = cblk->y0 - band->y0;
1082                                                         z = cblk->z0 - band->z0;
1083                                                 } else if (band->bandno == 2) {
1084                                                         opj_tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
1085                                                         x = cblk->x0 - band->x0;
1086                                                         y = pres->y1 - pres->y0 + cblk->y0 - band->y0;
1087                                                         z = cblk->z0 - band->z0;
1088                                                 } else if (band->bandno == 3) {
1089                                                         opj_tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
1090                                                         x = pres->x1 - pres->x0 + cblk->x0 - band->x0;
1091                                                         y = pres->y1 - pres->y0 + cblk->y0 - band->y0;
1092                                                         z = cblk->z0 - band->z0;
1093                                                 } else if (band->bandno == 4) {
1094                                                         opj_tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
1095                                                         x = cblk->x0 - band->x0;
1096                                                         y = cblk->y0 - band->y0;
1097                                                         z = pres->z1 - pres->z0 + cblk->z0 - band->z0;
1098                                                 } else if (band->bandno == 5) {
1099                                                         opj_tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
1100                                                         x = pres->x1 - pres->x0 + cblk->x0 - band->x0;
1101                                                         y = cblk->y0 - band->y0;
1102                                                         z = pres->z1 - pres->z0 + cblk->z0 - band->z0;
1103                                                 } else if (band->bandno == 6) {
1104                                                         opj_tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
1105                                                         x = cblk->x0 - band->x0;
1106                                                         y = pres->y1 - pres->y0 + cblk->y0 - band->y0;
1107                                                         z = pres->z1 - pres->z0 + cblk->z0 - band->z0;
1108                                                 } else if (band->bandno == 7) {
1109                                                         opj_tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
1110                                                         x = pres->x1 - pres->x0 + cblk->x0 - band->x0;
1111                                                         y = pres->y1 - pres->y0 + cblk->y0 - band->y0;
1112                                                         z = pres->z1 - pres->z0 + cblk->z0 - band->z0;
1113                                                 }
1114
1115                                                 if (tcp->tccps[compno].roishift) {
1116                                                         int thresh, val, mag;
1117                                                         thresh = 1 << tcp->tccps[compno].roishift;
1118                                                         for (k = 0; k < cblk->z1 - cblk->z0; k++) {
1119                                                                 for (j = 0; j < cblk->y1 - cblk->y0; j++) {
1120                                                                         for (i = 0; i < cblk->x1 - cblk->x0; i++) {
1121                                                                                 val = t1->data[k][j][i];
1122                                                                                 mag = int_abs(val);
1123                                                                                 if (mag >= thresh) {
1124                                                                                         mag >>= tcp->tccps[compno].roishift;
1125                                                                                         t1->data[k][j][i] = val < 0 ? -mag : mag;
1126                                                                                 }
1127                                                                         }
1128                                                                 }
1129                                                         }
1130                                                 }
1131
1132                                                 if (tcp->tccps[compno].reversible == 1) {
1133                                                         for (k = 0; k < cblk->z1 - cblk->z0; k++) {
1134                                                                 for (j = 0; j < cblk->y1 - cblk->y0; j++) {
1135                                                                         for (i = 0; i < cblk->x1 - cblk->x0; i++) {
1136                                                                                 int tmp = t1->data[k][j][i];
1137                                                                                 tilec->data[x + i + (y + j) * (tilec->x1 - tilec->x0) + (z + k) * (tilec->x1 - tilec->x0) * (tilec->y1 - tilec->y0)] = tmp/2;
1138                                                                         }
1139                                                                 }
1140                                                         }
1141                                                 } else {                /* if (tcp->tccps[compno].reversible == 0) */
1142                                                         for (k = 0; k < cblk->z1 - cblk->z0; k++) {
1143                                                                 for (j = 0; j < cblk->y1 - cblk->y0; j++) {
1144                                                                         for (i = 0; i < cblk->x1 - cblk->x0; i++) {
1145                                                                                 double tmp = (double)(t1->data[k][j][i] * band->stepsize * 4096.0);
1146                                                                                 if (t1->data[k][j][i] >> 1 == 0) {
1147                                                                                         tilec->data[x + i + (y + j) * (tilec->x1 - tilec->x0) + (z + k) * (tilec->x1 - tilec->x0) * (tilec->y1 - tilec->y0)] = 0;
1148                                                                                 } else {
1149                                                                                         int tmp2 = ((int) (floor(fabs(tmp)))) + ((int) floor(fabs(tmp*2))%2);
1150                                                                                         tilec->data[x + i + (y + j) * (tilec->x1 - tilec->x0) + (z + k) * (tilec->x1 - tilec->x0) * (tilec->y1 - tilec->y0)] = ((tmp<0)?-tmp2:tmp2);
1151                                                                                 }
1152                                                                         }
1153                                                                 }
1154                                                         }
1155                                                 }
1156                                         } /* cblkno */
1157                                 } /* precno */
1158                         } /* bandno */
1159                 } /* resno */
1160         } /* compno */
1161 }
1162
1163
1164 /** mod fixed_quality */
1165 double t1_getwmsedec(int nmsedec, int compno, int level[3], int orient, int bpno, double stepsize, int numcomps, int dwtid[3])  {
1166         double w1, w2, wmsedec;
1167
1168         if (dwtid[0] == 1 || dwtid[1] == 1 || dwtid[2] == 1) {
1169                 w1 = (numcomps > 1) ? mct_getnorm_real(compno) : 1;
1170         } else {
1171                 w1 = (numcomps > 1) ? mct_getnorm(compno) : 1;
1172         }
1173         w2 = dwt_getnorm(orient, level, dwtid);
1174
1175         /*fprintf(stdout,"nmsedec %d level %d %d %d orient %d bpno %d stepsize %f \n",nmsedec ,level[0],level[1],level[2],orient,bpno,stepsize);*/
1176         wmsedec = w1 * w2 * stepsize * (1 << bpno);
1177         wmsedec *= wmsedec * nmsedec / 8192.0;
1178
1179         return wmsedec;
1180 }
1181 /** mod fixed_quality */