enc.c

a small program by Yuval Fisher that has gone a long way to revealing some of the secrets of fracal

/**************************************************************************/
/* Encode a byte image using a fractal scheme with a quadtree partition   */
/*                                                                        */
/*       Copyright 1993,1994 Yuval Fisher. All rights reserved.           */
/*                                                                        */
/* Version 0.03 3/14/94                                                   */
/**************************************************************************/

#include <stdio.h>
#include <math.h>

#define DEBUG 0
#define GREY_LEVELS 255

#define bound(a)   ((a) < 0.0 ? 0 : ((a)>255.0? 255 : a))
#define IMAGE_TYPE unsigned char /* may be different in some applications */

/* various function declarations to keep compiler warnings away. ANSI     */
/* prototypes can go here, for the hearty.                                */
void fatal();
char *malloc();
char *strcpy(); 

/* The following #define allocates an hsize x vsize  matrix of type TYPE */
#define matrix_allocate(matrix, hsize, vsize, TYPE) {\
    TYPE *imptr; \
    int _i; \
    matrix = (TYPE **)malloc((vsize)*sizeof(TYPE *));\
    imptr = (TYPE*)malloc((long)(hsize)*(long)(vsize)*sizeof(TYPE));\
    if (imptr == NULL) \
        fatal("\nNo memory in matrix allocate."); \
    for (_i = 0; _i<vsize; ++_i, imptr += hsize) \
        matrix[_i] = imptr; \
 }

#define swap(a,b,TYPE)           {TYPE _temp; _temp=b; b=a; a= _temp;}

IMAGE_TYPE **image;          /* The input image data                      */ 
double     **domimage[4];    /* Decimated input image used for domains    */

double max_scale = 1.0;      /* Maximum allowable grey level scale factor */

int    s_bits = 5,           /* Number of bits used to store scale factor */
       o_bits = 7,           /* Number of bits used to store offset       */
       min_part = 4,         /* Min and max _part determine a range of    */
       max_part = 6,         /* Range sizes from hsize>>min to hsize>>max */
       dom_step = 1,         /* Density of domains relative to size       */
       dom_step_type = 0,    /* Flag for dom_step a multiplier or divisor */
       dom_type = 0,         /* Method of generating domain pool 0,1,2..  */
       only_positive = 0,    /* A flag specifying use of positive scaling */
       subclass_search = 0,  /* A flag specifying classes searched        */
       fullclass_search = 0, /* A flag specifying classes searched        */
       *bits_needed,         /* Number of bits to encode domain position. */
       zero_ialpha,          /* The const ialpha when alpha = 0           */
       max_exponent;         /* The max power of 2 side of square image   */
                             /* that fits in our input image.             */


                             /* The class_transform gives the transforms  */
                             /* between classification numbers for        */
                             /* negative scaling values, when brightest   */
                             /* becomes darkest, etc...                   */
int    class_transform[2][24] = {23,17,21,11,15,9,22,16,19,5,13,3,20,10,18,
                                 4,7,1,14,8,12,2,6,0,
                                 16,22,10,20,8,14,17,23,4,18,2,12,11,21,5,
                                 19,0,6,9,15,3,13,1,7};

                             /* rot_transform gives the rotations for     */
                             /* domains with negative scalings.           */
int    rot_transform[2][8] = {7,4,5,6,1,2,3,0, 2,3,0,1,6,7,4,5};

struct domain_data {
       int *no_h_domains,    /* The number of domains horizontally for    */
           *no_v_domains,    /* each size.                                */
           *domain_hsize,    /* The size of the domain.                   */
           *domain_vsize,    /* The size of the domain.                   */
           *domain_hstep,    /* The density of the domains.               */
           *domain_vstep;    /* The density of the domains.               */
struct domain_pixels {       /* This is a three (sigh) index array that   */
       int dom_x, dom_y;     /* dynamically allocated. The first index is */
       double sum,sum2;      /* the domain size, the other are two its    */
       int sym;              /* position. It contains the sum and sum^2   */
} ***pixel;                  /* of the pixel values in the domains, which */
} domain;                    /* are computed just once.                   */ 
                             

struct classified_domain {             /* This is a list which containes  */
       struct domain_pixels *the;      /* pointers to  the domain data    */
       struct classified_domain *next; /* in the structure above. There   */
} **the_domain[3][24];                 /* are three classes with 24 sub-  */
                                       /* classes. Using this array, only */
                                       /* domains and ranges in the same  */
                                       /* class are compared..            */
                                       /* The first pointer points to the */
                                       /* domain size the the second to   */
                                       /* list of domains.                */

FILE *output;                /* Output FILE containing compressed data    */

main(argc,argv)
int argc;
char **argv;
/* Usage: quadfrac [tol [inputfilename [outputfilename [hsize [vsize]]]]]  */
{
    /* Defaults are set initially */
    double         tol = 8.0;            /* Tolerance value for quadtree.  */
    char           inputfilename[200];
    char           outputfilename[200];
    int            i,j,k,
                   hsize = -1,           /* The horizontal and vertical    */
                   vsize = -1;           /* size of the input image.       */
    long           stripchar=0;          /* chars to ignore in input file. */
    FILE           *input; 

    inputfilename[0] = 1;  /* We initially set the input to this and */
    outputfilename[0] = 1; /* then check if the input/output names   */
                           /* have been set below.                   */

    /* scan through the input line and read in the arguments */
    for (i=1; i<argc; ++i) 
       if (argv[i][0] != '-' ) 
            if (inputfilename[0] == 1)
                strcpy(inputfilename, argv[i]);
            else if (outputfilename[0] == 1)
                strcpy(outputfilename, argv[i]);
            else;
       else { /* we have a flag */
            if (strlen(argv[i]) == 1) break;
            switch(argv[i][1]) {
               case 't': tol = atof(argv[++i]);
                         break;
               case 'S': stripchar = atoi(argv[++i]);
                         break;
               case 'x': 
               case 'w': hsize = atoi(argv[++i]);
                         break;
               case 'y': 
               case 'h': vsize = atoi(argv[++i]);
                         break;
               case 'D': dom_type = atoi(argv[++i]);
                         break;
               case 'd': dom_step = atoi(argv[++i]);
                         if (dom_step < 0 || dom_step > 15) 
                            fatal("\n Bad domain step.");
                         break;
               case 's': s_bits = atoi(argv[++i]);
                         break;
               case 'o': o_bits = atoi(argv[++i]);
                         break;
               case 'm': min_part = atoi(argv[++i]);
                         break;
               case 'M': max_part = atoi(argv[++i]);
                         break;
               case 'e': dom_step_type= 1;
                         break;
               case 'p': only_positive = 1;
                         break;
               case 'f': subclass_search = 1;
                         break;
               case 'F': fullclass_search = 1;
                         break;
               case 'N': max_scale = atof(argv[++i]);
                         break;
               case '?':
               case 'H':
               default:
		   printf("\nUsage: enc -[options] [inputfile [outputfile]]");
		   printf("\nOptions are: (# = number), defaults show in ()");
		   printf("\n  -t # tolerance criterion for fidelity. (%lf)", tol);
           printf("\n  -m # minimum quadtree partitions.  (%d)",min_part);
		   printf("\n  -M # maximum quadtree partitions. (%d)",max_part);
		   printf("\n  -S # number of input bytes to ignore. (%ld)",stripchar);
		   printf("\n  -w # width (horizontal size) of input data. (256)");
		   printf("\n  -h # height (vertical size) of input data. (256)");
		   printf("\n  -d # domain step size. (%d)", dom_step);
		   printf("\n  -D # method {0,1,2} for domain pool (%d)",dom_type);
		   printf("\n  -s # number of scaling quantizing bits. (%d)",s_bits);
		   printf("\n  -o # number of offset quantizing bits. (%d)",o_bits);
		   printf("\n  -N # maximum scaling in encoding. (%lf)",max_scale);
		   printf("\n  -e   domain step as multiplier not divisor. (off)");
		   printf("\n  -p   use only positive scaling (for speed). (off)");
		   printf("\n  -f   search 24 domain classes (for fidelity). (off)");
		   printf("\n  -F   search 3 domain classes (for fidelity). (off)");
		   fatal("\n       -F and -f can be used together.");
            }
       }
   
    if (inputfilename[0] == 1) strcpy(inputfilename, "lenna.dat");
    if (outputfilename[0] == 1) strcpy(outputfilename, "lenna.trn");

    if (hsize == -1) 
        if (vsize == -1) hsize = vsize = 256;
        else hsize = vsize;
    else 
        if (vsize == -1) vsize = hsize;

    /* allocate memory for the input image. Allocating one chunck saves  */
    /* work and time later.                                              */
    matrix_allocate(image, hsize, vsize, IMAGE_TYPE)
    matrix_allocate(domimage[0], hsize/2, vsize/2, double)
    matrix_allocate(domimage[1], hsize/2, vsize/2, double)
    matrix_allocate(domimage[2], hsize/2, vsize/2, double)
    matrix_allocate(domimage[3], hsize/2, vsize/2, double)

    /* max_ & min_ part are variable, so this must be run time allocated */
    bits_needed = (int *)malloc(sizeof(int)*(1+max_part-min_part));

    if ((input = fopen(inputfilename, "r")) == NULL)
        fatal("Can't open input file.");

    /* skip the first  stripchar  chars */ 
    fseek(input, stripchar, 0); 
    i = fread(image[0], sizeof(IMAGE_TYPE), hsize*vsize, input);
    fclose(input);

    if (i < hsize*vsize) 
        fatal("Not enough input data in the input file.");
    else
        printf("%dx%d=%d pixels read from %s.", hsize,vsize,i,inputfilename);

    /* allcate memory for domain data and initialize it */
    compute_sums(hsize,vsize);

    if ((output = fopen(outputfilename, "w")) == NULL)
        fatal("Can't open output file.");

    /* output some data into the outputfile.                       */
    pack(4,(long)min_part,output);
    pack(4,(long)max_part,output);
    pack(4,(long)dom_step,output);
    pack(1,(long)dom_step_type,output);
    pack(2,(long)dom_type,output);
    pack(12,(long)hsize,output);
    pack(12,(long)vsize,output);

	/* This is the quantized value of zero scaling.. needed later */
    zero_ialpha = 0.5 + (max_scale)/(2.0*max_scale)*(1<<s_bits);

    /* The following routine takes a rectangular image and calls the */
    /* quadtree routine to encode square sum-images in it.           */
    /* the tolerance is a parameter since in some applications different */
    /* regions of the image may need to be compressed to different tol's */
    printf("\nEncoding Image.....");
    fflush(stdout);
    partition_image(0, 0, hsize,vsize, tol);
    printf("Done.");
    fflush(stdout);

    /* stuff the last byte if needed */
    pack(-1,(long)0,output);
    
    fclose(output);
    i = pack(-2,(long)0,output);
    printf("\n Compression = %lf from %d bytes written in %s.\n",
            (double)(hsize*vsize)/(double)i, i, outputfilename);

    /* Free allocated memory*/
    free(bits_needed);
    free(domimage[0]);
    free(domimage[1]);
    free(domimage[2]);
    free(domimage[3]);
    free(domain.no_h_domains);
    free(domain.no_v_domains);
    free(domain.domain_hsize);
    free(domain.domain_vsize);
    free(domain.domain_hstep);
    free(domain.domain_vstep);
    for (i=0; i <= max_part-min_part; ++i) 
        free(domain.pixel[i]);
    free(domain.pixel);
    free(image[0]);
    for (i=0; i <= max_part-min_part; ++i) 
    for (k=0; k<3; ++k) 
     for (j=0; j<24; ++j) list_free(the_domain[k][j][i]);
    return(0);
}

/* ************************************************************** */
/* free memory allocated in the list structure the_domain         */
/* ************************************************************** */
list_free(node)
struct classified_domain *node;
{
     if (node->next != NULL)
       list_free(node->next);
     free(node);
}

/* ************************************************************** */
/* return the average pixel value of a region of the image.       */
/* ************************************************************** */
void average(x,y,xsize,ysize, psum, psum2)
int x,y,xsize,ysize;
double *psum, *psum2;
{
    register int i,j,k;
    register double pixel;
    *psum = *psum2 = 0.0;
    k = ((x%2)<<1) + y%2;
    x >>= 1; y >>= 1;
    xsize >>= 1; ysize >>= 1;
    for (i=x; i<x+xsize; ++i)
    for (j=y; j<y+ysize; ++j) {
	   pixel = domimage[k][j][i];
       *psum += pixel;
       *psum2 += pixel*pixel;
    }
}

/* ************************************************************** */
/* return the average pixel value of a region of the image. This  */
/* routine differs from the previous in one slight way. It does   */
/* not average 2x2 sub-images to pixels. This is needed for clas- */
/* sifying ranges rather than domain where decimation is needed.  */
/* ************************************************************** */