From 22f703cab05b7cd368f4de9e03991b7664dc5022 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B6rg=20Frings-F=C3=BCrst?= Date: Mon, 1 Sep 2014 13:56:46 +0200 Subject: Initial import of argyll version 1.5.1-8 --- imdi/imdi.c | 605 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 605 insertions(+) create mode 100644 imdi/imdi.c (limited to 'imdi/imdi.c') 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 +#include +#include +#include +#include + +#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); +} + + + + + + + + + + + + + + + + + + -- cgit v1.2.3