libopenjp3d/t1.c

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