fractalencode.c

分形图像压缩源代码

    /* 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.                                            */
/* ************************************************************ */
void compute_sums(int hsize, int vsize)
{
    int i,j,k,l,
         domain_x, 
         domain_y, 
         first_class, 
         second_class,
         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)
      ERROROUT("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) ERROROUT("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] = (int)(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(int size, long 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;
        }
     }
	return 0;
}

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