cgen.c

这是一个高速多维插值算法。当我们建模以后

/* Integer Multi-Dimensional Interpolation */

/*
 * Copyright 2000 - 2002 Graeme W. Gill
 * All rights reserved.
 *
 * This material is licenced under the GNU GENERAL PUBLIC LICENCE :-
 * see the Licence.txt file for licencing details.
 */

/* 'C' code color transform kernel code generator. */

/*
   This module generates C code routines which implement
   an integer multi-channel transform. The input values
   are read, passed through per channel lookup tables,
   a multi-dimentional interpolation table, and then
   a per channel output lookup table, before being written.
*/

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

#include "imdi_imp.h"
#include "imdi_gen.h"
#include "imdi_tab.h"

#undef VERBOSE
#undef FORCESORT		/* Use sort algorithm allways */

/*
 * TTBD:
 *		Need to implement g->dir
 *      Haven't used t->it_map[] or t->im_map[].
 *
 *
 */

/* ------------------------------------ */
/* Context */
typedef struct {
	FILE *of;			/* Output file */
	int indt;			/* Indent */

	/* Other info */
	genspec *g;			/* Generation specifications */
	tabspec *t;			/* Table setup data */
	mach_arch *a;		/* Machine architecture and tuning data */

	/* Code generation information */
	/* if() conditions are for entry usage */

	/* Pixel read information */
	int ipt[IXDI];		/* Input pointer types */
	int nip;			/* Actual number of input pointers, accounting for pint */
	int chv_bits;		/* Bits in chv temp variable ?? */

	/* Input table entry */
	int itet;			/* Input table entry type */
	int itvt;			/* Input table variable type */
	int itmnb;			/* Input table minimum bits (actual is it_ab) */

	/* Interpolation index */
	int ixet;			/* Interpolation index entry type */
	int ixvt;			/* Interpolation index variable type */
	int ixmnb;			/* Interpolation index minimum bits (actual is ix_ab) */
	int ixmxres;		/* Interpolation table maximum resolution */

	/* Simplex index: if(!sort && it_xs) */
	int sxet;			/* Simplex index entry type  */
	int sxvt;			/* Simplex index variable type */
	int sxmnb;			/* Simplex index bits minimum (actual is sx_ab) */
	int sxmxres;		/* Simplex table maximum resolution (0 if sort) */

	/* Combination Weighting + Vertex offset values: if(it_xs && !wo_xs) */
	int woet;			/* Weighting+offset entry type  */
	int wovt;			/* Weighting+offset variable type */
	int womnb;			/* Weighting+offset index bits minimum (actual is wo_ab) */

	/* Weighting value: if(it_xs && wo_xs) */
	int weet;			/* Weighting entry type  */
	int wevt;			/* Weighting variable type */
	int wemnb;			/* Weighting index bits minimum (actual is we_ab) */

	/* Vertex offset value: if(it_xs && wo_xs) */
	int voet;			/* Vertex offset entry type  */
	int vovt;			/* Vertex offset variable type */
	int vomnb;			/* Vertex offset index bits minimum (actual is vo_ab) */

	/* Interpolation table entry: */
	int imovb;			/* Interpolation output value bits per channel required */
	int imfvt;			/* Interpolation full entry & variable type */
	int impvt;			/* Interpolation partial entry variable type */

	/* Interpolation accumulators: */
	int iaovb;			/* Interpolation output value bits per channel required */
	int iafvt;			/* Interpolation full entry & variable type */
	int iapvt;			/* Interpolation partial entry variable type */
	int ian;			/* Total number of accumulators */

	/* Output table lookup */
	int otit;			/* Output table index type */
	int otvt;			/* Output table value type (size is ot_ts bytes) */

	/* Write information */
	int opt[IXDO];		/* Output pointer types */
	int nop;			/* Actual number of output pointers, accounting for pint */

} fileo;

void line(fileo *f, char *fmt, ...);	/* Output one line */
void sline(fileo *f, char *fmt, ...);	/* Output start of line line */
void mline(fileo *f, char *fmt, ...);	/* Output middle of line */
void eline(fileo *f, char *fmt, ...);	/* Output end of line */
void cr(fileo *f) { line(f,""); }		/* Output a blank line */
void inc(fileo *f) { f->indt++; }		/* Increment the indent level */
void dec(fileo *f) { f->indt--; }		/* Decrement the indent level */
/* ------------------------------------ */

int findord(fileo *f, int bits);		/* Find ordinal with bits or more */
int nord(fileo *f, int ov);				/* Round ordinal type up to natural size */
int findnord(fileo *f, int bits);		/* Find ordinal with bits, or natural larger */
int findint(fileo *f, int bits);		/* Find integer with bits or more */
int nint(fileo *f, int iv);				/* Round integer type up to natural size */
int findnint(fileo *f, int bits);		/* Find integer with bits, or natural larger */
static void doheader(fileo *f);

static int calc_bits(int dim, int res);
static int calc_res(int dim, int bits);
static int calc_obits(int dim, int res, int esize);
static int calc_ores(int dim, int bits, int esize);


/* return a hexadecimal mask string */
/* take care of the case when bits >= 32 */
char *hmask(int bits) {
	static char buf[20];

	if (bits < 32) {
		sprintf(buf, "0x%x",(1 << bits)-1);
	} else if (bits == 32) {
		return "0xffffffff";
	} else if (bits == 64) {
		return "0xffffffffffffffff";
	} else {	/* Bits > 32 */
		sprintf(buf, "0x%xffffffff",(1 << (bits-32))-1);
	}
	return buf;
}

/* Generate a source file to implement the specified */
/* interpolation kernel. Fill in return values and return 0 if OK. */
/* Return non-zero on error. */
int gen_c_kernel(
	genspec *g,				/* Specification of what to generate */
	mach_arch *a,
	FILE *fp,				/* File to write to */
	int index				/* Identification index, 1 = first */
) {
	unsigned char kk[] = { 0x43, 0x6F, 0x70, 0x79, 0x72, 0x69, 0x67, 0x68,
	                       0x74, 0x20, 0x32, 0x30, 0x30, 0x34, 0x20, 0x47,
	                       0x72, 0x61, 0x65, 0x6D, 0x65, 0x20, 0x57, 0x2E,
	                       0x20, 0x47, 0x69, 0x6C, 0x6C, 0x00 };
	fileo f[1];
	int e, i;
	tabspec tabsp, *t = &tabsp;
	int timp = 0;		/* Flag to use temporary imp pointer. */
						/* Seem to make x86 MSVC++ slower */
						/* Has no effect on x86 IBMCC */

	sprintf(g->kname, "imdi_k%d",index); /* Kernel routine base name */
	strcpy(g->kkeys, (char*)kk);				 /* Kernel keys for this session */
	
	/* Setup the file output context */
	f->of = fp;
	f->indt = 0;			/* Start with no indentation */
	f->g = g;
	f->t = t;
	f->a = a;

	if (g->prec == 8) {
		if (g->id <= 4)
			t->sort = 0;		/* Implicit sort using simplex table lookup */
		else
			t->sort = 1;		/* Explicit sort */

	} else  if (g->prec == 16) {
		t->sort = 1;			/* Explit sort, no simplex table */

	} else {
		fprintf(stderr,"Can't cope with requested precision of %d bits\n",g->prec);
		exit(-1);
	}
#ifdef FORCESORT
	t->sort = 1;
#endif

	/* Compute input read and input table lookup stuff */

	/* Compute number of input pointers */
	if (g->in.pint != 0)	/* Pixel interleaved */
		f->nip = 1;
	else
		f->nip = g->id;

	/* Figure out the input pointer types */
	for (e = 0; e < f->nip; e++) {
		if ((f->ipt[e] = findord(f, g->in.bpch[e])) < 0) {
			fprintf(stderr,"Input channel size can't be handled\n");
			exit(-1);
		}
	}

	/* Set a default input channel mapping */
	for (e = 0; e < g->id; e++)
		t->it_map[e] = e;

	/* Do the rest of the input table size calculations after figuring */
	/* out simplex and interpolation table sizes. */


	/* Figure out the interpolation multi-dimentional table structure */
	/* and output accumulation variable sizes. */

	if (g->prec == 8
	 || g->prec == 16 && a->ords[a->nords-1].bits >= (g->prec * 4)) {
		int tiby;		/* Total interpolation bytes needed */

		/* We assume that we can normally compute more than one */
		/* output value at a time, so we need to hold the interpolation */
		/* output data in the expanded fixed point format in both the */
		/* table and accumulator. */
		t->im_cd = 1;
		f->imovb = g->prec * 2;		/* 16 bits needed for 8 bit precision, */
		f->iaovb = g->prec * 2;		/* 32 bits needed for 16 bit precision */
		f->imfvt = a->nords-1;		/* Full variable entry type is biggest available */
		f->iafvt = a->nords-1;		/* Full variable accum. type is same */

		if (a->ords[f->imfvt].bits < f->imovb) {
			fprintf(stderr,"Interpolation table entry size can't be handled\n");
			exit(-1);
		}

		/* Compute details of table entry sizes, number */
		tiby = (f->imovb * g->od)/8;				/* Total table bytes needed */
		t->im_fs = a->ords[f->imfvt].bits/8;		/* Full entry bytes */
		t->im_fv = (t->im_fs * 8)/f->imovb;			/* output values per full entry . */
		t->im_fn = tiby/t->im_fs;					/* Number of full entries (may be 0) */
		t->im_ts = t->im_fn * t->im_fs;				/* Structure size so far */
		tiby -= t->im_fn * t->im_fs;				/* Remaining bytes */

		if (tiby <= 0) {
			t->im_pn = 0;		/* No partials */
			t->im_ps = 0;
			t->im_pv = 0;
			f->impvt = 0;
			f->iapvt = 0;

		} else {
			t->im_pn = 1;					/* Must be just 1 partial */
			t->im_pv = (tiby * 8)/f->imovb;	/* Partial holds remaining entries */ 

			if ((f->impvt = findnord(f, tiby * 8)) < 0) {
				fprintf(stderr,"Can't find partial interp table entry variable size\n");
				exit(-1);
			}
			f->iapvt = f->impvt;
			t->im_ps = a->ords[f->impvt].bits/8;/* Partial entry bytes */

			if (a->ords[f->imfvt].align)		/* If full entry's need to be aligned */
				t->im_ts += t->im_fs;			/* Round out struct size by full entry */
			else
				t->im_ts += t->im_ps;			/* Round out to natural size */
		}
	
	} else {
		/* One 16 bit output value per entry + 32 bit accumulator. */
		/* We can conserve table space by not holding the table data in expanded */
		/* fixed point format, but expanding it when it is read. */
		/* Without resorting to compicated code, this restricts us */
		/* to only computing one output value per accumulator. */
		t->im_cd = 0;
		f->imovb = g->prec;			/* Table holds 16 bit entries with no fractions */
		f->iaovb = g->prec * 2;		/* 32 bits needed for 16 bit precision in comp. */

		if ((f->imfvt = findord(f, f->imovb)) < 0) {
			fprintf(stderr,"Interpolation table entry size can't be handled\n");
			exit(-1);
		}
		if ((f->iafvt = findord(f, f->iaovb)) < 0) {
			fprintf(stderr,"Interpolation accumulator size can't be handled\n");
			exit(-1);
		}

		/* Compute details of table entry sizes, number */
		t->im_fs = a->ords[f->imfvt].bits/8;		/* Full entry bytes */
		t->im_fv = 1;								/* output values per full entry . */
		t->im_fn = g->od;							/* Number of full entries */
		t->im_ts = t->im_fn * t->im_fs;				/* Total structure size */

		t->im_pn = 0;		/* No partials */
		t->im_ps = 0;
		t->im_pv = 0;
		f->impvt = 0;
		f->iapvt = 0;
	}
	f->ian = t->im_fn + t->im_pn;			/* Total number of output accumulators */

	/* Figure out how much of the interpolation entry offset to put in the */
	/* vertex offset value, and how much to make explicit in accessing the */
	/* interpolation table enty. */
	if (a->oscale > 0) {		/* We have a scaled index mode */
		/* Use as much of the scaled index mode as possible */
		/* and then do the balance by scaling the simplex index entry. */
		for (t->im_oc = a->oscale; ; t->im_oc >>= 1) {
			t->vo_om = t->im_ts/t->im_oc;		/* Simplex index multiplier */
			if ((t->vo_om * t->im_oc) == t->im_ts)
				break;				/* Got appropriate offset scale */
		}
	} else if (a->smmul) {		/* Architecure supports fast small multiply */
		t->im_oc = t->im_ts;	/* Do scale by structure size explicitly */
		t->vo_om = 1;			/* Do none in the Simplex index */
	} else {					/* We have no fast tricks */
		t->im_oc = 1;			/* Do none explicitly */
		t->vo_om = t->im_ts;	/* Do all in Simplex index */
	}

	/* Compute the number of bits needed to hold an index into */
	/* the interpolation table (index is in terms of table entry size). */
	/* This value is used to figure out the room needed in the input */
	/* table to accumulate the interpolation cube base offset value. (IM_O macro) */
	f->ixmnb = calc_bits(g->id, g->itres);

	/* Set a default output channel mapping */
	for (e = 0; e < g->od; e++)
		t->im_map[e] = e;

#ifdef VERBOSE
	/* Summarise the interpolation table arrangements */
	printf("\n");
	printf("Interpolation table structure:\n");
	printf("  Minimum bits needed to index table %d\n", f->ixmnb);
	printf("  Entry total size %d bytes\n", t->im_ts);
	printf("  Simplex entry offset scale %d\n", t->vo_om);
	printf("  Explicit entry offset scale %d\n", t->im_oc);
	printf("  %d full entries, size %d bytes\n", t->im_fn, t->im_fs);
	printf("  %d partial entries, size %d bytes\n", t->im_pn, t->im_ps);
	printf("  to hold %d output values of %d bits\n", g->od, f->imovb);

#endif /* VERBOSE */

	/* Number of bits needed for the weighting value */
	f->wemnb = g->prec+1;	/* Need to hold a weighting factor of 0 - 256 for 8 bits */
							/* Need to hold a weighting factor of 0 - 65536 for 16 bits */

	/* Variable that would be used to hold it */
	if ((f->wevt = findnord(f, f->wemnb)) < 0) {
		fprintf(stderr,"Can't find entry size to hold weighting variable\n");
		exit(-1);
	}

	/* Number of bits needed for vertex offset value */
	f->vomnb = calc_obits(g->id, g->itres, t->vo_om);

	/* Variable that would be used to hold it */
	if ((f->vovt = findnord(f, f->vomnb)) < 0) {
		fprintf(stderr,"Can't find entry size to hold vertex offset variable\n");
		exit(-1);
	}

	if (t->sort) {
		/* If we are using an explicit sort, we need to figure how many */
		/* separate entries we need to use to hold the interpolation index, */
		/* weighting factor and vertex offset values in the input table. */

		/* First try all three in one entry */
		if ((f->itet = findord(f, f->ixmnb + f->wemnb + f->vomnb)) >= 0) {/* size to read */
			int rem;						/* Remainder bits */

			t->it_xs = 0;					/* Combined interp+weight+offset */
			t->wo_xs = 0;
			t->it_ab = a->ords[f->itet].bits;	/* Bits in combined input entry */
			rem = t->it_ab - f->ixmnb - f->wemnb - f->vomnb; /* Spair bits */
			t->we_ab = f->wemnb;				/* Get minimum weight bits */
			t->vo_ab = f->vomnb + rem/2;		/* vertex offset index bits actually available */
			t->ix_ab = t->it_ab - t->vo_ab - t->we_ab;	/* interp index bits actually available */
			t->wo_ab = t->we_ab + t->vo_ab;		/* Weight & offset total bits */
			t->it_ts = a->ords[f->itet].bits/8;	/* total size in bytes */
			f->itvt = nord(f, f->itet);			/* Variable type */

			if ((f->wovt = findnord(f, t->we_ab + t->vo_ab)) < 0) {
				fprintf(stderr,"Can't find variable size to hold weight/offset\n");
				exit(-1);
			}
			if ((f->wevt = findnord(f, t->we_ab)) < 0) {
				fprintf(stderr,"Can't find variable size to hold weighting factor\n");
				exit(-1);
			}
			if ((f->vovt = findnord(f, t->vo_ab)) < 0) {
				fprintf(stderr,"Can't find variable size to hold vertex offset index\n");
				exit(-1);
			}
			if ((f->ixvt = findnord(f, t->ix_ab)) < 0) {
				fprintf(stderr,"Interp index variable size can't be handled\n");
				exit(-1);
			}
		} else {	/* Interp index will be a separate entry */
			int wit, oft, bigt;		/* weighting type, offset type, biggest type */
			int combt;				/* Combined type */
			int sepbits, combits;	/* Total separate, combined bits */

			t->it_xs = 1;				/* Separate interp index and weighting+offset */
			if ((f->ixet = findord(f, f->ixmnb)) < 0) {
				fprintf(stderr,"Interp index entry size can't be handled\n");
				exit(-1);
			}
			f->ixvt = nord(f, f->ixet);		/* Variable type */
			t->ix_ab = a->ords[f->ixet].bits;
			t->ix_es = t->ix_ab/8;
			t->ix_eo = 0;
			t->it_ts = t->ix_es;			/* Input table size so far */

			/* Now figure weighting and vertex offset */

			/* See if we can fit them into separately readable entries, or whether */
			/* they should be combined to minimise overall table size. */

			if ((wit = findord(f, f->wemnb)) < 0) {
				fprintf(stderr,"Can't find entry size to hold weighting factor\n");
				exit(-1);
			}
			if ((oft = findord(f, f->vomnb)) < 0) {
				fprintf(stderr,"Can't find entry size to hold vertex offset index\n");
				exit(-1);
			}
			bigt = wit > oft ? wit : oft;			/* Bigest separate type */

			if ((combt = findord(f, f->wemnb + f->vomnb)) < 0) {/* Combined isn't possible */
				sepbits = 2 * a->ords[bigt].bits;		/* Total separate bits */
				combits = sepbits;						/* Force separate entries */
			} else {
				sepbits = 2 * a->ords[bigt].bits;		/* Total separate bits */
				combits = a->ords[combt].bits;			/* Total combined bits */
			}

			if (sepbits <= combits) {				/* We will use separate entries */
				t->wo_xs = 1;
				t->we_es = a->ords[bigt].bits/8;	/* size in bytes for weighting entry */
				t->we_ab = a->ords[bigt].bits;		/* bits available for weighting */
				t->we_eo = t->ix_es;				/* Entry offset in input table */
				t->vo_es = a->ords[bigt].bits/8;	/* size in bytes for vertex offset entry */
				t->vo_ab = a->ords[bigt].bits;		/* bits available for vertex offset */
				t->vo_eo = t->ix_es + t->we_es;		/* Entry offset in input table */
				t->wo_es = t->we_es + t->vo_es;		/* Total entry size for each vertex */
				t->it_ts += t->we_es + t->vo_es;	/* Total input entry size in bytes */

				f->weet = bigt;				/* Variable type for accessing weighting entry */
				f->voet = bigt;				/* Variable type for accessing vertex offset entry */
				f->wevt = nord(f, wit);		/* Variable type for holding weight value */
				f->vovt = nord(f, oft);		/* Variable type for holding offset value */

			} else {								/* We will combine the two entries */