Someone requested this code from me some time ago. I got it with about 20 compression routines last year from a book I bought. But this is the most interested (IMHO). The book talks about a compression ratio of about 97 - 99% on images.
The bad thing with this code is, that I couldn't compile it :(
frac.c
bitio.h
errhand.h
errhand.c
bitio.c
The bad thing with this code is, that I couldn't compile it :(
frac.c
/************************** 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;
}
bitio.h
/************************** Start of BITIO.H *************************/
#ifndef _BITIO_H
#define _BITIO_H
#include <stdio.h>
typedef struct bit_file {
FILE *file;
unsigned char mask;
int rack;
int pacifier_counter;
} BIT_FILE;
#ifdef __STDC__
BIT_FILE *OpenInputBitFile( char *name );
BIT_FILE *OpenOutputBitFile( char *name );
void OutputBit( BIT_FILE *bit_file, int bit );
void OutputBits( BIT_FILE *bit_file,
unsigned long code, int count );
int InputBit( BIT_FILE *bit_file );
unsigned long InputBits( BIT_FILE *bit_file, int bit_count );
void CloseInputBitFile( BIT_FILE *bit_file );
void CloseOutputBitFile( BIT_FILE *bit_file );
void FilePrintBinary( FILE *file, unsigned int code, int bits );
#else /* __STDC__ */
BIT_FILE *OpenInputBitFile();
BIT_FILE *OpenOutputBitFile();
void OutputBit();
void OutputBits();
int InputBit();
unsigned long InputBits();
void CloseInputBitFile();
void CloseOutputBitFile();
void FilePrintBinary();
#endif /* __STDC__ */
#endif /* _BITIO_H */
/*************************** End of BITIO.H **************************/
errhand.h
/************************* Start of ERRHAND.H ************************/
#ifndef _ERRHAND_H
#define _ERRHAND_H
#ifdef __STDC__
void fatal_error( char *fmt, ... );
#else /* __STDC__ */
void fatal_error();
#endif /* __STDC__ */
#endif /* _ERRHAND_H */
/************************** End of ERRHAND.H *************************/
errhand.c
/************************* Start of ERRHAND.C ************************
*
* This is a general purpose error handler used with every program in
* the book.
*/
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include "errhand.h"
#ifdef __STDC__
void fatal_error( char *fmt, ... )
#else
#ifdef __UNIX__
void fatal_error( fmt, va_alist )
char *fmt;
va_dcl
#else
void fatal_error( fmt )
char *fmt;
#endif
#endif
{
va_list argptr;
va_start( argptr, fmt );
printf( "Fatal error: " );
vprintf( fmt, argptr );
va_end( argptr );
exit( -1 );
}
/************************** End of ERRHAND.C *************************/
bitio.c
/************************** Start of BITIO.C *************************
*
* This utility file contains all of the routines needed to impement
* bit oriented routines under either ANSI or K&R C. It needs to be
* linked with every program used in the entire book.
*
*/
#include <stdio.h>
#include <stdlib.h>
#include "bitio.h"
#include "errhand.h"
#define PACIFIER_COUNT 2047
BIT_FILE *OpenOutputBitFile( name )
char *name;
{
BIT_FILE *bit_file;
bit_file = (BIT_FILE *) calloc( 1, sizeof( BIT_FILE ) );
if ( bit_file == NULL )
return( bit_file );
bit_file->file = fopen( name, "wb" );
bit_file->rack = 0;
bit_file->mask = 0x80;
bit_file->pacifier_counter = 0;
return( bit_file );
}
BIT_FILE *OpenInputBitFile( name )
char *name;
{
BIT_FILE *bit_file;
bit_file = (BIT_FILE *) calloc( 1, sizeof( BIT_FILE ) );
if ( bit_file == NULL )
return( bit_file );
bit_file->file = fopen( name, "rb" );
bit_file->rack = 0;
bit_file->mask = 0x80;
bit_file->pacifier_counter = 0;
return( bit_file );
}
void CloseOutputBitFile( bit_file )
BIT_FILE *bit_file;
{
if ( bit_file->mask != 0x80 )
if ( putc( bit_file->rack, bit_file->file ) != bit_file->rack )
fatal_error( "Fatal error in CloseBitFile!\n" );
fclose( bit_file->file );
free( (char *) bit_file );
}
void CloseInputBitFile( bit_file )
BIT_FILE *bit_file;
{
fclose( bit_file->file );
free( (char *) bit_file );
}
void OutputBit( bit_file, bit )
BIT_FILE *bit_file;
int bit;
{
if ( bit )
bit_file->rack |= bit_file->mask;
bit_file->mask >>= 1;
if ( bit_file->mask == 0 ) {
if ( putc( bit_file->rack, bit_file->file ) != bit_file->rack )
fatal_error( "Fatal error in OutputBit!\n" );
else
if ( ( bit_file->pacifier_counter++ & PACIFIER_COUNT ) == 0 )
putc( '.', stdout );
bit_file->rack = 0;
bit_file->mask = 0x80;
}
}
void OutputBits( bit_file, code, count )
BIT_FILE *bit_file;
unsigned long code;
int count;
{
unsigned long mask;
mask = 1L << ( count - 1 );
while ( mask != 0) {
if ( mask & code )
bit_file->rack |= bit_file->mask;
bit_file->mask >>= 1;
if ( bit_file->mask == 0 ) {
if ( putc( bit_file->rack, bit_file->file ) != bit_file->rack )
fatal_error( "Fatal error in OutputBit!\n" );
else if ( ( bit_file->pacifier_counter++ & PACIFIER_COUNT ) == 0 )
putc( '.', stdout );
bit_file->rack = 0;
bit_file->mask = 0x80;
}
mask >>= 1;
}
}
int InputBit( bit_file )
BIT_FILE *bit_file;
{
int value;
if ( bit_file->mask == 0x80 ) {
bit_file->rack = getc( bit_file->file );
if ( bit_file->rack == EOF )
fatal_error( "Fatal error in InputBit!\n" );
if ( ( bit_file->pacifier_counter++ & PACIFIER_COUNT ) == 0 )
putc( '.', stdout );
}
value = bit_file->rack & bit_file->mask;
bit_file->mask >>= 1;
if ( bit_file->mask == 0 )
bit_file->mask = 0x80;
return( value ? 1 : 0 );
}
unsigned long InputBits( bit_file, bit_count )
BIT_FILE *bit_file;
int bit_count;
{
unsigned long mask;
unsigned long return_value;
mask = 1L << ( bit_count - 1 );
return_value = 0;
while ( mask != 0) {
if ( bit_file->mask == 0x80 ) {
bit_file->rack = getc( bit_file->file );
if ( bit_file->rack == EOF )
fatal_error( "Fatal error in InputBit!\n" );
if ( ( bit_file->pacifier_counter++ & PACIFIER_COUNT ) == 0 )
putc( '.', stdout );
}
if ( bit_file->rack & bit_file->mask )
return_value |= mask;
mask >>= 1;
bit_file->mask >>= 1;
if ( bit_file->mask == 0 )
bit_file->mask = 0x80;
}
return( return_value );
}
void FilePrintBinary( file, code, bits )
FILE *file;
unsigned int code;
int bits;
{
unsigned int mask;
mask = 1 << ( bits - 1 );
while ( mask != 0 ) {
if ( code & mask )
fputc( '1', file );
else
fputc( '0', file );
mask >>= 1;
}
}
/*************************** End of BITIO.C **************************/
I have that book. :) I like that you can enlarge the image and you don't get the typical pixelization. I think a zip file might have been a better choice. Someone on the hugi size coding compo said they had an assembly version that was very small in size! What compiler are you using?
>I have that book. :)
"The Data Compression Book" by Nelson Gailly? ;)
>I think a zip file might have been a better choice.
Sure. But it would be interested to see it working.
>What compiler are you using?
I tried it with VC++, LCC-Win32 and GCC under Linux.
"The Data Compression Book" by Nelson Gailly? ;)
>I think a zip file might have been a better choice.
Sure. But it would be interested to see it working.
>What compiler are you using?
I tried it with VC++, LCC-Win32 and GCC under Linux.
Hi
bitrake:
what do you mean :
can you zoom on the image and it will still keep on the quality?
if i understood correctly
can you publish or send me the source for that or the program the that does it ?
bye
eko
bitrake:
what do you mean :
I like that you can enlarge the image and you don't get the typical pixelization.
can you zoom on the image and it will still keep on the quality?
if i understood correctly
can you publish or send me the source for that or the program the that does it ?
bye
eko
This is by the vary nature of the algorithm. :) I don't know if I have the sources, but I'll look. There is C code HERE. If I remember right, all I changed was the decode resolution and it worked. I've been meaning to do an asm version. I started it before, but ran up against time constrains of a contest and had to use another compression algo.
thanks
eko
eko
Try THIS. VC++ project for fractal compression program.