enc.c

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

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);