477-505.html

The primary purpose of this book is to explain various data-compression techniques using the C progr

 * - the domain density factor, ranging from 0 (fastest compression)
     to 2
 *   (best but very slow compression).
 * - horizontal and vertical images sizes (default 320 x 200). Both
     sizes
 *   must be multiple of 4 in this implementation (this restriction
     could be
 *   removed with slightly more complex code).
 */
void CompressFile(input, output, argc, argv)
    FILE *input;
    BIT_FILE *output;
    int argc;
    char *argv[];
{
    int x_size = 320;     /* horizontal image size */
    int y_size = 200;     /* vertical image size */
    double quality = 2.0; /* quality factor */
    int s;                /* size index for domains; their size is
                             1<<(s+1) */

    /* Check the command line parameters: */
    for ( ; argc != 0; argv++, argc--) {
        if (argv[0][0] != `-' || argc == 1) {
            fatal_error("Incorrect argument: %s\n", *argv);
        }
        switch(argv[0][1]) {
            case `q': quality     = atof(*++argv); argc--; break;
            case `d': dom_density = atoi(*++argv); argc--; break;
            case `h': x_size      = atoi(*++argv); argc--; break;
            case `v': y_size      = atoi(*++argv); argc--; break;
            default:  fatal_error("Incorrect argument: %s\n", *argv);
        }
    }
    if (dom_density < 0 || dom_density > 2) {
        fatal_error("Incorrect domain density.\n");
    }
    if (x_size % 4 != 0 || y_size % 4 != 0) {
        fatal_error ("Image sizes must be multiple of 4\n");
    }

    /* Allocate and initialize the image data and cumulative image
       data: */
    compress_init(x_size, y_size, input);

    /* Initialize the domain size information as in the
       decompressor: */
    dominfo_init(x_size, y_size, dom_density);

    /* Classify all domains: */
    for (s = MIN_BITS; s <= MAX_BITS; s++) {
        classify_domains(x_size, y_size, s);
    }

    /* Output the header of the compressed file. The first byte
       (`F' for
     * fractal') is just for a consistency check in the decompressor.
     */
    frac_file = output;
    OutputBits(frac_file, (uns_long)'F',    8);
    OutputBits(frac_file, (uns_long)x_size, 16);
    OutputBits(frac_file, (uns_long)y_size, 16);
    OutputBits(frac_file, (uns_long)dom_density, 2);

    /* Compress the whole image recursively, stopping when the image
     * quality is good enough:
     */
    max_error2 = quality*quality;
    traverse_image(0, 0, x_size, y_size, compress_range);

    /* Free all dynamically allocated memory: */
    compress_cleanup(y_size);
}
/* =================================================================
 * Allocate and initialize the image data and cumulative image data.
 */
void compress_init(x_size, y_size, image_file)
    int x_size;       /* horizontal image size */
    int y_size;       /* vertical image size */
    FILE *image_file; /* the input image file */
{
    int x, y;             /* horizontal and vertical indices */
    uns_long r_sum;       /* cumulative range and domain data */
    double r_sum2;        /* cumulative squared range data */
    double d_sum2;        /* cumulative squared domain data */

    range =   (image_data**)allocate(y_size,     x_size,
                                     sizeof(image_data));
    domain    = (unsigned**)allocate(y_size/2,   x_size/2,
                                     sizeof(unsigned));
    cum_range = (uns_long**)allocate(y_size/2+1, x_size/2+1,
                                     sizeof(uns_long));
    cum_range2  =  (float**)allocate(y_size/2+1, x_size/2+1,
                                     sizeof(float));
    cum_domain2 =  (float**)allocate(y_size/2+1, x_size/2+1,
                                     sizeof(float));
    /* Read the input image: */
    for (y = 0; y < y_size; y++) {
       if (fread(range[y], sizeof(image_data), x_size, image_file)
       != x_size) {
           fatal_error("error reading the image data\n");
       }
    }
    /* Compute the `domain' image from the `range' image. Each
       pixel in
     * the domain image is the sum of 4 pixels in the range image.
       We
     * don't average (divide by 4) to avoid losing precision.
     */
     for (y=0; y < y_size/2; y++)
     for (x=0; x < x_size/2; x++) {
         domain[y][x] = (unsigned)range[y<<1][x<<1] + range[y<<1]
         [(x<<1)+1] + range[(y<<1)+1][x<<1] + range[(y<<1)+1][(x<<1)+1];
}

    /* Compute the cumulative data, which will avoid repeated
       computations
     * later (see the region_sum() macro below).
     */
    for (x=0; x <= x_size/2; x++) {
        cum_range[0][x] = 0;
        cum_range2[0][x] = cum_domain2[0][x] = 0.0;
    }
    for (y=0; y < y_size/2; y++) {
        d_sum2 = r_sum2 = 0.0;
        r_sum = cum_range[y+1][0] = 0;
        cum_range2[y+1][0] = cum_domain2[y+1][0] = 0.0;

        for (x=0; x < x_size/2; x++) {
            r_sum += domain[y][x];
            cum_range[y+1][x+1] = cum_range[y][x+1] + r_sum;

            d_sum2 += (double)square(domain[y][x]);
            cum_domain2[y+1][x+1] = cum_domain2[y][x+1] + d_sum2;
            r_sum2 += (double) (square(range[y<<1][x<<1])
                    + square(range[y<<1][(x<<1)+1])
                    + square(range[(y<<1)+1][x<<1])
                    + square(range[(y<<1)+1][(x<<1)+1]));
            cum_range2[y+1][x+1] = cum_range2[y][x+1] + r_sum2;
        }
    }
}
/* ================================================================
 * Free all dynamically allocated data structures for compression.
 */
void compress_cleanup(y_size)
    int y_size;              /* vertical image size */
{
    int s;                   /* size index for domains */
    int class;               /* class number */
    domain_data *dom, *next; /* domain pointers */

    free_array((void**)range,       y_size);
    free_array((void**)domain,      y_size/2);
    free_array((void**)cum_range,   y_size/2 + 1);
    free_array((void**)cum_range2,  y_size/2 + 1);
    free_array((void**)cum_domain2, y_size/2 + 1);

    for (s = MIN_BITS; s <= MAX_BITS; s++)
    for (class = 0; class < NCLASSES; class++)
    for (dom = domain_head[class][s]; dom != NULL; dom = next) {
        next = dom->next;
        free(dom);
    }
}

/* ==================================================================
 * Compute the sum of pixel values or squared pixel values in a range
 * or domain from (x,y) to (x+size-1, y+size-1) included.
 * For a domain, the returned value is scaled by 4 or 16.0
   respectively.
 * x, y and size must all be even.
 */
#define region_sum(cum,x,y,size) \
   (cum[((y)+(size))>>1][((x)+(size))>>1] - cum[(y)>>1]
   [((x)+(size))>>1] \
  - cum[((y)+(size))>>1][(x)>>1]          + cum[(y)>>1][(x)>>1])

/* =================================================================
 * Classify all domains of a certain size. This is done only once to
   save
 * computations later. Each domain is inserted in a linked list
    according
 * to its class and size.
 */
void classify_domains(x_size, y_size, s)
    int x_size;  /* horizontal size of the complete image */
    int y_size;  /* vertical size   of the complete image */
    int s;       /* size index of the domains; their size is
                    1<<(s+1) */
{
    domain_data *dom = NULL; /* pointer to new domain */
    int x, y;                /* horizontal and vertical domain
                                position */
    int class;               /* domain class */
    int dom_size = 1<<(s+1); /* domain size */
    int dom_dist = dom_size >> dom_density; /* distance between
                                      domains */

    /* Initialize all domain lists to be empty: */
    for (class = 0; class < NCLASSES; class++) {
        domain_head[class][s] = NULL;
    }

    /* Classify all domains of this size: */
    for (y = 0; y <= y_size - dom_size; y += dom_dist)
    for (x = 0; x <= x_size - dom_size; x += dom_dist) {

       dom = (domain_data *)xalloc(sizeof(domain_data));
       dom->x = x;
       dom->y = y;
       dom->d_sum  = 0.25  *(double)region_sum(cum_range, x, y,
       dom_size);
       dom->d_sum2 = 0.0625*(double)region_sum(cum_domain2, x, y,
       dom_size);

       class = find_class(x, y, dom_size);

       dom->next = domain_head[class][s];
       domain_head[class][s] = dom;
    }

    /* Check that each domain class contains at least one domain.
     * If a class is empty, we do as if it contains the last created
     * domain (which is actually of a different class).
     */
    for (class = 0; class < NCLASSES; class++) {
        if (domain_head[class][s] == NULL) {
            domain_data *dom2 = (domain_data *)
                                 xalloc(sizeof(domain_data));
            *dom2 = *dom;
            dom2->next = NULL;
            domain_head[class][s] = dom2;
        }
    }
}

/* =================================================================
 * Classify a range or domain.  The class is determined by the
 * ordering of the image brightness in the four quadrants of the
   range
 * or domain. For each quadrant we compute the number of brighter
 * quadrants; this is sufficient to uniquely determine the
 * class. class 0 has quadrants in order of decreasing brightness;
 * class 23 has quadrants in order of increasing brightness.
 *
 * IN assertion: x, y and size are all multiple of 4.
 */
int find_class(x, y, size)
    int x, y;  /* horizontal and vertical position of the range or
                  domain */
    int size;  /* size of the range or domain */
{
    int class = 0;               /* the result class */
    int i,j;                     /* quadrant indices */
    uns_long sum[4];             /* sums for each quadrant */
    static delta[3] = {6, 2, 1}; /* table used to compute the class
                                    number */
    int size1 = size >> 1;

    /* Get the cumulative values of each quadrant. By the IN
       assertion,
     * size1, x+size1 and y+size1 are all even.
     */
    sum[0] = region_sum(cum_range, x,       y,       size1);
    sum[1] = region_sum(cum_range, x,       y+size1, size1);
    sum[2] = region_sum(cum_range, x+size1, y+size1, size1);
    sum[3] = region_sum(cum_range, x+size1, y,       size1);

    /* Compute the class from the ordering of these values */
    for (i = 0;   i <= 2; i++)
    for (j = i+1; j <= 3; j++) {
        if (sum[i] < sum[j]) class += delta[i];
    }
    return class;
}

/* =================================================================
 * Compress a range by searching a match with all domains of the same
   class.
 * Split the range if the mean square error with the best domain
   is larger
 * than max_error2.
 * IN assertion: MIN_BITS <= s_log <= MAX_BITS
 */
void compress_range(x, y, s_log)
    int x, y;    /* horizontal and vertical position of the range */
    int s_log;   /* log base 2 of the range size */
{
    int r_size = 1<<s_log; /* size of the range */
    int class;             /* range class */
    domain_data *dom;     /* used to iterate over all domains of this
    class */
    domain_data *best_dom = NULL; /* pointer to the best domain */
    range_data range;      /* range information for this range */
    affine_map map;        /* affine map for current domain  */
    affine_map best_map;   /* best map for this range */
    uns_long dom_number;   /* domain number */

    /* Compute the range class and cumulative sums: */
    class = find_class(x, y, r_size);
    range.r_sum =  (double)region_sum(cum_range,  x, y, r_size);
    range.r_sum2 = (double)region_sum(cum_range2, x, y, r_size);
    range.x = x;
    range.y = y;
    range.s_log = s_log;

    /* Searching all classes can improve image quality but
       significantly
       slows
     * down compression. Compile with -DCOMPLETE_SEARCH if you can
       wait...
     */
#ifdef COMPLETE_SEARCH
    for (class = 0; class < NCLASSES; class++)
#endif
    for (dom = domain_head[class][s_log];  dom != NULL; dom =
    dom->next) {

        /* Find the optimal map from the range to the domain:
         */
        find_map(&range, dom, &map);

        if (best_dom == NULL || map.error2 < best_map.error2) {
                best_map = map;
                best_dom = dom;
        }
    }

    /* Output the best affine map if the mean square error with the
     * best domain is smaller than max_error2, or if it not possible
     * to split the range because it is too small:
     */
    if (s_log == MIN_BITS ||