477-505.html

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

        best_map.error2 <= max_error2*((long)r_size*r_size)) {

        /* If the range is too small to be split, the decompressor
           knows
         * this, otherwise we must indicate that the range has not
           been
           split:
         */
        if (s_log != MIN_BITS) {
            OutputBit(frac_file, 1);  /* affine map follows */
        }
        OutputBits(frac_file, (uns_long)best_map.contrast,
        CONTRAST_BITS);
        OutputBits(frac_file, (uns_long)best_map.brightness,
        BRIGHTNESS_BITS);

        /* When the contrast is null, the decompressor does not
           need to
           know
         * which domain was selected:
         */
        if (best_map.contrast == 0) return;
        dom_number = (uns_long)best_dom->y * dom_info[s_log]
        .x_domains
                      + (uns_long)best_dom->x;

        /* The distance between two domains is the domain size
           1<<(s_log+1)
         * shifted right by the domain_density, so it is a power of two.
         * The domain x and y positions have (s_log + 1 - dom_density)
           zero
         * bits each, which we don't have to transmit.
         */
        OutputBits(frac_file, dom_number >> (s_log + 1 - dom_density),
                   dom_info[s_log].pos_bits);
    } else {
        /* Tell the decompressor that no affine map follows because
         * this range has been split:
         */
        OutputBit(frac_file, 0);

        /* Split the range into 4 squares and process them
           recursively: */
        compress_range(x,          y,          s_log-1);
        compress_range(x+r_size/2, y,          s_log-1);
        compress_range(x,          y+r_size/2, s_log-1);
        compress_range(x+r_size/2, y+r_size/2, s_log-1);
    }
}

/* ==================================================================
 * Find the best affine mapping from a range to a domain. This is
   done
 * by minimizing the sum of squared errors as a function of the
   contrast
 * and brightness:  sum on all range pixels ri and domain pixels
   di of
 *      square(contrast*domain[di] + brightness - range[ri])
 * and solving the resulting equations to get contrast and
   brightness.
 */
void find_map(rangep, dom, map)
    range_data  *rangep; /* range information (input parameter) */
    domain_data *dom;    /* domain information (input parameter) */
    affine_map  *map;    /* resulting map (output parameter) */
{
    int ry;            /* vertical position inside the range */
    int dy = dom->y >> 1; /* vertical position inside the
                                      domain */
    uns_long rd = 0;   /* sum of range*domain values (scaled
                          by 4) */
    double rd_sum;     /* sum of range*domain values
                         (normalized) */
    double contrast;   /* optimal contrast between range and
                          domain */
    double brightness; /* optimal brightness offset between range
                          and domain */
    double qbrightness;/* brightness after quantization */
    double max_scaled; /* maximum scaled value =
                          contrast*MAX_GREY */
    int r_size = 1 << rangep->s_log;                 /* the range
                                               size */
    double pixels = (double)((long)r_size*r_size); /* total number of
    pixels */

    for (ry = rangep->y; ry < rangep->y + r_size; ry++, dy++) {

        register image_data *r = &range[ry][rangep->x];
        register unsigned   *d = &domain[dy][dom->x >> 1];
        int i = r_size >> 2;

        /* The following loop is the most time consuming part of the
           whole
         * program, so it is unrolled a little. We rely on r_size
           being a
         * multiple of 4 (ranges smaller than 4 don't make sense
           because
         * of the very bad compression). rd cannot overflow with
           unsigned
         * 32-bit arithmetic since MAX_BITS <= 7 implies r_size
           <= 128.
         */
        do {
            rd += (uns_long)(*r++)*(*d++);
            rd += (uns_long)(*r++)*(*d++);
            rd += (uns_long)(*r++)*(*d++);
            rd += (uns_long)(*r++)*(*d++);
        } while (--i != 0);
    }
    rd_sum = 0.25*rd;

    /* Compute and quantize the contrast: */
    contrast = pixels * dom->d_sum2 - dom->d_sum * dom->d_sum;
    if (contrast != 0.0) {
        contrast = (pixels*rd_sum - rangep->r_sum*dom->d_sum)/contrast;
    }
    map->contrast = quantize(contrast, MAX_CONTRAST, MAX_QCONTRAST);

    /* Recompute the contrast as in the decompressor: */
    contrast = dequantize(map->contrast, MAX_CONTRAST, MAX_QCONTRAST);
    /* Compute and quantize the brightness. We actually quantize the
       value
     * (brightness + 255*contrast) to get a positive value:
     *    -contrast*255 <= brightness <= 255
     * so 0 <= brightness + 255*contrast <= 255 + contrast*255
     */
    brightness = (rangep->r_sum - contrast*dom->d_sum)/pixels;
    max_scaled = contrast*MAX_GREY;
    map->brightness = quantize(brightness + max_scaled,
                         max_scaled + MAX_GREY, MAX_QBRIGHTNESS);

    /* Recompute the quantized brightness as in the decompressor: */
    qbrightness = dequantize(map->brightness, max_scaled + MAX_GREY,
                             MAX_QBRIGHTNESS) - max_scaled;

    /* Compute the sum of squared errors, which is the quantity we
       are
     * trying to minimize:
     */
    map->error2 = contrast*(contrast*dom->d_sum2 -
    2.0*rd_sum) + 
    rangep->r_sum2 + qbrightness*pixels*(qbrightness - 2.0*brightness);
}

                /************************************/
                /* Functions used for decompression */
                /************************************/

/*
 * Scale factor for decompression (decompressed size divided by
   original size).
 * Only integer values are supported to simplify the implementation.
 */
int image_scale = 1;

/*
 * An affine map is described by a contrast, a brightness offset, a range
 * and a domain. The contrast and brightness are kept as integer values
 * to speed up the decompression on machines with slow floating point.
 */
typedef struct map_info_struct {
    int contrast;   /* contrast scaled by 16384 (to maintain
                       precision) */
    int brightness; /* brightness offset scaled by 128 */
    int x;          /* horizontal position of the range */
    int y;          /* vertical position of the range */
    int size;       /* range size */
    int dom_x;      /* horizontal position of the domain */
    int dom_y;      /* vertical position of the domain */
    struct map_info_struct *next; /* next map */
} map_info;

map_info *map_head = NULL; /* head of the linked list of all affine
                              maps */

/* ==================================================================
 * This is the main decompression routine.  By the time it gets called,
 * the input and output files have been properly opened, so all it has to
 * do is the decompression.  Note that the decompression routine optionally
 * accepts additional parameters:
 * - the number of iterations, ranging from 1 to 15. The image quality
 *   does not improve much after 8 to 10 iterations. The default is 8.
 * - the scale factor (decompressed size divided by original size)
 */
void ExpandFile(input, output, argc, argv)
    BIT_FILE *input;
    FILE *output;
    int argc;
    char *argv[];
{
    int x_size;         /* horizontal image size */
    int y_size;         /* vertical image size */
    int x_dsize;        /* horizontal size of decompressed image */
    int y_dsize;        /* vertical size of decompressed image */
    int iterations = 8; /* number of iterations */
    int y;              /* current row being written to disk */

    /* 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 `i': iterations  = atoi(*++argv); argc--; break;
            case `s': image_scale = atoi(*++argv); argc--; break;
            default:  fatal_error("Incorrect argument: %s\n", *argv);
        }
    }
    if (image_scale < 1) {
        fatal_error("Incorrect image scale\n");
    }
    /* Read the header of the fractal file: */
    frac_file = input;
    if (InputBits(frac_file, 8) != `F') {
        fatal_error("Bad fractal file format\n");
    }
    x_size = (int)InputBits(frac_file, 16);
    y_size = (int)InputBits(frac_file, 16);
    dom_density = (int)InputBits(frac_file, 2);

    /* Allocate the scaled image: */
    x_dsize = x_size * image_scale;
    y_dsize = y_size * image_scale;
    range = (image_data**)allocate(y_dsize, x_dsize, sizeof
    (image_data));

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

    /* Read all the affine maps, by using the same recursive
       traversal
     * of the image as the compressor:
     */
    traverse_image(0, 0, x_size, y_size, decompress_range);

    /* Iterate all affine maps over an initially random image.
       Since the
     * affine maps are contractive, this process converges.
     */
    while (iterations-- > 0) refine_image();

    /* Smooth the transition between adjacent ranges: */
    average_boundaries();

    /* Write the uncompressed file: */
    for (y = 0; y < y_dsize; y++) {
        if (fwrite(range[y], sizeof(image_data), x_dsize, output)
        != x_dsize) {
            fatal_error("Error writing uncompressed image\n");
        }
    }
    /* Cleanup: */
    free_array((void**)range, y_dsize);
}

/* =================================================================
 * Read the affine map for a range, or split the range if the
   compressor
 * did so in the function compress_range().
 */
void decompress_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; /* range size */
    map_info *map;         /* pointer to affine map information */
    double contrast;       /* contrast between range and domain */
    double brightness;     /* brightness offset between range and
                              domain */
    double max_scaled;     /* maximum scaled value =
                              contrast*MAX_GREY */
    uns_long dom_number;   /* domain number */

    /* Read an affine map if the compressor has written one at this
       point: */
    if (s_log == MIN_BITS || InputBit(frac_file)) {

        map = (map_info *)xalloc(sizeof(map_info));
        map->next = map_head;
        map_head = map;
        map->x = x;
        map->y = y;
        map->size = r_size;
        map->contrast   = (int)InputBits(frac_file, CONTRAST_BITS);
        map->brightness = (int)InputBits(frac_file, BRIGHTNESS_BITS);

        contrast = dequantize(map->contrast, MAX_CONTRAST,
        MAX_QCONTRAST);
        max_scaled = contrast*MAX_GREY;
        brightness = dequantize(map->brightness, max_scaled +
        MAX_GREY,
                                MAX_QBRIGHTNESS) - max_scaled;

        /* Scale the brightness by 128 to maintain precision later,
           while
         * avoiding overflow with 16-bit arithmetic:
         *     -255 <= -contrast*255 <= brightness <= 255
         * so -32767 < brightness*128 < 32767
         */
        map->brightness = (int)(brightness*128.0);

        /* When the contrast is null, the compressor did not encode
           the
         * domain number:
         */
        if (map->contrast != 0) {

            /* Scale the contrast by 16384 to maintain precision
               later.
             *   0.0 <= contrast <= 1.0 so 0 <= contrast*16384
                 <= 16384
             */
            map->contrast = (int)(contrast*16384.0);

            /* Read the domain number, and add the zero bits that
               the
             * compressor did not transmit:
             */
            dom_number = InputBits(frac_file, dom_info[s_log]
            .pos_bits);
            map->dom_x = (int)(dom_number % dom_info[s_log]
            .x_domains)
                          << (s_log + 1 - dom_density);
            map->dom_y = (int)(dom_number / dom_info[s_log]
            .x_domains)
                          << (s_log + 1 - dom_density);
        } else {
            /* For a null contrast, use an arbitrary domain: */
            map->dom_x = map->dom_y = 0;
        }

        /* Scale the range and domain if necessary. This
           implementation
         * uses only an integer scale to make sure that the union