| Previous | Table of Contents | Next |
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 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 |