From: lusardi@sybil.cs.Buffalo.EDU (Christopher Lusardi)
Newsgroups: comp.compression
Subject: Some Code from MN's Book for Image Compression
Date: 10 Apr 92 02:04:41 GMT
Organization: State University of New York at Buffalo/Comp Sci
Originator: lusardi@sybil.cs.Buffalo.EDU

Here is part of Mark Nelson's  code (from his book "The Data Compression
Book") from the files dct.c, and bitio.c for dct image compression. 
I've left out these 3 files : 
  1. main-c.h for main-c.c, and 2. main-e.c with decompression main routine. 
  3. main-e.h for main-e.c which are similar to main-c.c except it calls
  the routine ExpandFile.
  (3.errhand.h, and errhand.c in the book aren't required for image compression
   so I left them out.)

This program requires a raw image 320 rows by 200 cols (the code is included 
for generating it).

I modified the code from the book with an standard DCT (and inverse DCT) 
formula and a quantized array; so, this code doesn't use the quality value 
from the command line. I modified it because I couldn't find the bug 
--output from decompression can't be put on screen--in MN's version, but
this same bug is still in the program.

Suggests on why the code is wrong would be appreciated.
I have to apologize for the lack of comments.

                                                      Enjoy,
                                                      cml

/*-------------------- uncommented file dct.c -----------------------------*/

#include <stdio.h>
#include <math.h>
#include "bitio.h"

#define ROWS 200
#define COLS 320
#define N    8
/* Rounding an integer */
#define ROUND(a) (((a) < 0) ? (int) ((a) - 0.5) : \
                                       (int) ((a) + 0.5))
/* For user information */
char *CompressionName = "DCT compression";
char *Usage           = "infile outfile [quality]\nQuality from 1-25";

void Initialize ();
void ReadPixelStrip ();
int InputCode ();
void ReadDCTData ();
void OutputCode();
void WriteDCTData ();
void WritePixelStrip();
void ForwardDCT ();
void InverseDCT ();
void CompressFile ();
void ExpandFile ();

unsigned char PixelStrip [N][COLS];
double C[N][N];
double Ct[N][N];
int InputRunLength;
int OutputRunLength;
double Coefficients [N][N];
int Quantum [N][N] = {                      /* Quantization matrix: I added 
                                               this, no quality on com line. */
  {16, 11, 10, 16,  24,  40,  51,  61},
  {12, 12, 14, 19,  26,  58,  60,  55},
  {14, 13, 16, 24,  40,  57,  69,  56},
  {14, 17, 22, 29,  51,  87,  80,  62},
  {18, 22, 37, 56,  68, 109, 103,  77},
  {24, 35, 55, 64,  81, 104, 113,  92},
  {49, 64, 78, 87, 103, 121, 120, 101},
  {72, 92, 95, 98, 112, 100, 103,  99}
};

struct zigzag {                            /* path of quantification steps */
  int row;
  int col;
} ZigZag [N * N] =
{
  {0,0},
  {0,1},{1,0},
  {2,0},{1,1},{0,2},
  {0,3},{1,2},{2,1},{3,0},
  {4,0},{3,1},{2,2},{1,3},{0,4},
  {0,5},{1,4},{2,3},{3,2},{4,1},{5,0},
  {6,0},{5,1},{4,2},{3,3},{2,4},{1,5},{0,6},
  {0,7},{1,6},{2,5},{3,4},{4,3},{5,2},{6,1},{7,0},
  {7,1},{6,2},{5,3},{4,4},{3,5},{2,6},{1,7},
  {2,7},{3,6},{4,5},{5,4},{6,3},{7,2},
  {7,3},{6,4},{5,5},{4,6},{3,7},
  {4,7},{5,6},{6,5},{7,4},
  {7,5},{6,6},{5,7},
  {6,7},{7,6},
  {7,7}
};
  
void Initialize (quality) /* Called by main routines */
int quality;
{
  int i;
  int j;
  int k;
  double temp;
  double pi = atan (1.0) * 4.0;

  OutputRunLength = 0;
  InputRunLength = 0;
  
  for (j = 0; j < N; j++) 
    for (k = 0; k < N; k++) {
      C [j][k] = cos(((2.0 * (double) k + 1.0) * 
                      (double) j * pi) / (2.0 * (double) N));
     if (j == 0)
       temp = 1.0 / sqrt (2.0);
     else
       temp = 1;
     if (k == 0)
       temp *= 1.0 / sqrt (2.0);
    Coefficients [j][k] = (double) temp;
   }
}

void ReadPixelStrip (input, strip)
FILE *input;
unsigned char strip [N][COLS];
{
  int row;
  int col;
  int c;
  for (row = 0; row < N; row++)
    for (col = 0; col < COLS; col++) {
      c = getc (input);
      if (c == EOF) {
        printf ("Error reading input gray scale file\n");
        exit (1);
      }
      strip[row][col] = (unsigned char) c;
    }
}

int InputCode (input_file)
BIT_FILE *input_file;
{
  int bit_count;
  int result;
   
  if (InputRunLength > 0) {
    InputRunLength --;
    return (0);
  }
  bit_count = (int) InputBits (input_file, 2);
  if (bit_count == 0) {
    InputRunLength = (int) InputBits (input_file, 4);
    return (0);
  }
  if (bit_count == 1)
    bit_count = (int) InputBits (input_file, 1) + 1;
  else
    bit_count = (int) InputBits (input_file, 2) + 
                     (bit_count << 2) - 5;
  result = (int) InputBits (input_file, bit_count);
  if (result & (1 << (bit_count - 1)))
    return (result);
  return (result - (1 << bit_count) + 1);
}

void ReadDCTData (input_file, input_data)
BIT_FILE *input_file;
int input_data [N][N];
{
  int i;
  int row;
  int col;
  
  for (i = 0;i < (N * N); i ++) {
    row = ZigZag [i].row;
    col = ZigZag[i].col;
    input_data[row][col] = InputCode (input_file) *
                           Quantum[row][col]; 
   }
}

void OutputCode (output_file, code)
BIT_FILE *output_file;
int code;
{
  int top_of_range;
  int abs_code;
  int bit_count;

  if (code == 0) {
    OutputRunLength ++;
    return;
  }
  if (OutputRunLength != 0) {
    while (OutputRunLength > 0) {
      OutputBits (output_file, 0L, 2);
      if (OutputRunLength <= 16) {
        OutputBits (output_file, 
                    (unsigned long) (OutputRunLength - 1), 4);
        OutputRunLength = 0;
      }
      else {
        OutputBits (output_file, 15L, 4);
        OutputRunLength -= 16;
      }
   }
  }
  if (code < 0)
    abs_code = -code;
        else
           abs_code = code;
        top_of_range = 1;
        bit_count = 1;
        while (abs_code > top_of_range) {
          bit_count++;
          top_of_range = ((top_of_range + 1) * 2) - 1;
        }
        if (bit_count < 3)
           OutputBits (output_file, (unsigned long) (bit_count + 1), 3);
        else
           OutputBits (output_file, (unsigned long) (bit_count + 5), 4);
        if (code > 0)
	  OutputBits (output_file, (unsigned long) code, bit_count);
        else
           OutputBits (output_file, (unsigned long) (code + top_of_range),
                          bit_count);
}

void WriteDCTData (output_file, output_data)
BIT_FILE *output_file;
int output_data [N][N];
{
  int i;
  int row;
  int col;
  double result;
  
  for (i = 0; i < (N * N); i++) {
    row = ZigZag [i].row;
    col = ZigZag [i].col;

    result = (double) output_data [row][col] / (double) Quantum[row][col];
    OutputCode(output_file, result);
  }
}

void WritePixelStrip (output, strip)
FILE *output;
unsigned char strip [N][COLS];
{
  int row;
  int col;

  for (row = 0; row < N; row++)
    for (col = 0; col < COLS; col++)
      putc (strip [row][col], output);  
}

/* For image processing applications, the forward 2-D DCT of an n x n block
of pixels is often defined as:
                                             _          _     _          _
F(u,v) = 4C(u)C(v)   n - 1  n - 1           |            |   |            |
         ---------   Sigma  Sigma f(j,k) cos| (2j+1u)pi  |cos| (2k+1)v pi |
              2      j = 0  k = 0           |------------|   |------------|
            n                               |_    2n    _|   |_    2n    _|
 
where
                   C(w) =  1 / sqrt(2) for w = 0 
                           1           for w = 1,2,...,n-1
              and pi is the mathematical pi

*/
void ForwardDCT (input, output)
unsigned char *input [N];
int output [N][N];
{
  double temp;
  int u;
  int v;
  int j,k;

  for (u = 0; u < N; u++)
    for (v = 0; v < N; v++) {
      temp = 0.0;
      for (j = 0; j < N; j++) 
        for (k = 0; k < N; k++) 
	  temp  += (double) input [j] [k] * C [u] [j] * C[v] [k];
      temp *= 4.0 * Coefficients [u][v] / (double) (N * N);
      output [u][v] = ROUND(temp);
    }
}

/* The inverse 2-D DCT is defined as:

                                         _          _     _          _
          n-1    n-1                    |            |   |            |
f(j,k) = Sigma  Sigma C(u)C(v)F(j,k) cos| (2j+1u)pi  |cos| (2k+1)v pi |
          u=0    v=0                    |------------|   |------------|
                                        |_    2n    _|   |_    2n    _|

where
                   C(w) =  1 / sqrt(2) for w = 0 
                           1           for w = 1,2,...,n-1

*/
void InverseDCT (input, output)
int input[N][N];
unsigned char *output [N];
{ 
  double temp;
  int i;
  int j;
  int k;
  int u,v;

  for (j = 0; j < N; j++) 
    for (k = 0; k < N; k++) {
       temp = 0.0;
       for (u = 0; u < N; u++)
	 for (v = 0; v < N; v++) {
	  temp += input [u] [v] * C [u][j] * C [v] [k];
          temp *= Coefficients [u][v];
         }
       output[j][k] =  (unsigned char) ROUND (temp);
    }
}

void CompressFile (input, output, argc, argv) /* Called by main-c.c */
FILE *input;
BIT_FILE *output;
int argc;
char *argv [];
{
int r,c;
  int row;
  int col;
  int i;
  unsigned char *input_array [N];
  int output_array [N][N];
  int quality;

  if (argc-- > 0)
    quality = atoi(argv [0]);
  else
    quality = 3;
  if (quality < 0 || quality > 50) {
    printf ("Illegal quality factor of %d\n", quality);
    exit (1);
  }
  printf ("Using quality factor of %d\n", quality);
  Initialize (quality);
  OutputBits (output, (unsigned long) quality, 8);
  for (row = 0; row < ROWS; row += N) {
    ReadPixelStrip (input, PixelStrip);
    for (col = 0; col < COLS;col += N) {
      for (i = 0; i < N;i ++) 
	input_array [i] = PixelStrip [i] + col;
      ForwardDCT (input_array, output_array); 
      WriteDCTData (output, output_array);
    }
  }
  OutputCode (output, 1);
   while (argc-- > 0)
     printf ("Unused argument: %s\n", *argv++);
}

void ExpandFile (input, output, argc, argv)  /* Called by main-e.c */
BIT_FILE *input;
FILE *output;
int argc;
char *argv[];
{
  int u,v;
  int row;
  int col;
  int i;
  int input_array [N][N];
  unsigned char *output_array [N];
  int quality;

  quality = (int) InputBits (input, 8);
  printf ("\rUsing quality factor of %d\n", quality);
  Initialize (quality);
  for (row = 0; row < ROWS; row += N) {
    for (col = 0; col < COLS; col += N) {
      for (i = 0; i < N; i ++)
        output_array [i] = PixelStrip [i] + col;
      ReadDCTData (input, input_array);
      InverseDCT (input_array, output_array);
    }
    WritePixelStrip (output, PixelStrip);
  }
  while (argc-- > 0)
    printf ("Unused argument: %s\n", *argv++);
}
/*-------------------- uncommented file bitio.c -----------------------------*/
#include <stdio.h>
#include "bitio.h"

BIT_FILE *OpenOutputBitFile (name)
char *name;
{
  BIT_FILE *bit_file;
  
  bit_file = (BIT_FILE *) calloc (1,sizeof(BIT_FILE));
  if (bit_file == NULL)
       return (bit_file);
  bit_file->file = fopen(name,"wb");
  bit_file->rack = 0;
  bit_file->mask = 0x80;
  bit_file->pacifier_counter = 0;
  return(bit_file);
}

BIT_FILE *OpenInputBitFile(name)    
char *name;
{
  BIT_FILE *bit_file;
 
  bit_file = (BIT_FILE *) calloc(1,sizeof(BIT_FILE));
  if (bit_file == NULL)
    return(bit_file);
  bit_file->file = fopen(name,"rb");
  bit_file->rack = 0;
  bit_file->mask = 0x80;
  bit_file->pacifier_counter = 0;
  return(bit_file);
}

void CloseOutputBitFile(bit_file)
BIT_FILE *bit_file;
{
  if (bit_file->mask != 0x80)
    if (putc(bit_file->rack,bit_file->file) != bit_file->rack) {
       printf("Fatal error in CloseBitFile!\n");
       exit (1);
     }
  fclose(bit_file->file);
  free((char *) bit_file);
}

void CloseInputBitFile (bit_file)
BIT_FILE *bit_file;
{
   fclose (bit_file->file);
   free((char*) bit_file);
}

void OutputBits(bit_file,code,count)  /* Record bits, output when have byte */
BIT_FILE *bit_file;
unsigned long code;
int count;
{
  unsigned long mask;

  mask = 1L << (count - 1);
  while (mask != 0)  {
    if (mask & code) 
       bit_file->rack |= bit_file->mask;
     bit_file->mask >>= 1;
     if (bit_file->mask == 0) {
       if (putc(bit_file->rack, bit_file->file) != bit_file->rack) {
	printf("Fatal error in OutputBit!\n");
        exit (1);
      }
      else if ((bit_file->pacifier_counter++ & 2047) == 0)
	putc ('.', stdout);
      bit_file->rack = 0;
      bit_file->mask = 0x80;
   }
    mask >>= 1;     
  }
}

unsigned long InputBits (bit_file, bit_count)
BIT_FILE *bit_file;
int bit_count;
{
  unsigned long mask;
  unsigned long return_value;
  
  mask = 1L << (bit_count - 1);
  return_value = 0;
  while (mask != 0)  {
    if (bit_file->mask == 0x80) {
      bit_file->rack = getc(bit_file->file);
      if (bit_file->rack == EOF) {
	printf ("Fatal error \n InputBits!\n");
        exit (1);
      }
      if ((bit_file->pacifier_counter++ & 2047) == 0)
        putc('.',stdout);
    }
    if (bit_file->rack & bit_file->mask)
      return_value |= mask;
    mask >>= 1;
    bit_file->mask >>= 1;
    if (bit_file->mask == 0)
        bit_file->mask = 0x80;
  }
  return (return_value);
}

/*-------------------- uncommented file bitio.h -----------------------------*/
/*-------------------- contains data structure for bitio.c-------------------*/
#ifndef _BITIO_H
#define _BITIO_H

#include <stdio.h>

typedef struct bib_file {
    FILE *file;
    unsigned char mask; 
    int rack;
    int pacifier_counter;
} BIT_FILE;

BIT_FILE       *OpenInputBitFile ();
BIT_FILE       *OpenOutputBitFile ();
void           OutputBit ();
void           OutputBits ();
int            InputBit ();
unsigned long  InputBits ();
void           CloseInputBitFile ();
void           CloseOutputBitFile ();

#endif

/*-------------------- uncommented file main-c.c ----------------------------*/

#include <stdio.h>
#include <string.h>
#include "bitio.h"
#include "main.h"

void usage_exit ();
void print_ratios ();
long file_size ();
extern int OutputRunLength;

int main (argc, argv)
int argc;
char *argv[];
{
  char *temp_argv;
  BIT_FILE *output;
  FILE *input;
 
  setbuf(stdout, NULL);
  if (argc < 3) 
    usage_exit (argv [0]);
  input = fopen (argv [1],"rb");
  if (input == NULL)  {
    printf ("Error opening %s for input\n",argv[1]);
    exit (1);
  }
  output = OpenOutputBitFile (argv[2]);
  if (output == NULL)  {
    printf ("Error opening %s for output\n",argv[2]);
    exit (1);
  }
  printf("\nCompresing %s to %s\n",argv[1],argv[2]);
  printf("Using %s\n",CompressionName);
  argc -= 3;                        /* Small error probably here (?)*/
temp_argv = argv;
  argv += 3;
  CompressFile (input, output, argc, argv);
  CloseOutputBitFile (output);
  fclose (input);
argv = temp_argv;
  print_ratios(argv[1],argv[2]);
  exit (0);  
}

void usage_exit (prog_name)
char *prog_name;
{
  char *short_name;
  char *extension;

  short_name = strrchr(prog_name, '\\');
  if (short_name == NULL)
     short_name = strrchr (prog_name, '/');
  if (short_name = NULL)
     short_name = strrchr(prog_name, ':');
  if (short_name != NULL)
     short_name++;
  else
     short_name = prog_name;
  extension = strrchr (short_name, '.');
  if (extension != NULL)
    *extension = '\0';
  printf ("\nUsage: %s %s\n", short_name, Usage);
  exit (0);
}  
  
#ifndef SEEK_END
#define SEEK_END 2
#endif

long file_size (name)
char *name;
{
  long eof_ftell;
  FILE *file;
  
  file = fopen (name, "rb");
  if (file == NULL)
    return (0L);
  fseek (file,0L,SEEK_END);
  eof_ftell = ftell (file);
  fclose (file);
  return (eof_ftell);
}

void print_ratios (input, output)
char *input;
char *output;
{
  long input_size;
  long output_size;
  int ratio;

  input_size = file_size (input);
  if (input_size == 0)
      input_size = 1;
  output_size = file_size (output);  /* code for this size is mine */
  if (output_size == 0)
      output_size = 1;
  ratio = (output_size * 100)/ input_size;
  printf ("\nInput bytes:      %ld\n",input_size);
  printf ("Output bytes:       %ld\n",  output_size);
  printf("Compression ratio: %d%%\n", ratio);
}
/*-------------------- code to generate your raw image ---------------------*/

#include <stdio.h>

main(argc,argv)
int argc;
char *argv[];
{
	int i,j,width,height;

	if (argc < 3) {
		fprintf(stderr,"Usage: %s <width> <height>\n",argv[0]);
	} else {
		width = atol(argv[1]);
		height = atol(argv[2]);
		for (i=0;i<height;i++)
			for (j=0;j<width;j++)
				putchar((char)((i+j)%256));
	}
}

-- 
Christopher M. Lusardi     University At Buffalo     lusardi@acsu.buffalo.edu