libopenjpeg/t1.c

   1 /*
   2  * Copyright (c) 2001-2002, David Janssens
   3  * Copyright (c) 2002-2003, Yannick Verschueren
   4  * Copyright (c) 2002-2003,  Communications and remote sensing Laboratory, Universite catholique de Louvain, Belgium
   5  * All rights reserved.
   6  *
   7  * Redistribution and use in source and binary forms, with or without
   8  * modification, are permitted provided that the following conditions
   9  * are met:
  10  * 1. Redistributions of source code must retain the above copyright
  11  *    notice, this list of conditions and the following disclaimer.
  12  * 2. Redistributions in binary form must reproduce the above copyright
  13  *    notice, this list of conditions and the following disclaimer in the
  14  *    documentation and/or other materials provided with the distribution.
  15  *
  16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS `AS IS'
  17  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  18  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  19  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  20  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  21  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  22  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  23  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  24  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  25  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  26  * POSSIBILITY OF SUCH DAMAGE.
  27  */
  28
  29 #include "t1.h"
  30 #include "j2k.h"
  31 #include "mqc.h"
  32 #include "raw.h"                                                                /* Antonin */
  33 #include "int.h"
  34 #include "mct.h"
  35 #include "dwt.h"
  36 #include "fix.h"
  37 #include <stdio.h>
  38 #include <math.h>
  39 #include <stdlib.h>
  40
  41 #define T1_MAXCBLKW 1024
  42 #define T1_MAXCBLKH 1024
  43
  44 #define T1_SIG_NE 0x0001
  45 #define T1_SIG_SE 0x0002
  46 #define T1_SIG_SW 0x0004
  47 #define T1_SIG_NW 0x0008
  48 #define T1_SIG_N 0x0010
  49 #define T1_SIG_E 0x0020
  50 #define T1_SIG_S 0x0040
  51 #define T1_SIG_W 0x0080
  52 #define T1_SIG_OTH (T1_SIG_N|T1_SIG_NE|T1_SIG_E|T1_SIG_SE|T1_SIG_S|T1_SIG_SW|T1_SIG_W|T1_SIG_NW)
  53 #define T1_SIG_PRIM (T1_SIG_N|T1_SIG_E|T1_SIG_S|T1_SIG_W)
  54
  55 #define T1_SGN_N 0x0100
  56 #define T1_SGN_E 0x0200
  57 #define T1_SGN_S 0x0400
  58 #define T1_SGN_W 0x0800
  59 #define T1_SGN (T1_SGN_N|T1_SGN_E|T1_SGN_S|T1_SGN_W)
  60
  61 #define T1_SIG 0x1000
  62 #define T1_REFINE 0x2000
  63 #define T1_VISIT 0x4000
  64
  65 #define T1_NUMCTXS_AGG 1
  66 #define T1_NUMCTXS_ZC 9
  67 #define T1_NUMCTXS_MAG 3
  68 #define T1_NUMCTXS_SC 5
  69 #define T1_NUMCTXS_UNI 1
  70
  71 #define T1_CTXNO_AGG 0
  72 #define T1_CTXNO_ZC (T1_CTXNO_AGG+T1_NUMCTXS_AGG)
  73 #define T1_CTXNO_MAG (T1_CTXNO_ZC+T1_NUMCTXS_ZC)
  74 #define T1_CTXNO_SC (T1_CTXNO_MAG+T1_NUMCTXS_MAG)
  75 #define T1_CTXNO_UNI (T1_CTXNO_SC+T1_NUMCTXS_SC)
  76 #define T1_NUMCTXS (T1_CTXNO_UNI+T1_NUMCTXS_UNI)
  77
  78 #define T1_NMSEDEC_BITS 7
  79 #define T1_NMSEDEC_FRACBITS (T1_NMSEDEC_BITS-1)
  80
  81 /* add TONY */
  82 #define T1_TYPE_MQ 0
  83 #define T1_TYPE_RAW 1
  84 /* dda */
  85
  86 static int t1_lut_ctxno_zc[1024];
  87 static int t1_lut_ctxno_sc[256];
  88 static int t1_lut_ctxno_mag[4096];
  89 static int t1_lut_spb[256];
  90 static int t1_lut_nmsedec_sig[1 << T1_NMSEDEC_BITS];
  91 static int t1_lut_nmsedec_sig0[1 << T1_NMSEDEC_BITS];
  92 static int t1_lut_nmsedec_ref[1 << T1_NMSEDEC_BITS];
  93 static int t1_lut_nmsedec_ref0[1 << T1_NMSEDEC_BITS];
  94
  95 static int t1_data[T1_MAXCBLKH][T1_MAXCBLKH];
  96 static int t1_flags[T1_MAXCBLKH + 2][T1_MAXCBLKH + 2];
  97
  98 int t1_getctxno_zc(int f, int orient)
  99 {
 100         return t1_lut_ctxno_zc[(orient << 8) | (f & T1_SIG_OTH)];
 101 }
 102
 103 int t1_getctxno_sc(int f)
 104 {
 105         return t1_lut_ctxno_sc[(f & (T1_SIG_PRIM | T1_SGN)) >> 4];
 106 }
 107
 108 int t1_getctxno_mag(int f)
 109 {
 110         return t1_lut_ctxno_mag[(f & T1_SIG_OTH) | (((f & T1_REFINE) != 0) << 11)];
 111 }
 112
 113 int t1_getspb(int f)
 114 {
 115         return t1_lut_spb[(f & (T1_SIG_PRIM | T1_SGN)) >> 4];
 116 }
 117
 118 int t1_getnmsedec_sig(int x, int bitpos)
 119 {
 120         if (bitpos > T1_NMSEDEC_FRACBITS)
 121                 return t1_lut_nmsedec_sig[(x >> (bitpos - T1_NMSEDEC_FRACBITS)) & ((1 << T1_NMSEDEC_BITS) - 1)];
 122         else
 123                 return t1_lut_nmsedec_sig0[x & ((1 << T1_NMSEDEC_BITS) - 1)];
 124 }
 125
 126 int t1_getnmsedec_ref(int x, int bitpos)
 127 {
 128         if (bitpos > T1_NMSEDEC_FRACBITS)
 129                 return t1_lut_nmsedec_ref[(x >> (bitpos - T1_NMSEDEC_FRACBITS)) & ((1 << T1_NMSEDEC_BITS) - 1)];
 130         else
 131                 return t1_lut_nmsedec_ref0[x & ((1 << T1_NMSEDEC_BITS) - 1)];
 132 }
 133
 134 void t1_updateflags(int *fp, int s)
 135 {
 136         int *np = fp - (T1_MAXCBLKW + 2);
 137         int *sp = fp + (T1_MAXCBLKW + 2);
 138         np[-1] |= T1_SIG_SE;
 139         np[1] |= T1_SIG_SW;
 140         sp[-1] |= T1_SIG_NE;
 141         sp[1] |= T1_SIG_NW;
 142         *np |= T1_SIG_S;
 143         *sp |= T1_SIG_N;
 144         fp[-1] |= T1_SIG_E;
 145         fp[1] |= T1_SIG_W;
 146         if (s) {
 147                 *np |= T1_SGN_S;
 148                 *sp |= T1_SGN_N;
 149                 fp[-1] |= T1_SGN_E;
 150                 fp[1] |= T1_SGN_W;
 151         }
 152 }
 153
 154 void t1_enc_sigpass_step(int *fp, int *dp, int orient, int bpno, int one, int *nmsedec, char type, int vsc)
 155 {
 156         int v, flag;
 157         flag = vsc ? ((*fp) & (~(T1_SIG_S | T1_SIG_SE | T1_SIG_SW | T1_SGN_S))) : (*fp);
 158         if ((flag & T1_SIG_OTH) && !(flag & (T1_SIG | T1_VISIT))) {
 159                 v = int_abs(*dp) & one ? 1 : 0;
 160                 if (type == T1_TYPE_RAW) {      /* BYPASS/LAZY MODE */
 161                         mqc_setcurctx(t1_getctxno_zc(flag, orient));    /* ESSAI */
 162                         mqc_bypass_enc(v);
 163                 } else {
 164                         mqc_setcurctx(t1_getctxno_zc(flag, orient));
 165                         mqc_encode(v);
 166                 }
 167                 if (v) {
 168                         v = *dp < 0 ? 1 : 0;
 169                         *nmsedec += t1_getnmsedec_sig(int_abs(*dp), bpno + T1_NMSEDEC_FRACBITS);
 170                         if (type == T1_TYPE_RAW) {      /* BYPASS/LAZY MODE */
 171                                 mqc_setcurctx(t1_getctxno_sc(flag));    /* ESSAI */
 172                                 mqc_bypass_enc(v);
 173                         } else {
 174                                 mqc_setcurctx(t1_getctxno_sc(flag));
 175                                 mqc_encode(v ^ t1_getspb(flag));
 176                         }
 177                         t1_updateflags(fp, v);
 178                         *fp |= T1_SIG;
 179                 }
 180                 *fp |= T1_VISIT;
 181         }
 182 }
 183
 184 void t1_dec_sigpass_step(int *fp, int *dp, int orient, int oneplushalf, char type, int vsc)
 185 {
 186         int v, flag;
 187         flag =
 188                 vsc ? ((*fp) & (~(T1_SIG_S | T1_SIG_SE | T1_SIG_SW | T1_SGN_S)))
 189                 : (*fp);
 190         if ((flag & T1_SIG_OTH) && !(flag & (T1_SIG | T1_VISIT))) {
 191                 if (type == T1_TYPE_RAW) {
 192                         if (raw_decode()) {
 193                                 v = raw_decode();                               /* ESSAI */
 194                                 *dp = v ? -oneplushalf : oneplushalf;
 195                                 t1_updateflags(fp, v);
 196                                 *fp |= T1_SIG;
 197                         }
 198                 } else {
 199                         mqc_setcurctx(t1_getctxno_zc(flag, orient));
 200                         if (mqc_decode()) {
 201                                 mqc_setcurctx(t1_getctxno_sc(flag));
 202                                 v = mqc_decode() ^ t1_getspb(flag);
 203                                 *dp = v ? -oneplushalf : oneplushalf;
 204                                 t1_updateflags(fp, v);
 205                                 *fp |= T1_SIG;
 206                         }
 207                 }
 208                 *fp |= T1_VISIT;
 209         }
 210 }               /* VSC and  BYPASS by Antonin */
 211
 212 void t1_enc_sigpass(int w, int h, int bpno, int orient, int *nmsedec, char type, int cblksty)
 213 {
 214         int i, j, k, one, vsc;
 215         *nmsedec = 0;
 216         one = 1 << (bpno + T1_NMSEDEC_FRACBITS);
 217         for (k = 0; k < h; k += 4) {
 218                 for (i = 0; i < w; i++) {
 219                         for (j = k; j < k + 4 && j < h; j++) {
 220                                 vsc = ((cblksty & J2K_CCP_CBLKSTY_VSC) && (j == k + 3 || j == h - 1)) ? 1 : 0;
 221                                 t1_enc_sigpass_step(&t1_flags[1 + j][1 + i], &t1_data[j][i], orient, bpno, one, nmsedec, type, vsc);
 222                         }
 223                 }
 224         }
 225 }
 226
 227 void t1_dec_sigpass(int w, int h, int bpno, int orient, char type, int cblksty)
 228 {
 229         int i, j, k, one, half, oneplushalf, vsc;
 230         one = 1 << bpno;
 231         half = one >> 1;
 232         oneplushalf = one | half;
 233         for (k = 0; k < h; k += 4) {
 234           for (i = 0; i < w; i++) {
 235             for (j = k; j < k + 4 && j < h; j++) {
 236               vsc = ((cblksty & J2K_CCP_CBLKSTY_VSC)
 237                      && (j == k + 3 || j == h - 1)) ? 1 : 0;
 238               t1_dec_sigpass_step(&t1_flags[1 + j][1 + i], &t1_data[j][i], orient, oneplushalf, type, vsc);
 239             }
 240           }
 241         }
 242 }               /* VSC and  BYPASS by Antonin */
 243
 244 void t1_enc_refpass_step(int *fp, int *dp, int bpno, int one, int *nmsedec, char type, int vsc)
 245 {
 246         int v, flag;
 247         flag = vsc ? ((*fp) & (~(T1_SIG_S | T1_SIG_SE | T1_SIG_SW | T1_SGN_S))) : (*fp);
 248         if ((flag & (T1_SIG | T1_VISIT)) == T1_SIG) {
 249                 *nmsedec += t1_getnmsedec_ref(int_abs(*dp), bpno + T1_NMSEDEC_FRACBITS);
 250                 v = int_abs(*dp) & one ? 1 : 0;
 251                 if (type == T1_TYPE_RAW) {      /* BYPASS/LAZY MODE */
 252                         mqc_setcurctx(t1_getctxno_mag(flag));   /* ESSAI */
 253                         mqc_bypass_enc(v);
 254                 } else {
 255                         mqc_setcurctx(t1_getctxno_mag(flag));
 256                         mqc_encode(v);
 257                 }
 258                 *fp |= T1_REFINE;
 259         }
 260 }
 261
 262 void t1_dec_refpass_step(int *fp, int *dp, int poshalf, int neghalf, char type, int vsc)
 263 {
 264         int v, t, flag;
 265         flag = vsc ? ((*fp) & (~(T1_SIG_S | T1_SIG_SE | T1_SIG_SW | T1_SGN_S))) : (*fp);
 266         if ((flag & (T1_SIG | T1_VISIT)) == T1_SIG) {
 267                 if (type == T1_TYPE_RAW) {
 268                         mqc_setcurctx(t1_getctxno_mag(flag));   /* ESSAI */
 269                         v = raw_decode();
 270                 } else {
 271                         mqc_setcurctx(t1_getctxno_mag(flag));
 272                         v = mqc_decode();
 273                 }
 274                 t = v ? poshalf : neghalf;
 275                 *dp += *dp < 0 ? -t : t;
 276                 *fp |= T1_REFINE;
 277         }
 278 }       /* VSC and  BYPASS by Antonin  */
 279
 280 void t1_enc_refpass(int w, int h, int bpno, int *nmsedec, char type, int cblksty)
 281 {
 282         int i, j, k, one, vsc;
 283         *nmsedec = 0;
 284         one = 1 << (bpno + T1_NMSEDEC_FRACBITS);
 285         for (k = 0; k < h; k += 4) {
 286                 for (i = 0; i < w; i++) {
 287                         for (j = k; j < k + 4 && j < h; j++) {
 288                                 vsc = ((cblksty & J2K_CCP_CBLKSTY_VSC) && (j == k + 3 || j == h - 1)) ? 1 : 0;
 289                                 t1_enc_refpass_step(&t1_flags[1 + j][1 + i], &t1_data[j][i], bpno, one, nmsedec, type, vsc);
 290                         }
 291                 }
 292         }
 293 }
 294
 295 void t1_dec_refpass(int w, int h, int bpno, char type, int cblksty)
 296 {
 297         int i, j, k, one, poshalf, neghalf;
 298         int vsc;
 299         one = 1 << bpno;
 300         poshalf = one >> 1;
 301         neghalf = bpno > 0 ? -poshalf : -1;
 302         for (k = 0; k < h; k += 4) {
 303                 for (i = 0; i < w; i++) {
 304                         for (j = k; j < k + 4 && j < h; j++) {
 305                                 vsc = ((cblksty & J2K_CCP_CBLKSTY_VSC)  && (j == k + 3 || j == h - 1)) ? 1 : 0;
 306                                 t1_dec_refpass_step(&t1_flags[1 + j][1 + i], &t1_data[j][i], poshalf, neghalf, type, vsc);
 307                         }
 308                 }
 309         }
 310 }       /* VSC and  BYPASS by Antonin */
 311
 312 void t1_enc_clnpass_step(int *fp, int *dp, int orient, int bpno, int one, int *nmsedec, int partial, int vsc)
 313 {
 314         int v, flag;
 315         flag = vsc ? ((*fp) & (~(T1_SIG_S | T1_SIG_SE | T1_SIG_SW | T1_SGN_S))) : (*fp);
 316         if (partial)
 317                 goto label_partial;
 318         if (!(*fp & (T1_SIG | T1_VISIT))) {
 319                 mqc_setcurctx(t1_getctxno_zc(flag, orient));
 320                 v = int_abs(*dp) & one ? 1 : 0;
 321                 mqc_encode(v);
 322                 if (v) {
 323                 label_partial:
 324                         *nmsedec += t1_getnmsedec_sig(int_abs(*dp), bpno + T1_NMSEDEC_FRACBITS);
 325                         mqc_setcurctx(t1_getctxno_sc(flag));
 326                         v = *dp < 0 ? 1 : 0;
 327                         mqc_encode(v ^ t1_getspb(flag));
 328                         t1_updateflags(fp, v);
 329                         *fp |= T1_SIG;
 330                 }
 331         }
 332         *fp &= ~T1_VISIT;
 333 }
 334
 335 void t1_dec_clnpass_step(int *fp, int *dp, int orient, int oneplushalf, int partial, int vsc)
 336 {
 337         int v, flag;
 338         flag = vsc ? ((*fp) & (~(T1_SIG_S | T1_SIG_SE | T1_SIG_SW | T1_SGN_S))) : (*fp);
 339         if (partial)
 340                 goto label_partial;
 341         if (!(flag & (T1_SIG | T1_VISIT))) {
 342                 mqc_setcurctx(t1_getctxno_zc(flag, orient));
 343                 if (mqc_decode()) {
 344                 label_partial:
 345                         mqc_setcurctx(t1_getctxno_sc(flag));
 346                         v = mqc_decode() ^ t1_getspb(flag);
 347                         *dp = v ? -oneplushalf : oneplushalf;
 348                         t1_updateflags(fp, v);
 349                         *fp |= T1_SIG;
 350                 }
 351         }
 352         *fp &= ~T1_VISIT;
 353 }       /* VSC and  BYPASS by Antonin */
 354
 355 void t1_enc_clnpass(int w, int h, int bpno, int orient, int *nmsedec, int cblksty)
 356 {
 357         int i, j, k, one, agg, runlen, vsc;
 358         *nmsedec = 0;
 359         one = 1 << (bpno + T1_NMSEDEC_FRACBITS);
 360         for (k = 0; k < h; k += 4) {
 361                 for (i = 0; i < w; i++) {
 362                         if (k + 3 < h) {
 363                                 if (cblksty & J2K_CCP_CBLKSTY_VSC) {
 364                                         agg = !(t1_flags[1 + k][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH)
 365                                                 || t1_flags[1 + k + 1][1 + i] & (T1_SIG | T1_VISIT |
 366                                                 T1_SIG_OTH) || t1_flags[1 + k + 2][1 + i] & (T1_SIG | T1_VISIT |
 367                                                 T1_SIG_OTH) || (t1_flags[1 + k + 3][1 + i] & (~(T1_SIG_S | T1_SIG_SE
 368                                                 | T1_SIG_SW | T1_SGN_S))) & (T1_SIG | T1_VISIT | T1_SIG_OTH));
 369                                 } else {
 370                                         agg = !(t1_flags[1 + k][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH)
 371                                                 || t1_flags[1 + k + 1][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH)
 372                                                 || t1_flags[1 + k + 2][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH)
 373                                                 || t1_flags[1 + k + 3][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH));
 374                                 }
 375                         } else {
 376                                 agg = 0;
 377                         }
 378                         if (agg) {
 379                                 for (runlen = 0; runlen < 4; runlen++) {
 380                                         if (int_abs(t1_data[k + runlen][i]) & one)
 381                                                 break;
 382                                 }
 383                                 mqc_setcurctx(T1_CTXNO_AGG);
 384                                 mqc_encode(runlen != 4);
 385                                 if (runlen == 4) {
 386                                         continue;
 387                                 }
 388                                 mqc_setcurctx(T1_CTXNO_UNI);
 389                                 mqc_encode(runlen >> 1);
 390                                 mqc_encode(runlen & 1);
 391                         } else {
 392                                 runlen = 0;
 393                         }
 394                         for (j = k + runlen; j < k + 4 && j < h; j++) {
 395                                 vsc = ((cblksty & J2K_CCP_CBLKSTY_VSC) && (j == k + 3 || j == h - 1)) ? 1 : 0;
 396                                 t1_enc_clnpass_step(&t1_flags[1 + j][1 + i], &t1_data[j][i], orient, bpno, one, nmsedec, agg
 397                                                     && (j == k + runlen), vsc);
 398                         }
 399                 }
 400         }
 401 }
 402
 403 void t1_dec_clnpass(int w, int h, int bpno, int orient, int cblksty)
 404 {
 405         int i, j, k, one, half, oneplushalf, agg, runlen, vsc;
 406         int segsym = cblksty & J2K_CCP_CBLKSTY_SEGSYM;
 407         one = 1 << bpno;
 408         half = one >> 1;
 409         oneplushalf = one | half;
 410         for (k = 0; k < h; k += 4) {
 411                 for (i = 0; i < w; i++) {
 412                         if (k + 3 < h) {
 413                                 if (cblksty & J2K_CCP_CBLKSTY_VSC) {
 414                                         agg = !(t1_flags[1 + k][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH)
 415                                                 || t1_flags[1 + k + 1][1 + i] & (T1_SIG | T1_VISIT |
 416                                                 T1_SIG_OTH) || t1_flags[1 + k + 2][1 + i] & (T1_SIG | T1_VISIT |
 417                                                 T1_SIG_OTH) || (t1_flags[1 + k + 3][1 + i] & (~(T1_SIG_S | T1_SIG_SE |
 418                                                 T1_SIG_SW | T1_SGN_S))) & (T1_SIG | T1_VISIT | T1_SIG_OTH));
 419                                 } else {
 420                                         agg = !(t1_flags[1 + k][1 + i] & (T1_SIG | T1_VISIT | T1_SIG_OTH)
 421                                                 || t1_flags[1 + k + 1][1 + i] & (T1_SIG | T1_VISIT |
 422                                                 T1_SIG_OTH) || t1_flags[1 + k + 2][1 + i] & (T1_SIG | T1_VISIT |
 423                                                 T1_SIG_OTH) || t1_flags[1 + k + 3][1 + i] & (T1_SIG | T1_VISIT |
 424                                                 T1_SIG_OTH));
 425                                 }
 426                         } else {
 427                                 agg = 0;
 428                         }
 429                         if (agg) {
 430                                 mqc_setcurctx(T1_CTXNO_AGG);
 431                                 if (!mqc_decode()) {
 432                                         continue;
 433                                 }
 434                                 mqc_setcurctx(T1_CTXNO_UNI);
 435                                 runlen = mqc_decode();
 436                                 runlen = (runlen << 1) | mqc_decode();
 437                         } else {
 438                                 runlen = 0;
 439                         }
 440                         for (j = k + runlen; j < k + 4 && j < h; j++) {
 441                                 vsc = ((cblksty & J2K_CCP_CBLKSTY_VSC) && (j == k + 3 || j == h - 1)) ? 1 : 0;
 442                                 t1_dec_clnpass_step(&t1_flags[1 + j][1 + i], &t1_data[j][i], orient, oneplushalf,
 443                                                     agg && (j == k + runlen), vsc);
 444                         }
 445                 }
 446         }
 447         if (segsym) {
 448                 int v = 0;
 449                 mqc_setcurctx(T1_CTXNO_UNI);
 450                 v = mqc_decode();
 451                 v = (v << 1) | mqc_decode();
 452                 v = (v << 1) | mqc_decode();
 453                 v = (v << 1) | mqc_decode();
 454                 /*  if (v!=0xa)
 455                    {
 456                    fprintf(stderr, "warning: bad segmentation symbol %x\n",v);
 457                    } */
 458         }
 459 }       /* VSC and  BYPASS by Antonin */
 460
 461 double t1_getwmsedec(int nmsedec, int compno, int level, int orient, int bpno, int qmfbid, double stepsize)
 462 {
 463         double w1, w2, wmsedec;
 464         if (qmfbid == 1) {
 465                 w1 = mct_getnorm(compno);
 466                 w2 = dwt_getnorm(level, orient);
 467         } else {  /* if (qmfbid == 0) */
 468                 w1 = mct_getnorm_real(compno);
 469                 w2 = dwt_getnorm_real(level, orient);
 470         }
 471         wmsedec = w1 * w2 * stepsize * (1 << bpno);
 472         wmsedec *= wmsedec * nmsedec / 8192.0;
 473         return wmsedec;
 474 }
 475
 476 void t1_encode_cblk(tcd_cblk_t * cblk, int orient, int compno, int level, int qmfbid, double stepsize, int cblksty)
 477 {
 478         int i, j;
 479         int w, h;
 480         int passno;
 481         int bpno, passtype;
 482         int max;
 483         int nmsedec;
 484         double cumwmsedec = 0;
 485         char type = T1_TYPE_MQ;
 486
 487         w = cblk->x1 - cblk->x0;
 488         h = cblk->y1 - cblk->y0;
 489
 490         max = 0;
 491         for (j = 0; j < h; j++) {
 492           for (i = 0; i < w; i++) {
 493             max = int_max(max, int_abs(t1_data[j][i]));
 494           }
 495         }
 496
 497         cblk->numbps = max ? (int_floorlog2(max) + 1) - T1_NMSEDEC_FRACBITS : 0;
 498
 499         for (i = 0; i < sizeof(t1_flags) / sizeof(int); i++)
 500                 ((int *) t1_flags)[i] = 0;
 501         bpno = cblk->numbps - 1;
 502         passtype = 2;
 503
 504         mqc_resetstates();
 505         mqc_setstate(T1_CTXNO_UNI, 0, 46);
 506         mqc_setstate(T1_CTXNO_AGG, 0, 3);
 507         mqc_setstate(T1_CTXNO_ZC, 0, 4);
 508         mqc_init_enc(cblk->data);
 509
 510         for (passno = 0; bpno >= 0; passno++) {
 511                 tcd_pass_t *pass = &cblk->passes[passno];
 512                 int correction = 3;
 513                 type = ((bpno < (cblk->numbps - 4)) && (passtype < 2)
 514                         && (cblksty & J2K_CCP_CBLKSTY_LAZY)) ? T1_TYPE_RAW : T1_TYPE_MQ;
 515
 516                 switch (passtype) {
 517                 case 0:
 518                         t1_enc_sigpass(w, h, bpno, orient, &nmsedec, type, cblksty);
 519                         break;
 520                 case 1:
 521                         t1_enc_refpass(w, h, bpno, &nmsedec, type, cblksty);
 522                         break;
 523                 case 2:
 524                         t1_enc_clnpass(w, h, bpno, orient, &nmsedec, cblksty);
 525                         /* code switch SEGMARK (i.e. SEGSYM) */
 526                         if (cblksty & J2K_CCP_CBLKSTY_SEGSYM)
 527                                 mqc_segmark_enc();
 528                         break;
 529                 }
 530
 531                 cumwmsedec += t1_getwmsedec(nmsedec, compno, level, orient, bpno, qmfbid, stepsize);
 532
 533                 /* Code switch "RESTART" (i.e. TERMALL) */
 534                 if ((cblksty & J2K_CCP_CBLKSTY_TERMALL) && !((passtype == 2) && (bpno - 1 < 0)))
 535                   {
 536                     if (type == T1_TYPE_RAW)
 537                       {
 538                         mqc_flush();
 539                         correction=1;
 540                         /* correction = mqc_bypass_flush_enc();*/
 541                       }
 542                     else
 543                       {         /* correction = mqc_restart_enc(); */
 544                         mqc_flush();
 545                         correction = 1;
 546                       }
 547                     pass->term = 1;
 548                   } else
 549                     {
 550                       if (((bpno < (cblk->numbps - 4) && (passtype > 0)) || ((bpno == (cblk->numbps - 4)) && (passtype == 2)))
 551                           && (cblksty & J2K_CCP_CBLKSTY_LAZY))
 552                         {
 553                           if (type == T1_TYPE_RAW)
 554                             {
 555                               mqc_flush();
 556                               correction=1;
 557                               /* correction = mqc_bypass_flush_enc();*/
 558                             }
 559                           else
 560                             {  /* correction = mqc_restart_enc();*/
 561                               mqc_flush();
 562                               correction=1;
 563                             }
 564                           pass->term = 1;
 565                         } else
 566                           {
 567                             pass->term = 0;
 568                           }
 569                     }
 570
 571                 if (++passtype == 3)
 572                   {
 573                     passtype = 0;
 574                     bpno--;
 575                   }
 576
 577                 if (pass->term && bpno > 0) {
 578                         type = ((bpno < (cblk->numbps - 4)) && (passtype < 2)
 579                                 && (cblksty & J2K_CCP_CBLKSTY_LAZY)) ? T1_TYPE_RAW : T1_TYPE_MQ;
 580                         if (type == T1_TYPE_RAW)
 581                                 mqc_bypass_init_enc();
 582                         else
 583                                 mqc_restart_init_enc();
 584                 }
 585
 586                 pass->distortiondec = cumwmsedec;
 587                 pass->rate = mqc_numbytes() + correction;       /* FIXME */
 588                 pass->len = pass->rate - (passno == 0 ? 0 : cblk->passes[passno - 1].rate);
 589
 590                 /* Code-switch "RESET" */
 591                 if (cblksty & J2K_CCP_CBLKSTY_RESET)
 592                         mqc_reset_enc();
 593         }
 594
 595         /* Code switch "ERTERM" (i.e. PTERM) */
 596         if (cblksty & J2K_CCP_CBLKSTY_PTERM)
 597                 mqc_erterm_enc();
 598         else /* Default coding */ if (!(cblksty & J2K_CCP_CBLKSTY_LAZY))
 599                 mqc_flush();
 600
 601         cblk->totalpasses = passno;
 602 }
 603
 604 void t1_decode_cblk(tcd_cblk_t * cblk, int orient, int roishift, int cblksty)
 605 {
 606         int i;
 607         int w, h;
 608         int bpno, passtype;
 609         int segno, passno;
 610         /* add TONY */
 611         char type = T1_TYPE_MQ;
 612         /* dda */
 613
 614         for (i = 0; i < sizeof(t1_data) / sizeof(int); i++)
 615                 ((int *) t1_data)[i] = 0;
 616         for (i = 0; i < sizeof(t1_flags) / sizeof(int); i++)
 617                 ((int *) t1_flags)[i] = 0;
 618
 619         w = cblk->x1 - cblk->x0;
 620         h = cblk->y1 - cblk->y0;
 621         bpno = roishift + cblk->numbps - 1;
 622         passtype = 2;
 623
 624         mqc_resetstates();
 625         mqc_setstate(T1_CTXNO_UNI, 0, 46);
 626         mqc_setstate(T1_CTXNO_AGG, 0, 3);
 627         mqc_setstate(T1_CTXNO_ZC, 0, 4);
 628
 629         for (segno = 0; segno < cblk->numsegs; segno++) {
 630                 tcd_seg_t *seg = &cblk->segs[segno];
 631
 632                 /* add TONY */
 633                 type = ((bpno <= (cblk->numbps - 1) - 4) && (passtype < 2)
 634                         && (cblksty & J2K_CCP_CBLKSTY_LAZY)) ? T1_TYPE_RAW : T1_TYPE_MQ;
 635                 if (type == T1_TYPE_RAW)
 636                         raw_init_dec(seg->data, seg->len);
 637                 else
 638                         mqc_init_dec(seg->data, seg->len);
 639                 /* dda */
 640
 641                 for (passno = 0; passno < seg->numpasses; passno++) {
 642                         switch (passtype) {
 643                         case 0:
 644                                 t1_dec_sigpass(w, h, bpno, orient, type, cblksty);
 645                                 break;
 646                         case 1:
 647                                 t1_dec_refpass(w, h, bpno, type, cblksty);
 648                                 break;
 649                         case 2:
 650                                 t1_dec_clnpass(w, h, bpno, orient, cblksty);
 651                                 break;
 652                         }
 653
 654                         if ((cblksty & J2K_CCP_CBLKSTY_RESET) && type == T1_TYPE_MQ)
 655                                 mqc_reset_enc();
 656
 657                         if (++passtype == 3) {
 658                                 passtype = 0;
 659                                 bpno--;
 660                         }
 661                 }
 662         }
 663 }
 664
 665 void t1_encode_cblks(tcd_tile_t * tile, j2k_tcp_t * tcp)
 666 {
 667   int compno, resno, bandno, precno, cblkno;
 668   int x, y, i, j, orient;
 669   tcd_tilecomp_t *tilec;
 670   tcd_resolution_t *res;
 671   tcd_band_t *band;
 672   tcd_precinct_t *prc;
 673   tcd_cblk_t *cblk;
 674
 675   for (compno = 0; compno < tile->numcomps; compno++) {
 676     tilec = &tile->comps[compno];
 677     for (resno = 0; resno < tilec->numresolutions; resno++) {
 678       res = &tilec->resolutions[resno];
 679       for (bandno = 0; bandno < res->numbands; bandno++) {
 680         band = &res->bands[bandno];
 681         for (precno = 0; precno < res->pw * res->ph; precno++) {
 682           prc = &band->precincts[precno];
 683           for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++) {
 684             cblk = &prc->cblks[cblkno];
 685
 686             if (band->bandno == 0) {
 687               x = cblk->x0 - band->x0;
 688               y = cblk->y0 - band->y0;
 689             } else if (band->bandno == 1) {
 690               tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
 691               x = pres->x1 - pres->x0 + cblk->x0 - band->x0;
 692               y = cblk->y0 - band->y0;
 693             } else if (band->bandno == 2) {
 694               tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
 695               x = cblk->x0 - band->x0;
 696               y = pres->y1 - pres->y0 + cblk->y0 - band->y0;
 697             } else {  /* if (band->bandno == 3) */
 698               tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
 699               x = pres->x1 - pres->x0 + cblk->x0 - band->x0;
 700               y = pres->y1 - pres->y0 + cblk->y0 - band->y0;
 701             }
 702
 703             if (tcp->tccps[compno].qmfbid == 1) {
 704
 705               for (j = 0; j < cblk->y1 - cblk->y0; j++) {
 706                 for (i = 0; i < cblk->x1 - cblk->x0; i++) {
 707                   t1_data[j][i] = tilec->data[(x + i) + (y + j) * (tilec->x1 - tilec->x0)] << T1_NMSEDEC_FRACBITS;
 708                 }
 709               }
 710             } else if (tcp->tccps[compno].qmfbid == 0) {
 711               for (j = 0; j < cblk->y1 - cblk->y0; j++) {
 712                 for (i = 0; i < cblk->x1 - cblk->x0; i++) {
 713                   t1_data[j][i] = fix_mul(tilec->data[x + i + (y + j) * (tilec->x1 - tilec->x0)],
 714                                           8192 * 8192 / band->stepsize) >> (13 - T1_NMSEDEC_FRACBITS);
 715                 }
 716               }
 717             }
 718             orient = band->bandno;      /* FIXME */
 719             if (orient == 2) {
 720               orient = 1;
 721             } else if (orient == 1) {
 722               orient = 2;
 723             }
 724             t1_encode_cblk(cblk, orient, compno,
 725                            tilec->numresolutions - 1 - resno,
 726                            tcp->tccps[compno].qmfbid, band->stepsize,
 727                            tcp->tccps[compno].cblksty);
 728           }     /* cblkno */
 729         }       /* precno */
 730       }         /* bandno */
 731     }           /* resno  */
 732   }             /* compo  */
 733 }
 734
 735
 736 void t1_decode_cblks(tcd_tile_t * tile, j2k_tcp_t * tcp)
 737 {
 738   int compno, resno, bandno, precno, cblkno;
 739
 740   for (compno = 0; compno < tile->numcomps; compno++)
 741     {
 742       tcd_tilecomp_t *tilec = &tile->comps[compno];
 743       for (resno = 0; resno < tilec->numresolutions; resno++)
 744         {
 745           tcd_resolution_t *res = &tilec->resolutions[resno];
 746           for (bandno = 0; bandno < res->numbands; bandno++)
 747             {
 748               tcd_band_t *band = &res->bands[bandno];
 749               for (precno = 0; precno < res->pw * res->ph; precno++)
 750                 {
 751                   tcd_precinct_t *prc = &band->precincts[precno];
 752                   for (cblkno = 0; cblkno < prc->cw * prc->ch; cblkno++)
 753                     {
 754                       int x, y, i, j, orient;
 755                       tcd_cblk_t *cblk = &prc->cblks[cblkno];
 756                       orient = band->bandno;    /* FIXME */
 757                       if (orient == 2)
 758                         orient = 1;
 759                       else if (orient == 1)
 760                         orient = 2;
 761                       t1_decode_cblk(cblk, orient, tcp->tccps[compno].roishift, tcp->tccps[compno].cblksty);
 762                       if (band->bandno == 0) {
 763                         x = cblk->x0 - band->x0;
 764                         y = cblk->y0 - band->y0;
 765                       } else if (band->bandno == 1) {
 766                         tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
 767                         x = pres->x1 - pres->x0 + cblk->x0 - band->x0;
 768                         y = cblk->y0 - band->y0;
 769                       } else if (band->bandno == 2) {
 770                         tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
 771                         x = cblk->x0 - band->x0;
 772                         y = pres->y1 - pres->y0 + cblk->y0 - band->y0;
 773                       } else {  /* if (band->bandno == 3) */
 774                         tcd_resolution_t *pres = &tilec->resolutions[resno - 1];
 775                         x = pres->x1 - pres->x0 + cblk->x0 - band->x0;
 776                         y = pres->y1 - pres->y0 + cblk->y0 - band->y0;
 777                       }
 778
 779                       if (tcp->tccps[compno].roishift)
 780                         {
 781                           int thresh, val, mag;
 782                           thresh = 1 << tcp->tccps[compno].roishift;
 783                           for (j = 0; j < cblk->y1 - cblk->y0; j++) {
 784                             for (i = 0; i < cblk->x1 - cblk->x0; i++) {
 785                               val = t1_data[j][i];
 786                               mag = int_abs(val);
 787                               if (mag >= thresh) {
 788                                 mag >>= tcp->tccps[compno].roishift;
 789                                 t1_data[j][i] = val < 0 ? -mag : mag;
 790                               }
 791                             }
 792                           }
 793                         }
 794
 795                       if (tcp->tccps[compno].qmfbid == 1) {
 796                         for (j = 0; j < cblk->y1 - cblk->y0; j++) {
 797                           for (i = 0; i < cblk->x1 - cblk->x0; i++) {
 798                             tilec->data[x + i + (y + j) * (tilec->x1 - tilec->x0)] = t1_data[j][i];
 799                           }
 800                         }
 801                       } else /* if (tcp->tccps[compno].qmfbid == 0)*/ {
 802                         for (j = 0; j < cblk->y1 - cblk->y0; j++) {
 803                           for (i = 0; i < cblk->x1 - cblk->x0; i++) {
 804                             if (t1_data[j][i] == 0) {
 805                               tilec->data[x + i + (y + j) * (tilec->x1 - tilec->x0)] = 0;
 806                             } else {
 807                               tilec->data[x + i + (y + j) * (tilec->x1 - tilec->x0)] = fix_mul(t1_data[j][i] << 13, band->stepsize);
 808                             }
 809                           }
 810                         }
 811                       }
 812                     }
 813                 }
 814             }
 815         }
 816     }
 817 }
 818
 819 int t1_init_ctxno_zc(int f, int orient)
 820 {
 821         int h, v, d, n, t, hv;
 822         n = 0;
 823         h = ((f & T1_SIG_W) != 0) + ((f & T1_SIG_E) != 0);
 824         v = ((f & T1_SIG_N) != 0) + ((f & T1_SIG_S) != 0);
 825         d = ((f & T1_SIG_NW) != 0) + ((f & T1_SIG_NE) != 0) + ((f & T1_SIG_SE) != 0) + ((f & T1_SIG_SW) != 0);
 826         switch (orient) {
 827         case 2:
 828                 t = h;
 829                 h = v;
 830                 v = t;
 831         case 0:
 832         case 1:
 833                 if (!h) {
 834                         if (!v) {
 835                                 if (!d)
 836                                         n = 0;
 837                                 else if (d == 1)
 838                                         n = 1;
 839                                 else
 840                                         n = 2;
 841                         } else if (v == 1)
 842                                 n = 3;
 843                         else
 844                                 n = 4;
 845                 } else if (h == 1) {
 846                         if (!v) {
 847                                 if (!d)
 848                                         n = 5;
 849                                 else
 850                                         n = 6;
 851                         } else
 852                                 n = 7;
 853                 } else
 854                         n = 8;
 855                 break;
 856         case 3:
 857                 hv = h + v;
 858                 if (!d) {
 859                         if (!hv)
 860                                 n = 0;
 861                         else if (hv == 1)
 862                                 n = 1;
 863                         else
 864                                 n = 2;
 865                 } else if (d == 1) {
 866                         if (!hv)
 867                                 n = 3;
 868                         else if (hv == 1)
 869                                 n = 4;
 870                         else
 871                                 n = 5;
 872                 } else if (d == 2) {
 873                         if (!hv)
 874                                 n = 6;
 875                         else
 876                                 n = 7;
 877                 } else
 878                         n = 8;
 879                 break;
 880         }
 881         return T1_CTXNO_ZC + n;
 882 }
 883
 884 int t1_init_ctxno_sc(int f)
 885 {
 886         int hc, vc, n;
 887         n = 0;
 888         hc = int_min(((f & (T1_SIG_E | T1_SGN_E)) == T1_SIG_E) + ((f & (T1_SIG_W | T1_SGN_W)) == T1_SIG_W),
 889                         1) - int_min(((f & (T1_SIG_E | T1_SGN_E)) == (T1_SIG_E | T1_SGN_E)) + ((f & (T1_SIG_W | T1_SGN_W)) ==
 890                         (T1_SIG_W | T1_SGN_W)), 1);
 891         vc = int_min(((f & (T1_SIG_N | T1_SGN_N)) == T1_SIG_N) + ((f & (T1_SIG_S | T1_SGN_S)) == T1_SIG_S),
 892                         1) - int_min(((f & (T1_SIG_N | T1_SGN_N)) == (T1_SIG_N | T1_SGN_N)) + ((f & (T1_SIG_S | T1_SGN_S)) ==
 893                         (T1_SIG_S | T1_SGN_S)), 1);
 894         if (hc < 0) {
 895                 hc = -hc;
 896                 vc = -vc;
 897         }
 898         if (!hc) {
 899                 if (vc == -1)
 900                         n = 1;
 901                 else if (!vc)
 902                         n = 0;
 903                 else
 904                         n = 1;
 905         } else if (hc == 1) {
 906                 if (vc == -1)
 907                         n = 2;
 908                 else if (!vc)
 909                         n = 3;
 910                 else
 911                         n = 4;
 912         }
 913         return T1_CTXNO_SC + n;
 914 }
 915
 916 int t1_init_ctxno_mag(int f)
 917 {
 918         int n;
 919         if (!(f & T1_REFINE))
 920                 n = (f & (T1_SIG_OTH)) ? 1 : 0;
 921         else
 922                 n = 2;
 923         return T1_CTXNO_MAG + n;
 924 }
 925
 926 int t1_init_spb(int f)
 927 {
 928         int hc, vc, n;
 929         hc = int_min(((f & (T1_SIG_E | T1_SGN_E)) == T1_SIG_E) + ((f & (T1_SIG_W | T1_SGN_W)) == T1_SIG_W),
 930                         1) - int_min(((f & (T1_SIG_E | T1_SGN_E)) == (T1_SIG_E | T1_SGN_E)) + ((f & (T1_SIG_W | T1_SGN_W)) ==
 931                         (T1_SIG_W | T1_SGN_W)), 1);
 932         vc = int_min(((f & (T1_SIG_N | T1_SGN_N)) == T1_SIG_N) + ((f & (T1_SIG_S | T1_SGN_S)) == T1_SIG_S),
 933                         1) - int_min(((f & (T1_SIG_N | T1_SGN_N)) == (T1_SIG_N | T1_SGN_N)) + ((f & (T1_SIG_S | T1_SGN_S)) ==
 934                         (T1_SIG_S | T1_SGN_S)), 1);
 935         if (!hc && !vc)
 936                 n = 0;
 937         else
 938                 n = (!(hc > 0 || (!hc && vc > 0)));
 939         return n;
 940 }
 941
 942 void t1_init_luts()
 943 {
 944         int i, j;
 945         double u, v, t;
 946         for (j = 0; j < 4; j++) {
 947                 for (i = 0; i < 256; ++i) {
 948                         t1_lut_ctxno_zc[(j << 8) | i] = t1_init_ctxno_zc(i, j);
 949                 }
 950         }
 951         for (i = 0; i < 256; i++) {
 952                 t1_lut_ctxno_sc[i] = t1_init_ctxno_sc(i << 4);
 953         }
 954         for (j = 0; j < 2; j++) {
 955                 for (i = 0; i < 2048; ++i) {
 956                         t1_lut_ctxno_mag[(j << 11) + i] =
 957                                 t1_init_ctxno_mag((j ? T1_REFINE : 0) | i);
 958                 }
 959         }
 960         for (i = 0; i < 256; ++i) {
 961                 t1_lut_spb[i] = t1_init_spb(i << 4);
 962         }
 963         /* FIXME FIXME FIXME */
 964         /* printf("nmsedec luts:\n"); */
 965         for (i = 0; i < (1 << T1_NMSEDEC_BITS); i++) {
 966                 t = i / pow(2, T1_NMSEDEC_FRACBITS);
 967                 u = t;
 968                 v = t - 1.5;
 969                 t1_lut_nmsedec_sig[i] = int_max(0, (int) (floor((u * u - v * v) * pow(2,T1_NMSEDEC_FRACBITS) +
 970                                                 0.5) / pow(2, T1_NMSEDEC_FRACBITS) * 8192.0));
 971                 t1_lut_nmsedec_sig0[i] = int_max(0, (int) (floor((u * u) * pow(2, T1_NMSEDEC_FRACBITS) + 0.5) /
 972                                                 pow(2, T1_NMSEDEC_FRACBITS) * 8192.0));
 973                 u = t - 1.0;
 974                 if (i & (1 << (T1_NMSEDEC_BITS - 1))) {
 975                         v = t - 1.5;
 976                 } else {
 977                         v = t - 0.5;
 978                 }
 979                 t1_lut_nmsedec_ref[i] = int_max(0, (int) (floor((u * u - v * v) * pow(2,T1_NMSEDEC_FRACBITS) +
 980                                                 0.5) / pow(2, T1_NMSEDEC_FRACBITS) * 8192.0));
 981                 t1_lut_nmsedec_ref0[i] = int_max(0, (int) (floor((u * u) * pow(2, T1_NMSEDEC_FRACBITS) + 0.5) /
 982                                                  pow(2, T1_NMSEDEC_FRACBITS) * 8192.0));
 983         }
 984 }