libopenjpeg/t1.c

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