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