477-505.html

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

           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;
}
</PRE>
<!--  END CODE //-->
<P><BR></P>
<CENTER>
<TABLE BORDER>
<TR>
<TD><A HREF="474-477.html">Previous</A></TD>
<TD><A HREF="../ewtoc.html">Table of Contents</A></TD>
<TD><A HREF="506-509.html">Next</A></TD>
</TR>
</TABLE>
</CENTER>


</TD>
</TR>
</TABLE>

</BODY>
</HTML>