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_d.c --- the main pipeline which processes a scanline by doing
43 * prediction, context computation, context quantization,
44 * and statistics gathering. (for lossless mode)
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
63 /* Do Golomb-Rice statistics and DECODING for LOSSLESS images */
64 inline int lossless_regular_mode_d(int Q, int SIGN, int Px)
66 int At, Bt, Nt, Errval, absErrval;
69 /* This function is called only for regular contexts.
70 End_of_run context is treated separately */
76 register int nst = Nt;
77 for(k=0; nst < At; nst *=2, k++);
80 /* Get the number of leading zeros */
85 temp = zeroLUT[reg >> 24];
94 if ( absErrval < limit ) {
95 /* now add the binary part of the Rice code */
97 register unsigned long temp;
104 /* the original unary would have been too long:
105 (mapped value)-1 was sent verbatim */
106 GETBITS(absErrval, qbpp);
110 /* Do the Rice mapping */
111 if ( absErrval & 1 ) { /* negative */
112 absErrval = (absErrval + 1) / 2;
121 if ( k==0 && (2*Bt <= -Nt) )
123 /* special case: see encoder side */
124 Errval = -(Errval+1);
125 absErrval = (Errval<0)? (-Errval):Errval;
128 /* center, clip if necessary, and mask final error */
134 /* this is valid if alpha is a power of 2 */
135 current = (Px - Errval)&(alpha-1);
137 current = Px - Errval;
145 /* valid if alpha is a power of 2 */
146 current = (Px + Errval)&(alpha-1);
148 current = Px + Errval;
154 /* reduce mod alpha, for arbitrary alpha */
157 else if (current >= alpha)
162 /* update bias stats */
163 B[Q] = (Bt += Errval);
165 /* update Golomb-Rice stats */
169 /* check reset (joint for Rice-Golomb and bias cancelation) */
177 /* Do bias estimation for NEXT pixel */
189 } else if ( Bt > 0 ) {
208 /* Do end of run DECODING for LOSSLESS images */
209 inline pixel lossless_end_of_run_d(pixel Ra, pixel Rb, int RItype)
230 for(k=0; Nt < At; Nt *=2, k++);
232 /* read and decode the Golomb code */
233 /* Get the number of leading zeros */
238 temp = zeroLUT[reg >> 24];
241 FILLBUFFER(temp + 1);
247 eor_limit = limit - limit_reduce;
249 if ( MErrval < eor_limit ) {
250 /* now add the binary part of the Golomb code */
252 register unsigned long temp;
259 /* the original unary would have been too long:
260 (mapped value)-1 was sent verbatim */
261 GETBITS(MErrval, qbpp);
265 oldmap = ( k==0 && (RItype||MErrval) && (2*B[Q]<Nt));
267 Note: the Boolean variable 'oldmap' is not
268 identical to the variable 'map' in the
269 JPEG-LS draft. We have
270 oldmap = (qdiff<0) ? (1-map) : map;
273 MErrval += ( RItype + oldmap );
275 if ( MErrval & 1 ) { /* negative */
276 Errval = oldmap - (MErrval+1)/2;
277 absErrval = -Errval-RItype;
280 else { /* nonnegative */
282 absErrval = Errval-RItype;
288 Ix = ( Rb - Errval ) & (alpha-1);
292 else /* includes case a==b */
294 Ix = ( Rb + Errval ) & (alpha-1);
301 /* reduce mod alpha, for arbitrary alpha */
304 else if (Ix >= alpha)
316 N[Q]++; /* for next pixel */
326 /* For line and plane interleaved mode in LOSSLESS mode */
328 int lossless_undoscanline( pixel *psl, /* previous scanline */
329 pixel *sl, /* current scanline */
330 int no, int color) /* number of values in it */
331 /*** watch it! actual pixels in the scan line are numbered 1 to no .
332 pixels with indices < 1 or > no are dummy "border" pixels */
335 pixel Ra, Rb, Rc, Rd;
342 /**********************************************/
343 /* Do for all pixels in the row in 8-bit mode */
344 /**********************************************/
358 /* Quantize the gradient */
359 cont = vLUT[0][Rd - Rb + LUTMAX8] +
360 vLUT[1][Rb - Rc + LUTMAX8] +
361 vLUT[2][Rc - Ra + LUTMAX8];
365 /*********** RUN STATE **********/
369 /* get length of the run */
370 /* arg is # of pixels left */
371 m = n = process_run_dec(no-i+1, color);
373 if ( m > 0 ) { /* run of nonzero length, otherwise
374 we go directly to the end-of-run
384 /* update context pixels */
389 /* Do end of run encoding for LOSSLESS images */
390 Ra = lossless_end_of_run_d(Ra, Rb, (Ra==Rb));
392 } /* Run state block */
395 /************ REGULAR CONTEXT **********/
399 /* map symmetric contexts */
400 cont = classmap[cont];
410 /* decode a Rice code of a given context */
411 Ra = lossless_regular_mode_d(cont, SIGN, Px);
424 /***********************************************/
425 /* Do for all pixels in the row in 16-bit mode */
426 /***********************************************/
429 Rc = ENDIAN16(psl[0]);
430 Rb = ENDIAN16(psl[1]);
431 Ra = ENDIAN16(sl[0]);
438 Rd = ENDIAN16(psl[i + 1]);
440 /* Quantize the gradient */
444 /* Following segment assumes that T3 <= LUTMAX16 */
445 /* This condition should have been checked when the
446 lookup tables were built */
449 cont = (diff > -LUTMAX16) ? vLUT[0][diff + LUTMAX16] : 7*CREGIONS*CREGIONS;
451 cont = (diff < LUTMAX16) ? vLUT[0][diff + LUTMAX16] : 8*CREGIONS*CREGIONS;
455 cont += (diff > -LUTMAX16) ? vLUT[1][diff + LUTMAX16] : 7*CREGIONS;
457 cont += (diff < LUTMAX16) ? vLUT[1][diff + LUTMAX16] : 8*CREGIONS;
461 cont += (diff > -LUTMAX16) ? vLUT[2][diff + LUTMAX16] : 7;
463 cont += (diff < LUTMAX16) ? vLUT[2][diff + LUTMAX16] : 8;
468 /********** RUN STATE **********/
472 /* get length of the run */
473 /* arg is # of pixels left */
474 m = n = process_run_dec(no-i+1, color);
476 if ( m > 0 ) { /* run of nonzero length, otherwise
477 we go directly to the end-of-run
480 sl[i++] = ENDIAN16(Ra);
487 /* update context pixels */
488 Rb = ENDIAN16(psl[i]);
489 Rd = ENDIAN16(psl[i + 1]);
493 /* Do end of run encoding for LOSSLESS images */
494 Ra = lossless_end_of_run_d(Ra, Rb, (Ra==Rb));
499 /********** REGULAR CONTEXT **********/
503 /* map symmetric contexts */
504 cont = classmap[cont];
514 /* decode a Rice code of a given context */
515 Ra = lossless_regular_mode_d(cont, SIGN, Px);
519 sl[i] = ENDIAN16(Ra);
526 } /* End "if 8/16 bit" */
537 /* For DEOCODING pixel interleavde mode for LOSSLESS images */
539 int lossless_undoscanline_pixel(pixel *psl, /* previous scanline */
540 pixel *sl, /* current scanline */
541 int no) /* number of values in it */
542 /*** watch it! actual pixels in the scan line are numbered 1 to no .
543 pixels with indices < 1 or > no are dummy "border" pixels */
545 int i, psfix, n_c, color, enter_run=0, was_in_run = 0,
547 pixel Ra, Rb, Rc, Rd;
548 pixel c_aa[MAX_COMPONENTS],
549 c_bb[MAX_COMPONENTS],
550 c_cc[MAX_COMPONENTS],
551 c_dd[MAX_COMPONENTS],
552 c_xx[MAX_COMPONENTS];
556 int cont,c_cont[MAX_COMPONENTS];
560 /**********************************************/
561 /* Do for all pixels in the row in 8-bit mode */
562 /**********************************************/
566 for (n_c=0; n_c<components; n_c++) {
567 c_cc[n_c] = psl[n_c];
568 c_bb[n_c] = psl[components+n_c];
578 if (!was_in_run) color = (color+1)%components;
583 for (n_c=0;n_c<components;n_c++) {
584 c_dd[n_c] = psl[i + components + n_c];
586 /* Quantize the gradient */
587 c_cont[n_c] = vLUT[0][c_dd[n_c] - c_bb[n_c] + LUTMAX8] +
588 vLUT[1][c_bb[n_c] - c_cc[n_c] + LUTMAX8] +
589 vLUT[2][c_cc[n_c] - c_aa[n_c] + LUTMAX8];
598 enter_run = was_in_run = test_run = 0;
602 for (n_c=0;n_c<components;n_c++)
603 if (c_cont[n_c]!=0) {
611 /********** RUN STATE *********/
615 enter_run = was_in_run = 1;
617 /* get length of the run */
618 /* arg is # of pixels left */
619 m = n = process_run_dec((no+components-1-i+1)/components, 0);
621 if ( m > 0 ) { /* run of nonzero length, otherwise
622 we go directly to the end-of-run
625 for (n_c=0;n_c<components;n_c++) {
630 if (i > no+components-1)
634 /* update context pixels */
635 for (n_c=0;n_c<components;n_c++) {
636 c_bb[n_c] = psl[i+n_c];
637 c_dd[n_c] = psl[i+components+n_c];
641 /* here we handle the "end-of-run" stat */
643 for (n_c=0;n_c<components;n_c++) {
644 /* The end of run is processed for each component */
647 c_aa[n_c] = c_xx[n_c] = lossless_end_of_run_d(Ra, Rb, 0);
648 } /* Components loop */
650 } /* Run state block */
653 /******* REGULAR CONTEXT *******/
657 cont = classmap[cont];
666 /* decode a Rice code of a given context */
667 c_aa[color] = Ra = lossless_regular_mode_d(cont, SIGN, Px);
677 for (n_c=0;n_c<components;n_c++) {
678 sl[i+n_c] = c_aa[n_c];
679 c_cc[n_c] = c_bb[n_c];
680 c_bb[n_c] = c_dd[n_c];
685 } while (i <= (no+components-1));
689 /***********************************************/
690 /* Do for all pixels in the row in 16-bit mode */
691 /***********************************************/
694 for (n_c=0; n_c<components; n_c++) {
695 c_cc[n_c] = ENDIAN16(psl[n_c]);
696 c_bb[n_c] = ENDIAN16(psl[components+n_c]);
697 c_aa[n_c] = ENDIAN16(sl[n_c]);
706 if (!was_in_run) color = (color+1)%components;
710 for (n_c=0;n_c<components;n_c++) {
712 c_dd[n_c] = ENDIAN16(psl[i + components + n_c]);
714 /* Quantize the gradient */
718 /* Following segment assumes that T3 <= LUTMAX16 */
719 /* This condition should have been checked when the
720 lookup tables were built */
721 diff = c_dd[n_c] - c_bb[n_c];
723 c_cont[n_c] = (diff > -LUTMAX16) ? vLUT[0][diff + LUTMAX16] : 7*CREGIONS*CREGIONS;
725 c_cont[n_c] = (diff < LUTMAX16) ? vLUT[0][diff + LUTMAX16] : 8*CREGIONS*CREGIONS;
727 diff = c_bb[n_c] - c_cc[n_c];
729 c_cont[n_c] += (diff > -LUTMAX16) ? vLUT[1][diff + LUTMAX16] : 7*CREGIONS;
731 c_cont[n_c] += (diff < LUTMAX16) ? vLUT[1][diff + LUTMAX16] : 8*CREGIONS;
733 diff = c_cc[n_c] - c_aa[n_c];
735 c_cont[n_c] += (diff > -LUTMAX16) ? vLUT[2][diff + LUTMAX16] : 7;
737 c_cont[n_c] += (diff < LUTMAX16) ? vLUT[2][diff + LUTMAX16] : 8;
747 enter_run = was_in_run = test_run = 0;
751 for (n_c=0;n_c<components;n_c++)
752 if (c_cont[n_c]!=0) {
760 /********* RUN STATE *********/
764 enter_run = was_in_run = 1;
766 /* get length of the run */
767 /* arg is # of pixels left */
768 m = n = process_run_dec((no+components-1-i+1)/components, 0);
770 if ( m > 0 ) { /* run of nonzero length, otherwise
771 we go directly to the end-of-run
774 for (n_c=0;n_c<components;n_c++) {
775 sl[i++] = ENDIAN16(c_aa[n_c]);
779 if (i > no+components-1)
783 /* update context pixels */
784 for (n_c=0;n_c<components;n_c++) {
785 c_bb[n_c] = ENDIAN16(psl[i+n_c]);
786 c_dd[n_c] = ENDIAN16(psl[i+components+n_c]);
790 /* here we handle the "end-of-run" state */
791 for (n_c=0;n_c<components;n_c++) {
792 /* The end of run is processed for each component */
796 c_aa[n_c] = c_xx[n_c] = lossless_end_of_run_d(Ra, Rb, 0);
797 } /* Components loop */
799 } /* Run state block */
802 /******** REGULAR CONTEXT ********/
806 cont = classmap[cont];
816 /* decode a Rice code of a given context */
817 c_aa[color] = Ra = lossless_regular_mode_d(cont, SIGN, Px);
822 sl[i] = ENDIAN16(Ra);
828 for (n_c=0;n_c<components;n_c++) {
829 sl[i+n_c] = ENDIAN16(c_aa[n_c]);
830 c_cc[n_c] = c_bb[n_c];
831 c_bb[n_c] = c_dd[n_c];
836 } while (i <= (no+components-1));
838 } /* ends "if 8/16 bit */