fractalencode.c

分形图像压缩源代码

/*   FractalEncode.c is a modified version of enc.c by Yuval Fisher  */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#include <png.h>

#include "FractalSqueeze.h"

static IMAGE_TYPE **image;          /* The input image data                      */ 
static double     **domimage[4];    /* Decimated input image used for domains    */

static double max_scale = 1.0;      /* Maximum allowable grey level scale factor */

static int    s_bits = 5,           /* Number of bits used to store scale factor */
       o_bits = 7,           /* Number of bits used to store offset       */
       min_part = 4,         /* Min and max _part determine a range of    */
       max_part = 6,         /* Range sizes from hsize>>min to hsize>>max */
       dom_step = 1,         /* 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..  */
       only_positive = 0,    /* A flag specifying use of positive scaling */
       subclass_search = 0,  /* A flag specifying classes searched        */
       fullclass_search = 0, /* A flag specifying classes searched        */
       *bits_needed,         /* Number of bits to encode domain position. */
       zero_ialpha,          /* The const ialpha when alpha = 0           */
       max_exponent;         /* The max power of 2 side of square image   */
                             /* that fits in our input image.             */


                             /* The class_transform gives the transforms  */
                             /* between classification numbers for        */
                             /* negative scaling values, when brightest   */
                             /* becomes darkest, etc...                   */
static int    class_transform[2][24] = {23,17,21,11,15,9,22,16,19,5,13,3,20,10,18,
                                 4,7,1,14,8,12,2,6,0,
                                 16,22,10,20,8,14,17,23,4,18,2,12,11,21,5,
                                 19,0,6,9,15,3,13,1,7};

                             /* rot_transform gives the rotations for     */
                             /* domains with negative scalings.           */
static int    rot_transform[2][8] = {7,4,5,6,1,2,3,0, 2,3,0,1,6,7,4,5};

struct domain_data {
       int *no_h_domains,    /* The number of domains horizontally for    */
           *no_v_domains,    /* each size.                                */
           *domain_hsize,    /* The size of the domain.                   */
           *domain_vsize,    /* The size of the domain.                   */
           *domain_hstep,    /* The density of the domains.               */
           *domain_vstep;    /* The density of the domains.               */
struct domain_pixels {       /* This is a three (sigh) index array that   */
       int dom_x, dom_y;     /* dynamically allocated. The first index is */
       double sum,sum2;      /* the domain size, the other are two its    */
       int sym;              /* position. It contains the sum and sum^2   */
} ***pixel;                  /* of the pixel values in the domains, which */
} domain;                    /* are computed just once.                   */ 
                             

struct classified_domain {             /* This is a list which containes  */
       struct domain_pixels *the;      /* pointers to  the domain data    */
       struct classified_domain *next; /* in the structure above. There   */
} **the_domain[3][24];                 /* are three classes with 24 sub-  */
                                       /* classes. Using this array, only */
                                       /* domains and ranges in the same  */
                                       /* class are compared..            */
                                       /* The first pointer points to the */
                                       /* domain size the the second to   */
                                       /* list of domains.                */

static FILE *output;                /* Output FILE containing compressed data    */

// Prototypes:
int pack(int size, long value, FILE *foutf);
void average(int x, int y, int xsize, int ysize, double *psum, double *psum2);
void compute_sums(int hsize, int vsize);
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);
void quadtree(int atx, int aty, int xsize, int ysize, double tol, int depth);
void partition_image(int atx, int aty, int hsize, int vsize, double tol);
void list_free(struct classified_domain *node);

void fractal_encode(png_infop image_data, char* outputfilename)
{
    /* Defaults are set initially */
    double         tol = 8.0;            /* Tolerance value for quadtree.  */
    int            i,j,k,
                   hsize = image_data->height,
                   vsize = image_data->width;
 
    /* allocate memory for the input image. Allocating one chunck saves  */
    /* work and time later.                                              */
    matrix_allocate(image, hsize, vsize, IMAGE_TYPE)
    matrix_allocate(domimage[0], hsize/2, vsize/2, double)
    matrix_allocate(domimage[1], hsize/2, vsize/2, double)
    matrix_allocate(domimage[2], hsize/2, vsize/2, double)
    matrix_allocate(domimage[3], hsize/2, vsize/2, double)

    /* max_ & min_ part are variable, so this must be run time allocated */
    bits_needed = (int *)malloc(sizeof(int)*(1+max_part-min_part));

	{
		unsigned int i,j;

		for (i = 0; i < image_data->width; ++i) {
			for (j = 0; j < image_data->height; ++j) {
				image[i][j] = image_data->row_pointers[j][i];
			}
		}
	}

    /* allcate memory for domain data and initialize it */
    compute_sums(hsize,vsize);

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

    /* output some data into the outputfile.                       */
    pack(4,(long)min_part,output);
    pack(4,(long)max_part,output);
    pack(4,(long)dom_step,output);
    pack(1,(long)dom_step_type,output);
    pack(2,(long)dom_type,output);
    pack(12,(long)hsize,output);
    pack(12,(long)vsize,output);

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

    /* The following routine takes a rectangular image and calls the */
    /* quadtree routine to encode square sum-images in it.           */
    /* the tolerance is a parameter since in some applications different */
    /* regions of the image may need to be compressed to different tol's */
    printf("\nEncoding Image.....");
    fflush(stdout);
    partition_image(0, 0, hsize,vsize, tol);
    printf("Done.");
    fflush(stdout);

    /* stuff the last byte if needed */
    pack(-1,(long)0,output);
    
    fclose(output);
    i = pack(-2,(long)0,output);
    printf("\n Compression = %lf from %d bytes written in %s.\n",
            (double)(hsize*vsize)/(double)i, i, outputfilename);

    /* Free allocated memory*/
    free(bits_needed);
    free(domimage[0]);
    free(domimage[1]);
    free(domimage[2]);
    free(domimage[3]);
    free(domain.no_h_domains);
    free(domain.no_v_domains);
    free(domain.domain_hsize);
    free(domain.domain_vsize);
    free(domain.domain_hstep);
    free(domain.domain_vstep);
    for (i=0; i <= max_part-min_part; ++i) 
        free(domain.pixel[i]);
    free(domain.pixel);
    free(image[0]);
    for (i=0; i <= max_part-min_part; ++i) 
    for (k=0; k<3; ++k) 
     for (j=0; j<24; ++j) list_free(the_domain[k][j][i]);
    
	return;
}

/* ************************************************************** */
/* free memory allocated in the list structure the_domain         */
/* ************************************************************** */
void list_free(struct classified_domain *node)
{
     if (node->next != NULL)
       list_free(node->next);
     free(node);
}

/* ************************************************************** */
/* return the average pixel value of a region of the image.       */
/* ************************************************************** */
void average(int x, int y, int xsize, int ysize, double *psum, double *psum2)
{
    register int i,j,k;
    register double pixel;
    *psum = *psum2 = 0.0;
    k = ((x%2)<<1) + y%2;
    x >>= 1; y >>= 1;
    xsize >>= 1; ysize >>= 1;
    for (i=x; i<x+xsize; ++i)
    for (j=y; j<y+ysize; ++j) {
	   pixel = domimage[k][j][i];
       *psum += pixel;
       *psum2 += pixel*pixel;
    }
}

/* ************************************************************** */
/* return the average pixel value of a region of the image. This  */
/* routine differs from the previous in one slight way. It does   */
/* not average 2x2 sub-images to pixels. This is needed for clas- */
/* sifying ranges rather than domain where decimation is needed.  */
/* ************************************************************** */
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)
    }