enc.c

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

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