Previous | Table of Contents | Next

#### The decompression module

The decompression code uses the same traverse_image() function as the compressor, to ensure that both use exactly the same partitioning strategy. The only difference is that during the traversal, the decoder only reads the encoded affine maps; it does not actually generates the output image at this step.

Once the affine maps have been read in, the decoder starts with a random initial image and calls refine_image() a number of times selected by the user. For each iteration, refine_image() goes through all the maps and applies them to the current image. The “pure” method would compute a separate new image and then swap the roles of the old and new image. However the convergence towards the final image happens to be quicker if we overwrite the same image while applying the affine maps; for the same quality of reconstructed image we need fewer iterations. Overwriting the same image also reduces the memory requirements since we need only a buffer for one image instead of two.

After the image has been reconstructed, the decoder calls the function average_boundaries() to smooth the transition between adjacent ranges. This is only done for large ranges (8x8 and 16x16), since averaging the small 4x4 ranges can degrade the image quality.

### The complete code listing

The complete listing of FRAC.C follows.

    /*********************** Start of FRAC.C  ***********************
     *
     * This is the FRAC module, which implements a graphics fractal
       compression
     * program.  It needs to be linked with the standard support routines.
     * Copyright 1995 Jean-loup Gailly
     */

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

    #ifdef unix
    #  define float double /* better accuracy but more memory usage */
    #endif

    char *CompressionName = "Fractal compression";
    char *Usage =
    "infile outfile [-q quality] [-d density] [-h h_size] [-v v_size]\n\
       quality from 1..20, domain density from 0..2\n\
    Decompression parameters:\n\
        infile outfile [-i iterations] [-s scale]\n";

    typedef unsigned char image_data;
    typedef unsigned long uns_long;

    /*
     * Maximum gray level in an image
     */
    #define MAX_GREY 255

    /*
     * Number of classes. Each class corresponds to one specific ordering
     * of the image brightness in the four quadrants of a range or domain.
     * There are 4*3*2 = 24 classes.
     */
    #define NCLASSES 24

    /*
     * Minimum and maximum number of bits for the side of a range. The actual
     * range sizes are between 1<<MIN_BITS and 1<<MAX_BITS. To simplify the
     * implementation and avoid ridiculously small ranges, MIN_BITS must be
       >= 2.
     * This implementation also requires MAX_BITS <= 7.
     */
    #define MIN_BITS 2
    #define MAX_BITS 4

    /*
     * Maximum contrast factor in a range to domain mapping.
     */
    #define MAX_CONTRAST 1.0

    /*
     * Bit sizes for encodings of contrast and brightness in an affine map.
     * Using smaller sizes increases compression and degrades image quality.
     */
    #define CONTRAST_BITS    4
    #define BRIGHTNESS_BITS  6

    #define MAX_QCONTRAST   ((1<<CONTRAST_BITS)-1)   /* max quantized
                                                        contrast */
    #define MAX_QBRIGHTNESS ((1<<BRIGHTNESS_BITS)-1) /* max quantized
                                                        brightness */

    /*
     * De-quantize an integer value in the range 0 .. imax to the range
       0.0 .. max
     * while preserving the mapping 0 -> 0.0 and imax -> max.
     */
    #define dequantize(value, max, imax) ((double)(value)*(max)/(double)imax)

    /*
     * Compute the square of a pixel value and return the result as unsigned
       long
     */
    #define square(pixel) (uns_long)(pixel)*(pixel)

    /*
     * Range data: range[i][j] is the brightness at row i and column j
     */
    image_data **range;

    /*
     * Domain data, summed over 4 pixels: domain[i][j] is the sum of the
     * pixel values at (2j, 2i), (2j+1, 2i), (2j, 2i+1) and (2j+1, 2i+1)
     */
    unsigned **domain;

    /*
     * Cumulative range data, kept only for pixels of even coordinates.
     * cum_range[i][j] is the sum of all pixel values strictly above and to
       the
     * left of pixel (2j, 2i). In particular, cum_range[y_size/2][x_size/2]
       is
     * the sum of all pixel values in the image. This table is also used for
     * the cumulative domain data.
     */
    uns_long **cum_range;

    /*
     * Cumulative squared range data, kept only for pixels of even
       coordinates.
     * cum_range2[i][j] is the sum of the squares of all pixel values
       strictly
     * above and to the left of pixel (2j, 2i). In particular,
     * cum_range2[y_size/2][x_size/2] is the sum of all the squared pixel
       values
     * in the image.
     */
    float **cum_range2;

    /*
     * Cumulative squared domain data. cum_domain2[i][j] is the sum of the
       squares
     * of all domain values strictly above and to the left of domain (j,i),
       which
     * corresponds to pixel (2j, 2i). The values in cum_domain2 are scaled
       by
     * a factor of 16.0 to avoid some multiplications.
     */
    float **cum_domain2;

    /*
     * Domain density: domains of size s*s are located every
       (s>>dom_density)
     * pixels. The density factor can range from 0 to 2 (smallest domains
     * have a size of 8 and must start on even pixels boundaries). Density
     * factors 1 and 2 get better image quality but significantly slow
     * down compression.
     */
    int dom_density = 0;

    /*
     * Maximum tolerated mean square error between original image and
     * reconstructed uncompressed image.
     */
    double max_error2;

    /*
     * The fractal (compressed) file
     */
    BIT_FILE *frac_file;

    /*
     * Information common to all domains of a certain size: info[s]
       describes
     * domains of size 1<<(s+1), corresponding to ranges of size 1<<s
     */
    struct domain_info {
        int pos_bits;   /* Number of bits required to encode a domain
                           position */
        int x_domains;  /* Number of domains in x (horizontal) dimension */
    } dom_info[MAX_BITS+1];

    /*
     * Each domain is described by a `domain_data' structure.
     * domain_head[c][s] is the head of the list of domains of class c
     * and size 1<<(s+1) (corresponding to ranges of size 1<<s).
     */
    typedef struct domain_struct {
        int x;                      /* horizontal position */
        int y;                      /* vertical position */
        float d_sum;                /* sum of all values in the domain */
        float d_sum2;               /* sum of all squared values in the
                                       domain */
        struct domain_struct *next; /* next domain in same class */
    } domain_data;

    domain_data *domain_head[NCLASSES][MAX_BITS+1];

    /*
     * Ranges are described by a `range_data' structure. This structure
     * is computed on the fly for each range as it is compressed.
     */
    typedef struct range_struct {
        int x;         /* horizontal position */
        int y;         /* vertical position */
        int s_log;     /* log base 2 of the range size */
        double r_sum;  /* sum of all values in the range */
        double r_sum2; /* sum of all squared values in the range */
    } range_data;

    /*
     * Range to domain mappings are described by an `affine_map' structure.
     */
    typedef struct map_struct {
        int    contrast;   /* quantized best contrast between range and
                              domain */
        int    brightness; /* quantized best brightness offset */
        double error2;    /* sum of squared differences between range
                             and domain */
    } affine_map;

    /*
     * Function prototypes for both ANSI and K&R.
     */
    #ifdef __STDC__
    #  define OF(args)  args
    #else
    #  define OF(args)  ()
    #endif

    /*
     * Functions used for compression
     */
    void CompressFile  OF((FILE *input, BIT_FILE *output, int argc, char
    *argv[]));
    void compress_init OF((int x_size, int y_size, FILE *image_file));
    void compress_cleanup OF((int y_size));
    void classify_domains OF((int x_size, int y_size, int s));
    int  find_class       OF((int x, int y, int size));
    void compress_range   OF((int x, int y, int s_log));
    void find_map  OF((range_data *rangep, domain_data *dom, affine_map
    *map));

    /*
     * Functions used for decompression
     */
    void ExpandFile OF((BIT_FILE *input, FILE *output, int argc, char
    *argv[]));
    void decompress_range   OF((int x, int y, int s_log));
    void refine_image       OF((void));
    void average_boundaries OF((void));

    /*
     * Functions common to compression and decompression
     */
    typedef void (*process_func) OF((int x, int y, int s_log));

    void traverse_image OF((int x, int y, int x_size, int y_size,
                            process_func process));
    int  quantize       OF((double value, double max, int imax));
    void dominfo_init   OF((int x_size, int y_size, int density));
    void *xalloc        OF((unsigned size));
    void **allocate     OF((int rows, int columns, int elem_size));
    void free_array     OF((void **array, int rows));
    int  bitlength      OF((uns_long val));

                    /**********************************/
                    /* Functions used for compression */
                    /**********************************/

    /* =====================================================================
     * This is the main compression routine.  By the time it gets called,
     * the input and output files have been properly opened, so all it
       has to
     * do is the compression.  Note that the compression routine optionally
     * accepts additional parameters:
     * - the quality value, ranging from 0 to 20. It is used as average
       tolerated
     *   error between the original image and its uncompressed version.
         (Non
     *   integer values are also accepted.)
     * - 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 ||
            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
               of all
             * ranges is exactly the scaled image, that ranges never
               overlap,
             * and that all range sizes are even.
             */
            if (image_scale != 1) {
                map->x *= image_scale;
                map->y *= image_scale;
                map->size *= image_scale;
                map->dom_x *= image_scale;
                map->dom_y *= image_scale;
            }
        } else {
            /* Split the range into 4 squares and process them
               recursively
             * as in the compressor:
             */
            decompress_range(x,          y,          s_log-1);
            decompress_range(x+r_size/2, y,          s_log-1);
            decompress_range(x,          y+r_size/2, s_log-1);
            decompress_range(x+r_size/2, y+r_size/2, s_log-1);
        }
    }

    /* ===================================================================
     * Refine the image by applying one round of all affine maps on
       the
     * image. The "pure" method would compute a separate new image and
       then
     * copy it to the original image. However the convergence towards
       the
     * final image happens to be quicker if we overwrite the same image
     * while applying the affine maps; for the same quality of
       reconstructed
     * image we need fewer iterations. Overwriting the same image also
     * reduces the memory requirements.
     */
    void refine_image()
    {
        map_info *map;   /* pointer to current affine map */
        long brightness; /* brightness offset of the map, scaled by
                            65536 */
        long val;        /* new pixel value */
        int y;           /* vertical position in range */
        int dom_y;       /* vertical position in domain */
        int j;

        for (map = map_head; map != NULL; map = map->next) {

            /* map->brightness is scaled by 128, so scale it again
               by 512 to
             * get a total scale factor of 65536:
             */
            brightness = (long)map->brightness << 9;

            dom_y = map->dom_y;
            for (y = map->y; y < map->y + map->size; y++) {

                /* The following loop is the most time consuming, so
                   we move
                 * some address calculations outside the loop:
                 */
                image_data *r  = &range[y][map->x];
                image_data *d  = &range[dom_y++][map->dom_x];
                image_data *d1 = &range[dom_y++][map->dom_x];
                j = map->size;
                do {
                    val  = *d++ + *d1++;
                    val += *d++ + *d1++;
                    /* val is now scaled by 4 and map->contrast is
                       scaled by
                       16384,
                     * so val * map->contrast will be scaled by
                       65536.
                     */
                    val = val * map->contrast + brightness;
                    if (val < 0) val = 0;
                    val >>= 16; /* get rid of the 65536 scaling */
                    if (val >= MAX_GREY) val = MAX_GREY;

                    *r++ = (image_data)val;
                } while (--j != 0);
            }
        }
    }

    /* =================================================================
     * Go through all ranges to smooth the transition between adjacent
     * ranges, except those of minimal size.
     */
    void average_boundaries()
    {
        map_info *map;   /* pointer to current affine map */
        unsigned val;    /* sum of pixel value for current and adjacent
                            ranges */
        int x;           /* horizontal position in current range */
        int y;           /* vertical position in current range */

        for (map = map_head; map != NULL; map = map->next) {

            if (map->size == (1<<MIN_BITS)) continue; /* range
                                                        too small */
            if (map->x > 1) {
                /* Smooth the left boundary of the range and the
                   right boundary
                 * of the adjacent range(s) to the left:
                 */
                for (y = map->y; y < map->y + map->size; y++) {
                     val  = range[y][map->x - 1] + range[y][map->x];
                     range[y][map->x - 1] =
                         (image_data)((range[y][map->x - 2] + val)/3);
                     range[y][map->x] =
                         (image_data)((val + range[y][map->x + 1])/3);
                }
            }
            if (map->y > 1)  {
                /* Smooth the top boundary of the range and the bottom
                   boundary
                 * of the range(s) above:
                 */
                for (x = map->x; x < map->x + map->size; x++) {
                     val  = range[map->y - 1][x] + range[map->y][x];
                     range[map->y - 1][x] =
                         (image_data)((range[map->y - 2][x] + val)/3);
                     range[map->y][x] =
                         (image_data)((val + range[map->y + 1][x])/3);
                }
            }
        }
    }

            /*****************************************************/
            /* Functions common to compression and decompression */
            /*****************************************************/

    /* ===================================================================
     * Split a rectangle sub-image into a square and potentially two
       rectangles,
     * then split the square and rectangles recursively if necessary.
       To simplify
     * the algorithm, the size of the square is chosen as a power of
       two.
     * If the square if small enough as a range, call the appropriate
       compression
     * or decompression function for this range.
     * IN assertions: x, y, x_size and y_size are multiple of 4.
     */
    void traverse_image(x, y, x_size, y_size, process)
        int x, y;             /* sub-image horizontal and vertical
                                 position */
        int x_size, y_size;   /* sub-image horizontal and vertical
                                 sizes */
        process_func process; /* the compression or decompression
                                 function */
    {
        int s_size;  /* size of the square; s_size = 1<<s_log */
        int s_log;   /* log base 2 of this size */

        s_log = bitlength(x_size < y_size ? (uns_long)x_size :
        (uns_long)y_size)-1;
        s_size = 1 << s_log;
        /* Since x_size and y_size are >= 4, s_log >= MIN_BITS */

        /* Split the square recursively if it is too large for a
           range: */
        if (s_log > MAX_BITS) {
            traverse_image(x,          y,          s_size/2, s_size/2,
            process);
            traverse_image(x+s_size/2, y,          s_size/2, s_size/2,
            process);
            traverse_image(x,          y+s_size/2, s_size/2, s_size/2,
            process);
            traverse_image(x+s_size/2, y+s_size/2, s_size/2, s_size/2,
            process);
        } else {
            /* Compress or decompress the square as a range: */
            (*process)(x, y, s_log);
        }

        /* Traverse the rectangle on the right of the square: */
        if (x_size > s_size) {
            traverse_image(x + s_size, y, x_size - s_size, y_size,
             process);

            /* Since x_size and s_size are multiple of 4, x + s_size and
             * x_size - s_size are also multiple of 4.
             */
        }
        /* Traverse the rectangle below the square: */
        if (y_size > s_size) {
            traverse_image(x, y + s_size, s_size, y_size - s_size,
            process);
        }
    }

    /* =================================================================
     * Initialize the domain information dom_info. This must be done
       in the
     * same manner in the compressor and the decompressor.
     */
    void dominfo_init(x_size, y_size, density)
        int x_size;       /* horizontal size of original image */
        int y_size;       /* vertical size of original image */
        int density;      /* domain density (0 to 2) */
    {
        int s;            /* size index for domains; their size is
                             1<<(s+1) */

        for (s = MIN_BITS; s <= MAX_BITS; s++) {
            int y_domains;            /* number of domains vertically */
            int dom_size = 1<<(s+1);  /* domain size */

            /* The distance between two domains is the domain size
               1<<(s+1)
             * shifted right by the domain density, so it is a power
               of two.
             */
            dom_info[s].x_domains = ((x_size - dom_size)>>(s + 1 -
            density)) + 1;
            y_domains             = ((y_size - dom_size)>>(s + 1 -
           density)) + 1;

            /* Number of bits required to encode a domain position: */
            dom_info[s].pos_bits =  bitlength
                ((uns_long)dom_info[s].x_domains * y_domains - 1);
        }
    }

    /* ==============================================================
     * Quantize a value in the range 0.0 .. max to the range 0..imax
     * ensuring that 0.0 is encoded as 0 and max as imax.
     */
    int quantize(value, max, imax)
        double value, max;
        int imax;
    {
        int ival = (int) floor((value/max)*(double)(imax+1));

        if (ival < 0) return 0;
        if (ival > imax) return imax;
        return ival;
    }

    /* ==============================================================
     * Allocate memory and check that the allocation was successful.
     */
    void *xalloc(size)
        unsigned size;
    {
        void *p = malloc(size);

        if (p == NULL) {
            fatal_error("insufficient memory\n");
        }
        return p;
    }

    /* ==============================================================
     * Allocate a two dimensional array. For portability to 16-bit
     * architectures with segments limited to 64K, we allocate one
     * array per row, so the two dimensional array is allocated
     * as an array of arrays.
     */
    void **allocate(rows, columns, elem_size)
        int rows;      /* number of rows */
        int columns;   /* number of columns */
        int elem_size; /* element size */
    {
        int row;
        void **array = (void**)xalloc(rows * sizeof(void *));

        for (row = 0; row < rows; row++) {
            array[row] = (void*)xalloc(columns * elem_size);
        }
        return array;
    }

    /* ==========================================================
     * Free a two dimensional array allocated as a set of rows.
     */
    void free_array(array, rows)
        void **array;  /* the two-dimensional array */
        int rows;      /* number of rows */
    {
        int row;
        for (row = 0; row < rows; row++) {
            free(array[row]);
        }
    }

    /* =============================================================
     * Return the number of bits needed to represent an integer:
     * 0 to 1 -> 1,
     * 2 to 3 -> 2,
     * 3 to 7 -> 3, etc...
     * This function could be made faster with a lookup table.
     */
    int  bitlength(val)
        uns_long val;
    {

        int bits = 1;

        if (val > 0xffff) bits += 16, val >>= 16;
        if (val > 0xff)   bits += 8,  val >>= 8;
        if (val > 0xf)    bits += 4,  val >>= 4;
        if (val > 0x3)    bits += 2,  val >>= 2;
        if (val > 0x1)    bits += 1;
        return bits;
    }

Previous | Table of Contents | Next