/**************************************************************************/
/* Decode an image encoded with a quadtree partition based fractal scheme */
/*                                                                        */
/*       Copyright 1993,1994 Yuval Fisher. All rights reserved.           */
/*                                                                        */
/* Version 0.05 10/10/94                                                  */
/**************************************************************************/

/* The following belong in a encdec.h file, but nevermind...              */
/* -----------------------------------------------------------------------*/
#include <stdio.h>
#include <math.h>

#define DEBUG 0
#define GREY_LEVELS 255

#define IMAGE_TYPE unsigned char /* may be different in some applications */

/* 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 bound(a)   ((a) < 0.0 ? 0 : ((a)>255.0? 255 : a))

/* various function declarations to keep compiler warnings away           */
void fatal();
char *malloc();
char *strcpy(); 

/* ---------------------------------------------------------------------- */

IMAGE_TYPE **image,*imptr,**image1; 
                             /* The input image data and a dummy         */ 

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       */
       hsize = -1,           /* The horizontal size of the input image    */
       vsize = -1,           /* The vertical size                         */
       scaledhsize,          /* hsize*scalefactor                         */
       scaledvsize,          /* vsize*scalefactor                         */
       size,                 /* largest square image that fits in image   */
       min_part = 3,         /* min and max _part determine a range of    */
       max_part = 4,         /* Range sizes from hsize>>min to hsize>>max */
       dom_step = 4,         /* 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..  */
       post_process = 1,     /* Flag for postprocessing.                  */
       num_iterations= 10,   /* Number of decoding iterations used.       */
       *no_h_domains,        /* Number of horizontal domains.             */
       *domain_hstep,        /* Domain density step size.                 */
       *domain_vstep,        /* Domain density step size.                 */
       *bits_needed,         /* Number of bits to encode domain position. */
       zero_ialpha,          /* the const ialpha when alpha = 0           */
       output_partition=0,   /* A flag for outputing the partition        */
       max_exponent;         /* The max power of 2 side of square image   */
                             /* that fits in our input image.             */

struct transformation_node {
       int rx,ry,            /* The range position and size in a trans.   */
              xsize, ysize,  
              rrx,rry,
              dx,dy;         /* The domain position.                      */
       int sym_op;           /* The symmetry operation used in the trans. */
       int depth;            /* The depth in the quadtree partition.      */
       double scale, offset; /* scalling and offset values.               */
       struct transformation_node *next;       /* The next trans. in list */
} transformations, *trans; 

FILE *input,*output,*other_input; 


/* fatal is for when there is a fatal error... print a message and exit */
void fatal(s)
char *s;
{
     printf("\n%s\n",s);
     exit(-1);
}

/* ************************************************************ */
/* unpack value using size bits read from fin.                  */
/* ************************************************************ */
long unpack(size, fin)
int size;
FILE *fin;
{
    int i;
    int value = 0;
    static int ptr = 1; /* how many bits are packed in sum so far */
    static int sum; 


    /* size == -2 means we initialize things */
    if (size == -2) {
        sum = fgetc(fin);
        sum <<= 1;
        return((long)0);
    }

    /* size == -1 means we want to peek at the next bit without */
    /* advancing the pointer */
    if (size == -1)
        return((long)((sum&256)>>8));

     for (i=0; i<size; ++i, ++ptr,  sum <<= 1) {
        if (sum & 256) value |= 1<<i;

         if (ptr == 8) {
            sum = getc(fin);
            ptr=0;
        }
     }
    return((long)value);
}

main(argc,argv)
int argc;
char **argv;
{
    /* Defaults are set initially */
    double         scalefactor = 1.0;           /* Scale factor for output */
    char           inputfilename[200];
    char           outputfilename[200];
    char           other_input_file[200];
    int            i,j, x_exponent, y_exponent;
    int            domain_size, no_domains;

 
    inputfilename[0] = 1;    /* We initially set the input to this and */
    outputfilename[0] = 1;   /* then check if the input/output names   */
    other_input_file[0] = 1; /* 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 'i': strcpy(other_input_file,argv[++i]); 
                         break;
               case 'n': num_iterations = atoi(argv[++i]);
                         break;
               case 'f': scalefactor = atof(argv[++i]);
                         break;
               case 'P': output_partition = 1;
                         break;
               case 'p': post_process = 0;
                         break;
               case 's': s_bits = atoi(argv[++i]);
                         break;
               case 'o': o_bits = atoi(argv[++i]);
                         break;
               case 'N': max_scale = atof(argv[++i]);
                         break;
               case '?':
               case 'H':
               default:
		   printf("\nUsage: dec -[options] [inputfile [outputfile]]");
		   printf("\nOptions are: (# = number), defaults show in ()");
		   printf("\n  -f # scale factor of output size. (%lf)", scalefactor);
		   printf("\n  -i file name. An initial image to iteration from.");
		   printf("\n  -n # no. of decoding iterations. (%d)", num_iterations);
		   printf("\n  -N # maximum allowed scaling. (%lf)",max_scale);
		   printf("\n  -s # number of scaling quantizing bits. (%d)",s_bits);
		   printf("\n  -o # number of offset quantizing bits. (%d)",o_bits);
		   printf("\n  -P   output the partition into the output file. (off)");
		   printf("\n  -p   supress artifact postprocessing. (off)");
           fatal(" ");
            }
       }

    if (inputfilename[0] == 1) strcpy(inputfilename, "lenna.trn");
    if (outputfilename[0] == 1) strcpy(outputfilename, "lenna.out");

    if ((input = fopen(inputfilename, "r")) == NULL)
        fatal("Can't open input file.");
        
    unpack(-2,input); /* initialize the unpacking routine */

	/* read the header data from the input file. This should probably */
    /* be put into one read which reads a structure with the info     */
    min_part = (int)unpack(4,input);
    max_part = (int)unpack(4,input);
    dom_step = (int)unpack(4,input);
    dom_step_type = (int)unpack(1,input);
    dom_type = (int)unpack(2,input);
    hsize = (int)unpack(12,input);
    vsize = (int)unpack(12,input);

    /* we now compute size */
    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; 

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

    /* allocate memory for the output image. Allocating one chunck saves  */
    /* work and time later.                                              */
	scaledhsize = (int)(scalefactor*hsize);
	scaledvsize = (int)(scalefactor*vsize);
    matrix_allocate(image, scaledhsize,scaledvsize, IMAGE_TYPE);
    matrix_allocate(image1, scaledhsize, scaledvsize, IMAGE_TYPE);

    if (other_input_file[0] != 1) {
        other_input = fopen(other_input_file, "r");
        i = fread(image[0], sizeof(IMAGE_TYPE), 
				scaledhsize*scaledvsize, other_input);
        if (i < scaledhsize*scaledvsize) 
           fatal("Couldn't read input... not enough data.");
        else
            printf("\n%d pixels read from %s.\n", i,other_input_file);
        fclose(other_input);
    }
          
    /* since max_ and min_ part are variable, these must be allocated */
    i = max_part - min_part + 1;
    bits_needed = (int *)malloc(sizeof(int)*i);
    no_h_domains = (int *)malloc(sizeof(int)*i);
    domain_hstep = (int *)malloc(sizeof(int)*i);
    domain_vstep = (int *)malloc(sizeof(int)*i);
    
    /* compute bits needed to read each domain type */
    for (i=0; i <= max_part-min_part; ++i) {
       /* first compute how many domains there are horizontally */
       domain_size = size >> (min_part+i-1);
       if (dom_type == 2) 
             domain_hstep[i] = dom_step; 
       else if (dom_type == 1)
             if (dom_step_type ==1)
                domain_hstep[i] = (size >> (max_part - i-1))*dom_step;
             else 
                domain_hstep[i] = (size >> (max_part - i-1))/dom_step;
       else 
             if (dom_step_type ==1)
                domain_hstep[i] = domain_size*dom_step;
             else 
                domain_hstep[i] = domain_size/dom_step;

       no_h_domains[i] = 1+(hsize-domain_size)/domain_hstep[i];
       /* bits_needed[i][0] = ceil(log(no_domains)/log(2));  */

       /* now compute how many domains there are vertically */
       if (dom_type == 2) 
             domain_vstep[i] = dom_step; 
       else if (dom_type == 1)
             if (dom_step_type ==1)
             domain_vstep[i] = (size >> (max_part - i-1))*dom_step;
             else 
             domain_vstep[i] = (size >> (max_part - i-1))/dom_step;
       else 
             if (dom_step_type ==1)
                domain_vstep[i] = domain_size*dom_step;
             else 
                domain_vstep[i] = domain_size/dom_step;

       no_domains = 1+(vsize-domain_size)/domain_vstep[i];
       bits_needed[i] = ceil(log((double)no_domains*(double)no_h_domains[i])/
                           log(2.0)); 
     }

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

    /* Read in the transformation data */
    trans = &transformations;
    printf("\nReading transformations.....");
    fflush(stdout);
    partition_image(0, 0, hsize,vsize );
    fclose(input);
    printf("Done.");
    fflush(stdout);
   
	if (scalefactor != 1.0) {
		printf("\nScaling image to %d x %d.", scaledhsize,scaledvsize);
		scale_transformations(scalefactor); 
	}
 
    /* when we output the partition, we just read the transformations */
    /* in and write them to the outputfile                            */ 
    if (output_partition) {
          fprintf(output,"\n%d %d\n %d %d\n%d %d\n\n", 
                  0, 0, scaledhsize, 0, scaledhsize, scaledvsize);
          printf("\nOutputed partition data in %s\n",outputfilename);
          fclose(output);
          return;
    }

    for (i=0; i<num_iterations; ++i) 
      apply_transformations();

    if (post_process) 	
      smooth_image();

    i = fwrite(image[0], sizeof(IMAGE_TYPE), scaledhsize*scaledvsize, output);
    if (i < scaledhsize*scaledvsize) 
        fatal("Couldn't write output... not enough disk space ?.");
    else
        printf("\n%d pixels written to output file.\n", i);

    fclose(output);
}

/* ************************************************************ */
/* Read in the transformation data from *input.                 */
/* This is a recursive routine whose recursion tree follows the */
/* recursion done by the encoding program.                      */
/* ************************************************************ */
read_transformations(atx,aty,xsize,ysize,depth)
int atx,aty,xsize,ysize,depth; 
{
    /* Having all these locals in a recursive procedure is hard on the  */
    /* stack.. but it is more readable.                                 */
    int i,j,
        sym_op,                   /* A symmetry operation of the square */
        ialpha,                   /* Intgerized scalling factor         */
        ibeta;                    /* Intgerized offset                  */

    long domain_ref; 

    double alpha, beta;

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

    if (depth < max_part && unpack(1,input)) {
         /* A 1 means we subdivided.. so quadtree */
         read_transformations(atx,aty, xsize/2, ysize/2, depth+1);
         read_transformations(atx+xsize/2,aty, xsize/2, ysize/2, depth+1);
         read_transformations(atx,aty+ysize/2, xsize/2, ysize/2, depth+1);
         read_transformations(atx+xsize/2,aty+ysize/2,xsize/2,ysize/2,depth+1);
    } else {  
         /* we have a transformation to read */
         trans->next = (struct transformation_node *)
                  malloc(sizeof(struct transformation_node ));
         trans = trans->next; trans->next = NULL;
         ialpha = (int)unpack(s_bits,  input);
         ibeta = (int)unpack(o_bits,  input);
         alpha = (double)ialpha/(double)(1<<s_bits)*(2.0*max_scale)-max_scale;

         beta = (double)ibeta/(double)((1<<o_bits)-1)*
               ((1.0+fabs(alpha))*GREY_LEVELS);
         if (alpha > 0.0) beta -= alpha*GREY_LEVELS;

         trans->scale = alpha;
         trans->offset = beta;
         if (ialpha != zero_ialpha) {
            trans-> sym_op = (int)unpack(3, input);
            domain_ref = unpack(bits_needed[depth-min_part], input);
            trans->dx = (double)(domain_ref % no_h_domains[depth-min_part])
                                              * domain_hstep[depth-min_part];
            trans->dy = (double)(domain_ref / no_h_domains[depth-min_part])
                                              * domain_vstep[depth-min_part];
         } else {
            trans-> sym_op = 0;
            trans-> dx  = 0;
            trans-> dy = 0;
         }
         trans->rx = atx;
         trans->ry = aty;
         trans->depth = depth;
        
         trans->rrx = atx + xsize;
         trans->rry = aty + ysize;

         if (output_partition) 
                fprintf(output,"\n%d %d\n %d %d\n%d %d\n\n", 
                  atx,       vsize-aty-ysize, 
                  atx,       vsize-aty, 
                  atx+xsize, vsize-aty);
  
    }
}
     
/* ************************************************************ */
/* Apply the transformations once to an initially black image.  */
/* ************************************************************ */
apply_transformations()
{
    IMAGE_TYPE **tempimage;
	int	i,j,i1,j1,count=0;
    double pixel;

    trans = &transformations;
    while (trans->next != NULL) {
        trans = trans->next;
		++count;

/* Since the inner loop is the same in each case of the switch below */
/* we just define it once for easy modification.                     */
#define COMPUTE_LOOP              {                                  \
        pixel = (image[j1][i1]+image[j1][i1+1]+image[j1+1][i1]+      \
                image[j1+1][i1+1])/4.0;                              \
        image1[j][i] = bound(0.5 + trans->scale*pixel+trans->offset);\
        }

        switch(trans->sym_op) {
           case 0: for (i=trans->rx, i1 = trans->dx; 
                         i<trans->rrx; ++i, i1 += 2)
                   for (j=trans->ry, j1 = trans->dy; 
                         j<trans->rry; ++j, j1 += 2) 
                       COMPUTE_LOOP
                   break;
           case 1: for (j=trans->ry, i1 = trans->dx;
                           j<trans->rry; ++j, i1 += 2)
                   for (i=trans->rrx-1,
                           j1 = trans->dy; i>=(int)trans->rx; --i, j1 += 2) 
                       COMPUTE_LOOP
                   break;
           case 2: for (i=trans->rrx-1, 
                          i1 = trans->dx; i>=(int)trans->rx; --i, i1 += 2)
                   for (j=trans->rry-1,
                           j1 = trans->dy; j>=(int)trans->ry; --j, j1 += 2) 
                       COMPUTE_LOOP
                   break;
           case 3: for (j=trans->rry-1,
                           i1 = trans->dx; j>=(int)trans->ry; --j, i1 += 2)
                   for (i=trans->rx, j1 = trans->dy;
                           i<trans->rrx; ++i, j1 += 2) 
                       COMPUTE_LOOP
                   break;
           case 4: for (j=trans->ry, i1 = trans->dx;
                           j<trans->rry; ++j, i1 += 2)
                   for (i=trans->rx, j1 = trans->dy;
                           i<trans->rrx; ++i, j1 += 2) 
                       COMPUTE_LOOP
                   break;
           case 5: for (i=trans->rx, i1 = trans->dx;
                           i<trans->rrx; ++i, i1 += 2)
                   for (j=trans->rry-1,
                           j1 = trans->dy; j>=(int)trans->ry; --j, j1 += 2) 
                       COMPUTE_LOOP
                   break;
           case 6: for (j=trans->rry-1,
                           i1 = trans->dx; j>=(int)trans->ry; --j, i1 += 2)
                   for (i=trans->rrx-1,
                           j1 = trans->dy; i>=(int)trans->rx; --i, j1 += 2) 
                       COMPUTE_LOOP
                   break;
           case 7: for (i=trans->rrx-1, 
                          i1 = trans->dx; i>=(int)trans->rx; --i, i1 += 2)
                   for (j=trans->ry, j1 = trans->dy; 
                          j<trans->rry; ++j, j1 += 2) 
                       COMPUTE_LOOP
                   break;
        }  
  
    }
    tempimage = image;
    image = image1;
    image1 = tempimage;

    printf("\n%d transformations applied.",count);
}

/*        This should really be done when they are read in.     */
/* ************************************************************ */
scale_transformations(scalefactor)
double scalefactor;
{
    trans = &transformations;
    while (trans->next != NULL) {
        trans = trans->next;

		trans->rrx *= scalefactor;
		trans->rry *= scalefactor;
		trans->rx *= scalefactor;
		trans->ry *= scalefactor;
		trans->dx *= scalefactor;
		trans->dy *= scalefactor;
    }
}

/* ************************************************************* */
/* Recursively partition an image, finding the largest contained */
/* square and call read_transformations .                        */
/* ************************************************************* */
partition_image(atx, aty, hsize,vsize )
int atx, aty, hsize,vsize;
{
   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;

   read_transformations(atx,aty,size,size,depth);

   if (size != hsize) 
      partition_image(atx+size, aty, hsize-size,vsize );

   if (size != vsize) 
      partition_image(atx, aty+size, size,vsize-size );
}        
  
/* ************************************************************* */
/* Scan the image and average the transformation boundaries.     */
/* ************************************************************* */
smooth_image()
{
    IMAGE_TYPE pixel1, pixel2;
	int i,j;
	int w1,w2;

	printf("\nPostprocessing Image.");
    trans = &transformations;
    while (trans->next != NULL) {
        trans = trans->next;
		if (trans->rx == 0 || trans->ry == 0) 
            continue;
		
		if (trans->depth == max_part) {
			w1 = 5;
			w2 = 1;
		} else {
			w1 = 2;
			w2 = 1;
		}

        for (i=trans->rx; i<trans->rrx; ++i) {
             pixel1 = image[(int)trans->ry][i];
             pixel2 = image[(int)trans->ry-1][i];
             image[(int)trans->ry][i] = (w1*pixel1 + w2*pixel2)/(w1+w2);
             image[(int)trans->ry-1][i] = (w2*pixel1 + w1*pixel2)/(w1+w2);
        }

        for (j=trans->ry; j<trans->rry; ++j) {
             pixel1 = image[j][(int)trans->rx];
             pixel2 = image[j][(int)trans->rx-1];
             image[j][(int)trans->rx] = (w1*pixel1 + w2*pixel2)/(w1+w2);
             image[j][(int)trans->rx-1] = (w2*pixel1 + w1*pixel2)/(w1+w2);
        }
    }
}