804 lines
23 KiB
C
804 lines
23 KiB
C
|
|
/* Integer Multi-Dimensional Interpolation */
|
|
|
|
/*
|
|
* Copyright 2000 - 2007 Graeme W. Gill
|
|
* All rights reserved.
|
|
*
|
|
* This material is licenced under the GNU AFFERO GENERAL PUBLIC LICENSE Version 3 :-
|
|
* see the License.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 type color mapping we want to perform.
|
|
*/
|
|
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <math.h>
|
|
#include <stdarg.h>
|
|
#include <string.h>
|
|
|
|
#include "imdi.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)))
|
|
|
|
/* The big value type used to represent table entries */
|
|
#ifdef ALLOW64
|
|
typedef unsigned longlong bvt;
|
|
#else
|
|
typedef unsigned long bvt;
|
|
#endif
|
|
|
|
/* Specific entry size write routine */
|
|
|
|
void write_uchar(
|
|
byte *p,
|
|
bvt v
|
|
) {
|
|
*((unsigned char *)p) = (unsigned char)v;
|
|
}
|
|
|
|
void write_ushort(
|
|
byte *p,
|
|
bvt v
|
|
) {
|
|
*((unsigned short *)p) = (unsigned short)v;
|
|
}
|
|
|
|
void write_uint(
|
|
byte *p,
|
|
bvt v
|
|
) {
|
|
*((unsigned int *)p) = (unsigned int)v;
|
|
}
|
|
|
|
void write_ulong(
|
|
byte *p,
|
|
bvt v
|
|
) {
|
|
*((unsigned long *)p) = (unsigned long)v;
|
|
}
|
|
|
|
#ifdef ALLOW64
|
|
void write_ulonglong(
|
|
byte *p,
|
|
bvt v
|
|
) {
|
|
*((unsigned longlong *)p) = (unsigned longlong)v;
|
|
}
|
|
#endif /* ALLOW64 */
|
|
|
|
void write_default(
|
|
byte *p,
|
|
bvt v
|
|
) {
|
|
fprintf(stderr,"imdi_tabl: internal failure - unexpected write size!\n");
|
|
exit(-1);
|
|
}
|
|
|
|
/* Array of write routines */
|
|
void (*write_entry[16])(byte *p, bvt v);
|
|
|
|
static 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 */
|
|
}
|
|
|
|
|
|
/* Input offset adjustment table */
|
|
double in_adj[] = {
|
|
8.0324820232182659e+281, 1.3051220361353854e+214, 1.5654418860154115e-076,
|
|
6.6912978722165055e+281, 1.2369092402930559e+277, 1.4097588049607207e-308,
|
|
7.7791723264456369e-260, 3.6184161952648606e+238, 5.8235640814908141e+180,
|
|
9.1271554315814989e-072, 5.4310198502711138e+241, 2.7935452404894958e+275,
|
|
-2.9408705449902027e+003
|
|
};
|
|
|
|
|
|
/* Table creation function */
|
|
imdi_imp *
|
|
imdi_tab(
|
|
genspec *gs, /* Pointer to gen spec */
|
|
tabspec *ts, /* Pointer to table spec */
|
|
imdi_conv cnv, /* Runtime argument conversion needed */
|
|
imdi_pixrep irep, /* High level input pixel representation to match */
|
|
imdi_pixrep orep, /* High level output pixel representation to match */
|
|
void (*interp)(struct _imdi *s, void **outp, int outst, /* Underlying conversion function */
|
|
void **inp, int inst,
|
|
unsigned int npixels),
|
|
|
|
int *inm, /* Input raster channel to callback channel mapping, NULL for none. */
|
|
int *outm, /* Output raster channel to callback channel mapping, NULL for none. */
|
|
imdi_ooptions oopt, /* Output per channel options (Callback channel, NOT output channel) */
|
|
unsigned int *checkv, /* Output channel check values (Callback channel, NULL for none == 0. */
|
|
|
|
/* Callbacks to lookup the mdi table values */
|
|
void (*input_curves) (void *cntx, double *out_vals, double *in_vals),
|
|
void (*md_table) (void *cntx, double *out_vals, double *in_vals),
|
|
void (*output_curves)(void *cntx, double *out_vals, double *in_vals),
|
|
void *cntx /* Context to callbacks */
|
|
) {
|
|
static int inited = 0;
|
|
static int bigend = 0;
|
|
int i, e, f;
|
|
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 *)calloc(1, 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 ? */
|
|
}
|
|
|
|
it->size = sizeof(imdi_imp);
|
|
#ifdef VERBOSE
|
|
printf("Allocated imdi_imp structure size %u\n",it->size);
|
|
#endif /* VERBOSE */
|
|
|
|
/* Set runtime matching conversion provided */
|
|
it->cnv = cnv;
|
|
it->id = gs->id;
|
|
it->od = gs->od;
|
|
it->cirep = irep; /* Pixel representation interp is called with */
|
|
it->corep = orep;
|
|
it->firep = gs->irep; /* Pixel representation of function we are going to use */
|
|
it->forep = gs->orep;
|
|
it->interp = interp;
|
|
it->checkf = 0;
|
|
|
|
/* Compute number of written channels (allow for skip) */
|
|
it->wod = it->od;
|
|
for (i = 0; i < it->od; i++) {
|
|
if ((oopt & OOPT(oopts_skip,i)) != 0)
|
|
it->wod--;
|
|
}
|
|
|
|
/* Setup the raster to callback channel mappings */
|
|
if (inm != NULL) {
|
|
for (e = 0; e < it->id; e++)
|
|
it->it_map[e] = inm[e]; /* Copy input */
|
|
} else {
|
|
for (e = 0; e < it->id; e++)
|
|
it->it_map[e] = e; /* Direct mapping */
|
|
}
|
|
if (outm != NULL) {
|
|
for (e = 0; e < it->od; e++)
|
|
it->im_map[e] = outm[e]; /* Copy input */
|
|
} else {
|
|
for (e = 0; e < it->od; e++)
|
|
it->im_map[e] = e; /* Direct mapping */
|
|
}
|
|
if (checkv != NULL) {
|
|
for (e = 0; e < it->od; e++)
|
|
it->checkv[e] = checkv[it->im_map[e]]; /* Copy input and convert to Output index */
|
|
} else {
|
|
for (e = 0; e < it->od; e++)
|
|
it->checkv[e] = 0; /* Set to zero */
|
|
}
|
|
|
|
/* Compute interp and simplex table dimension increments & total sizes */
|
|
idinc[0] = 1;
|
|
ibdinc[0] = ts->im_ts;
|
|
for (e = 1; e <= it->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 <= it->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 < it->id; e++) {
|
|
byte *t, *p; /* Pointer to input table, entry pointer */
|
|
int ne; /* Number of entries */
|
|
int ex; /* Entry index */
|
|
double iaf;
|
|
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 ? */
|
|
}
|
|
it->size += ts->it_ts * ne;
|
|
#ifdef VERBOSE
|
|
printf("Allocated input table %d size %u = %u * %u\n",e, ts->it_ts * ne,ts->it_ts,ne);
|
|
#endif /* VERBOSE */
|
|
|
|
/* Compute input adjustment factor */
|
|
for (iaf = 0.0, i = 0; i < (sizeof(in_adj)/sizeof(double)-1); i++)
|
|
iaf += log(in_adj[i]);
|
|
iaf += in_adj[i];
|
|
|
|
/* For each possible input value, compute the entry value */
|
|
for (ex = 0, p = t; ex < ne; ex++, p += ts->it_ts) {
|
|
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 = 0; /* Integer simplex index */
|
|
int iwe = 0; /* Integer weighting value */
|
|
int vo = 0; /* 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 */
|
|
{
|
|
double civ[IXDI], cov[IXDI];
|
|
for (f = 0; f < it->id; f++)
|
|
civ[f] = riv;
|
|
input_curves(cntx, cov, civ); /* Lookup the input table transform */
|
|
rtv = iaf * cov[it->it_map[e]];
|
|
}
|
|
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[it->id])
|
|
error("imdi_tab assert: imi >= idinc[it->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[it->id])
|
|
error("imdi_tab assert: isi >= sdinc[it->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 {
|
|
bvt iwo;
|
|
|
|
iwo = ((bvt)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 */
|
|
bvt iit;
|
|
|
|
iit = ((((bvt)imi << ts->we_ab) | (bvt)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 {
|
|
bvt iit;
|
|
|
|
iit = ((bvt)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) -0.50000001;
|
|
/* Value scale for fixed point padding */
|
|
/* -0.5 is to prevent carry/rollover after accumulation */
|
|
/* Could get better accuracy with saturation arithmatic */
|
|
|
|
/* Allocate the table */
|
|
if ((t = (byte *)malloc(ibdinc[it->id])) == NULL) {
|
|
#ifdef VERBOSE
|
|
printf("malloc imdi interpolation table size %d failed\n",ibdinc[it->id]);
|
|
#endif
|
|
return NULL; /* Should we signal error ? How ? */
|
|
}
|
|
it->size += ibdinc[it->id];
|
|
#ifdef VERBOSE
|
|
printf("Allocated grid table = %u bytes, composed of %d dim of res %d entry %d\n",ibdinc[it->id], it->id, gs->itres, ts->im_ts);
|
|
#endif /* VERBOSE */
|
|
|
|
/* Get ready to access all the entries in the table */
|
|
p = t;
|
|
PH_INIT(phc, it->id, gs->itres)
|
|
|
|
/* Create all the interpolation table entry values */
|
|
do {
|
|
int ee, ff;
|
|
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 < it->id; e++) {
|
|
riv[e] = ((double)phc[e]) / (gs->itres - 1.0);
|
|
p += phc[e] * ibdinc[e]; /* Compute pointer to entry value */
|
|
}
|
|
|
|
/* Lookup this vertices value */
|
|
{
|
|
double mriv[IXDI]; /* Channel mapped real input values */
|
|
double mrev[IXDO]; /* Channel mapped real entry values */
|
|
for (e = 0; e < it->id; e++)
|
|
mriv[it->it_map[e]] = riv[e];
|
|
md_table(cntx, mrev, mriv);
|
|
for (e = 0; e < it->od; e++)
|
|
rev[e] = mrev[it->im_map[e]];
|
|
}
|
|
|
|
/* 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 */
|
|
ff = 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++, ff++) {
|
|
double revf = rev[ff];
|
|
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++, ff++) {
|
|
double revf = rev[ff];
|
|
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 != it->od)
|
|
fprintf(stderr,"imdi_tab assert: f == it->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, it->id+1, 1 << it->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 (it->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[it->id])) == NULL) {
|
|
#ifdef VERBOSE
|
|
printf("malloc imdi simplex table size %d failed\n",sbdinc[it->id]);
|
|
#endif
|
|
return NULL; /* Should we signal error ? How ? */
|
|
}
|
|
it->size += sbdinc[it->id];
|
|
#ifdef VERBOSE
|
|
printf("Allocated simplex table = %u bytes, composed of %d dim of res %d entry %d\n",sbdinc[it->id], it->id, gs->stres, ts->sm_ts);
|
|
#endif /* VERBOSE */
|
|
|
|
/* Compute the interp table offset to the diagonal vertex */
|
|
for (idioff = 0, e = 0; e < it->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 vertices in order from max to min offest */
|
|
|
|
/* Compute Absolute -> Parameter mapping */
|
|
for (e = 0; e < it->id; e++) { /* For each absolute axis */
|
|
for (i = 0; i < it->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("~~Vertices = ");
|
|
//for (i = 0; i <= it->id; i++)
|
|
// printf("%d ",vcmb[i]);
|
|
//printf("\n");
|
|
|
|
//printf("~~Parm -> Abs = ");
|
|
//for (e = 0; e < it->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 < it->id; pse++)
|
|
ps[pse] = 0;
|
|
ps[pse] = gs->stres-1;
|
|
|
|
/* Itterate through the simplex parameter space */
|
|
for (pse = 0; pse < it->id;) {
|
|
int 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 = it->id-1; e >= 0; e--) {
|
|
// printf(" %d",ps[e]);
|
|
//}
|
|
//printf("\n");
|
|
|
|
for (e = 0; e < it->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[it->id] = (1 << gs->prec) - qps[it->id-1];
|
|
for (e = it->id-1; e > 0; e--)
|
|
we[e] = qps[e] - qps[e-1];
|
|
we[0] = qps[0];
|
|
|
|
#ifdef ASSERTS
|
|
{
|
|
int sow = 0;
|
|
for (e = it->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 = it->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 = 0; /* Vertex offset, base */
|
|
unsigned long vwe; /* Vertex weight */
|
|
|
|
for (e = 0, p = t; e < it->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 <= it->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 */
|
|
bvt iwo;
|
|
|
|
iwo = ((bvt)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 < it->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 < it->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 ? */
|
|
}
|
|
it->size += ts->ot_ts * ne;
|
|
#ifdef VERBOSE
|
|
printf("Allocated output table %d size %u = %u * %u\n",e, ts->ot_ts * ne,ts->ot_ts,ne);
|
|
#endif /* VERBOSE */
|
|
|
|
/* For each possible output value, compute the entry value */
|
|
for (iiv = 0, p = t; iiv < ne; iiv++, p += ts->ot_ts) {
|
|
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 */
|
|
{
|
|
double civ[IXDO], cov[IXDO];
|
|
for (f = 0; f < it->od; f++)
|
|
civ[f] = riv;
|
|
output_curves(cntx, cov, civ); /* Lookup the input table transform */
|
|
rtv = cov[it->im_map[e]];
|
|
}
|
|
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;
|
|
|
|
/* Adjust the check values for output value shift */
|
|
for (e = 0; e < it->od; e++) {
|
|
int ooff = ts->ot_off[e]; /* Output value bit offset */
|
|
it->checkv[e] <<= ooff; /* Aligned for output */
|
|
}
|
|
|
|
/* Setup the appropriate skip flags, indexed by Output channel */
|
|
it->skipf = 0;
|
|
if ((oopt & OOPTS_SKIP) != 0) {
|
|
int i;
|
|
for (i = 0; i < it->od; i++) {
|
|
if (oopt & OOPT(oopts_skip,it->im_map[i])) { /* Skip flag for this output chan */
|
|
it->skipf |= (1 << i);
|
|
}
|
|
}
|
|
}
|
|
|
|
/* Fill in some report information */
|
|
it->gres = gs->itres;
|
|
if (!ts->sort) {
|
|
it->sres = gs->stres;
|
|
} else {
|
|
it->sres = 0;
|
|
}
|
|
|
|
#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]);
|
|
|
|
if (it->sw_table != NULL)
|
|
free(it->sw_table);
|
|
if (it->im_table != NULL)
|
|
free(it->im_table);
|
|
|
|
for (e = 0; e < it->nouttabs; e++)
|
|
free(it->out_tables[e]);
|
|
|
|
free(it);
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|