fractaldecode.c

分形图像压缩源代码

    /* 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 = (int)(double)(domain_ref % no_h_domains[depth-min_part])
                                              * domain_hstep[depth-min_part];
            trans->dy = (int)(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.  */
/* ************************************************************ */
void apply_transformations(void)
{
    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 = (int)((image[j1][i1]+image[j1][i1+1]+image[j1+1][i1]+\
                image[j1+1][i1+1])/4.0);                             \
        image1[j][i] = (int)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.     */
/* ************************************************************ */
void scale_transformations(double scalefactor)
{
    trans = &transformations;
    while (trans->next != NULL) {
        trans = trans->next;

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

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

   if (size != vsize) 
      partition_image_read(atx, aty+size, size,vsize-size );
}        
  
/* ************************************************************* */
/* Scan the image and average the transformation boundaries.     */
/* ************************************************************* */
void smooth_image(void)
{
    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);
        }
    }
}