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;
341 /**********************************************/
342 /* Do for all pixels in the row in 8-bit mode */
343 /**********************************************/
357 /* Quantize the gradient */
358 cont = vLUT[0][Rd - Rb + LUTMAX8] +
359 vLUT[1][Rb - Rc + LUTMAX8] +
360 vLUT[2][Rc - Ra + LUTMAX8];
364 /*********** RUN STATE **********/
368 /* get length of the run */
369 /* arg is # of pixels left */
370 m = n = process_run_dec(no-i+1, color);
372 if ( m > 0 ) { /* run of nonzero length, otherwise
373 we go directly to the end-of-run
383 /* update context pixels */
388 /* Do end of run encoding for LOSSLESS images */
389 Ra = lossless_end_of_run_d(Ra, Rb, (Ra==Rb));
391 } /* Run state block */
394 /************ REGULAR CONTEXT **********/
398 /* map symmetric contexts */
399 cont = classmap[cont];
409 /* decode a Rice code of a given context */
410 Ra = lossless_regular_mode_d(cont, SIGN, Px);
423 /***********************************************/
424 /* Do for all pixels in the row in 16-bit mode */
425 /***********************************************/
428 Rc = ENDIAN16(psl[0]);
429 Rb = ENDIAN16(psl[1]);
430 Ra = ENDIAN16(sl[0]);
437 Rd = ENDIAN16(psl[i + 1]);
439 /* Quantize the gradient */
443 /* Following segment assumes that T3 <= LUTMAX16 */
444 /* This condition should have been checked when the
445 lookup tables were built */
448 cont = (diff > -LUTMAX16) ? vLUT[0][diff + LUTMAX16] : 7*CREGIONS*CREGIONS;
450 cont = (diff < LUTMAX16) ? vLUT[0][diff + LUTMAX16] : 8*CREGIONS*CREGIONS;
454 cont += (diff > -LUTMAX16) ? vLUT[1][diff + LUTMAX16] : 7*CREGIONS;
456 cont += (diff < LUTMAX16) ? vLUT[1][diff + LUTMAX16] : 8*CREGIONS;
460 cont += (diff > -LUTMAX16) ? vLUT[2][diff + LUTMAX16] : 7;
462 cont += (diff < LUTMAX16) ? vLUT[2][diff + LUTMAX16] : 8;
467 /********** RUN STATE **********/
471 /* get length of the run */
472 /* arg is # of pixels left */
473 m = n = process_run_dec(no-i+1, color);
475 if ( m > 0 ) { /* run of nonzero length, otherwise
476 we go directly to the end-of-run
479 sl[i++] = ENDIAN16(Ra);
486 /* update context pixels */
487 Rb = ENDIAN16(psl[i]);
488 Rd = ENDIAN16(psl[i + 1]);
492 /* Do end of run encoding for LOSSLESS images */
493 Ra = lossless_end_of_run_d(Ra, Rb, (Ra==Rb));
498 /********** REGULAR CONTEXT **********/
502 /* map symmetric contexts */
503 cont = classmap[cont];
513 /* decode a Rice code of a given context */
514 Ra = lossless_regular_mode_d(cont, SIGN, Px);
518 sl[i] = ENDIAN16(Ra);
525 } /* End "if 8/16 bit" */
536 /* For DEOCODING pixel interleavde mode for LOSSLESS images */
538 int lossless_undoscanline_pixel(pixel *psl, /* previous scanline */
539 pixel *sl, /* current scanline */
540 int no) /* number of values in it */
541 /*** watch it! actual pixels in the scan line are numbered 1 to no .
542 pixels with indices < 1 or > no are dummy "border" pixels */
544 int i, n_c, color, was_in_run = 0,
546 pixel Ra, Rb, Rc, Rd;
547 pixel c_aa[MAX_COMPONENTS],
548 c_bb[MAX_COMPONENTS],
549 c_cc[MAX_COMPONENTS],
550 c_dd[MAX_COMPONENTS],
551 c_xx[MAX_COMPONENTS];
555 int cont,c_cont[MAX_COMPONENTS];
558 /**********************************************/
559 /* Do for all pixels in the row in 8-bit mode */
560 /**********************************************/
564 for (n_c=0; n_c<components; n_c++) {
565 c_cc[n_c] = psl[n_c];
566 c_bb[n_c] = psl[components+n_c];
576 if (!was_in_run) color = (color+1)%components;
581 for (n_c=0;n_c<components;n_c++) {
582 c_dd[n_c] = psl[i + components + n_c];
584 /* Quantize the gradient */
585 c_cont[n_c] = vLUT[0][c_dd[n_c] - c_bb[n_c] + LUTMAX8] +
586 vLUT[1][c_bb[n_c] - c_cc[n_c] + LUTMAX8] +
587 vLUT[2][c_cc[n_c] - c_aa[n_c] + LUTMAX8];
596 was_in_run = test_run = 0;
600 for (n_c=0;n_c<components;n_c++)
601 if (c_cont[n_c]!=0) {
609 /********** RUN STATE *********/
615 /* get length of the run */
616 /* arg is # of pixels left */
617 m = n = process_run_dec((no+components-1-i+1)/components, 0);
619 if ( m > 0 ) { /* run of nonzero length, otherwise
620 we go directly to the end-of-run
623 for (n_c=0;n_c<components;n_c++) {
628 if (i > no+components-1)
632 /* update context pixels */
633 for (n_c=0;n_c<components;n_c++) {
634 c_bb[n_c] = psl[i+n_c];
635 c_dd[n_c] = psl[i+components+n_c];
639 /* here we handle the "end-of-run" stat */
641 for (n_c=0;n_c<components;n_c++) {
642 /* The end of run is processed for each component */
645 c_aa[n_c] = c_xx[n_c] = lossless_end_of_run_d(Ra, Rb, 0);
646 } /* Components loop */
648 } /* Run state block */
651 /******* REGULAR CONTEXT *******/
655 cont = classmap[cont];
664 /* decode a Rice code of a given context */
665 c_aa[color] = Ra = lossless_regular_mode_d(cont, SIGN, Px);
675 for (n_c=0;n_c<components;n_c++) {
676 sl[i+n_c] = c_aa[n_c];
677 c_cc[n_c] = c_bb[n_c];
678 c_bb[n_c] = c_dd[n_c];
683 } while (i <= (no+components-1));
687 /***********************************************/
688 /* Do for all pixels in the row in 16-bit mode */
689 /***********************************************/
692 for (n_c=0; n_c<components; n_c++) {
693 c_cc[n_c] = ENDIAN16(psl[n_c]);
694 c_bb[n_c] = ENDIAN16(psl[components+n_c]);
695 c_aa[n_c] = ENDIAN16(sl[n_c]);
704 if (!was_in_run) color = (color+1)%components;
708 for (n_c=0;n_c<components;n_c++) {
710 c_dd[n_c] = ENDIAN16(psl[i + components + n_c]);
712 /* Quantize the gradient */
716 /* Following segment assumes that T3 <= LUTMAX16 */
717 /* This condition should have been checked when the
718 lookup tables were built */
719 diff = c_dd[n_c] - c_bb[n_c];
721 c_cont[n_c] = (diff > -LUTMAX16) ? vLUT[0][diff + LUTMAX16] : 7*CREGIONS*CREGIONS;
723 c_cont[n_c] = (diff < LUTMAX16) ? vLUT[0][diff + LUTMAX16] : 8*CREGIONS*CREGIONS;
725 diff = c_bb[n_c] - c_cc[n_c];
727 c_cont[n_c] += (diff > -LUTMAX16) ? vLUT[1][diff + LUTMAX16] : 7*CREGIONS;
729 c_cont[n_c] += (diff < LUTMAX16) ? vLUT[1][diff + LUTMAX16] : 8*CREGIONS;
731 diff = c_cc[n_c] - c_aa[n_c];
733 c_cont[n_c] += (diff > -LUTMAX16) ? vLUT[2][diff + LUTMAX16] : 7;
735 c_cont[n_c] += (diff < LUTMAX16) ? vLUT[2][diff + LUTMAX16] : 8;
745 was_in_run = test_run = 0;
749 for (n_c=0;n_c<components;n_c++)
750 if (c_cont[n_c]!=0) {
758 /********* RUN STATE *********/
764 /* get length of the run */
765 /* arg is # of pixels left */
766 m = n = process_run_dec((no+components-1-i+1)/components, 0);
768 if ( m > 0 ) { /* run of nonzero length, otherwise
769 we go directly to the end-of-run
772 for (n_c=0;n_c<components;n_c++) {
773 sl[i++] = ENDIAN16(c_aa[n_c]);
777 if (i > no+components-1)
781 /* update context pixels */
782 for (n_c=0;n_c<components;n_c++) {
783 c_bb[n_c] = ENDIAN16(psl[i+n_c]);
784 c_dd[n_c] = ENDIAN16(psl[i+components+n_c]);
788 /* here we handle the "end-of-run" state */
789 for (n_c=0;n_c<components;n_c++) {
790 /* The end of run is processed for each component */
794 c_aa[n_c] = c_xx[n_c] = lossless_end_of_run_d(Ra, Rb, 0);
795 } /* Components loop */
797 } /* Run state block */
800 /******** REGULAR CONTEXT ********/
804 cont = classmap[cont];
814 /* decode a Rice code of a given context */
815 c_aa[color] = Ra = lossless_regular_mode_d(cont, SIGN, Px);
820 sl[i] = ENDIAN16(Ra);
826 for (n_c=0;n_c<components;n_c++) {
827 sl[i+n_c] = ENDIAN16(c_aa[n_c]);
828 c_cc[n_c] = c_bb[n_c];
829 c_bb[n_c] = c_dd[n_c];
834 } while (i <= (no+components-1));
836 } /* ends "if 8/16 bit */