imdi_tab.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.
 */

/* Run time table allocater and initialiser */

/*
 * The function here that knows how to create the
 * appropriate run time tables for our chosen kernel,
 * and the color mapping we want to perform.
 */

#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 ASSERTS			/* Check asserts */

#ifdef ASSERTS
#include <numlib.h>
#endif

typedef unsigned char byte;


/* Left shift, handles >= 32 properly */
#define LSHIFT(aa, bb)  ((bb) <= 31 ? ((aa) << (bb)) : (((aa) << 31) << ((bb)-31)))

/* Specific entry size write routine */

void write_uchar(
byte *p,
#ifdef ALLOW64
unsigned longlong v
#else
unsigned long v
#endif
) {
	*((unsigned char *)p) = (unsigned char)v;
}

void write_ushort(
byte *p,
#ifdef ALLOW64
unsigned longlong v
#else
unsigned long v
#endif
) {
	*((unsigned short *)p) = (unsigned short)v;
}

void write_uint(
byte *p,
#ifdef ALLOW64
unsigned longlong v
#else
unsigned long v
#endif
) {
	*((unsigned int *)p) = (unsigned int)v;
}

void write_ulong(
byte *p,
#ifdef ALLOW64
unsigned longlong v
#else
unsigned long v
#endif
) {
	*((unsigned long *)p) = (unsigned long)v;
}

#ifdef ALLOW64
void write_ulonglong(
byte *p,
unsigned longlong v
) {
	*((unsigned longlong *)p) = (unsigned longlong)v;
}
#endif /* ALLOW64 */

void write_default(
byte *p,
#ifdef ALLOW64
unsigned longlong v
#else
unsigned long v
#endif
) {
	fprintf(stderr,"imdi_tabl: internal failure - unexpected write size!\n");
*((char *)NULL) = 0;	// ~~999
	exit(-1);
}

/* Array of write routines */
#ifdef ALLOW64
void (*write_entry[16])(byte *p, unsigned longlong v);
#else
void (*write_entry[16])(byte *p, unsigned long v);
#endif

void
init_write_tab(void) {
	int i;

	for (i = 0; i < 16; i++)
		write_entry[i] = write_default;		/* Make sure any un-inited access bombs */

	write_entry[sizeof(unsigned char)]  = write_uchar;
	write_entry[sizeof(unsigned short)] = write_ushort;
	write_entry[sizeof(unsigned int)]   = write_uint;
	write_entry[sizeof(unsigned long)]  = write_ulong;
#ifdef ALLOW64
	write_entry[sizeof(unsigned longlong)]  = write_ulonglong;
#endif /* ALLOW64 */
}

imdi_imp *
imdi_tab(
	genspec *gs,	/* Pointer to gen spec */
	tabspec *ts,	/* Pointer to table spec */

	/* Callbacks to lookup the mdi table values */
	double (*input_curve) (void *cntx, int ch, double in_val),
	void   (*md_table)    (void *cntx, double *out_vals, double *in_vals),
	double (*output_curve)(void *cntx, int ch, double in_val),
	void *cntx		/* Context to callbacks */
) {
	static inited = 0;
	static bigend = 0;
	int e;
	imdi_imp *it;
	unsigned long etest = 0xff;
	int idinc[IXDI+1];		/* Increment for each dimension of interp table. */
	int ibdinc[IXDI+1];		/* idinc[] in bytes */
	int sdinc[IXDI+1];		/* Increment for each dimension of simplex table. */
	int sbdinc[IXDI+1];		/* sdinc[] in bytes */

#ifdef VERBOSE
	printf("imdi_tab called\n");
#endif

	if (inited == 0) {
		init_write_tab();
		if (*((unsigned char *)&etest) == 0xff)
			bigend = 0;		/* Little endian */
		else
			bigend = 1;		/* Big endian */
		inited = 1;
	}

	if ((it = (imdi_imp *)malloc(sizeof(imdi_imp))) == NULL) {
#ifdef VERBOSE
		printf("malloc imdi_imp size %d failed\n",sizeof(imdi_imp));
#endif
		return NULL;	/* Should we signal error ? How ? */
	}

	/* Compute interp and simplex table dimension increments & total sizes */
	idinc[0]  = 1;
	ibdinc[0] = ts->im_ts;
	for (e = 1; e <= gs->id; e++) {
		idinc[e]  = idinc[e-1]  * gs->itres;
		ibdinc[e] = ibdinc[e-1] * gs->itres;
	}

	if (!ts->sort) {
		sdinc[0]  = 1;
		sbdinc[0] = ts->sm_ts;
		for (e = 1; e <= gs->id; e++) {
			sdinc[e]  = sdinc[e-1]  * gs->stres;
			sbdinc[e] = sbdinc[e-1] * gs->stres;
		}
	}

	/* First we setup the input tables */
	for (e = 0; e < gs->id; e++) {
		byte *t, *p;	/* Pointer to input table, entry pointer */
		int ne;			/* Number of entries */
		int ex;			/* Entry index */
		int ix = 0;		/* Extract flag */

		/* Compute number of entries */
		if (ts->it_ix && !gs->in.packed) {	/* Input is the whole bpch[] size */
			ix = 1;					/* Need to do extraction in lookup */
			if (gs->in.pint) {
				ne = (1 << (gs->in.bpch[0]));		/* Same size used for all input tables */
			} else {
				ne = (1 << (gs->in.bpch[e]));		/* This input channels size */
			}
		} else {				/* Input is the value size */
			ne = (1 << (gs->in.bpv[e]));		/* This input values size */
		}

		/* Allocate the table */
		if ((t = (byte *)malloc(ts->it_ts * ne)) == NULL) {
#ifdef VERBOSE
			printf("malloc imdi input table size %d failed\n",ts->it_ts * ne);
#endif
			return NULL;	/* Should we signal error ? How ? */
		}

		/* For each possible input value, compute the entry value */
		for (ex = 0, p = t; ex < ne; ex++, p += ts->it_ts) {
			int ee;
			int iiv;		/* Integer input value */
			int ivr;		/* Input value range */
			int isb;		/* Input sign bit/signed to offset displacement */
			double riv;		/* Real input value, 0.0 - 1.0 */
			double rtv;		/* Real transformed value, 0.0 - 1.0 */
			double rmi;		/* Real interpolation table index */
			double rsi;		/* Real simplex index */
			int imi;		/* Interpiolation table index */
			int isi;		/* Integer simplex index */
			int iwe;		/* Integer weighting value */
			int vo;			/* Vertex offset value */

			if (ix) {		/* Extract value from index */
				ivr = ((1 << (gs->in.bpv[e])) -1);
				iiv = (ex >> gs->in.bov[e]) & ((1 << (gs->in.bpv[e])) -1);
			} else {
				ivr = (ne - 1);	/* (Should be bpv[e], but take no chances!) */
				iiv = ex;	/* Input value is simply index */
			}
			isb = ivr & ~(((unsigned int)ivr) >> 1);			/* Top bit */
			if (gs->in_signed & (1 << e))		/* Treat input as signed */
				iiv = (iiv & isb) ? iiv - isb : iiv + isb;	/* Convert to offset from signed */
			riv = (double) iiv / (double)ivr;	/* Compute floating point */
			rtv = input_curve(cntx, e, riv);	/* Lookup the input table transform */
			if (rtv < 0.0)						/* Guard against sillies */
				rtv = 0.0;
			else if (rtv > 1.0)
				rtv = 1.0;

			/* divide into interp base and cube sub index */
			rmi = rtv * (gs->itres - 1);
			imi = (int)floor(rmi);		/* Interp. entry coordinate */
			if (imi >= (gs->itres-1))	/* Keep cube base one row back from far edge */
				imi = gs->itres-2;
			rsi = rmi - (double)imi;	/* offset into entry cube */
			if (ts->sort) {
				iwe = (int)((rsi * (1 << gs->prec)) + 0.5);	/* Weighting scale */
				vo = idinc[e] * ts->vo_om;	/* Vertex offset */
			} else {
				isi = (int)((rsi * gs->stres) + 0.5);
				if (isi == gs->stres) {		/* Keep simplex index within table */
					isi = 0;
					imi++;					/* Move to next interp. lattice */
				}
				isi *= sdinc[e]; 		/* Convert the raw indexes into offset in this dim */
			}
			imi *= idinc[e]; 			/* Convert the raw indexes into offset in this dim */

#ifdef ASSERTS
			/* ~~~ needs fixing for sort ~~~~ */
			if ((imi & (LSHIFT(1,ts->it_ab)-1)) != imi)
				error("imdi_tab assert: (imi & ((1 << ts->it_ab)-1)) != imi, imi = 0x%x, it_ab = 0x%x\n",imi,ts->it_ab);
			if (imi >= idinc[gs->id])
				error("imdi_tab assert: imi >= idinc[gs->id]\n");
			if ((isi & (LSHIFT(1,ts->sx_ab)-1)) != isi) 
				error("imdi_tab assert: (isi & ((1 << ts->sx_ab)-1)) != isi, isi = 0x%x, sx_ab = 0x%x\n",isi,ts->sx_ab);
			if (!ts->sort && isi >= sdinc[gs->id])
				error("imdi_tab assert: isi >= sdinc[gs->id]\n");
#endif

			/* Now stuff them into the table entry */
			if (ts->sort) {
				if (ts->it_xs) {		/* Separate interp index and weight/offset*/
					if (ts->wo_xs) {	/* All 3 are separate */
						write_entry[ts->ix_es](p + ts->ix_eo, imi);
						write_entry[ts->we_es](p + ts->we_eo, iwe);
						write_entry[ts->vo_es](p + ts->vo_eo, vo);
					} else {
#ifdef ALLOW64
						unsigned longlong iwo;
#else
						unsigned long iwo;
#endif
	
						iwo = (iwe << ts->vo_ab) | vo; 	/* Combined weight+vertex offset */
						write_entry[ts->ix_es](p + ts->ix_eo, imi);
						write_entry[ts->wo_es](p + ts->wo_eo, iwo);
					}
				} else {			/* All 3 are combined  */
#ifdef ALLOW64
					unsigned longlong iit;
#else
					unsigned long iit;
#endif
					iit = (((imi << ts->we_ab) | iwe) << ts->vo_ab) | vo;
					write_entry[ts->it_ts](p, iit);
				}
			} else {
				if (ts->it_xs) {		/* Separate interp index and weight/offset*/
					write_entry[ts->ix_es](p + ts->ix_eo, imi);
					write_entry[ts->sx_es](p + ts->sx_eo, isi);
				} else {
#ifdef ALLOW64
					unsigned longlong iit;
#else
					unsigned long iit;
#endif
					iit = (imi << ts->sx_ab) | isi;	/* Combine interp and simplex indexes */
					write_entry[ts->it_ts](p, iit);
				}
			}
		}

		/* Put table into place */
		it->in_tables[e] = (void *)t;
	}
	it->nintabs = e;

	/* Setup the interpolation table */
	{
		byte *t, *p;	/* Pointer to interp table, pointer to total entry */
		PHILBERT(phc)	/* Pseudo Hilbert counter */
		double vscale;	/* Value scale for fixed point */
		int vsize;		/* Fixed point storage size */

		if (ts->im_cd)
			vsize = (gs->prec * 2)/8;	/* Fixed point entry & computation size */
		else
			vsize = gs->prec/8;			/* Fixed point entry size */
		vscale = (1 << gs->prec) -1.0;	/* Value scale for fixed point padding */
										/* -1.0 is to prevent carry after accumulation */
		/* Allocate the table */
		if ((t = (byte *)malloc(ibdinc[gs->id])) == NULL) {
#ifdef VERBOSE
			printf("malloc imdi interpolation table size %d failed\n",ibdinc[gs->id]);
#endif
			return NULL;	/* Should we signal error ? How ? */
		}
#ifdef VERBOSE
		printf("Allocated grid table = %d bytes\n",ibdinc[gs->id]);