]> Creatis software - gdcm.git/blob - src/gdcmjpegls/Encoder/lossy_e.c
Reverse order sorting now works, even with user supplied function.
[gdcm.git] / 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