src/gdcmjpegls/Encoder/lossy_e.c

   1 /* SPMG/JPEG-LS IMPLEMENTATION V.2.1
   2    =====================================
   3    These programs are Copyright (c) University of British Columbia. All rights reserved.
   4    They may be freely redistributed in their entirety provided that this copyright
   5    notice is not removed. THEY MAY NOT BE SOLD FOR PROFIT OR INCORPORATED IN
   6    COMMERCIAL PROGRAMS WITHOUT THE WRITTEN PERMISSION OF THE COPYRIGHT HOLDER.
   7    Each program is provided as is, without any express or implied warranty,
   8    without even the warranty of fitness for a particular purpose.
   9
  10    =========================================================
  11    THIS SOFTWARE IS BASED ON HP's implementation of jpeg-ls:
  12    =========================================================
  13
  14    LOCO-I/JPEG-LS IMPLEMENTATION V.0.90
  15    -------------------------------------------------------------------------------
  16    (c) COPYRIGHT HEWLETT-PACKARD COMPANY, 1995-1999.
  17         HEWLETT-PACKARD COMPANY ("HP") DOES NOT WARRANT THE ACCURACY OR
  18    COMPLETENESS OF THE INFORMATION GIVEN HERE.  ANY USE MADE OF, OR
  19    RELIANCE ON, SUCH INFORMATION IS ENTIRELY AT USER'S OWN RISK.
  20         BY DOWNLOADING THE LOCO-I/JPEG-LS COMPRESSORS/DECOMPRESSORS
  21    ("THE SOFTWARE") YOU AGREE TO BE BOUND BY THE TERMS AND CONDITIONS
  22    OF THIS LICENSING AGREEMENT.
  23         YOU MAY DOWNLOAD AND USE THE SOFTWARE FOR NON-COMMERCIAL PURPOSES
  24    FREE OF CHARGE OR FURTHER OBLIGATION.  YOU MAY NOT, DIRECTLY OR
  25    INDIRECTLY, DISTRIBUTE THE SOFTWARE FOR A FEE, INCORPORATE THIS
  26    SOFTWARE INTO ANY PRODUCT OFFERED FOR SALE, OR USE THE SOFTWARE
  27    TO PROVIDE A SERVICE FOR WHICH A FEE IS CHARGED.
  28         YOU MAY MAKE COPIES OF THE SOFTWARE AND DISTRIBUTE SUCH COPIES TO
  29    OTHER PERSONS PROVIDED THAT SUCH COPIES ARE ACCOMPANIED BY
  30    HEWLETT-PACKARD'S COPYRIGHT NOTICE AND THIS AGREEMENT AND THAT
  31    SUCH OTHER PERSONS AGREE TO BE BOUND BY THE TERMS OF THIS AGREEMENT.
  32         THE SOFTWARE IS NOT OF PRODUCT QUALITY AND MAY HAVE ERRORS OR DEFECTS.
  33    THE JPEG-LS STANDARD IS STILL UNDER DEVELOPMENT. THE SOFTWARE IS NOT A
  34    FINAL OR FULL IMPLEMENTATION OF THE STANDARD.  HP GIVES NO EXPRESS OR
  35    IMPLIED WARRANTY OF ANY KIND AND ANY IMPLIED WARRANTIES OF
  36    MERCHANTABILITY AND FITNESS FOR PURPOSE ARE DISCLAIMED.
  37         HP SHALL NOT BE LIABLE FOR ANY DIRECT, INDIRECT, SPECIAL, INCIDENTAL,
  38    OR CONSEQUENTIAL DAMAGES ARISING OUT OF ANY USE OF THE SOFTWARE.
  39    -------------------------------------------------------------------------------
  40 */
  41
  42 /* lossy_e.c ---  the main pipeline which processes a scanline by doing
  43  *              prediction, context computation, context quantization,
  44  *              and statistics gathering. (for LOSSY compression)
  45  *
  46  * Initial code by Alex Jakulin,  Aug. 1995
  47  *
  48  * Modified and optimized: Gadiel Seroussi, October 1995
  49  *
  50  * Modified and added Restart marker and input tables by:
  51  * David Cheng-Hsiu Chu, and Ismail R. Ismail march 1999
  52  *
  53  */
  54
  55 #include <stdio.h>
  56 #include <math.h>
  57
  58 #include "global.h"
  59 #include "bitio.h"
  60
  61
  62 /*byte getk[65][3000];
  63 int clipPx[510];*/
  64
  65 static int eor_limit;
  66
  67
  68 /* Do Golomb statistics and ENCODING for LOSSY images */
  69 inline void lossy_regular_mode(int Q, int SIGN, int Px, pixel *xp)
  70 {
  71   int At, Bt, Nt, absErrval, Errval, MErrval,
  72       qErrval, iqErrval, Rx,
  73       Ix = *xp;                          /* current pixel */
  74   int unary;
  75   int temp;
  76   byte k;
  77
  78
  79   Nt = N[Q];
  80   At = A[Q];
  81   /* Estimate k - Golomb coding variable computation (A.5.1) */
  82   {
  83       register nst = Nt;
  84     for(k=0; nst < At; nst<<=1, k++);
  85   }
  86
  87
  88   /* Prediction correction (A.4.2), compute prediction error (A.4.3)
  89      , and error quantization (A.4.4) */
  90
  91   Px = Px + (SIGN) * C[Q];
  92   clip(Px,alpha);
  93   Errval = SIGN * (Ix - Px);
  94   qErrval = qdiv[Errval];
  95   iqErrval = qmul[qErrval];
  96   Rx = Px + SIGN * iqErrval;
  97
  98   clip(Rx, alpha);
  99   *xp = Rx;          /* store reconstructed pixel in scan line */
 100
 101
 102   /* Modulo reduction of predication error (A.4.5) */
 103   if ( qErrval < 0 )
 104     qErrval += qbeta;      /* qErrval is now in [0..qbeta-1] */
 105
 106
 107   /* Do Rice mapping and compute magnitude of diff */
 108   Bt = B[Q];
 109
 110
 111   /* Error Mapping (A.5.2) */
 112   temp = ( k==0 && NEAR==0 && ((Bt<<1) <= -Nt) );
 113   if (qErrval >= ceil_half_qbeta) {
 114     qErrval -= qbeta;
 115     absErrval = -qErrval;
 116     MErrval = 2*absErrval - 1 - temp;
 117   } else {
 118     absErrval = qErrval;
 119     MErrval = 2*qErrval + temp;
 120   }
 121
 122
 123   /* update bias stats (after correction of the difference) (A.6.1) */
 124
 125   Errval = qmul[qErrval];         /* convert back to alphabet space */
 126
 127   B[Q] = (Bt += Errval);
 128
 129   /* update Golomb stats */
 130   A[Q] += absErrval;
 131
 132   /* check for reset */
 133   if (Nt == reset) {
 134   /* reset for Golomb and bias cancelation at the same time */
 135     N[Q] = (Nt >>= 1);
 136     A[Q] >>= 1;
 137     B[Q] = (Bt >>= 1);
 138   }
 139
 140
 141   /* Do bias estimation for NEXT pixel */
 142   /* Bias cancelation tries to put error in (-1,0] (A.6.2)*/
 143   N[Q] = (++Nt);
 144   if  ( Bt <= -Nt ) {
 145
 146       if (C[Q] > MIN_C)
 147       --C[Q];
 148
 149       Bt = (B[Q] += Nt);
 150
 151       if ( Bt <= -Nt )
 152       B[Q] = -Nt+1;
 153
 154   } else if ( Bt > 0 ) {
 155
 156     if (C[Q] < MAX_C)
 157       ++C[Q];
 158
 159       Bt = (B[Q] -= Nt);
 160
 161       if ( Bt > 0 )
 162       B[Q] = 0;
 163
 164   }
 165
 166
 167   /* Actually output the code: Mapped Error Encoding (A.5.3) */
 168   unary = MErrval >> k;
 169   if ( unary < limit ) {
 170       put_zeros(unary);
 171       putbits((1 << k) + (MErrval & ((1 << k) - 1)), k + 1);
 172   }
 173   else {
 174       put_zeros(limit);
 175       putbits((1<<qbpp) + MErrval - 1, qbpp+1);
 176   }
 177
 178 }
 179
 180
 181
 182
 183 /* Do end of run encoding for LOSSY images -  returns reconstructed value of Ix */
 184 pixel lossy_end_of_run(pixel Ra, pixel Rb, pixel Ix, int RItype)
 185 {
 186   int qErrval, iqErrval, xpr,
 187     MErrval,
 188     Q,
 189     absErrval,
 190     oldmap,
 191     k,
 192     Nt,
 193     At,
 194     Errval;
 195   int Rx;      /* reconstructed pixel */
 196   int unary;
 197
 198   Q = EOR_0 + RItype;
 199   Nt = N[Q];
 200   At = A[Q];
 201
 202   if (RItype) {
 203     Errval = Ix - (xpr = Ra);
 204     At += Nt/2;
 205   }
 206   else {
 207     Errval = Ix - (xpr = Rb);
 208     if ( Rb < Ra )
 209       Errval = -Errval;
 210   }
 211
 212   qErrval = qdiv[Errval];
 213   iqErrval = qmul[qErrval];
 214
 215   if ( RItype || (Rb >= Ra) )
 216     Rx = xpr + iqErrval;
 217   else
 218     Rx = xpr - iqErrval;
 219
 220   clip(Rx,alpha);    /* reconstructed pixel */
 221   Ix = Rx;
 222
 223   /* Estimate k */
 224   for(k=0; Nt < At; Nt *=2, k++);
 225
 226   if (qErrval < 0)
 227     qErrval += qbeta;
 228
 229   if( qErrval >= ceil_half_qbeta )
 230     qErrval -= qbeta;
 231
 232   oldmap = ( k==0 && qErrval && 2*B[Q]<Nt );
 233
 234   /*
 235     Note: the Boolean variable 'oldmap' is not
 236     identical to the variable 'map' in the
 237     JPEG-LS draft. We have
 238     oldmap = (qErrval<0) ? (1-map) : map;
 239   */
 240
 241   /* Error mapping for run-interrupted sample (Figure A.22) */
 242   if( qErrval < 0) {
 243     MErrval = -2*qErrval-1-RItype+oldmap;
 244     B[Q]++;
 245   }
 246   else
 247     MErrval = 2*qErrval-RItype-oldmap;
 248
 249   absErrval = (MErrval+1-RItype)/2;
 250
 251   /* update stats */
 252   A[Q] += absErrval;
 253   if (N[Q] == reset) {
 254     N[Q] >>= 1;
 255     A[Q] >>= 1;
 256     B[Q] >>= 1;
 257   }
 258
 259   N[Q]++; /* for next pixel */
 260
 261   /* Do the actual Golomb encoding: */
 262   eor_limit = limit - limit_reduce;
 263   unary = MErrval >> k;
 264   if ( unary < eor_limit ) {
 265     put_zeros(unary);
 266     putbits((1 << k) + (MErrval & ((1 << k) - 1)), k + 1);
 267   }
 268   else {
 269     put_zeros(eor_limit);
 270     putbits((1<<qbpp) + MErrval-1, qbpp+1);
 271   }
 272
 273   return Ix;
 274 }
 275
 276
 277
 278
 279
 280
 281 /* For line and plane interleaved mode in LOSSY mode */
 282
 283 void lossy_doscanline(pixel *psl,          /* previous scanline */
 284           pixel *sl,           /* current scanline */
 285           int no, int color)   /* number of values in it */
 286
 287 /*** watch it! actual pixels in the scan line are numbered 1 to no .
 288      pixels with indices < 1 or > no are dummy "border" pixels  */
 289 {
 290   int i;
 291   pixel Ra, Rb, Rc, Rd,           /* context pixels */
 292         Ix,                  /* current pixel */
 293         Px;             /* predicted current pixel */
 294   int Rx;              /* reconstructed current pixel */
 295
 296   int SIGN;        /* sign of current context */
 297   int cont;        /* context */
 298   int unary;
 299   int RItype;
 300
 301   i = 1;    /* pixel indices in a scan line go from 1 to no */
 302
 303   /**********************************************/
 304   /* Do for all pixels in the row in 8-bit mode */
 305   /**********************************************/
 306
 307   if (bpp16==FALSE) {
 308
 309     Rc = psl[0];
 310     Rb = psl[1];
 311     Ra = sl[0];
 312
 313     /*  For 8-bit Image */
 314
 315     do {
 316       int RUNcnt;
 317
 318       Ix = sl[i];
 319       Rd = psl[i + 1];
 320
 321       /* Context determination */
 322
 323       /* Quantize the gradient */
 324       /* partial context number: if (b-e) is used then its
 325          contribution is added after determination of the run state.
 326          Also, sign flipping, if any, occurs after run
 327          state determination */
 328
 329       cont =  vLUT[0][Rd - Rb + LUTMAX8] +
 330           vLUT[1][Rb - Rc + LUTMAX8] +
 331           vLUT[2][Rc - Ra + LUTMAX8];
 332
 333
 334       if ( cont == 0 ) {      /* Run state? */
 335       /*************** RUN STATE ***************************/
 336
 337         register delta = Ix - Ra;
 338         RUNcnt = 0;
 339
 340         if ( delta <= NEAR && delta >= negNEAR )
 341         {
 342           while ( 1 )
 343           {
 344             ++RUNcnt;
 345
 346             sl[i] = Ra;
 347
 348             if (++i > no) {
 349               /* Run-lenght coding when reach end of line (A.7.1.2) */
 350               process_run(RUNcnt, EOLINE, color);
 351               return;   /* end of line */
 352             }
 353
 354             Ix = sl[i];
 355
 356             delta = Ix-Ra;
 357             if ( delta > NEAR || delta < negNEAR )  /*  Run is broken */
 358             {
 359               Rd = psl[i + 1];
 360               Rb = psl[i];
 361               break;  /* out of while loop */
 362             }
 363             /* Run continues */
 364           }
 365         }
 366
 367         /* we only get here if the run is broken by
 368            a non-matching symbol */
 369
 370         /* Run-lenght coding when end of line not reached (A.7.1.2) */
 371         process_run(RUNcnt,NOEOLINE, color);
 372
 373
 374         /* This is the END_OF_RUN state */
 375         RItype = ((Rb-Ra)<=NEAR && (Rb-Ra)>=negNEAR);
 376         Ix = lossy_end_of_run(Ra, Rb, Ix, RItype);
 377
 378       }
 379       else {
 380
 381       /*************** REGULAR CONTEXT *******************/
 382
 383         predict(Rb, Ra, Rc);
 384
 385         /* do symmetric context merging */
 386         cont = classmap[cont];
 387
 388         if (cont<0) {
 389           SIGN=-1;
 390           cont = -cont;
 391         } else
 392           SIGN=+1;
 393
 394         /* output a rice code */
 395         lossy_regular_mode(cont, SIGN, Px, &Ix);
 396       }
 397
 398       /* context for next pixel: */
 399       sl[i] = Ix;
 400       Ra = Ix;
 401       Rc = Rb;
 402       Rb = Rd;
 403     } while (++i <= no);
 404
 405   } else { /* 16 bit mode instead of 8 */
 406
 407     /***********************************************/
 408     /* Do for all pixels in the row in 16-bit mode */
 409     /***********************************************/
 410
 411     Rc = ENDIAN16(psl[0]);
 412     Rb = ENDIAN16(psl[1]);
 413     Ra = ENDIAN16(sl[0]);
 414
 415     /*  For 16-bit Image */
 416
 417     do {
 418       int RUNcnt;
 419
 420       Ix = ENDIAN16(sl[i]);
 421       Rd = ENDIAN16(psl[i + 1]);
 422
 423       /* Context determination */
 424
 425       /* Quantize the gradient */
 426       /* partial context number: if (b-e) is used then its
 427         contribution is added after determination of the run state.
 428         Also, sign flipping, if any, occurs after run
 429         state determination */
 430
 431       {
 432         register int diff;
 433
 434         /* Following segment assumes that Sc <= LUTMAX16 */
 435         /* This condition should have been checked when the
 436           lookup tables were built */
 437         diff = Rd - Rb;
 438         if (diff < 0)
 439           cont = (diff > -LUTMAX16) ? vLUT[0][diff + LUTMAX16] : 7*CREGIONS*CREGIONS;
 440         else
 441           cont = (diff < LUTMAX16) ? vLUT[0][diff + LUTMAX16] : 8*CREGIONS*CREGIONS;
 442
 443         diff = Rb - Rc;
 444         if (diff < 0)
 445           cont += (diff > -LUTMAX16) ? vLUT[1][diff + LUTMAX16] : 7*CREGIONS;
 446         else
 447           cont += (diff < LUTMAX16) ? vLUT[1][diff + LUTMAX16] : 8*CREGIONS;
 448
 449         diff = Rc - Ra;
 450         if (diff < 0)
 451           cont += (diff > -LUTMAX16) ? vLUT[2][diff + LUTMAX16] : 7;
 452         else
 453           cont += (diff < LUTMAX16) ? vLUT[2][diff + LUTMAX16] : 8;
 454       }
 455
 456       if ( cont == 0 ) {      /* Run state? */
 457     /*************** RUN STATE ***************************/
 458
 459         register delta = Ix - Ra;
 460         RUNcnt = 0;
 461
 462         if ( delta <= NEAR && delta >= negNEAR )
 463         {
 464           while ( 1 )
 465           {
 466             ++RUNcnt;
 467
 468             sl[i] = ENDIAN16(Ra);
 469
 470             if (++i > no) {
 471               /* Run-lenght coding when reach end of line (A.7.1.2) */
 472               process_run(RUNcnt, EOLINE, color);
 473               return;   /* end of line */
 474             }
 475
 476             Ix = ENDIAN16(sl[i]);
 477
 478             delta = Ix-Ra;
 479             if ( delta > NEAR || delta < negNEAR )  /*  Run is broken */
 480             {
 481               Rd = ENDIAN16(psl[i + 1]);
 482               Rb = ENDIAN16(psl[i]);
 483               break;  /* out of while loop */
 484             }
 485             /* Run continues */
 486           }
 487         }
 488
 489         /* we only get here if the run is broken by
 490           a non-matching symbol */
 491
 492         /* Run-lenght coding when end of line not reached (A.7.1.2) */
 493         process_run(RUNcnt,NOEOLINE, color);
 494
 495         /* This is the END_OF_RUN state */
 496         RItype = ((Rb-Ra)<=NEAR && (Rb-Ra)>=negNEAR);
 497         Ix = lossy_end_of_run(Ra, Rb, Ix, RItype);
 498
 499       }
 500       else {
 501
 502       /*************** REGULAR CONTEXT *******************/
 503
 504         predict(Rb, Ra, Rc);
 505
 506         /* do symmetric context merging */
 507         cont = classmap[cont];
 508
 509         if (cont<0) {
 510           SIGN=-1;
 511           cont = -cont;
 512         } else
 513           SIGN=+1;
 514
 515         /* output a rice code */
 516         lossy_regular_mode(cont, SIGN, Px, &Ix);
 517       }
 518
 519       /* context for next pixel: */
 520       sl[i] = ENDIAN16(Ix);
 521       Ra = Ix;
 522       Rc = Rb;
 523       Rb = Rd;
 524     } while (++i <= no);
 525
 526   } /* for "if" 16 or 8 bit mode */
 527
 528 }
 529
 530
 531
 532
 533
 534
 535 /* For pixel interleavde mode for LOSSY encoding */
 536
 537 void lossy_doscanline_pixel(pixel *psl,            /* previous scanline */
 538               pixel *sl,             /* current scanline */
 539               int no)                /* number of values in it */
 540 /*** watch it! actual pixels in the scan line are numbered 1 to no .
 541      pixels with indices < 1 or > no are dummy "border" pixels  */
 542 {
 543   int i,n_c, enter_run=0, break_run, was_in_run=0, test_run;
 544   int color;    /* Index to the component, 0..COMPONENTS-1 */
 545   pixel c_aa[MAX_COMPONENTS],
 546         c_bb[MAX_COMPONENTS],
 547         c_cc[MAX_COMPONENTS],
 548         c_dd[MAX_COMPONENTS],
 549         c_xx[MAX_COMPONENTS],
 550         Ra, Rb, Rc, Rd,  /* context pixels */
 551         Ix,        /* current pixel */
 552         Px;        /* predicted current pixel */
 553   int SIGN;        /* sign of current context */
 554   int cont,c_cont[MAX_COMPONENTS];        /* context */
 555
 556
 557   if (bpp16==FALSE)
 558   {
 559   /**********************************************/
 560   /* Do for all pixels in the row in 8-bit mode */
 561   /**********************************************/
 562
 563     for (n_c=0; n_c<components; n_c++) {
 564       c_cc[n_c] = psl[n_c];
 565       c_bb[n_c] = psl[components+n_c];
 566       c_aa[n_c] = sl[n_c];
 567     }
 568
 569     i = components;    /* pixel indices in a scan line go from COMPONENTS to no */
 570     color = -1;
 571
 572     do {
 573       int RUNcnt;
 574
 575       if (!was_in_run) color = (color+1)%components;
 576       else color = 0;
 577       Ix = sl[i];
 578
 579       for (n_c=0;n_c<components;n_c++)
 580         c_xx[n_c] = sl[i+n_c];
 581
 582       if (color == 0)
 583       for (n_c=0;n_c<components;n_c++) {
 584
 585         c_dd[n_c] = psl[i+components+n_c];
 586
 587         /* Context determination */
 588
 589         /* Quantize the gradient */
 590         /* partial context number: if (b-e) is used
 591           then its contribution is added after
 592           determination of the run state.
 593           Also, sign flipping, if any, occurs after run
 594           state determination */
 595
 596           c_cont[n_c] =  vLUT[0][c_dd[n_c] - c_bb[n_c] + LUTMAX8] +
 597                 vLUT[1][c_bb[n_c] - c_cc[n_c] + LUTMAX8] +
 598                 vLUT[2][c_cc[n_c] - c_aa[n_c] + LUTMAX8];
 599       }
 600
 601       Ra=c_aa[color];
 602       Rb=c_bb[color];
 603       Rc=c_cc[color];
 604       Rd=c_dd[color];
 605       cont=c_cont[color];
 606
 607       enter_run = was_in_run = test_run = 0;
 608
 609       if (color == 0) {
 610         test_run = 1;
 611         for (n_c=0;n_c<components;n_c++)
 612           if (c_cont[n_c]!=0) {
 613           test_run=0;
 614           break;
 615         }
 616       }
 617
 618       if ( test_run ) {
 619       /*************** RUN STATE ***************************/
 620
 621         int delta[MAX_COMPONENTS];
 622
 623         enter_run = was_in_run = 1;
 624         for (n_c=0;n_c<components;n_c++) {
 625           delta[n_c] = sl[i+n_c] - c_aa[n_c];
 626           if (delta[n_c]>NEAR || delta[n_c]<negNEAR) enter_run=0;
 627         }
 628         RUNcnt = 0;
 629
 630         if (enter_run)
 631         {
 632           while ( 1 ) {
 633             ++RUNcnt;
 634
 635             for (n_c=0; n_c<components; n_c++)
 636               sl[i+n_c] = c_aa[n_c];
 637
 638             if((i=i+components)>(no+components-1)){
 639               process_run(RUNcnt, EOLINE, 0);
 640               return;   /* end of line */
 641             }
 642
 643             for (n_c=0;n_c<components;n_c++)
 644               c_xx[n_c] = sl[i+n_c];
 645
 646             break_run=0;
 647             for (n_c=0;n_c<components;n_c++) {
 648               delta[n_c] = c_xx[n_c] - c_aa[n_c];
 649               if(delta[n_c]>NEAR || delta[n_c]<negNEAR) break_run=1;
 650             }
 651
 652             if ( break_run )
 653             {
 654               for(n_c=0; n_c<components; n_c++){
 655                 c_dd[n_c] = psl[i+components+n_c];
 656                 c_bb[n_c] = psl[i+n_c];
 657               }
 658               break;  /* out of while loop */
 659             }
 660             /* Run continues */
 661           }
 662         }
 663
 664         /* we only get here if the run is broken by
 665           a non-matching symbol */
 666
 667         process_run(RUNcnt, NOEOLINE, 0);
 668
 669         /* This is the END_OF_RUN state */
 670
 671         for (n_c=0;n_c<components;n_c++) {
 672           /* The end of run is done for each component */
 673           Ix = c_xx[n_c];
 674           Rb = c_bb[n_c];
 675           Ra = c_aa[n_c];
 676
 677           /* Handle two special EOR states */
 678           c_xx[n_c] = Ix = lossy_end_of_run(Ra, Rb, Ix, 0);
 679
 680         }   /* loop for components */
 681
 682       }           /* Run state block */
 683       else {
 684
 685       /*************** REGULAR CONTEXT *******************/
 686
 687         predict(Rb, Ra, Rc);
 688         cont = classmap[cont];
 689
 690         if (cont<0)
 691         {  SIGN=-1;
 692           cont=-cont;
 693         }
 694         else
 695           SIGN=1;
 696
 697         /* output a rice code */
 698         lossy_regular_mode(cont, SIGN, Px, &Ix);
 699       }
 700
 701       /* context for next pixel: */
 702       if (!was_in_run) {
 703
 704         c_aa[color] = Ix;
 705         sl[i] = Ix;   /* store reconstructed x */
 706
 707         c_cc[color] = Rb;
 708         c_bb[color] = Rd;
 709         i++;
 710       }
 711       else {
 712         for(n_c=0;n_c<components;n_c++) {
 713           c_aa[n_c] = c_xx[n_c];
 714           sl[i+n_c] = c_xx[n_c];
 715
 716           c_cc[n_c] = c_bb[n_c];
 717           c_bb[n_c] = c_dd[n_c];
 718         }
 719             i += components;
 720       }
 721
 722     } while (i <= (no+components-1));
 723
 724   } else
 725
 726   {
 727
 728   /***********************************************/
 729   /* Do for all pixels in the row in 16-bit mode */
 730   /***********************************************/
 731
 732     for (n_c=0; n_c<components; n_c++) {
 733       c_cc[n_c] = ENDIAN16(psl[n_c]);
 734       c_bb[n_c] = ENDIAN16(psl[components+n_c]);
 735       c_aa[n_c] = ENDIAN16(sl[n_c]);
 736     }
 737
 738     i = components;    /* pixel indices in a scan line go from COMPONENTS to no */
 739     color = -1;
 740
 741     do {
 742       int RUNcnt;
 743
 744       if (!was_in_run) color = (color+1)%components;
 745       else color = 0;
 746       Ix = ENDIAN16(sl[i]);
 747
 748       for (n_c=0;n_c<components;n_c++)
 749         c_xx[n_c] = ENDIAN16(sl[i+n_c]);
 750
 751       if (color == 0)
 752       for (n_c=0;n_c<components;n_c++) {
 753
 754         c_dd[n_c] = ENDIAN16(psl[i+components+n_c]);
 755
 756         /* Context determination */
 757
 758         /* Quantize the gradient */
 759         /* partial context number: if (b-e) is used
 760           then its contribution is added after
 761           determination of the run state.
 762           Also, sign flipping, if any, occurs after run
 763           state determination */
 764         {
 765         register int diff;
 766
 767         /* Following segment assumes that Sc <= LUTMAX16 */
 768         /* This condition should have been checked when the
 769          l  ookup tables were built */
 770         diff = c_dd[n_c] - c_bb[n_c];
 771         if (diff < 0)
 772           c_cont[n_c] = (diff > -LUTMAX16) ? vLUT[0][diff + LUTMAX16] : 7*CREGIONS*CREGIONS;
 773         else
 774           c_cont[n_c] = (diff < LUTMAX16) ? vLUT[0][diff + LUTMAX16] : 8*CREGIONS*CREGIONS;
 775         diff = c_bb[n_c] - c_cc[n_c];
 776         if (diff < 0)
 777           c_cont[n_c] += (diff > -LUTMAX16) ? vLUT[1][diff + LUTMAX16] : 7*CREGIONS;
 778         else
 779           c_cont[n_c] += (diff < LUTMAX16) ? vLUT[1][diff + LUTMAX16] : 8*CREGIONS;
 780         diff = c_cc[n_c] - c_aa[n_c];
 781         if (diff < 0)
 782           c_cont[n_c] += (diff > -LUTMAX16) ? vLUT[2][diff + LUTMAX16] : 7;
 783         else
 784           c_cont[n_c] += (diff < LUTMAX16) ? vLUT[2][diff + LUTMAX16] : 8;
 785         }
 786
 787       }
 788
 789       Ra=c_aa[color];
 790       Rb=c_bb[color];
 791       Rc=c_cc[color];
 792       Rd=c_dd[color];
 793       cont=c_cont[color];
 794
 795       enter_run = was_in_run = test_run = 0;
 796
 797       if (color == 0) {
 798         test_run = 1;
 799         for (n_c=0;n_c<components;n_c++)
 800           if (c_cont[n_c]!=0) {
 801           test_run=0;
 802           break;
 803         }
 804       }
 805
 806       if ( test_run ) {
 807       /*************** RUN STATE ***************************/
 808
 809         int delta[MAX_COMPONENTS];
 810
 811         enter_run = was_in_run = 1;
 812         for (n_c=0;n_c<components;n_c++) {
 813           delta[n_c] = ENDIAN16(sl[i+n_c]) - c_aa[n_c];
 814           if (delta[n_c]>NEAR || delta[n_c]<negNEAR) enter_run=0;
 815         }
 816         RUNcnt = 0;
 817
 818         if (enter_run)
 819         {
 820           while ( 1 ) {
 821             ++RUNcnt;
 822
 823             for (n_c=0; n_c<components; n_c++)
 824               sl[i+n_c] = ENDIAN16(c_aa[n_c]);
 825
 826             if((i=i+components)>(no+components-1)){
 827               process_run(RUNcnt, EOLINE, 0);
 828               return;   /* end of line */
 829             }
 830
 831             for (n_c=0;n_c<components;n_c++)
 832               c_xx[n_c] = ENDIAN16(sl[i+n_c]);
 833
 834             break_run=0;
 835             for (n_c=0;n_c<components;n_c++) {
 836               delta[n_c] = c_xx[n_c] - c_aa[n_c];
 837               if(delta[n_c]>NEAR || delta[n_c]<negNEAR) break_run=1;
 838             }
 839
 840             if ( break_run )
 841             {
 842               for(n_c=0; n_c<components; n_c++){
 843                 c_dd[n_c] = ENDIAN16(psl[i+components+n_c]);
 844                 c_bb[n_c] = ENDIAN16(psl[i+n_c]);
 845               }
 846               break;  /* out of while loop */
 847             }
 848             /* Run continues */
 849           }
 850         }
 851
 852         /* we only get here if the run is broken by
 853           a non-matching symbol */
 854
 855         process_run(RUNcnt, NOEOLINE, 0);
 856
 857         /* This is the END_OF_RUN state */
 858
 859         for (n_c=0;n_c<components;n_c++) {
 860           /* The end of run is done for each component */
 861           Ix = c_xx[n_c];
 862           Rb = c_bb[n_c];
 863           Ra = c_aa[n_c];
 864
 865           /* Handle two special EOR states */
 866           c_xx[n_c] = Ix = lossy_end_of_run(Ra, Rb, Ix, 0);
 867
 868         }   /* loop for components */
 869
 870       }           /* Run state block */
 871       else {
 872
 873       /*************** REGULAR CONTEXT *******************/
 874
 875         predict(Rb, Ra, Rc);
 876         cont = classmap[cont];
 877
 878         if (cont<0)
 879         {  SIGN=-1;
 880           cont=-cont;
 881         }
 882         else
 883           SIGN=1;
 884
 885         /* output a rice code */
 886         lossy_regular_mode(cont, SIGN, Px, &Ix);
 887       }
 888
 889       /* context for next pixel: */
 890       if (!was_in_run) {
 891
 892         c_aa[color] = Ix;
 893         sl[i] = ENDIAN16(Ix);   /* store reconstructed x */
 894
 895         c_cc[color] = Rb;
 896         c_bb[color] = Rd;
 897         i++;
 898       }
 899       else {
 900         for(n_c=0;n_c<components;n_c++) {
 901           c_aa[n_c] = c_xx[n_c];
 902           sl[i+n_c] = ENDIAN16(c_xx[n_c]);
 903
 904           c_cc[n_c] = c_bb[n_c];
 905           c_bb[n_c] = c_dd[n_c];
 906         }
 907             i += components;
 908       }
 909
 910     } while (i <= (no+components-1));
 911
 912   }  /* ends 'if' for 8 or 16 bit */
 913 }
 914
 915