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