imdi_tab.c

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

#endif

		/* Get ready to access all the entries in the table */
		p = t;
		PH_INIT(phc, gs->id, gs->itres)

		/* Create all the interpolation table entry values */
		do {
			int ee, f;
			double riv[IXDI];	/* Real input values */
			double rev[IXDO];	/* Real entry values */
			unsigned long iev; 
			byte *pp;			/* Pointer to sub-entry */

			for (e = 0, p = t; e < gs->id; e++) {
				riv[e] = ((double)phc[e]) / (gs->itres - 1.0);
				p += phc[e] * ibdinc[e];		/* Compute pointer to entry value */
			}

			/* Lookup this verticies value */
			md_table(cntx, rev, riv);

			/* Create all the output values */

			/* I'm trying to avoid having to declare the actual entry sized */
			/* variables, since it is difficult dynamically. */

			/* For all the full entries */
			f = 0;
			pp = p;
			for (e = 0; e < ts->im_fn; e++, pp += ts->im_fs) {
				/* For all channels within full entry */
				for (ee = 0; ee < ts->im_fv; ee++, f++) {
					double revf = rev[f];
					if (revf < 0.0)						/* Guard against sillies */
						revf = 0.0;
					else if (revf > 1.0)
						revf = 1.0;
					iev = (unsigned long)(revf * vscale + 0.5);

					if (bigend) {
						write_entry[vsize](pp + (ts->im_fs - (ee+1) * vsize), iev);
					} else {
						write_entry[vsize](pp + ee * vsize, iev);
					}
				}
			}

			/* For all the 0 or 1 partial entry */
			for (e = 0; e < ts->im_pn; e++) {
				/* For all channels within partial entry */
				for (ee = 0; ee < ts->im_pv; ee++, f++) {
					double revf = rev[f];
					if (revf < 0.0)						/* Guard against sillies */
						revf = 0.0;
					else if (revf > 1.0)
						revf = 1.0;
					iev = (unsigned long)(revf * vscale + 0.5);

					if (bigend) {
						write_entry[vsize](pp + (ts->im_ps - (ee+1) * vsize), iev);
					} else {
						write_entry[vsize](pp + ee * vsize, iev);
					}
				}
			}
#ifdef ASSERTS
			if (f != gs->od)
				fprintf(stderr,"imdi_tab assert: f == gs->od\n");
#endif

			PH_INC(phc)
		} while (!PH_LOOPED(phc));

		/* Put table into place */
		it->im_table = (void *)t;
	}

	/* Setup the simplex table */
	if (ts->sort) {
		it->sw_table = (void *)NULL;

	} else {
		byte *t, *p;		/* Pointer to input table, pointer to total entry */
		int nsplx;			/* Total number of simplexes */
		XCOMBO(vcmb, gs->id+1, 1 << gs->id);/* Simplex dimension id out of cube dimention id */
		int comb[24][IXDI];	/* Parameter[id]->Absolute[id] coordinate index */
		int ps[IXDI+1];		/* Base simplex parameter space counter */
		int pse;			/* Base simplex parameter space counter index */
		int idioff;			/* Interpolation table diagonal offset value */

		if (gs->id > 4) {
			fprintf(stderr,"imdi_tabl: internal failure - trying to create simplex table with di > 4!\n");
			exit(-1);
		}

		/* Allocate the table */
		if ((t = (byte *)malloc(sbdinc[gs->id])) == NULL) {
#ifdef VERBOSE
			printf("malloc imdi simplex table size %d failed\n",sbdinc[gs->id]);
#endif
			return NULL;	/* Should we signal error ? How ? */
		}
#ifdef VERBOSE
		printf("Allocated simplex table = %d bytes\n",sbdinc[gs->id]);
#endif

		/* Compute the interp table offset to the diagonal vertex */
		for (idioff = 0, e = 0; e < gs->id; e++)
			idioff += idinc[e];		/* Sum one offset in each dimension */

		/* Figure out how many simplexes fit into this dimension cube, */
		/* and how to map from the base simplex to each actual simplex. */
		XCB_INIT(vcmb);
		for (nsplx = 0; ;) {
			int i;

			/* XCOMB generates verticies in order from max to min offest */
	
			/* Compute Absolute -> Parameter mapping */
			for (e = 0; e < gs->id; e++) {		/* For each absolute axis */
				for (i = 0; i < gs->id; i++) {	/* For each verticy, order large to small */
					if ((vcmb[i]   & (1<<e)) != 0 && 
					    (vcmb[i+1] & (1<<e)) == 0) {/* Transition from offset 1 to 0 */
						comb[nsplx][i] = e;
						break;
					}
				}
			}
	
//printf("~~Verticies   = ");
//for (i = 0; i <= gs->id; i++)
//	printf("%d ",vcmb[i]);
//printf("\n");

//printf("~~Parm -> Abs = ");
//for (e = 0; e < gs->id; e++)
//	printf("%d ",comb[nsplx][e]);
//printf("\n");

			/* Increment the counter value */
			XCB_INC(vcmb);
			nsplx++;
			if (XCB_DONE(vcmb))
				break;
		}

		/* Now generate the contents of the base simplex, */
		/* and map it to all the symetrical simplexes */

		/* Init parameter space counter. */
		/* Note that ps[id-1] >= ps[id-2] >=  ... >= ps[1] >= ps[0] */
		for (pse = 0; pse < gs->id; pse++)
			ps[pse] = 0;
		ps[pse] = gs->stres-1;

		/* Itterate through the simplex parameter space */
		for (pse = 0; pse < gs->id;) {
			double qps[IXDI];	/* Quantized parameter values */
			int    we[IXDI+1];	/* Baricentric coords/vertex weighting */
			double wvscale = (1 << gs->prec);	/* Weighting value scale */
			int sx;				/* Simplex */

//printf("Param coord =");
//for (e = gs->id-1; e >= 0; e--) {
//	printf(" %d",ps[e]);
//}
//printf("\n");
			for (e = 0; e < gs->id; e++) {
				/* (Should try wvscale + 0.49999999, or something ?) */
				double tt = (wvscale * (double)ps[e])/((double)gs->stres);
				qps[e] = (int)(tt + 0.5);
			}
	
			/* Convert quantized parameter values into weighting values */
			we[gs->id] = (1 << gs->prec) - qps[gs->id-1];
			for (e = gs->id-1; e > 0; e--)
				we[e] = qps[e] - qps[e-1];
			we[0] = qps[0];

#ifdef ASSERTS
			{
				int sow = 0;
				for (e = gs->id; e >= 0; e--)
					sow += we[e];
				
				if (sow != (1 << gs->prec))
					fprintf(stderr,"imdi_tab assert: sum weights == (1 << gs->prec)\n");
			}
#endif

//printf("Baricentric coord =");
//for (e = gs->id; e >= 0; e--) {
//	printf(" %d",we[e]);
//}
//printf("\n");

			/* For each simplex, compute the interp. and */
			/* and entry offsets, and write the entry. */
			for (sx = 0; sx < nsplx; sx++ ) {
				int v;					/* Vertex index */
				byte *pp;				/* Pointer to sub-entry */
				unsigned long vofb;		/* Vertex offset, base */
				unsigned long vwe;		/* Vertex weight */

				for (e = 0, p = t; e < gs->id; e++) {
					int ee = comb[sx][e];		/* Absolute coord index */
					p += ps[e] * sbdinc[ee];	/* Pointer to entry */
				}

				/* For each vertex entry */
				for (v = 0, pp = p; v <= gs->id; v++) {
					unsigned long vof;
					if (v == 0) {
						vofb = idioff;		/* Start at diagonal offset */
					} else {
						vofb -= idinc[comb[sx][v-1]];/* Move to next vertex */
					}
					vwe = we[v];					/* Weight for this vertex */

					if (vwe == 0)
						vof = 0;					/* Use zero offset if weight is zero */
					else
						vof = vofb * ts->vo_om;		/* Strength reduce kernel scaling */

					/* Write vwe and vof to entry */
					if (ts->wo_xs) {	/* Separate entries */
						write_entry[ts->we_es](pp + ts->we_eo, vwe);
						write_entry[ts->vo_es](pp + ts->vo_eo, vof);
						pp += ts->wo_es;
					} else {			/* Combined entries */
#ifdef ALLOW64
						unsigned longlong iwo;
#else
						unsigned long iwo;
#endif
						iwo = (vwe << ts->vo_ab) | vof; 	/* Combined weight+vertex offset */
						write_entry[ts->wo_es](pp + ts->wo_eo, iwo);
						pp += ts->wo_es;
					}
				}

				/* Assert vofb == 0 */
#ifdef ASSERTS
				if (vofb != 0)
					fprintf(stderr,"imdi_tab assert: vofb == 0\n");
#endif
			}	/* Next simplex */

			/* Increment the parameter coords */
			for (pse = 0; pse < gs->id; pse++) {
				ps[pse]++;
				if (ps[pse] <= ps[pse+1])
					break;	/* No carry */
				ps[pse] = 0;
			}
		}

		/* Put table into place */
		it->sw_table = (void *)t;
	}
	
	/* Last, setup the output tables */
	for (e = 0; e < gs->od; e++) {
		byte *t, *p;	/* Pointer to output table, entry pointer */
		int ne;			/* Number of entries */
		int iiv;		/* Integer input value */
		double ivr = (double)((1 << gs->prec)-1); /* Input value range */
		double ovr = (double)((1 << ts->ot_bits[e])-1); /* Output value range */
		int	osb = (1 << (ts->ot_bits[e]-1)); /* Output offset to signed displacement */
		int ooff = ts->ot_off[e];	/* Output value bit offset */

		ne = (1 << gs->prec);		/* Output of clut is prec bits */

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

		/* For each possible output value, compute the entry value */
		for (iiv = 0, p = t; iiv < ne; iiv++, p += ts->ot_ts) {
			int ee;
			double riv;		/* Real input value, 0.0 - 1.0 */
			double rtv;		/* Real transformed value, 0.0 - 1.0 */
			unsigned long iov;	/* Integer output value */

			riv = (double) iiv / ivr;			/* Compute floating point */
			rtv = output_curve(cntx, e, riv);	/* Lookup the output table transform */
			if (rtv < 0.0)						/* Guard against sillies */
				rtv = 0.0;
			else if (rtv > 1.0)
				rtv = 1.0;
			iov = (unsigned long)(rtv * ovr + 0.5);	/* output value */
			if (gs->out_signed & (1 << e))		/* Treat output as signed */
				iov = (iov >= osb) ? iov - osb : iov + osb; /* Convert to signed from offset */
			iov <<= ooff;						/* Aligned for output */

			write_entry[ts->ot_ts](p, iov);		/* Write entry */
		}

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

#ifdef VERBOSE
	printf("imdi_tabl returning OK\n");
#endif
	return it;
}

/* Free up the data allocated */
void
imdi_tab_free(
imdi_imp *it
) {
	int e;

	for (e = 0; e < it->nintabs; e++)
		free(it->in_tables[e]);

	free(it->sw_table);
	free(it->im_table);

	for (e = 0; e < it->nouttabs; e++)
		free(it->out_tables[e]);

	free(it);

}