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.
10 =========================================================
11 THIS SOFTWARE IS BASED ON HP's implementation of jpeg-ls:
12 =========================================================
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 -------------------------------------------------------------------------------
42 /* lossless_e.c --- the main pipeline which processes a scanline by doing
43 * prediction, context computation, context quantization,
44 * and statistics gathering.
46 * Initial code by Alex Jakulin, Aug. 1995
48 * Modified and optimized: Gadiel Seroussi, October 1995
50 * Modified and added Restart marker and input tables by:
51 * David Cheng-Hsiu Chu, and Ismail R. Ismail march 1999
61 /*byte getk[65][3000];*/
67 /* Do Golomb statistics and ENCODING for LOSS-LESS images */
68 inline void lossless_regular_mode(int Q, int SIGN, int Px, pixel *xp)
70 int At, Nt, Bt, absErrval, Errval, MErrval,
71 Ix = *xp; /* current pixel */
80 /* Prediction correction (A.4.2), compute prediction error (A.4.3)
81 , and error quantization (A.4.4) */
82 Px = Px + (SIGN) * C[Q];
83 /*Px = clipPx[Px+127];*/
85 Errval = SIGN * (Ix - Px);
88 /* Modulo reduction of predication error (A.4.5) */
90 Errval += alpha; /* Errval is now in [0.. alpha-1] */
93 /* Estimate k - Golomb coding variable computation (A.5.1) */
96 for(k=0; nst < At; nst<<=1, k++);
101 /* Do Rice mapping and compute magnitude of Errval */
104 /* Error Mapping (A.5.2) */
105 temp = ( k==0 && ((Bt<<1) <= -Nt) );
106 if (Errval >= ceil_half_alpha) {
109 MErrval = (absErrval<<1) - 1 - temp;
112 MErrval = (Errval<<1) + temp;
116 /* update bias stats (after correction of the difference) (A.6.1) */
117 B[Q] = (Bt += Errval);
120 /* update Golomb stats */
124 /* check for reset */
126 /* reset for Golomb and bias cancelation at the same time */
134 /* Do bias estimation for NEXT pixel */
135 /* Bias cancelation tries to put error in (-1,0] (A.6.2)*/
141 if ( (B[Q] += Nt) <= -Nt )
144 } else if ( Bt > 0 ) {
149 if ( (B[Q] -= Nt) > 0 )
154 /* Actually output the code: Mapped Error Encoding (Appendix G) */
155 unary = MErrval >> k;
156 if ( unary < limit ) {
158 putbits((1 << k) + (MErrval & ((1 << k) - 1)), k + 1);
162 putbits((1<<qbpp) + MErrval - 1, qbpp+1);
169 /* Do end of run encoding for LOSSLESS images */
170 inline void lossless_end_of_run(pixel Ra, pixel Rb, pixel Ix, int RItype)
197 for(k=0; Nt < At; Nt<<=1, k++);
201 if( Errval >= ceil_half_alpha )
205 oldmap = ( k==0 && Errval && (B[Q]<<1)<Nt );
206 /* Note: the Boolean variable 'oldmap' is not
207 identical to the variable 'map' in the
208 JPEG-LS draft. We have
209 oldmap = (Errval<0) ? (1-map) : map;
212 /* Error mapping for run-interrupted sample (Figure A.22) */
214 MErrval = -(Errval<<1)-1-RItype+oldmap;
217 MErrval = (Errval<<1)-RItype-oldmap;
219 absErrval = (MErrval+1-RItype)>>1;
221 /* Update variables for run-interruped sample (Figure A.23) */
229 N[Q]++; /* for next pixel */
231 /* Do the actual Golomb encoding: */
232 eor_limit = limit - limit_reduce;
233 unary = MErrval >> k;
234 if ( unary < eor_limit ) {
236 putbits((1 << k) + (MErrval & ((1 << k) - 1)), k + 1);
239 put_zeros(eor_limit);
240 putbits((1<<qbpp) + MErrval-1, qbpp+1);
249 /* For line and plane interleaved mode in LOSS-LESS mode */
251 void lossless_doscanline( pixel *psl, /* previous scanline */
252 pixel *sl, /* current scanline */
253 int no, int color) /* number of values in it */
255 /*** watch it! actual pixels in the scan line are numbered 1 to no .
256 pixels with indices < 1 or > no are dummy "border" pixels */
259 pixel Ra, Rb, Rc, Rd, /* context pixels */
260 Ix, /* current pixel */
261 Px; /* predicted current pixel */
263 int SIGN; /* sign of current context */
264 int cont; /* context */
266 i = 1; /* pixel indices in a scan line go from 1 to no */
268 /**********************************************/
269 /* Do for all pixels in the row in 8-bit mode */
270 /**********************************************/
277 /* For 8-bit Image */
285 /* Context determination */
287 /* Quantize the gradient */
288 /* partial context number: if (b-e) is used then its
289 contribution is added after determination of the run state.
290 Also, sign flipping, if any, occurs after run
291 state determination */
294 cont = vLUT[0][Rd - Rb + LUTMAX8] +
295 vLUT[1][Rb - Rc + LUTMAX8] +
296 vLUT[2][Rc - Ra + LUTMAX8];
300 /*************** RUN STATE ***************************/
310 /* Run-lenght coding when reach end of line (A.7.1.2) */
311 process_run(RUNcnt, EOLINE, color);
312 return; /* end of line */
317 if (Ix != Ra) /* Run is broken */
321 break; /* out of while loop */
327 /* we only get here if the run is broken by
328 a non-matching symbol */
330 /* Run-lenght coding when end of line not reached (A.7.1.2) */
331 process_run(RUNcnt,NOEOLINE, color);
334 /* This is the END_OF_RUN state */
335 lossless_end_of_run(Ra, Rb, Ix, (Ra==Rb));
340 /*************** REGULAR CONTEXT *******************/
344 /* do symmetric context merging */
345 cont = classmap[cont];
354 /* output a rice code */
355 lossless_regular_mode(cont, SIGN, Px, &Ix);
358 /* context for next pixel: */
369 /***********************************************/
370 /* Do for all pixels in the row in 16-bit mode */
371 /***********************************************/
373 Rc = ENDIAN16(psl[0]);
374 Rb = ENDIAN16(psl[1]);
375 Ra = ENDIAN16(sl[0]);
377 /* For 16-bit Image */
382 Ix = ENDIAN16(sl[i]);
383 Rd = ENDIAN16(psl[i + 1]);
385 /* Context determination */
387 /* Quantize the gradient */
388 /* partial context number: if (b-e) is used then its
389 contribution is added after determination of the run state.
390 Also, sign flipping, if any, occurs after run
391 state determination */
397 /* Following segment assumes that T3 <= LUTMAX16 */
398 /* This condition should have been checked when the
399 lookup tables were built */
402 cont = (diff > -LUTMAX16) ? vLUT[0][diff + LUTMAX16] : 7*CREGIONS*CREGIONS;
404 cont = (diff < LUTMAX16) ? vLUT[0][diff + LUTMAX16] : 8*CREGIONS*CREGIONS;
408 cont += (diff > -LUTMAX16) ? vLUT[1][diff + LUTMAX16] : 7*CREGIONS;
410 cont += (diff < LUTMAX16) ? vLUT[1][diff + LUTMAX16] : 8*CREGIONS;
414 cont += (diff > -LUTMAX16) ? vLUT[2][diff + LUTMAX16] : 7;
416 cont += (diff < LUTMAX16) ? vLUT[2][diff + LUTMAX16] : 8;
420 if ( cont == 0 ) { /* Run state? */
422 /*************** RUN STATE ***************************/
432 /* Run-lenght coding when reach end of line (A.7.1.2) */
433 process_run(RUNcnt, EOLINE, color);
434 return; /* end of line */
437 Ix = ENDIAN16(sl[i]);
439 if (Ix != Ra) /* Run is broken */
441 Rd = ENDIAN16(psl[i + 1]);
442 Rb = ENDIAN16(psl[i]);
443 break; /* out of while loop */
449 /* we only get here if the run is broken by
450 a non-matching symbol */
452 /* Run-lenght coding when end of line not reached (A.7.1.2) */
453 process_run(RUNcnt,NOEOLINE, color);
455 /* This is the END_OF_RUN state */
456 lossless_end_of_run(Ra, Rb, Ix, (Ra==Rb));
461 /*************** REGULAR CONTEXT *******************/
465 /* do symmetric context merging */
466 cont = classmap[cont];
475 /* output a rice code */
476 lossless_regular_mode(cont, SIGN, Px, &Ix);
479 /* context for next pixel: */
480 sl[i] = ENDIAN16(Ix);
493 /* For pixel interleaved mode for LOSSLESS encoding */
495 void lossless_doscanline_pixel( pixel *psl, /* previous scanline */
496 pixel *sl, /* current scanline */
497 int no) /* number of values in it */
499 /*** watch it! actual pixels in the scan line are numbered 1 to no .
500 pixels with indices < 1 or > no are dummy "border" pixels */
502 int i,n_c, enter_run=0, break_run, was_in_run=0, test_run;
503 int color; /* Index to the component, 0..COMPONENTS-1 */
504 pixel c_aa[MAX_COMPONENTS],
505 c_bb[MAX_COMPONENTS],
506 c_cc[MAX_COMPONENTS],
507 c_dd[MAX_COMPONENTS],
508 c_xx[MAX_COMPONENTS],
509 Ra, Rb, Rc, Rd, /* context pixels */
510 Ix, /* current pixel */
511 Px; /* predicted current pixel */
512 int SIGN; /* sign of current context */
513 int cont,c_cont[MAX_COMPONENTS]; /* context */
518 /**********************************************/
519 /* Do for all pixels in the row in 8-bit mode */
520 /**********************************************/
522 for (n_c=0; n_c<components; n_c++) {
523 c_cc[n_c] = ENDIAN8(psl[n_c]);
524 c_bb[n_c] = ENDIAN8(psl[components+n_c]);
525 c_aa[n_c] = ENDIAN8(sl[n_c]);
528 i = components; /* pixel indices in a scan line go from COMPONENTS to no */
534 if (!was_in_run) color = (color+1)%components;
539 for (n_c=0;n_c<components;n_c++)
540 c_xx[n_c]=ENDIAN8(sl[i+n_c]);
543 for (n_c=0;n_c<components;n_c++)
545 c_dd[n_c] = ENDIAN8(psl[i+components+n_c]);
547 /* Context determination */
549 /* Quantize the gradient */
550 /* partial context number: if (b-e) is used
551 then its contribution is added after
552 determination of the run state.
553 Also, sign flipping, if any, occurs after run
554 state determination */
556 c_cont[n_c] = vLUT[0][c_dd[n_c] - c_bb[n_c] + LUTMAX8] +
557 vLUT[1][c_bb[n_c] - c_cc[n_c] + LUTMAX8] +
558 vLUT[2][c_cc[n_c] - c_aa[n_c] + LUTMAX8];
567 enter_run = was_in_run = test_run = 0;
571 for (n_c=0;n_c<components;n_c++)
572 if (c_cont[n_c]!=0) {
580 /*************** RUN STATE ***************************/
582 enter_run = was_in_run = 1;
583 for (n_c=0;n_c<components;n_c++)
584 if (ENDIAN8(sl[i+n_c]) != c_bb[n_c]) enter_run=0;
591 if((i=i+components)>(no+components-1)){
592 process_run(RUNcnt, EOLINE, 0);
593 return; /* end of line */
596 for (n_c=0;n_c<components;n_c++)
597 c_xx[n_c] = ENDIAN8(sl[i+n_c]);
601 for (n_c=0;n_c<components;n_c++)
602 if (c_xx[n_c] != c_aa[n_c]) break_run=1;
604 if (break_run) /* Run is broken */
606 for(n_c=0;n_c<components;n_c++){
607 c_dd[n_c] = ENDIAN8(psl[i+components+n_c]);
608 c_bb[n_c] = ENDIAN8(psl[i+n_c]);
610 break; /* out of while loop */
616 /* we only get here if the run is broken by
617 a non-matching symbol */
619 process_run(RUNcnt,NOEOLINE, 0);
621 /* This is the END_OF_RUN state */
622 for (n_c=0;n_c<components;n_c++)
624 /* The end of run is done for each component */
629 lossless_end_of_run(Ra, Rb, Ix, 0);
631 } /* loop for components */
633 } /* Run state block */
636 /*************** REGULAR CONTEXT *******************/
639 cont = classmap[cont];
647 /* output a rice code */
648 lossless_regular_mode(cont, SIGN, Px, &Ix);
651 /* context for next pixel: */
659 for(n_c=0;n_c<components;n_c++) {
660 c_aa[n_c] = c_xx[n_c];
661 c_cc[n_c] = c_bb[n_c];
662 c_bb[n_c] = c_dd[n_c];
667 } while (i <= (no+components-1));
671 /**********************************************/
672 /* Do for all pixels in the row in 16-bit mode*/
673 /**********************************************/
675 for (n_c=0; n_c<components; n_c++) {
676 c_cc[n_c] = ENDIAN16(psl[n_c]);
677 c_bb[n_c] = ENDIAN16(psl[components+n_c]);
678 c_aa[n_c] = ENDIAN16(sl[n_c]);
681 i = components; /* pixel indices in a scan line go from COMPONENTS to no */
687 if (!was_in_run) color = (color+1)%components;
690 Ix = ENDIAN16(sl[i]);
692 for (n_c=0;n_c<components;n_c++)
693 c_xx[n_c]=ENDIAN16(sl[i+n_c]);
696 for (n_c=0;n_c<components;n_c++)
698 c_dd[n_c] = ENDIAN16(psl[i+components+n_c]);
700 /* Context determination */
702 /* Quantize the gradient */
703 /* partial context number: if (b-e) is used
704 then its contribution is added after
705 determination of the run state.
706 Also, sign flipping, if any, occurs after run
707 state determination */
712 /* Following segment assumes that Sc <= LUTMAX16 */
713 /* This condition should have been checked when the
714 lookup tables were built */
715 diff = c_dd[n_c] - c_bb[n_c];
717 c_cont[n_c] = (diff > -LUTMAX16) ? vLUT[0][diff + LUTMAX16] : 7*CREGIONS*CREGIONS;
719 c_cont[n_c] = (diff < LUTMAX16) ? vLUT[0][diff + LUTMAX16] : 8*CREGIONS*CREGIONS;
721 diff = c_bb[n_c] - c_cc[n_c];
723 c_cont[n_c] += (diff > -LUTMAX16) ? vLUT[1][diff + LUTMAX16] : 7*CREGIONS;
725 c_cont[n_c] += (diff < LUTMAX16) ? vLUT[1][diff + LUTMAX16] : 8*CREGIONS;
727 diff = c_cc[n_c] - c_aa[n_c];
729 c_cont[n_c] += (diff > -LUTMAX16) ? vLUT[2][diff + LUTMAX16] : 7;
731 c_cont[n_c] += (diff < LUTMAX16) ? vLUT[2][diff + LUTMAX16] : 8;
741 enter_run = was_in_run = test_run = 0;
745 for (n_c=0;n_c<components;n_c++)
746 if (c_cont[n_c]!=0) {
754 /*************** RUN STATE ***************************/
756 enter_run = was_in_run = 1;
757 for (n_c=0;n_c<components;n_c++)
758 if (ENDIAN16(sl[i+n_c]) != c_bb[n_c]) enter_run=0;
765 if((i=i+components)>(no+components-1)){
766 process_run(RUNcnt, EOLINE, 0);
767 return; /* end of line */
770 for (n_c=0;n_c<components;n_c++)
771 c_xx[n_c] = ENDIAN16(sl[i+n_c]);
775 for (n_c=0;n_c<components;n_c++)
776 if (c_xx[n_c] != c_aa[n_c]) break_run=1;
778 if (break_run) /* Run is broken */
780 for(n_c=0;n_c<components;n_c++){
781 c_dd[n_c] = ENDIAN16(psl[i+components+n_c]);
782 c_bb[n_c] = ENDIAN16(psl[i+n_c]);
784 break; /* out of while loop */
790 /* we only get here if the run is broken by
791 a non-matching symbol */
793 process_run(RUNcnt,NOEOLINE, 0);
795 /* This is the END_OF_RUN state */
796 for (n_c=0;n_c<components;n_c++)
798 /* The end of run is done for each component */
803 lossless_end_of_run(Ra, Rb, Ix, 0);
805 } /* loop for components */
807 } /* Run state block */
810 /*************** REGULAR CONTEXT *******************/
813 cont = classmap[cont];
821 /* output a rice code */
822 lossless_regular_mode(cont, SIGN, Px, &Ix);
825 /* context for next pixel: */
833 for(n_c=0;n_c<components;n_c++) {
834 c_aa[n_c] = c_xx[n_c];
835 c_cc[n_c] = c_bb[n_c];
836 c_bb[n_c] = c_dd[n_c];
841 } while (i <= (no+components-1));
843 } /* ends "if" for 8 or 16 bit */