summaryrefslogtreecommitdiff
path: root/imdi/imdi.c
diff options
context:
space:
mode:
Diffstat (limited to 'imdi/imdi.c')
-rw-r--r--imdi/imdi.c605
1 files changed, 605 insertions, 0 deletions
diff --git a/imdi/imdi.c b/imdi/imdi.c
new file mode 100644
index 0000000..8174c8e
--- /dev/null
+++ b/imdi/imdi.c
@@ -0,0 +1,605 @@
+
+/* 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.
+ */
+
+/*
+ * This is the implementation of the run time code.
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <stdarg.h>
+#include <string.h>
+
+#include "imdi.h"
+#include "imdi_tab.h"
+#include "imdi_k.h" /* Declaration of all the kernel functions */
+
+#undef VERBOSE
+#undef VVERBOSE
+
+static unsigned int imdi_get_check(imdi *im);
+static void imdi_reset_check(imdi *im);
+static void imdi_info(imdi *s, unsigned long *size, int *gres, int *sres);
+static void imdi_del(imdi *im);
+static void interp_match(imdi *s, void **outp, int outst, void **inp, int inst,
+ unsigned int npixels);
+
+
+/* Create a new imdi */
+/* Return NULL if request is not supported */
+/* Note that we use the high level pixel layout description to locate */
+/* a suitable run-time. */
+imdi *new_imdi(
+ int id, /* Number of input dimensions */
+ int od, /* Number of output lookup dimensions */
+ /* Number of output channels written = od - no. of oopt skip flags */
+ imdi_pixrep in, /* Input pixel representation */
+ int in_signed, /* Bit flag per channel, NZ if treat as signed */
+ int *inm, /* Input raster channel to callback channel mapping, NULL for none. */
+ imdi_iprec iprec, /* Internal processing precision */
+ imdi_pixrep out, /* Output pixel representation */
+ int out_signed, /* Bit flag per channel, NZ if treat as signed */
+ int *outm, /* Output raster channel to callback channel mapping, NULL for none. */
+ /* Mapping must include skipped channels. */
+ int res, /* Desired table resolution */
+ imdi_ooptions oopt, /* Output per channel options (by callback channel) */
+ unsigned int *checkv, /* Output channel check values (by callback channel, NULL == 0) */
+ imdi_options opt, /* Direction and stride options */
+
+ /* Callbacks to lookup the imdi table values. */
+ /* (Skip output channels are looked up) */
+ 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 */
+) {
+ int i;
+ int prec; /* Target internal precision */
+ int bk = -1; /* Best table */
+ int bfig = 0x7fffffff; /* Best tables figure of merit */
+ int bstres = 0; /* Best tables target stres */
+ genspec gs; /* Generation specification */
+ tabspec ts; /* Table specifications */
+ genspec bgs; /* Best gen spec */
+ tabspec bts; /* Best tab spec */
+ imdi_conv bcnv = conv_none; /* Best tables conversion flags */
+ imdi_ooptions Ooopt; /* oopt re-aranged to correspond to output channel index */
+
+ imdi *im;
+
+ /* Compute the Output channel index oopt mask */
+ if (outm == NULL)
+ Ooopt = oopt;
+ else {
+ int ff;
+ Ooopt = 0;
+ for (i = 0; i < od; i++) {
+ ff = OOPTX(oopt,outm[i]);
+ Ooopt |= OOPT(ff,i);
+ }
+ }
+
+#ifdef VERBOSE
+ printf("new_imdi called with id %d, od %d, res %d Ooopt 0x%x, opt 0x%x\n", id, od, res, Ooopt, opt);
+ printf("about to checking %d kernels\n", no_kfuncs);
+#endif
+
+ /* Figure out the internal precision requested */
+ COMPUTE_IPREC(prec, in, iprec, out)
+
+ /* Zero out the gen and tabspecs, to give diff a place to start */
+ memset((void *)&gs, 0, sizeof(genspec));
+ memset((void *)&ts, 0, sizeof(tabspec));
+
+ /* The first thing to do is see if there is an available kernel function */
+ for (i = 0; i < no_kfuncs; i++) {
+ int stres; /* Computed stres needed */
+ imdi_conv cnv = conv_none; /* Conversions needed for this choice */
+ int fig; /* Figure of merit - smaller is better */
+ ktable[i].gentab(&gs, &ts); /* Udate the kernel functions genspec and tabspec */
+
+#ifdef VERBOSE
+ printf("\n");
+ printf("kernel %d has id %d, od %d, irep %d orep %d oopt 0x%x, opt 0x%x\n",
+ i, gs.id, gs.od, gs.irep, gs.orep, gs.oopt, gs.opt);
+ printf("Input req is id %d, od %d, irep %d orep %d, oopt 0x%x, opt 0x%x\n",
+ id, od, in, out, Ooopt, opt);
+#endif
+ /* First check mandatory things */
+ if (!(
+ id == gs.id /* Input dimension */
+ && od == gs.od /* Output dimension */
+ )) {
+#ifdef VVERBOSE
+ printf(" Input or Output dimension mismatch\n");
+#endif
+ continue;
+ }
+
+ /* Check if we have an exact or conversion match for input stride */
+ if (!(
+ (opt & opts_istride) == (gs.opt & opts_istride)
+ || (!(opt & opts_istride) && (gs.opt & opts_istride))
+ )) {
+#ifdef VVERBOSE
+ printf(" Input stride mismatch\n");
+#endif
+ continue;
+ }
+
+ /* Check if we have an exact or conversion match for output stride */
+ if (!(
+ (opt & opts_ostride) == (gs.opt & opts_ostride)
+ || (!(opt & opts_ostride) && (gs.opt & opts_ostride))
+ )) {
+#ifdef VVERBOSE
+ printf(" Output stride mismatch\n");
+#endif
+ continue;
+ }
+
+ /* Check if we have an exact or conversion match for input representation */
+ if (!(
+ in == gs.irep
+ || (in == pixint8 && gs.irep == planeint8 && (gs.opt & opts_istride))
+ || (in == pixint16 && gs.irep == planeint16 && (gs.opt & opts_istride))
+ )) {
+#ifdef VVERBOSE
+ printf(" Input representation mismatch\n");
+#endif
+ continue;
+ }
+
+ /* Check if we have an exact or conversion match for output representation */
+ if (!(
+ out == gs.orep
+ || (out == pixint8 && gs.orep == planeint8 && (gs.opt & opts_ostride))
+ || (out == pixint16 && gs.orep == planeint16 && (gs.opt & opts_ostride))
+ )) {
+#ifdef VVERBOSE
+ printf(" Output representation mismatch\n");
+#endif
+ continue;
+ }
+
+ /* See if we have the output per channel options needed */
+ if (!((Ooopt & gs.oopt) == Ooopt)) {
+#ifdef VVERBOSE
+ printf(" Output per channel options mismatch\n");
+#endif
+ continue;
+ }
+
+ /* See if we have an exact or conversion match for reverse direction */
+ if (!(
+ ((!(opt & opts_bwd) && !(opt & opts_fwd))) /* Don't care */
+ || (((opt & opts_bwd) && (gs.opt & opts_bwd))) /* Match */
+ || (((opt & opts_fwd) && (gs.opt & opts_fwd))) /* Match */
+ || ((gs.opt & opts_istride) && (gs.opt & opts_ostride)) /* Can convert */
+ )) {
+#ifdef VVERBOSE
+ printf(" Direction mismatch\n");
+#endif
+ continue;
+ }
+
+#ifdef VERBOSE
+ printf(" found match\n");
+#endif
+ fig = 0;
+
+ /* Apply penalty if we will have to do a conversion match */
+ if ((opt & opts_istride) != (gs.opt & opts_istride)) {
+ cnv |= conv_istr;
+ fig += 1000;
+#ifdef VERBOSE
+ printf(" needs input stride conversion though\n");
+#endif
+ }
+
+ if ((opt & opts_ostride) != (gs.opt & opts_ostride)) {
+ cnv |= conv_ostr;
+ fig += 1000;
+#ifdef VERBOSE
+ printf(" needs output stride conversion though\n");
+#endif
+ }
+
+ if (in != gs.irep) {
+ cnv |= conv_irep;
+ fig += 5000;
+#ifdef VERBOSE
+ printf(" needs input representation conversion though\n");
+#endif
+ }
+
+ if (out != gs.orep) {
+ cnv |= conv_orep;
+ fig += 5000;
+#ifdef VERBOSE
+ printf(" needs output representation conversion though\n");
+#endif
+ }
+
+ /* If we don't need output checking or skipping, but we've got it */
+ if ((~Ooopt & gs.oopt) != 0) {
+ unsigned int tmp = (~Ooopt & gs.oopt);
+ /* Count number of bits set that we don't need */
+ /* (From HACKMEM 169) */
+ tmp = tmp - ((tmp >> 1) & 033333333333) - ((tmp >> 2) & 011111111111);
+ tmp = ((tmp + (tmp >> 3)) & 030707070707) % 63;
+ fig += tmp;
+#ifdef VERBOSE
+ printf(" has output options we don't need though\n");
+#endif
+ }
+
+ /* If we've got output skipping, we need to skip pointers in interp. */
+ if ((Ooopt & OOPTS_SKIP) != 0) {
+ cnv |= conv_skip;
+ }
+
+ if (((opt & opts_fwd) && (gs.opt & opts_bwd))
+ || ((opt & opts_bwd) && (gs.opt & opts_fwd))) {
+ cnv |= conv_rev;
+ fig += 1000;
+#ifdef VERBOSE
+ printf(" needs direction conversion though\n");
+#endif
+ }
+
+ if (prec != gs.prec) { /* Internal precision mismatch */
+ fig += 100000; /* Will have a major effect, so discourage using it */
+#ifdef VERBOSE
+ printf(" internal precision doesn't match though (want %d, got %d)\n",prec,gs.prec);
+#endif
+ }
+
+ /* If we've got a choice of algorithm possible */
+ if (gs.opt & (opts_splx_sort | opts_sort_splx)) {
+
+ /* If we've got a preference of algorithm, */
+ if ((opt & (opts_splx_sort | opts_sort_splx))) {
+
+ /* and there is a mismatch */
+ if (((opt & opts_splx_sort) && ts.sort != 0)
+ || ((opt & opts_sort_splx) && ts.sort == 0)) {
+ fig += 10000; /* Will have a great effect, so discourage using it */
+#ifdef VERBOSE
+ printf(" sort/simplex algorithm mismatch\n");
+#endif
+ }
+
+ } else { /* There is no preference chosen at run time, */
+ /* and there is a mismatch to the compile time preference */
+ if (((gs.opt & opts_splx_sort) && ts.sort != 0)
+ || ((gs.opt & opts_sort_splx) && ts.sort == 0)) {
+ fig += 10000; /* Will have a great effect, so discourage using it */
+#ifdef VERBOSE
+ printf(" sort/simplex algorithm mismatch\n");
+#endif
+ }
+ }
+ }
+
+ if (ts.sort) {
+ stres = 0;
+#ifdef VERBOSE
+ printf("gres = %d, gs.gres = %d\n",res,gs.itres);
+#endif
+ /* We want one that is equals or exceeds the desired */
+ /* resolution, but doesn't exceed it too much, or the code */
+ /* will be inefficient. */
+ /* If there are no routines that can meet the desired precision, */
+ /* then it is ok to use the one closest to the desired precision. */
+ if (gs.itres >= res) {
+ fig += 10 * (gs.itres - res);
+ } else {
+ fig += 10000 + 10 * (res - gs.itres);
+ }
+ } else {
+ /* compute the needed stres (Assuming not sort) */
+ stres = ((1 << gs.prec)-1 + res-2)/(res-1);
+
+#ifdef VERBOSE
+ printf("gres = %d, sres = %d, gs.gres = %d, gs.sres = %d\n",res,stres,gs.itres,gs.stres);
+#endif
+ /* We want one that is equals or exceeds the desired */
+ /* resolution, but doesn't exceed it too much, or the code */
+ /* will be inefficient. */
+ /* If there are no routines that can meet the desired precision, */
+ /* then it is ok to use the one closest to the desired precision. */
+ if (gs.itres >= res && gs.stres >= stres) {
+ fig += 10 * (gs.itres - res) + (gs.stres - stres);
+ } {
+ if (gs.itres < res) {
+ fig += 10000 + 10 * (res - gs.itres);
+ }
+ if (gs.stres < stres) {
+ fig += 1000 + 10 * (stres - gs.stres);
+ }
+ }
+ }
+
+#ifdef VERBOSE
+ printf(" figure of merit %d\n",fig);
+#endif
+ /* Is this the best one so far ? */
+ if (fig < bfig) {
+ bfig = fig;
+ bk = i;
+ bstres = stres;
+ bgs = gs; /* Structure copy */
+ bts = ts; /* Structure copy */
+ bcnv = cnv;
+#ifdef VERBOSE
+ printf(" best so far\n");
+#endif
+ }
+ }
+
+ if (bk < 0) {
+#ifdef VERBOSE
+ printf("new_imdi failed - dimensionality or representations couldn't be matched\n");
+#endif
+ return NULL; /* Nothing matches */
+ }
+
+ if ((im = (imdi *)calloc(1, sizeof(imdi))) == NULL) {
+#ifdef VERBOSE
+ printf("new_imdi malloc imdi failed\n");
+#endif
+ /* Should we return an error somehow ? */
+ return NULL;
+ }
+
+ /* We've decided kernel function bk is going to be the best, */
+ /* so now setup the appropriate tables to use with it. */
+ if (bgs.itres > res)
+ bgs.itres = res; /* Tell table create what the res is */
+ if (bgs.stres > bstres)
+ bgs.stres = bstres;
+
+ /* Tel table setup how to treat integer input in per channel lookups */
+ bgs.in_signed = in_signed;
+ bgs.out_signed = out_signed;
+
+#ifdef VERBOSE
+ if (!bts.sort) {
+ if ((bgs.stres * (bgs.itres-1)) < ((1 << bgs.prec)-1)) {
+ printf("Warning: table chosen doesn't reach desired precision!\n");
+ printf("Wanted %d, got %d\n", ((1 << bgs.prec)-1), (bgs.stres * (bgs.itres-1)));
+ }
+ }
+#endif
+
+ /* Allocate and initialise the appropriate tables */
+ im->impl = (void *)imdi_tab(&bgs, &bts, bcnv, in, out, ktable[bk].interp,
+ inm, outm, oopt, checkv, input_curves, md_table,
+ output_curves, cntx);
+
+ if (im->impl == NULL) {
+#ifdef VERBOSE
+ printf("imdi_tab failed\n");
+#endif
+ imdi_del(im);
+ return NULL;
+ }
+
+#ifdef VERBOSE
+ if (bcnv != conv_none)
+ printf("imdi_tab: using a runtime match, cnv flags 0x%x\n",bcnv);
+#endif
+
+ if (bcnv == conv_none) /* No runtime match conversion needed */
+ im->interp = ktable[bk].interp;
+ else
+ im->interp = interp_match;
+ im->get_check = imdi_get_check;
+ im->reset_check = imdi_reset_check;
+ im->info = imdi_info;
+ im->del = imdi_del;
+
+#ifdef VERBOSE
+ printf("new_imdi returning 0x%x\n", im);
+#endif
+ return im;
+}
+
+/* Runtime matching adapter */
+/* We can emulate many combinations of interpolation function */
+/* by converting to a routine that supports stride, or one that */
+/* supports plane interleaved and stride. */
+/* This also lets us emulate functions that can convert between */
+/* pixel and plane interleaved at the same time as doing the */
+/* color interpolation. Warp is in components (NOT bytes) */
+/* We cope with skipping writes to output channels by */
+/* only implementing that support in plane interleaved, */
+/* and depending on that being used for pixel interleaved. */
+static void interp_match(
+imdi *s,
+void **outp, int outst, /* Output pointers and stride */
+void **inp, int inst, /* Input pointers and stride */
+unsigned int npixels /* Number of pixels */
+) {
+ int j, i;
+ void *minp[IXDI];
+ void *moutp[IXDI];
+ imdi_imp *impl = (imdi_imp *)s->impl;
+ imdi_conv cnv = impl->cnv;
+
+ /* Need to supply default stride. */
+ /* This is simply the input dimension for pixel interleaved, */
+ /* and 1 for plane interleaved. */
+ if (cnv & conv_istr) {
+ if (impl->cirep == pixint8 || impl->cirep == pixint16)
+ inst = impl->id;
+ else
+ inst = 1;
+ }
+ if (cnv & conv_ostr) {
+ if (impl->corep == pixint8 || impl->corep == pixint16)
+ outst = impl->wod;
+ else
+ outst = 1;
+ }
+
+ /* Convert from pixel to plane by setting the plane pointers */
+ /* to each pixel component, and using the pixel interleaved stride */
+ if (cnv & conv_irep) {
+ if (impl->cirep == pixint8) { /* Convert from pix8 to plane8 */
+ for (j = 0; j < impl->id; j++)
+ minp[j] = (void *)((char *)inp[0] + j);
+ } else if (impl->cirep == pixint16) { /* Convert from pix16 to plane16 */
+ for (j = 0; j < impl->id; j++)
+ minp[j] = (void *)((char *)inp[0] + 2 * j);
+ }
+ } else { /* Copy pointers, because we call with them. */
+ if (impl->cirep == pixint8 || impl->cirep == pixint16) {
+ minp[0] = inp[0];
+ } else {
+ for (j = 0; j < impl->id; j++)
+ minp[j] = inp[j];
+ }
+ }
+ if (cnv & conv_orep) {
+ if (impl->corep == pixint8) { /* Convert from pix8 to plane8 */
+ for (j = i = 0; j < impl->od; j++) {
+ if ((impl->skipf & (1 << j)) == 0) /* If not skipped */
+ moutp[j] = (void *)((char *)outp[0] + i++);
+ else
+ moutp[j] = NULL;
+ }
+ } else if (impl->corep == pixint16) { /* Convert from pix16 to plane16 */
+ for (j = i = 0; j < impl->od; j++)
+ if ((impl->skipf & (1 << j)) == 0) /* If not skipped */
+ moutp[j] = (void *)((char *)outp[0] + 2 * i++);
+ else
+ moutp[j] = NULL;
+ }
+ } else { /* Copy pointers, because we call with them. */
+ if (impl->corep == pixint8 || impl->corep == pixint16) {
+ moutp[0] = outp[0];
+ } else { /* Plane interleaved */
+ for (j = i = 0; j < impl->od; j++) {
+ if ((impl->skipf & (1 << j)) == 0) /* If not skipped */
+ moutp[j] = outp[i++];
+ else
+ moutp[j] = NULL;
+ }
+ }
+ }
+
+ /* Convert fwd function to a backwards one, or visa-versa. */
+ /* Move the pointers to the last pixel, and reverse the stride */
+ if (cnv & conv_rev) {
+ if (impl->firep == pixint8) {
+ minp[0] = (void *)((char *)minp[0] + inst * (npixels - 1));
+ } else if (impl->firep == pixint16) {
+ minp[0] = (void *)((char *)minp[0] + inst * 2 * (npixels - 1));
+ } else if (impl->firep == planeint8) {
+ for (j = 0; j < impl->id; j++)
+ minp[j] = (void *)((char *)minp[j] + inst * (npixels - 1));
+ } else if (impl->firep == planeint16) {
+ for (j = 0; j < impl->id; j++)
+ minp[j] = (void *)((char *)minp[j] + inst * 2 * (npixels - 1));
+ }
+ inst = -inst;
+ if (impl->forep == pixint8) {
+ moutp[0] = (void *)((char *)moutp[0] + outst * (npixels - 1));
+ } else if (impl->forep == pixint16) {
+ moutp[0] = (void *)((char *)moutp[0] + outst * 2 * (npixels - 1));
+ } else if (impl->forep == planeint8) {
+ for (j = 0; j < impl->od; j++)
+ moutp[j] = (void *)((char *)moutp[j] + outst * (npixels - 1));
+ } else if (impl->forep == planeint16) {
+ for (j = 0; j < impl->od; j++)
+ moutp[j] = (void *)((char *)moutp[j] + outst * 2 * (npixels - 1));
+ }
+ outst = -outst;
+ }
+
+ /* Call the conversion routine we have */
+ impl->interp(s, moutp, outst, minp, inst, npixels);
+}
+
+/* Get the per output channel check flags - bit corresponds to output interpolation channel */
+static unsigned int imdi_get_check(imdi *im) {
+ imdi_imp *impl = (imdi_imp *)im->impl;
+ unsigned int rv;
+ int i;
+
+ /* Convert from output channel index to callback index */
+ for (rv = 0, i = 0; i < impl->od; i++) {
+ unsigned int b;
+ b = 1 & (impl->checkf >> i); /* Output index flag */
+ rv |= (b << impl->im_map[i]); /* to callback index */
+ }
+ return rv;
+}
+
+/* Reset the output check flags (not reset by interp) */
+static void imdi_reset_check(imdi *im) {
+ imdi_imp *impl = (imdi_imp *)im->impl;
+
+ impl->checkf = 0;
+}
+
+/* Return some information about the imdi */
+static void imdi_info(imdi *im, unsigned long *psize, int *pgres, int *psres) {
+ imdi_imp *impl = (imdi_imp *)im->impl;
+ unsigned long size;
+
+ size = sizeof(imdi);
+
+ if (impl != NULL) {
+ size += impl->size;
+ if (psize != NULL)
+ *psize = size;
+
+ if (pgres != NULL)
+ *pgres = impl->gres;
+
+ if (psres != NULL)
+ *psres = impl->sres;
+ }
+}
+
+
+/* Delete the object */
+static void imdi_del(imdi *im) {
+ /* Free all the allocated tables */
+ if (im->impl != NULL)
+ imdi_tab_free((imdi_imp *)im->impl);
+
+ /* Free this structure */
+ free(im);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-18 14:33:36 +0000