/**************************************************************************/
/* 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.  */
/* ************************************************************** */
void average1(x,y,xsize,ysize, psum, psum2)
int x,y,xsize,ysize;
double *psum, *psum2;
{
    register int i,j;
    register double pixel;
    *psum = *psum2 = 0.0;

    for (i=x; i<x+xsize; ++i)
    for (j=y; j<y+ysize; ++j) {
	   pixel = (double)image[j][i];
       *psum += pixel;
       *psum2 += pixel*pixel;
    }
}

/* ************************************************************** */
/* Take a region of the image at x,y and classify it.             */
/* The four quadrants of the region are ordered from brightest to */
/* least bright average value, then it is rotated into one of the */
/* three cannonical orientations possible with the brightest quad */
/* in the upper left corner.                                      */
/* The routine returns two indices that are class numbers: pfirst */
/* and psecond; the symmetry operation that bring the square into */
/* cannonical position; and the sum and sum^2 of the pixel values */
/* ************************************************************** */
classify(x, y, xsize, ysize, pfirst, psecond, psym, psum, psum2, type)
int x,y,xsize,ysize,   /* position, size of subimage to be classified */
    *pfirst, *psecond, /* returned first and second class numbers     */
    *psym;             /* returned symmetry operation that brings the */
                       /* subimage to cannonical position.            */
double *psum, *psum2;  /* returned sum and sum^2 of pixel values      */
int type;              /* flag for decimating (for domains) or not    */ 
{

    int order[4], i,j; 
    double a[4],a2[4];
    void (*average_func)();
    
    if (type == 2) average_func = average; else average_func = average1;

    /* get the average values of each quadrant                         */


    (*average_func)(x,y,                 xsize/2,ysize/2,   &a[0], &a2[0]);
    (*average_func)(x,y+ysize/2,         xsize/2,ysize/2,   &a[1], &a2[1]);
    (*average_func)(x+xsize/2,y+ysize/2, xsize/2,ysize/2,   &a[2], &a2[2]);
    (*average_func)(x+xsize/2,y,         xsize/2,ysize/2,   &a[3], &a2[3]);

    *psum = a[0] + a[1] + a[2] + a[3];
    *psum2 =  a2[0] + a2[1] + a2[2] + a2[3];

    for (i=0; i<4; ++i) {
         /* after the sorting below order[i] is the i-th brightest   */
         /* quadrant.                                                */
         order[i] = i; 
         /* convert a2[] to store the variance of each quadrant      */
         a2[i] -= (double)(1<<(2*type))*a[i]*a[i]/(double)(xsize*ysize);
    }

    /* Now order the average value and also in order[],   which will */
    /* then tell us the indices (in a[]) of the brightest to darkest */
    for (i=2; i>=0; --i)
    for (j=0; j<=i; ++j)
        if (a[j]<a[j+1]) {
            swap(order[j], order[j+1],int)
            swap(a[j], a[j+1],double)
    }

    /* because of the way we ordered the a[] the rotation can be */
    /* read right off of order[]. That will make the brightest   */
    /* quadrant be in the upper left corner. But we must still   */ 
    /* decide which cannonical class the image portion belogs    */
    /* to and whether to do a flip or just a rotation. This is   */
    /* the following table summarizes the horrid lines below     */
    /* order      class            do a rotation                 */
    /* 0,2,1,3      0                   0                        */
    /* 0,2,3,1      0                   1                        */
    /* 0,1,2,3      1                   0                        */
    /* 0,3,2,1      1                   1                        */
    /* 0,1,3,2      2                   0                        */
    /* 0,3,1,2      2                   1                        */

    *psym = order[0]; 
    /* rotate the values */
    for (i=0; i<4; ++i)
        order[i] = (order[i] - (*psym) + 4)%4; 

    for (i=0; order[i] != 2; ++i); 
    *pfirst = i-1;
    if (order[3] == 1 || (*pfirst == 2 && order[2] == 1)) *psym += 4;

    /* Now further classify the sub-image by the variance of its    */
    /* quadrants. This give 24 subclasses for each of the 3 classes */
    for (i=0; i<4; ++i) order[i] = i; 

    for (i=2; i>=0; --i)
    for (j=0; j<=i; ++j)
        if (a2[j]<a2[j+1]) {
            swap(order[j], order[j+1],int)
            swap(a2[j], a2[j+1],double)
    }

    /* Now do the symmetry operation */
    for (i=0; i<4; ++i)
        order[i] = (order[i] - (*psym%4) + 4)%4; 
    if (*psym > 3) 
        for (i=0; i<4; ++i)
           if (order[i]%2) order[i] = (2 + order[i])%4;

    /* We want to return a class number from 0 to 23 depending on */
    /* the ordering of the quadrants according to their variance  */
    *psecond = 0;
    for (i=2; i>=0; --i)
    for (j=0; j<=i; ++j)
        if (order[j] > order[j+1]) {
            swap(order[j],order[j+1], int);
            if (order[j] == 0 || order [j+1] == 0) 
                 *psecond += 6;
            else if (order[j] == 1 || order [j+1] == 1) 
                 *psecond += 2;
            else if (order[j] == 2 || order [j+1] == 2) 
                 *psecond += 1;
        }
}

/* ************************************************************ */
/* Compute sum and sum^2 of pixel values in domains for use in  */
/* the rms computation later. Since a domain is compared with   */
/* many ranges, doing this just once saves a lot of computation */
/* This routine also fills a list structure with the domains    */
/* as they are classified and creates the memory for the domain */
/* data in a matrix.                                            */
/* ************************************************************ */
compute_sums(hsize,vsize)
int hsize,vsize;
{
    int i,j,k,l,
         domain_x, 
         domain_y, 
         first_class, 
         second_class,
         domain_size, 
         domain_step_size, 
         size,  
         x_exponent, 
         y_exponent;

    struct classified_domain *node;

    printf("\nComputing domain sums... ");
    fflush(stdout); 

    /* pre-decimate the image into domimage to avoid having to      */
    /* do repeated averaging of 2x2 pixel groups.                   */
    /* There are 4 ways to decimate the image, depending on the     */
    /* location of the domain, odd or even address.                 */
    for (i=0; i<2; ++i)
    for (j=0; j<2; ++j)
    for (k=i; k<hsize-i; k += 2)
    for (l=j; l<vsize-j; l += 2) 
        domimage[(i<<1)+j][l>>1][k>>1] = 
                     ((double)image[l][k] + (double)image[l+1][k+1] +
                     (double)image[l][k+1] + (double)image[l+1][k])*0.25;


    /* Allocate memory for the sum and sum^2 of domain pixels        */
    /* We first compute the size of the largest square that fits in  */
    /* the image.                                                    */
    x_exponent = (int)floor(log((double)hsize)/log(2.0));
    y_exponent = (int)floor(log((double)vsize)/log(2.0));
   
    /* exponent is min of x_ and y_ exponent */
    max_exponent = (x_exponent > y_exponent ? y_exponent : x_exponent);

    /* size is the size of the largest square that fits in the image */
    /* It is used to compute the domain and range sizes.             */
    size =  1<<max_exponent; 

    if (max_exponent < max_part)
      fatal("Reduce maximum number of quadtree partitions.");
    if (max_exponent-2 < max_part)
      printf("\nWarning: so many quadtree partitions yield absurd ranges.");

    i = max_part - min_part + 1;
    domain.no_h_domains = (int *)malloc(sizeof(int)*i);
    domain.no_v_domains = (int *)malloc(sizeof(int)*i);
    domain.domain_hsize = (int *)malloc(sizeof(int)*i);
    domain.domain_vsize = (int *)malloc(sizeof(int)*i);
    domain.domain_hstep = (int *)malloc(sizeof(int)*i);
    domain.domain_vstep = (int *)malloc(sizeof(int)*i);

    domain.pixel= (struct domain_pixels ***)
              malloc(i*sizeof(struct domain_pixels **));
    if (domain.pixel == NULL) fatal("No memory for domain pixel sums.");

    for (i=0; i <= max_part-min_part; ++i) {
       /* first compute how many domains there are horizontally */
       domain.domain_hsize[i] = size >> (min_part+i-1);
       if (dom_type == 2) 
             domain.domain_hstep[i] = dom_step; 
       else if (dom_type == 1)
             if (dom_step_type == 1) 
               domain.domain_hstep[i] = (size >> (max_part - i-1))*dom_step;
             else 
               domain.domain_hstep[i] = (size >> (max_part - i-1))/dom_step;
       else 
             if (dom_step_type == 1) 
               domain.domain_hstep[i] = domain.domain_hsize[i]*dom_step;
             else 
               domain.domain_hstep[i] = domain.domain_hsize[i]/dom_step;

       domain.no_h_domains[i] = 1+(hsize-domain.domain_hsize[i])/
                                                  domain.domain_hstep[i];
       /* bits_needed[i][0] = ceil(log(domain.no_h_domains[i])/log(2.0));  */

       /* now compute how many domains there are vertically. The sizes  */
       /* are the same for square domains, but not for rectangular ones */
       domain.domain_vsize[i] = size >> (min_part+i-1);
       if (dom_type == 2) 
             domain.domain_vstep[i] = dom_step; 
       else if (dom_type == 1)
             if (dom_step_type == 1) 
               domain.domain_vstep[i] = (size >> (max_part - i-1))*dom_step;
             else
               domain.domain_vstep[i] = (size >> (max_part - i-1))/dom_step;
       else 
             if (dom_step_type == 1) 
               domain.domain_vstep[i] = domain.domain_vsize[i]*dom_step;
             else
               domain.domain_vstep[i] = domain.domain_vsize[i]/dom_step;

       domain.no_v_domains[i] = 1+(vsize-domain.domain_vsize[i])/
                                                  domain.domain_vstep[i];

       /* Now compute the number of bits needed to store the domain data */
       bits_needed[i] = ceil(log((double)domain.no_h_domains[i]*
                                 (double)domain.no_v_domains[i])/log(2.0)); 

       matrix_allocate(domain.pixel[i], domain.no_h_domains[i],
                     domain.no_v_domains[i], struct domain_pixels) 

    }

    /* allocate and zero the list containing the classified domain data */
    i = max_part - min_part + 1;
    for (first_class = 0; first_class < 3; ++first_class)
    for (second_class = 0; second_class < 24; ++second_class) {
       the_domain[first_class][second_class] = 
                           (struct classified_domain **) 
                           malloc(i*sizeof(struct classified_domain *));
       for (j=0; j<i; ++j)
                    the_domain[first_class][second_class][j] = NULL;
    }

    /* Precompute sum and sum of squares for domains                 */
    /* This part can get faster for overlapping domains if repeated  */
    /* sums are avoided                                              */
    for (i=0; i <= max_part-min_part; ++i) {
      for (j=0,domain_x=0; j<domain.no_h_domains[i]; ++j,
             domain_x+=domain.domain_hstep[i]) 
      for (k=0,domain_y=0; k<domain.no_v_domains[i]; ++k,
             domain_y+=domain.domain_vstep[i]) {
		classify(domain_x, domain_y, 
                 domain.domain_hsize[i], 
                 domain.domain_vsize[i], 
					 &first_class, &second_class,
					 &domain.pixel[i][k][j].sym, 
                     &domain.pixel[i][k][j].sum, 
                     &domain.pixel[i][k][j].sum2, 2);

        /* When the domain data is referenced from the list, we need to */
        /* know where the domain is.. so we have to store the position  */
        domain.pixel[i][k][j].dom_x = j;
        domain.pixel[i][k][j].dom_y = k;
        node = (struct classified_domain *)
                                malloc(sizeof(struct classified_domain));
        
        /* put this domain in the classified list structure */
        node->the = &domain.pixel[i][k][j];
        node->next = the_domain[first_class][second_class][i];
        the_domain[first_class][second_class][i] = node;
	  }
    }

    /* Now we make sure no domain class is actually empty.  */
    for (i=0; i <= max_part-min_part; ++i) 
    for (first_class = 0; first_class < 3; ++first_class)
    for (second_class = 0; second_class < 24; ++second_class) 
       if (the_domain[first_class][second_class][i] == NULL) {
           node = (struct classified_domain *)
                                malloc(sizeof(struct classified_domain));
           node->the = &domain.pixel[i][0][0];
           node->next = NULL;
           the_domain[first_class][second_class][i] = node;
       }

    printf("Done.");
    fflush(stdout); 
}

/* ************************************************************ */
/* pack value using size bits and output into foutf */
/* ************************************************************ */
int pack(size, value, foutf)
int size; long int value;
FILE *foutf;
{
     int i;
     static int ptr = 1, /* how many bits are packed in sum so far */
                sum = 0, /* packed bits */ 
                num_of_packed_bytes = 0; /* total bytes written out */

    /* size == -1 means we are at the end, so write out what is left */
    if (size == -1 && ptr != 1) {
        fputc(sum<<(8-ptr), foutf);
        ++num_of_packed_bytes;
        return(0);
    }
    
    /* size == -2 means we want to know how many bytes we have written */
    if (size == -2) 
            return(num_of_packed_bytes);
    
    for (i=0; i<size; ++i, ++ptr, value = value>>1, sum = sum<<1) {
        if (value & 1) sum |= 1;

         if (ptr == 8) {
            fputc(sum,foutf);
            ++num_of_packed_bytes;
            sum=0;
            ptr=0;
        }
     }
}

/* ************************************************************ */
/* Compare a range to a domain and return rms and the quantized */
/* scaling and offset values (pialhpa, pibeta).                 */
/* ************************************************************ */
double compare(atx,aty, xsize, ysize, depth, rsum, rsum2, dom_x,dom_y, 
                                  sym_op, pialpha,pibeta)
int atx, aty, xsize, ysize, depth, dom_x, dom_y, sym_op, *pialpha, *pibeta;
double rsum, rsum2;
{
    int i, j, i1, j1, k,
        domain_x, 
        domain_y;        /* The domain position                */

    double pixel,
           det,          /* determinant of solution */
           dsum,         /* sum of domain values */
           rdsum = 0,    /* sum of range*domain values   */
           dsum2,        /* sum of domain^2 values   */
           w2 = 0,       /* total number of values tested */
           rms = 0,      /* root means square difference */
           alpha,        /* the scale factor */
           beta;         /* The offset */


	
    /* xsize = hsize >> depth; */
    /* ysize = vsize >> depth; */
    w2 = xsize * ysize;

    dsum = domain.pixel[depth-min_part][dom_y][dom_x].sum;
    dsum2 = domain.pixel[depth-min_part][dom_y][dom_x].sum2;
    domain_x = (dom_x * domain.domain_hstep[depth-min_part]);
    domain_y = (dom_y * domain.domain_vstep[depth-min_part]);
    k = ((domain_x%2)<<1) + domain_y%2;
    domain_x >>= 1;
    domain_y >>= 1;

    /* The main statement in this routine is a switch statement which */
    /* scans through the domain and range to compare them. The loop   */
    /* center is the same so we #define it for easy modification      */
#define COMPUTE_LOOP              {                                 \
        pixel = domimage[k][j1][i1];                                \
        rdsum += image[j][i]*pixel;                                 \
        }
    
     switch(sym_op) {
         case 0: for (i=atx, i1 = domain_x; i<atx+xsize; ++i, ++i1)
                 for (j=aty, j1 = domain_y; j<aty+ysize; ++j, ++j1) 
                     COMPUTE_LOOP
                 break;
         case 1: for (j=aty, i1 = domain_x; j<aty+ysize; ++j, ++i1)
                 for (i=atx+xsize-1, j1 = domain_y; i>=atx; --i, ++j1) 
                     COMPUTE_LOOP
                 break;
         case 2: for (i=atx+xsize-1, i1 = domain_x; i>=atx; --i, ++i1)
                 for (j=aty+ysize-1, j1 = domain_y; j>=aty; --j, ++j1) 
                     COMPUTE_LOOP
                 break;
         case 3: for (j=aty+ysize-1, i1 = domain_x; j>=aty; --j, ++i1)
                 for (i=atx, j1 = domain_y; i<atx+xsize; ++i, ++j1) 
                     COMPUTE_LOOP
                 break;
         case 4: for (j=aty, i1 = domain_x; j<aty+ysize; ++j, ++i1)
                 for (i=atx, j1 = domain_y; i<atx+xsize; ++i, ++j1) 
                     COMPUTE_LOOP
                 break;
         case 5: for (i=atx, i1 = domain_x; i<atx+xsize; ++i, ++i1)
                 for (j=aty+ysize-1, j1 = domain_y; j>=aty; --j, ++j1) 
                     COMPUTE_LOOP
                 break;
         case 6: for (j=aty+ysize-1, i1 = domain_x; j>=aty; --j, ++i1)
                 for (i=atx+xsize-1, j1 = domain_y; i>=atx; --i, ++j1) 
                     COMPUTE_LOOP
                 break;
         case 7: for (i=atx+xsize-1, i1 = domain_x; i>=atx; --i, ++i1)
                 for (j=aty, j1 = domain_y; j<aty+ysize; ++j, ++j1) 
                     COMPUTE_LOOP
                 break;
     }
    
     det = (xsize*ysize)*dsum2 - dsum*dsum;
    
     if (det == 0.0) 
         alpha = 0.0; 
      else 
         alpha = (w2*rdsum - rsum*dsum)/det;
        
     if  (only_positive && alpha < 0.0) alpha = 0.0;
     *pialpha = 0.5 + (alpha + max_scale)/(2.0*max_scale)*(1<<s_bits);
     if (*pialpha < 0) *pialpha = 0;
     if (*pialpha >= 1<<s_bits) *pialpha = (1<<s_bits)-1;

     /* Now recompute alpha back */
     alpha = (double)*pialpha/(double)(1<<s_bits)*(2.0*max_scale)-max_scale;
        
     /* compute the offset */
     beta = (rsum - alpha*dsum)/w2;

     /* Convert beta to an integer */
     /* we use the sign information of alpha to pack efficiently */
     if (alpha > 0.0)  beta += alpha*GREY_LEVELS;
     *pibeta = 0.5 + beta/
             ((1.0+fabs(alpha))*GREY_LEVELS)*((1<<o_bits)-1);
     if (*pibeta< 0) *pibeta = 0;
     if (*pibeta>= 1<<o_bits) *pibeta = (1<<o_bits)-1;

     /* Recompute beta from the integer */
     beta = (double)*pibeta/(double)((1<<o_bits)-1)*
               ((1.0+fabs(alpha))*GREY_LEVELS);
     if (alpha > 0.0) beta  -= alpha*GREY_LEVELS;

     /* Compute the rms based on the quantized alpha and beta! */
     rms= sqrt((rsum2 + alpha*(alpha*dsum2 - 2.0*rdsum + 2.0*beta*dsum) +
                 beta*(beta*w2 - 2.0*rsum))/w2);
    
     return(rms);
}

/* ************************************************************ */
/* Recursively partition an image, computing the best transfoms */
/* ************************************************************ */
quadtree(atx,aty,xsize,ysize,tol,depth)
int atx, aty, xsize, ysize, depth; 
double tol;  /* the tolerance fit  */
{
    int i,
        sym_op,                   /* A symmetry operation of the square */
        ialpha,                   /* Intgerized scalling factor         */
        ibeta,                    /* Intgerized offset                  */
        best_ialpha,              /* best ialpha found                  */
        best_ibeta,
        best_sym_op,
        best_domain_x,
        best_domain_y, 
        first_class, 
        the_first_class, 
        first_class_start,        /* loop beginning and ending values   */
        first_class_end, 
        second_class[2],
        the_second_class, 
        second_class_start,       /* loop beginning and ending values   */
        second_class_end, 
        range_sym_op[2],          /* the operations to bring square to  */
        domain_sym_op;            /* cannonical position.               */

    struct classified_domain *node;  /* var for domains we scan through */

    double rms, best_rms,          /* rms value and min rms found so far */
           sum=0, sum2=0;          /* sum and sum^2 of range pixels      */


    /* keep breaking it down until we are small enough */
    if (depth < min_part) {
         quadtree(atx,aty, xsize/2, ysize/2, tol,depth+1);
         quadtree(atx+xsize/2,aty, xsize/2, ysize/2,tol,depth+1);
         quadtree(atx,aty+ysize/2, xsize/2, ysize/2,tol,depth+1);
         quadtree(atx+xsize/2,aty+ysize/2, xsize/2, ysize/2,tol,depth+1);
         return;
    }

    /* now search for the best domain-range match and write it out */
    best_rms = 10000000000.0;      /* just a big number */

    classify(atx, aty, xsize,ysize, 
              &the_first_class, &the_second_class, 
              &range_sym_op[0], &sum, &sum2, 1);


    /* sort out how much to search based on -f and -F input flags */
    if (fullclass_search) {
         first_class_start = 0; 
         first_class_end = 3;
    } else {
         first_class_start = the_first_class; 
         first_class_end = the_first_class+1;
    } 

    if (subclass_search) {
         second_class_start = 0; 
         second_class_end = 24;
    } else {
         second_class_start = the_second_class; 
         second_class_end = the_second_class+1;
    } 

    /* these for loops vary according to the optimization flags we set */
    /* for subclass_search and fullclass_search==1, we search all the  */
    /* domains (except not all rotations).                             */
    for (first_class = first_class_start; 
         first_class < first_class_end; ++first_class)
    for (second_class[0] = second_class_start; 
         second_class[0] < second_class_end; ++second_class[0]) {

       /* We must check each domain twice. Once for positive scaling,  */
       /* once for negative scaling. Each has its own class and sym_op */
       if (!only_positive) {
          second_class[1] = 
	         class_transform[(first_class == 2 ? 1 : 0)][second_class[0]];
          range_sym_op[1] = 
             rot_transform[(the_first_class == 2 ? 1 : 0)][range_sym_op[0]];
       } 

       /* only_positive is 0 or 1, so we may or not scan                */
       for (i=0; i<(2-only_positive); ++i) 
       for (node = the_domain[first_class][second_class[i]][depth-min_part]; 
           node != NULL; 
           node = node->next) {
           domain_sym_op = node->the->sym;
           /* The following if statement figures out how to transform  */
           /* the domain onto the range, given that we know how to get */
           /* each into cannonical position.                           */
           if (((domain_sym_op>3 ? 4: 0) + (range_sym_op[i]>3 ? 4: 0))%8 == 0)
             sym_op = (4 + domain_sym_op%4 - range_sym_op[i]%4)%4;
           else 
             sym_op = (4 + (domain_sym_op%4 + 3*(4-range_sym_op[i]%4))%4)%8;
            
           rms = compare(atx,aty, xsize, ysize,  depth, sum,sum2, 
                                  node->the->dom_x, 
                                  node->the->dom_y, 
                                  sym_op, &ialpha,&ibeta); 
    
           if (rms < best_rms) {
                   best_ialpha = ialpha;
                   best_ibeta = ibeta;
                   best_rms = rms;
                   best_sym_op = sym_op;
                   best_domain_x = node->the->dom_x;
                   best_domain_y = node->the->dom_y;
           }
       }
    }
        
    if (best_rms > tol && depth < max_part) {
       /* We didn't find a good enough fit so quadtree down */
       pack(1,(long)1,output);  /* This bit means we quadtree'd down */
       quadtree(atx,aty, xsize/2, ysize/2, tol,depth+1);
       quadtree(atx+xsize/2,aty, xsize/2, ysize/2, tol,depth+1);
       quadtree(atx,aty+ysize/2, xsize/2, ysize/2, tol,depth+1);
       quadtree(atx+xsize/2,aty+ysize/2, xsize/2, ysize/2, tol,depth+1);
    } else {
       /* The fit was good enough or we just can't get smaller ranges */
       /* So write out the data */
       if (depth < max_part)          /* if we are not at the smallest range */
               pack(1,(long)0,output);/* then we must tell the decoder we    */
                                      /* stopped quadtreeing                 */
       pack(s_bits, (long)best_ialpha, output);
       pack(o_bits, (long)best_ibeta, output);
       /* When the scaling is zero, there is no need to store the domain */
       if (best_ialpha != zero_ialpha) {
          pack(3, (long)best_sym_op, output);
          pack(bits_needed[depth-min_part], (long)(best_domain_y*
            domain.no_h_domains[depth-min_part]+best_domain_x), output);
      }
     }
}

/* ************************************************************* */
/* Recursively partition an image, finding the largest contained */
/* square and call the quadtree routine the encode that square.  */
/* This enables us to encode rectangular image easily.           */
/* ************************************************************* */
partition_image(atx, aty, hsize,vsize, tol)
int atx, aty, hsize,vsize;
double tol;
{
   int x_exponent,    /* the largest power of 2 image size that fits */
       y_exponent,    /* horizontally or vertically the rectangle.   */
       exponent,      /* The actual size of image that's encoded.    */
       size, 
       depth; 
   
   x_exponent = (int)floor(log((double)hsize)/log(2.0));
   y_exponent = (int)floor(log((double)vsize)/log(2.0));
   
   /* exponent is min of x_ and y_ exponent */
   exponent = (x_exponent > y_exponent ? y_exponent : x_exponent);
   size = 1<<exponent; 
   depth = max_exponent - exponent;
   quadtree(atx,aty,size,size,tol,depth);
   if (size != hsize) 
      partition_image(atx+size, aty, hsize-size,vsize, tol);

   if (size != vsize) 
      partition_image(atx, aty+size, size,vsize-size, tol);
}        
  
/* fatal is for when there is a fatal error... print a message and exit */
void fatal(s)
char *s;
{
     printf("%s\n",s);
     exit(-1);
}