diff options
author | Jörg Frings-Fürst <debian@jff-webhosting.net> | 2014-09-01 13:56:46 +0200 |
---|---|---|
committer | Jörg Frings-Fürst <debian@jff-webhosting.net> | 2014-09-01 13:56:46 +0200 |
commit | 22f703cab05b7cd368f4de9e03991b7664dc5022 (patch) | |
tree | 6f4d50beaa42328e24b1c6b56b6ec059e4ef21a5 /rspl/gam.c |
Initial import of argyll version 1.5.1-8debian/1.5.1-8
Diffstat (limited to 'rspl/gam.c')
-rw-r--r-- | rspl/gam.c | 1220 |
1 files changed, 1220 insertions, 0 deletions
diff --git a/rspl/gam.c b/rspl/gam.c new file mode 100644 index 0000000..f8d02d0 --- /dev/null +++ b/rspl/gam.c @@ -0,0 +1,1220 @@ + +/* + * Argyll Color Correction System + * Multi-dimensional regularized spline data structure + * + * Precice gamut surface, gamut pruning, ink limiting and K min/max + * support routine.s + * + * Author: Graeme W. Gill + * Date: 2008/11/21 + * + * Copyright 1999 - 2008 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. + * + * Latest simplex/linear equation version. + */ + +/* TTBD: + + Add ouutput curve lookup callback support. + Cache these values in the vertex structures ? + + Add ink limit support. This be done by breaking + a cell into a fixed geometry of smaller simplexes + by dividing the cell into two on each axis. + + Need to then add scan that detects areas to prune, + that then ties in with rev code to mark such + areas out of gamut. + + */ + +#include <stdio.h> +#include <stdlib.h> +#include <stdarg.h> +#include <math.h> +#include <memory.h> +#include <time.h> + +#include "rspl_imp.h" +#include "numlib.h" +#include "sort.h" /* Heap sort */ +#include "counters.h" /* Counter macros */ + +#include "vrml.h" /* If there is VRML stuff here */ + +/* Print a vectors value */ +#define DBGVI(text, dim, out, vec, end) \ +{ int pveci; \ + printf("%s",text); \ + for (pveci = 0 ; pveci < (dim); pveci++) \ + printf(out,(vec)[pveci]); \ + printf(end); \ +} + +/* Print a matrix value */ +#define DBGMI(text, rows, cols, out, mat, end) \ +{ int pveci, pvecr; \ + printf("%s",text); \ + for (pvecr = 0 ; pvecr < (rows); pvecr++) { \ + for (pveci = 0 ; pveci < (cols); pveci++) \ + printf(out,(mat)[pvecr][pveci]); \ + if ((pvecr+1) < (rows)) \ + printf("\n"); \ + } \ + printf(end); \ +} + +#undef VRML_TRACE /* Save a vrml at each step */ + +/* Do an arbitrary printf */ +#define DBGI(text) printf text ; + +#define DEBUG1 +#undef DBG +#undef DBGV +#undef DBGM + +#undef NEVER +#define ALWAYS + +#ifdef DEBUG1 +#undef DBGS +#undef DBG +#undef DBGV +#undef DBGM +#define DEBUG +#define DBGS(xxx) xxx +#define DBG(xxx) DBGI(xxx) +#define DBGV(xxx) DBGVI xxx +#define DBGM(xxx) DBGMI xxx +#else +#undef DEBUG +#undef DBGS +#undef DBG +#undef DBGV +#undef DBGM +#define DBGS(xxx) +#define DBG(xxx) +#define DBGV(xxx) +#define DBGM(xxx) +#endif + +/* Convention is to use: + i to index grid points u.a + n to index data points d.a + e to index position dimension di + f to index output function dimension fdi + j misc and cube corners + k misc + */ + +#define EPS (1e-10) /* Allowance for numeric error */ + +/* ====================================================== */ +/* Support functions */ + +/* Compute the norm (length) squared of a vector define by two points */ +static double norm33sq(double in1[3], double in0[3]) { + int j; + double rv; + for (rv = 0.0, j = 0; j < 3; j++) { + double tt = in1[j] - in0[j]; + rv += tt * tt; + } + return rv; +} + + +static rvert *get_vert(rspl *s, int gix); + +/* Given an output value, return the gamut radius */ +static double gvprad(rspl *s, double *v) { + int f, fdi = s->fdi; + double rr = 0.0; + + for (f = 0; f < fdi; f++) { + double tt; + tt = s->gam.scale[f] * (v[f] - s->gam.cent[f]); + rr += tt * tt; + } + rr = sqrt(rr); + return rr; +} + + +/* Given an output value, create the radial coordinate */ +static void radcoord(rspl *s, double *rad, double *v) { + int f, fdi = s->fdi; + double rr = 0.0; + + for (f = 0; f < fdi; f++) { + double tt; + tt = s->gam.scale[f] * (v[f] - s->gam.cent[f]); + rr += tt * tt; + } + rr = sqrt(rr); + + rad[0] = rr; +} + +/* Given an output value and an edge, return the side, */ +/* 0 = -vem 1 = +ve */ +static int eside(rspl *s, redge *ep, double *v) { + int f, fdi = s->fdi; + double tt; + + for (tt = 0.0, f = 0; f < fdi; f++) { + tt += v[f] * ep->pe[f]; + } + tt += ep->pe[f]; + return (tt >= 0.0 ? 1 : 0); +} + + +/* Given a list of nodes that form an sdi-1 sub-simplex, */ +/* return a list of nodes that could be added to make */ +/* an sdi sub-simplex. */ +/* Nodes are identified by their grid index. */ +/* All possible nodes are returned, even if they are already */ +/* in our surface triangulation. */ +/* NOTE that inodes will be sorted into sub-simplex order! */ +// ~~99 we aren't currently taking ink limit into account + +/* return nz on fatal error */ +static int get_ssimplex_nodes( +rspl *s, +int sdi, /* Dimensionality of target sub-simplex */ +rvert **inodes, /* sdi input nodes that form an sdi-1 dimensional input sub-simplex */ +int aonodes, /* Number of output nodes allowed for */ +int *nonodes, /* Number of output nodes set */ +rvert **onodes /* Space for up to 3^MXDI-1 output nodes */ +) { + int e, di = s->di; + int i, j, k; + + *nonodes = 0; + +//printf("\n~1 get_ssimplex_nodes called with sdi = %d\n",sdi); + /* Sort the input nodes into normal sub-simplex order */ + /* (We are assuming all the nodes are in the grid) */ + for (i = 0; i < sdi-1; i++) { + for (j = i+1; j < sdi; j++) { + if (inodes[i]->gix < inodes[j]->gix) { + rvert *tt; + tt = inodes[i]; inodes[i] = inodes[j]; inodes[j] = tt; + } + } + } +//printf("~2 input nodes = %d %d\n", inodes[0]-gix, inodes[1]->gix); + + /* For each sub-simplex */ + for (i = 0; i < s->gam.ssi[sdi].nospx; i++) { + int kk; +//printf("~1 sub simplex %d\n",i); + /* For leaving out one node of the ssimplex in turn, */ + /* check the inodes match the remaining ssimplex nodes */ + for (kk = 0; kk <= sdi; kk++) { /* ssimplex node being left out */ + int bix = 0; /* ssimplex base node */ + + if (kk == 0) + bix++; + +//printf("~1 anchor node %d = %d\n",j,s->gam.ssi[sdi].spxi[i].goffs[bix]); +//printf("~2 candidate node offsets = %d %d %d\n", +//s->gam.ssi[sdi].spxi[i].offs[0], +//s->gam.ssi[sdi].spxi[i].offs[1], +//s->gam.ssi[sdi].spxi[i].offs[2]); + + for (j = k = 0; j < sdi; j++, k++) { /* Check all inodes[] */ + if (k == kk) + k++; /* Skip the ssimplex node being left out */ + if (inodes[j]->gix != (inodes[0]->gix + s->gam.ssi[sdi].spxi[i].goffs[k] + - s->gam.ssi[sdi].spxi[i].goffs[bix])) { +//printf("~1 not a match\n"); + break; + } +//printf("~1 matched node %d\n",inodes[k]->gix); + } + if (j >= sdi) { /* they all match */ +//printf("~1 all match\n"); + /* Check if kk offset to the base is still in the grid */ + for (e = 0; e < di; e++) { + int doff = ((s->gam.ssi[sdi].spxi[i].offs[kk] >> e) & 1) + - ((s->gam.ssi[sdi].spxi[i].offs[bix] >> e) & 1); + int eflags = G_FL(inodes[0]->fg,e); +//printf("~1 checking dim %d, doff = %d, eflags = %d\n",e,doff,eflags); + if ((doff < 0 && (eflags & 4) && (eflags & 3) < 1) + || (doff > 0 && !(eflags & 4) && (eflags & 3) < 1)) { +//printf("~1 outside grid\n"); + break; /* Offset will take us past edge */ + } + } + if (e >= di) { /* Remaining point is within grid */ + int gix; + + /* Add the remaining node to the onodes */ + if (*nonodes >= aonodes) + return 1; /* Oops - ran out of return space! */ + + gix = inodes[0]->gix + s->gam.ssi[sdi].spxi[i].goffs[kk] + - s->gam.ssi[sdi].spxi[i].goffs[bix]; + onodes[*nonodes] = get_vert(s, gix); + +//printf("~1 Returning node %d\n",(onodes[*nonodes]->gix); +//printf("~1 Value %f %f %f\n", onodes[*nonodes]->p[0], onodes[*nonodes]->p[1], onodes[*nonodes]->p[2]); + (*nonodes)++; + } + } + } + } +//printf("~1 onodes = %d\n",*nonodes); + return 0; +} + +/* ====================================================== */ +/* Search for an existing vertex, and return it. */ +/* If there is no existing edge, create it. */ +/* (This is where we apply the output functions to the grid values) */ +static rvert *get_vert(rspl *s, int gix) { + int f, fdi = s->fdi; + int hash; + rvert *vp = NULL; + + if (gix < 0 || gix >= s->g.no) { /* Assert */ + error("rspl_gam: get_vert got out of range gix %d\n",gix); + } + + /* See if it is in our hash list */ + hash = gix % s->gam.vhsize; + + for (vp = s->gam.verts[hash]; vp != NULL; vp = vp->next) { + if (vp->gix == gix) + break; + } + if (vp == NULL) { /* No such vertex */ + float *fg; + + if ((vp = calloc(1, sizeof(rvert))) == NULL) + error("rspl_gam: get_vert calloc failed"); + + fg = s->g.a + gix * s->g.pss; + vp->n = s->gam.rvert_no++; /* serial number */ + vp->fg = fg; /* Pointer to node float data */ + vp->gix = gix; /* grid index */ + for (f = 0; f < fdi; f++) /* Node output value */ + vp->v[f] = (double)fg[f]; + if (s->gam.outf != NULL) + s->gam.outf(s->gam.cntxf, vp->v, vp->v); /* Apply output lookup */ + + radcoord(s, vp->r, vp->v); /* Compute radial coordinate */ + + /* Add vertex to hash */ + vp->next = s->gam.verts[hash]; + s->gam.verts[hash] = vp; + + /* Add the vertex to the bottom of the list of vertex */ + if (s->gam.vbot == NULL) { + s->gam.vtop = s->gam.vbot = vp; + } else { + s->gam.vbot->list = vp; + s->gam.vbot = vp; + } + } + + return vp; +} + +/* Search for an existing edge containing the given verticies, */ +/* and return it. If there is no existing edge, create it. */ +/* Note that edges have fdi-2 dimensions == fdi-1 verticies. */ +static redge *get_edge(rspl *s, rvert **_vv) { + int i, j, f, fdi = s->fdi; + rvert *vv[MXDO-1]; /* Sorted verticies */ + int gix, hash; + redge *ep = NULL; + + /* Sort the input nodes into normal sub-simplex order */ + /* (We are assuming all the nodes are in the grid) */ + for (i = 0; i < (fdi-1); i++) + vv[i] = _vv[i]; + for (i = 0; i < (fdi-2); i++) { + for (j = i+1; j < (fdi-1); j++) { + if (vv[i]->gix < vv[j]->gix) { + rvert *tt; + tt = vv[i]; vv[i] = vv[j]; vv[j] = tt; + } + } + } + for (gix = i = 0; i < (fdi-1); i++) + gix += vv[i]->gix; +//printf("~1 get edge with nodes = %d %d\n", vv[0]->gix, vv[1]->gix); + + /* See if it is in our hash list */ + hash = gix % s->gam.ehsize; +//printf("~1 get edge gix = %d, hash = %d\n", gix, hash); + + for (ep = s->gam.edges[hash]; ep != NULL; ep = ep->next) { + for (i = 0; i < (fdi-1); i++) { + if (ep->v[i] != vv[i]) { +//printf("~1 verticies don't match\n"); + break; /* No match */ + } + } + if (i >= (fdi-1)) { +//printf("~1 all verticies match\n"); + break; /* All match */ + } + } + if (ep == NULL) { /* No such edge */ + int sm, lg; + + if ((ep = calloc(1, sizeof(redge))) == NULL) + error("rspl_gam: get_edge calloc failed"); + + ep->n = s->gam.redge_no++; /* serial number */ + for (i = 0; i < (fdi-1); i++) + ep->v[i] = vv[i]; /* Vertex */ +printf("~1 new edge %d with nodes = %d %d\n", ep->n, ep->v[0]->gix, ep->v[1]->gix); + + /* Compute plane equation to center of gamut, so */ + /* that we can quickly determine which side a triangle lies on. */ + + /* Compute plane equation */ + if (fdi < 2 || fdi > 3) + error("rspl_gam: plane equation for out dimensions other than 2 or 3 not supported!"); + if (fdi == 2) { + + } else { + double v1[3], v2[3]; + + /* Compute two vectors from the three points */ + for (f = 0; f < fdi; f++) { + v1[f] = ep->v[0]->v[f] - s->gam.cent[f]; + v2[f] = ep->v[1]->v[f] - s->gam.cent[f]; + } + + /* Compute normal to the plane using the cross product */ + ep->pe[0] = ep->v[0]->v[1] * (ep->v[1]->v[2] - s->gam.cent[2]) + + ep->v[1]->v[1] * (s->gam.cent[2] - ep->v[0]->v[2]) + + s->gam.cent[1] * (ep->v[0]->v[2] - ep->v[1]->v[2]); + ep->pe[1] = ep->v[0]->v[2] * (ep->v[1]->v[0] - s->gam.cent[0]) + + ep->v[1]->v[2] * (s->gam.cent[0] - ep->v[0]->v[0]) + + s->gam.cent[2] * (ep->v[0]->v[0] - ep->v[1]->v[0]); + ep->pe[2] = ep->v[0]->v[0] * (ep->v[1]->v[1] - s->gam.cent[1]) + + ep->v[1]->v[0] * (s->gam.cent[1] - ep->v[0]->v[1]) + + s->gam.cent[0] * (ep->v[0]->v[1] - ep->v[1]->v[1]); + ep->pe[3] = - (ep->v[0]->v[0] * (ep->v[1]->v[1] * s->gam.cent[2] - s->gam.cent[1] * ep->v[1]->v[2]) + + ep->v[1]->v[0] * (s->gam.cent[1] * ep->v[0]->v[2] - ep->v[0]->v[1] * s->gam.cent[2]) + + s->gam.cent[0] * (ep->v[0]->v[1] * ep->v[1]->v[2] - ep->v[1]->v[1] * ep->v[0]->v[2])); + + } + + /* Add edge to hash */ + ep->next = s->gam.edges[hash]; + s->gam.edges[hash] = ep; + + /* Add the edge to the bottom of the list of edges */ + if (s->gam.ebot == NULL) { + s->gam.etop = s->gam.ebot = ep; + } else { + s->gam.ebot->list = ep; + s->gam.ebot = ep; + } + } + +printf("~1 returning edge no %d\n",ep->n); + return ep; +} + +/* Check whether a triangle like this already exists */ +/* Return NULL if it doesn't */ +static rtri *check_tri(rspl *s, redge *ep, rvert *_vv) { + int i, j, k, f, fdi = s->fdi; + rvert *vv[MXDO]; /* Sorted verticies */ + int gix, hash; + rtri *tp = NULL; + + /* Copy edge verticies from edge */ + for (i = 0; i < (fdi-1); i++) + vv[i] = ep->v[i]; + vv[i] = _vv; /* And new vertex */ + + /* Sort verticies */ + for (i = 0; i < (fdi-1); i++) { + for (j = i+1; j < fdi; j++) { + if (vv[i]->gix < vv[j]->gix) { + rvert *tv; + tv = vv[i]; vv[i] = vv[j]; vv[j] = tv; + } + } + } + + /* Create hash */ + for (gix = i = 0; i < fdi; i++) + gix += vv[i]->gix; + + /* See if it is in our hash list */ + hash = gix % s->gam.thsize; +//printf("~1 make tri gix = %d, hash = %d\n", gix, hash); + + for (tp = s->gam.tris[hash]; tp != NULL; tp = tp->next) { + for (i = 0; i < fdi; i++) { + if (tp->v[i] != vv[i]) { +//printf("~1 verticies don't match\n"); + break; /* No match */ + } + } + if (i >= fdi) { +//printf("~1 all verticies match\n"); + break; /* All match */ + } + } + + return tp; +} + +/* Create a triangle given an edge with fdi-1 verticies and a grid vertex. */ +/* if inv is set, triangle is upside down wrt to center point, */ +/* and so all the oppositive verticies should be placed on the */ +/* opposite to their natural side. */ +static rtri *make_tri(rspl *s, redge *ep, rvert *_vv, int inv) { + int i, j, k, f, fdi = s->fdi; + rvert *vv[MXDO]; /* Sorted verticies */ + int gix, hash; + rtri *tp = NULL; + +printf("~1 make_tri called\n"); + + /* Copy edge verticies from edge */ + for (i = 0; i < (fdi-1); i++) + vv[i] = ep->v[i]; + vv[i] = _vv; /* And new vertex */ + + /* Sort verticies */ + for (i = 0; i < (fdi-1); i++) { + for (j = i+1; j < fdi; j++) { + if (vv[i]->gix < vv[j]->gix) { + rvert *tv; + tv = vv[i]; vv[i] = vv[j]; vv[j] = tv; + } + } + } + + /* Create hash */ + for (gix = i = 0; i < fdi; i++) + gix += vv[i]->gix; + + /* See if it is in our hash list */ + hash = gix % s->gam.thsize; +printf("~1 make tri gix = %d, hash = %d\n", gix, hash); + + for (tp = s->gam.tris[hash]; tp != NULL; tp = tp->next) { + for (i = 0; i < fdi; i++) { + if (tp->v[i] != vv[i]) { +//printf("~1 verticies don't match\n"); + break; /* No match */ + } + } + if (i >= fdi) { +//printf("~1 all verticies match\n"); + break; /* All match */ + } + } + + if (tp == NULL) { /* No such triangle */ + +//printf("~1 creating a new triangle\n"); + /* Create triangle */ + if ((tp = calloc(1, sizeof(rtri))) == NULL) + error("rspl_gam: make_tri calloc failed"); + + tp->n = s->gam.rtri_no++; /* serial number */ + + /* Copy edge verticies to triangle */ + for (i = 0; i < fdi; i++) + tp->v[i] = vv[i]; + + /* Link all the triangles edges to this triangle */ + printf("~1 triangle nodes = %d %d %d\n", tp->v[0]->gix, tp->v[1]->gix, tp->v[2]->gix); +//printf("~2 triangle vert 0 = %f %f %f\n", tp->v[0]->v[0], tp->f[0]->v[1], tp->f[0]->v[2]); +//printf("~2 triangle vert 1 = %f %f %f\n", tp->v[1]->v[0], tp->f[1]->v[1], tp->f[1]->v[2]); +//printf("~2 triangle vert 2 = %f %f %f\n", tp->v[2]->v[0], tp->f[2]->v[1], tp->f[2]->v[2]); + for (i = 0; i < fdi; i++) { /* For each tri verticy being odd one out */ + rvert *ov, *rr[MXDO-1]; /* Odd verticy, remaining edge verticies */ + int ss; /* Side */ + + ov = tp->v[i]; /* Odd node */ + for (k = j = 0; j < fdi; j++) { + if (i == j) + continue; + rr[k++] = tp->v[j]; /* Remaining nodes */ + } +//printf("~1 edge nodes = %d %d, odd = %d\n", rr[0]->gix, rr[1]->gix, ov->gix); + ep = get_edge(s, rr); + + if (ep->nt >= MXNE) + error("rspl_gam: make_tri run out of triangle space %d in edge",MXNE); + ep->t[ep->nt++] = tp; + + /* See which side of the edge the remaining vertex is */ + ss = eside(s, ep, ov->v); +//printf("~1 node gix %d has side %d to edge %d %d\n", ov->gix, ss, ep->v[0]->gix, ep->v[1]->gix); + if (inv) + ss = 1-ss; + if (ss) { /* +ve side */ + ep->t[ep->npt++] = tp; + } else { + ep->t[ep->nnt++] = tp; + } + } + + /* Add triangle to hash */ + tp->next = s->gam.tris[hash]; + s->gam.tris[hash] = tp; + + /* Add them to the linked list */ + tp->list = s->gam.ttop; + s->gam.ttop = tp; + + } + + return tp; +} + + +/* ====================================================== */ + +/* Create a surface gamut representation. */ +/* Return NZ on error */ +/* Could add more flexibility with: + optional function instead of default radial distance function + option use sub set of output dimensions (ie allow for CMYK->LabK etc. ? ) +*/ +static int +gam_comp_gamut( +rspl *s, +double *cent, /* Optional center of gamut [fdi], default center of out range */ +double *scale, /* Optional Scale of output values in vector to center [fdi], def. 1.0 */ +void (*outf)(void *cntxf, double *out, double *in), /* Optional rspl val -> output value */ +void *cntxf, /* Context for function */ +void (*outb)(void *cntxb, double *out, double *in), /* Optional output value -> rspl val */ +void *cntxb /* Context for function */ +) { + int e, f, di = s->di, fdi = s->fdi; + int i, j, ssdi; + int maxp[MXDO]; /* Grid indexes of maxium function values */ + rvert *fedge[MXDO-1]; /* The first "edge" containing fdi-1 verticies */ + rvert *onodes[50]; /* float pointers of canditate nodes */ + int aonodes = 50; /* Allocated out nodes */ + int nonodes; /* Number of nodes returned */ + + rvert *cnodes[2][50]; /* Negative, positive candidate nodes */ + double rcnodes[2][50]; /* Radius of negative, positive candidate nodes */ + int ncnodes[2]; /* Number of negative, positive candidate nodes */ + + if (fdi < 2 || fdi > di) { + DBG(("gam: gam_comp_gamut called for di = %d, fdi = %d\n", di, fdi)); + return 2; + } + + /* Save output value conversion functions */ + s->gam.outf = outf; + s->gam.cntxf = cntxf; + s->gam.outb = outb; + s->gam.cntxb = cntxb; + + /* Deal with gamut center point, and scale */ + if (cent == NULL) { + double min[MXDO], max[MXDO]; + + s->get_out_range(s, min, max); + if (s->gam.outf != NULL) { + s->gam.outf(s->gam.cntxf, min, min); /* Apply output lookup */ + s->gam.outf(s->gam.cntxf, max, max); /* Apply output lookup */ + } + for (f = 0; f < fdi; f++) + s->gam.cent[f] = 0.5 * (min[f] + max[f]); + } else { + for (f = 0; f < fdi; f++) + s->gam.cent[f] = cent[f]; + } + DBGVI("Gamut center is ", fdi, "%f ", s->gam.cent, "\n") + + if (scale == NULL) { + for (f = 0; f < fdi; f++) + s->gam.scale[f] = 1.0; + } else { + for (f = 0; f < fdi; f++) + s->gam.scale[f] = scale[f]; + } + DBGVI("Gamut scale is ", fdi, "%f ", s->gam.scale, "\n") + + for (ssdi = 1; ssdi <= fdi-1; ssdi++) { + int i, j; + + /* Compute gamut surface sub-simplex geometry info */ + rspl_init_ssimplex_info(s, &s->gam.ssi[ssdi], ssdi); + + /* Filter the sub-simplex geometry to only include subsimplexes that */ + /* use the base vertex, so that there are no duplicates. */ + for (i = j = 0; i < s->gam.ssi[ssdi].nospx; i++) { + if (s->gam.ssi[ssdi].spxi[i].offs[s->gam.ssi[ssdi].sdi] == 0) { + s->gam.ssi[ssdi].spxi[j] = s->gam.ssi[ssdi].spxi[i]; /* Structure copy */ + j++; + } + } + s->gam.ssi[ssdi].nospx = j; + +#ifdef DEBUG + printf("Sub-simplex dim %d out of input %d\n",s->gam.ssi[ssdi].sdi,di); + printf("Number of subsimplex = %d\n",s->gam.ssi[ssdi].nospx); + for (i = 0; i < s->gam.ssi[ssdi].nospx; i++) { + printf("Cube Offset = "); + for (e = 0; e <= s->gam.ssi[ssdi].sdi; e++) + printf("%d ",s->gam.ssi[ssdi].spxi[i].offs[e]); + printf("\n"); + + printf("Grid Offset = "); + for (e = 0; e <= s->gam.ssi[ssdi].sdi; e++) + printf("%d ",s->gam.ssi[ssdi].spxi[i].foffs[e]); + printf("\n"); + printf("\n"); + } +#endif /* DEBUG */ + } + + /* Allocate the vertex hash array */ + if ((s->gam.verts = calloc(VHASHSIZE, sizeof(rvert *))) == NULL) { + DBG(("gam: allocating vertex hash array failed\n")); + return 1; + } + s->gam.vhsize = VHASHSIZE; + + /* Allocate the edge hash array */ + if ((s->gam.edges = calloc(EHASHSIZE, sizeof(redge *))) == NULL) { + DBG(("gam: allocating edge hash array failed\n")); + return 1; + } + s->gam.ehsize = EHASHSIZE; + + /* Allocate the triangle hash array */ + if ((s->gam.tris = calloc(THASHSIZE, sizeof(rtri *))) == NULL) { + DBG(("gam: allocating tris hash array failed\n")); + return 1; + } + s->gam.thsize = THASHSIZE; + + /* Get a starting gid point for surface */ + // ~~99 this isn't right if the point we get is over the ink limit!!!. */ + s->get_out_range_points(s, NULL, maxp); /* Maximum */ +// s->get_out_range_points(s, maxp, NULL); /* Minium */ + + DBG(("Starting point = gix %d/%d\n",maxp[0], s->g.no)); + + /* Now work it up to an fdi-2 sub-simplex, so that we have a */ + /* "triangle edge", and can enter the main loop. */ + // ~~~9 this needs fixing to switch to "maximum angle" calculation + fedge[0] = get_vert(s, maxp[0]); /* grid index of first vertex */ + + if (fdi == 2) { + /* We just need one point, and we've got it */ + + } else if (fdi == 3) { /* Usual case */ + double ba = -1.0; /* Bigest angle */ + /* We just need two points, so find the other points */ + /* that makes the greatest angle to the center */ + + if (get_ssimplex_nodes(s, 1, fedge, aonodes, &nonodes, onodes)) { + error("rspl_gam: get_ssimplex_nodes fatal error - too many nodes?"); + } + if (nonodes == 0) + error("rspl_gam: get_ssimplex_nodes fatal error - retrurned no nodes?"); + +printf("~1 get_ssimplex_nodes returned %d nodes\n",nonodes); +for(i = 0; i < nonodes; i++) +printf(" ~1 %d: %d\n",i,onodes[i]->gix); + + /* Evaluate the choice of verticies to choose the one that */ + /* will best enclose the gamut. (We use a really dumb criteria - maximum radius) */ + for (i = 0; i < nonodes; i++) { + double tt, a, b, c; /* Lenght of sides of triangle squared */ + a = norm33sq(onodes[i]->v, s->gam.cent); + b = norm33sq(onodes[i]->v, fedge[0]->v); + c = norm33sq(fedge[0]->v, s->gam.cent); + tt = acos((b + c - a)/(2.0 * sqrt(b * c))); +printf("~1 node %d angle = %f\n",onodes[i]->gix,tt); + if (tt > ba) { + ba = tt; + fedge[1] = onodes[i]; /* Use candidate with largest gamut as next base */ + } + } +printf("~1 chosen node %d\n",fedge[1]->gix); + + } else if (fdi > 3) { /* General case */ + /* This isn't a correct approach ... */ + + for (ssdi = 1; ssdi <= fdi-2; ssdi++) { + double br = -1.0; /* Bigest radius */ + + DBG(("Working up ssdim %d -> %d\n",ssdi-1, ssdi)); + + if (get_ssimplex_nodes(s, ssdi, fedge, aonodes, &nonodes, onodes)) { + error("rspl_gam: get_ssimplex_nodes fatal error - too many nodes?"); + } + if (nonodes == 0) + error("rspl_gam: get_ssimplex_nodes fatal error - retrurned no nodes?"); + +printf("~1 get_ssimplex_nodes returned %d nodes\n",nonodes); +for(i = 0; i < nonodes; i++) +printf(" ~1 %d: %d\n",i,onodes[i]->gix); + + /* Evaluate the choice of verticies to choose the one that */ + /* will best enclose the gamut. (We use a really dumb criteria - maximum radius) */ + for (i = 0; i < nonodes; i++) { + double tt; + tt = gvprad(s, onodes[i]->v); /* Compute radius */ + if (tt > br) { + br = tt; + fedge[ssdi] = onodes[i]; /* Use candidate with largest gamut as next base */ + } + } +printf("~1 chosen node %d\n",fedge[ssdi]->gix); + } + } + + /* Creat the initial "edge" */ + get_edge(s, fedge); + +printf("~1 Created initial edge\n"); + + { + redge *ep; + + double **A; /* lu value -> baricentric matrix */ + double *B; + int *pivx; + + pivx = ivector(0, fdi-1); /* pixv[fdi] */ + A = dmatrix(0, fdi-1, 0, fdi-1); /* A[fdi][fdi] */ + B = dvector(0, fdi-1); /* B[fdi] */ + + /* The main loop: */ + /* We start with any "edges" that have less than two associated "triangles", */ + /* and evaluate the vertex nodes that can make "triangles" with the "edge". */ + /* We choose the node that will give the greatest slope, and make a "triangle". */ + /* Loop until there are no "edges" with less than two associated "triangles". */ + for (ep = s->gam.etop; ep != NULL; ep = ep->list) { + +printf("~1 expanding from edge no %d\n",ep->n); +printf("~1 edge v1 = %d = %f %f %f\n", ep->v[0]->gix, ep->v[0]->v[0], ep->v[0]->v[1], ep->v[0]->v[2]); +printf("~1 edge v2 = %d = %f %f %f\n", ep->v[1]->gix, ep->v[1]->v[0], ep->v[1]->v[1], ep->v[1]->v[2]); + +// if (ep->npt < 1 || ep->nnt < 1) + { + int ss; + + if (get_ssimplex_nodes(s, fdi-1, ep->v, aonodes, &nonodes, onodes)) { + error("rspl_gam: get_ssimplex_nodes fatal error - too many nodes?"); + } + + /* Clasify the returned nodes as positive or negative side */ + ncnodes[0] = ncnodes[1] = 0; + for (i = 0; i < nonodes; i++) { + double tt; + tt = gvprad(s, onodes[i]->v); /* Compute radius */ + ss = eside(s, ep, onodes[i]->v); /* Side */ +printf("~1 node gix %d has rad %f side %d to edge %d %d\n", onodes[i]->gix, tt, ss, ep->v[0]->gix, ep->v[1]->gix); + +#ifdef NEVER /* This messes things up ? */ + /* Check if this node is already in a triangle with this edge, */ + /* to avoid costly matrix solution check below. */ + if (check_tri(s, ep, onodes[i]) == NULL) +#endif + { + cnodes[ss][ncnodes[ss]] = onodes[i]; + rcnodes[ss][ncnodes[ss]++] = tt; + } + } + + /* Sort them */ + for (ss = 0; ss < 2; ss++) { /* Negative side then positive */ + for (i = 0; i < (ncnodes[ss]-1); i++) { + for (j = i+1; j < ncnodes[ss]; j++) { + if (rcnodes[ss][i] < rcnodes[ss][j]) { + rvert *tt; + double tr; + tt = cnodes[ss][i]; cnodes[ss][i] = cnodes[ss][j]; + cnodes[ss][j] = tt; + tr = rcnodes[ss][i]; rcnodes[ss][i] = rcnodes[ss][j]; + rcnodes[ss][j] = tr; + } + } + } + } + +// ~~1 +for (ss = 0; ss < 2; ss++) { /* Negative side then positive */ + if (ss == 0) + printf("~1 -ve nodes:\n"); + else + printf("~1 +ve nodes:\n"); + for (i = 0; i < ncnodes[ss]; i++) { + printf("~1 node %d, rad %f\n",cnodes[ss][i]->gix,rcnodes[ss][i]); + } +} + +#ifdef VRML_TRACE +{ + vrml *wrl; + int i, j, k; + ECOUNT(gc, MXDIDO, di, 0, s->g.res, 0);/* coordinates */ + DCOUNT(cc, MXDIDO, s->di, 0, 0, 2); /* Surrounding cube counter */ + float *gp; /* Grid point pointer */ + rvert *vp; + rtri *tp; + double col[3], pos[3]; + + if ((wrl = new_vrml("gam_diag.wrl", 1)) == NULL) + error("new_vrml failed\n"); + + /* Display the grid */ + EC_INIT(gc); + for (gp = s->g.a; !EC_DONE(gc); gp += s->g.pss) { + /* Itterate cube from this base */ + DC_INIT(cc) + for (i = 0; !DC_DONE(cc); i++ ) { + float *sp = s->g.a; + + for (e = 0; e < s->di; e++) { /* Input tables */ + int j; + j = gc[e] + cc[e]; + if (j < 0 || j >= s->g.res[e]) { + sp = NULL; /* outside grid */ + break; + } + sp += s->g.fci[e] * j; /* Compute pointer to surrounder */ + } + if (i> 0 && e >= s->di) { + /* Create vector from base to surrounder */ + pos[0] = (double)gp[0]; + pos[1] = (double)gp[1]; + pos[2] = (double)gp[2]; + if (s->gam.outf != NULL) + s->gam.outf(s->gam.cntxf, pos, pos); /* Apply output lookup */ + wrl->add_vertex(wrl, pos); + pos[0] = (double)sp[0]; + pos[1] = (double)sp[1]; + pos[2] = (double)sp[2]; + if (s->gam.outf != NULL) + s->gam.outf(s->gam.cntxf, pos, pos); /* Apply output lookup */ + wrl->add_vertex(wrl, pos); + } + DC_INC(cc); + } + EC_INC(gc); + } + wrl->make_lines(wrl, 2); + wrl->clear(wrl); + + /* Display the current triangles transparently */ + if (s->gam.ttop != NULL) { + for (vp = s->gam.vtop; vp != NULL; vp = vp->list) { + wrl->add_vertex(wrl, vp->v); + } + + /* Set the triangles */ + for (tp = s->gam.ttop; tp != NULL; tp = tp->list) { + int ix[3]; + ix[0] = tp->v[0]->n; + ix[1] = tp->v[1]->n; + ix[2] = tp->v[2]->n; + wrl->add_triangle(wrl, ix); + } + +// wrl->make_triangles(wrl, 0.3, NULL); + wrl->make_triangles(wrl, 0.0, NULL); + wrl->clear(wrl); + } + + /* Show the active edge as a red cone */ + col[0] = 1.0; col[1] = 0.0; col[2] = 0.0; /* Red */ + wrl->add_cone(wrl, ep->v[0]->v, ep->v[1]->v, col, 1.0); +printf("~1 edge v1 = %d = %f %f %f\n", ep->v[0]->gix, ep->v[0]->v[0], ep->v[0]->v[1], ep->v[0]->v[2]); +printf("~1 edge v2 = %d = %f %f %f\n", ep->v[1]->gix, ep->v[1]->v[0], ep->v[1]->v[1], ep->v[1]->v[2]); + + /* Show the candidate verticies */ + for (ss = 0; ss < 2; ss++) { + for (i = 0; i < ncnodes[ss]; i++) { /* +ve */ + if (ss) { + col[0] = 0.0; col[1] = 1.0; col[2] = 0.0; /* Green +ve */ + } else { + col[0] = 0.0; col[1] = 0.0; col[2] = 1.0; /* Blue -ve */ + } + wrl->add_marker(wrl, cnodes[ss][i]->v, col, 1.0); + } + } + wrl->del(wrl); + getchar(); +} +#endif + /* See if we can make a triangle */ + for (ss = 0; ss < 2; ss++) { /* For -ve and +ve sides */ + int ii, si; + int sta, end, inc; + rvert **nods; + double rip; /* Row interchange parity */ + +printf("~1 direction = %d\n",ss); +#ifdef NEVER + if (( ss && ep->npt >= 1) /* Not looking for a positive node */ + || (!ss && ep->nnt >= 1)) { /* Not looking for a negative node */ +printf("~1 no need to look for node in this direction\n"); + continue; + } +#endif + + if (ncnodes[ss] > 0) { /* There are nodes in wanted direction */ + si = ss; /* use wanted direction nodes */ + sta = 0; /* Start at max radius node */ + end = ncnodes[si]; /* End and least radius node */ + inc = 1; /* Increment */ + nods = cnodes[si]; /* +ve nodes */ +printf("~1 Looking for biggest angle, inc = %d\n",inc); +#ifdef NEVER /* Convex tracing ? */ + } else if (ncnodes[1-ss] > 0) { /* There are opposited direction nodes */ + si = 1-ss; /* use opposite direction nodes */ + sta = ncnodes[si]-1; /* Start at min radius node */ + end = -1; /* end and max radius node */ + inc = -1; /* Decrement */ + nods = cnodes[si]; +printf("~1 Looking for smallest angle, inc = %d\n",-1); +#endif /* NEVER */ + } else { +printf("~1 No points to search\n"); + continue; + } + ii = 0; + if (ncnodes[si] > 1) { /* If there is more than one to choose from */ + /* Go through each candidate in the most likely order */ + for (ii = sta; ii != end; ii += inc) { + +printf("~1 Candidate %d: node %d\n",ii,nods[ii]->gix); + /* Create the baricentric conversion for this candidate */ + for (f = 0; f < fdi; f++) /* The center point */ + A[f][0] = s->gam.cent[f] - nods[ii]->v[f]; + for (j = 0; j < (fdi-1); j++) { /* The edge points */ + for (f = 0; f < fdi; f++) + A[f][j+1] = ep->v[j]->v[f] - nods[ii]->v[f]; + } + + if (lu_decomp(A, fdi, pivx, &rip)) { +printf("~1 lu_decomp failed\n"); +for (f = 0; f < fdi; f++) /* The center point */ + A[f][0] = s->gam.cent[f] - nods[ii]->v[f]; +for (j = 0; j < (fdi-1); j++) { /* The edge points */ + for (f = 0; f < fdi; f++) + A[f][j+1] = ep->v[j]->v[f] - nods[ii]->v[f]; +} +printf("~1 A = \n"); +printf("~1 %f %f %f\n", A[0][0], A[0][1], A[0][2]); +printf("~1 %f %f %f\n", A[1][0], A[1][1], A[1][2]); +printf("~1 %f %f %f\n", A[2][0], A[2][1], A[2][2]); + warning("lu_decomp failed"); + continue; + } + + /* Test the remainder candidates against this simplex */ + for (i = sta; i != end; i += inc) { + if (i == ii) + continue; + +printf("~1 Test %d: node %d\n",i,nods[i]->gix); + for (f = 0; f < fdi; f++) /* The candidate point */ + B[f] = nods[i]->v[f] - nods[ii]->v[f]; + lu_backsub(A, fdi, pivx, B); + +#ifdef NEVER +printf("~1 baricentric = %f %f %f %f\n",B[0],B[1],B[2],1.0-B[0]-B[1]-B[2]); +{ +double tt[3], B3; +/* Check baricentric */ +B3 = 1.0-B[0]-B[1]-B[2]; +for (f = 0; f < fdi; f++) + tt[f] = 0.0; +for (f = 0; f < fdi; f++) + tt[f] += B[0] * s->gam.cent[f]; +for (f = 0; f < fdi; f++) + tt[f] += B[1] * ep->v[0]->v[f]; +for (f = 0; f < fdi; f++) + tt[f] += B[2] * ep->v[1]->v[f]; +for (f = 0; f < fdi; f++) + tt[f] += B3 * nods[ii]->v[f]; +printf("~1 target point %f %f %f\n", (double)nods[i][0], (double)nods[i][1], (double)nods[i][2]); +printf("~1 barice check %f %f %f\n", tt[0],tt[1],tt[2]); +} +#endif /* NEVER */ + if ((inc == 1 && B[0] < -EPS) /* other point is at higher angle */ + || (inc == -1 && B[0] > EPS)) { /* other point is at lower angle */ +printf("~1 candidate isn't best\n"); + break; + } + } + if (i == end) { +printf("~1 candidate IS the best\n"); + break; /* Candidate is at highest angle */ + } + } + if (ii == end) { +printf("~1Inconsistent candidate ordering\n"); + error("Inconsistent candidate ordering"); + } + } else { +printf("~1 there are only %d nodes, so don't search them\n",ncnodes[si]); + } +printf("~1 Making triangle with %d: node %d\n",ii,nods[ii]->gix); + make_tri(s, ep, nods[ii], inc == -1); + } + if (ep->npt < 1 || ep->nnt < 1) { + if (ep->npt < 1) + warning("###### Unable to locate +ve triangles for edge %d\n",ep->n); + if (ep->nnt < 1) + warning("###### Unable to locate -ve triangles for edge %d\n",ep->n); + } + } + } + + free_ivector(pivx, 0, fdi-1); + free_dmatrix(A, 0, fdi-1, 0, fdi-1); + free_dvector(B, 0, fdi-1); + + /* Dump out the edges */ + for (ep = s->gam.etop; ep != NULL; ep = ep->list) { + +printf("~1 edge no %d, npt = %d, nnt = %d\n",ep->n, ep->npt, ep->nnt); + } + } + return 0; +} + +/* ====================================================== */ +/* Gamut rspl setup functions */ + +/* Called by rspl initialisation */ +void +init_gam(rspl *s) { + + /* Methods */ + s->comp_gamut = gam_comp_gamut; +} + +/* Free up all the gamut info */ +void free_gam( +rspl *s /* Pointer to rspl grid */ +) { + int i; + int ssdi; + rvert *vp, *nvp; + redge *ep, *nep; + rtri *tp, *ntp; + + for (ssdi = 1; ssdi <= s->fdi-1; ssdi++) + rspl_free_ssimplex_info(s, &s->gam.ssi[ssdi]); + + /* Free the verticies */ + for (vp = s->gam.vtop; vp != NULL; vp = nvp) { + nvp = vp->list; + free(vp); + } + free(s->gam.verts); + + /* Free the edges */ + for (ep = s->gam.etop; ep != NULL; ep = nep) { + nep = ep->list; + free(ep); + } + free(s->gam.edges); + + /* Free the triangles */ + for (tp = s->gam.ttop; tp != NULL; tp = ntp) { + ntp = tp->list; + free(tp); + } + free(s->gam.tris); +} + + +/* Create the gamut surface structure. */ +/* Current this is: + + not rationalized to be non-overlapping + culled to eliminate overlaps + Have ink limits applied + +*/ + +/* ~~~ we need space ing gam to: + + store radial coordinates of output values + mark nodes as being visited. + store surface triangle structures + hold tadial output bounding acceleration structures + + + There is a lot in common with rev here. + we need to know input channel significance (what's black) + plus ink limit info. + we're going to create black generation information used by rev. +*/ + +/* ---------------------- */ +/* rspl gam diagnostic */ + +/* Diagnostic */ +void rspl_gam_plot(rspl *s, char *name) { + int i; + double col[3] = { 0.7, 0.7, 0.7 }; + rvert *vp, *nvp; + rtri *tp, *ntp; + + vrml *wrl; + + if ((wrl = new_vrml(name, 1, 0)) == NULL) + error("new_vrml failed\n"); + + /* Set the verticies */ + for (vp = s->gam.vtop; vp != NULL; vp = vp->list) { + wrl->add_vertex(wrl, 0, vp->v); + } + + /* Set the triangles */ + for (tp = s->gam.ttop; tp != NULL; tp = ntp) { + int ix[3]; + ntp = tp->list; + ix[0] = tp->v[0]->n; + ix[1] = tp->v[1]->n; + ix[2] = tp->v[2]->n; + wrl->add_triangle(wrl, 0, ix); + } + + wrl->make_triangles(wrl, 0, 0.0, NULL); +// wrl->make_triangles(wrl, 0, 0.0, col); + + wrl->del(wrl); +} + +#undef DEBUG +#undef DBGV +#undef DBG +#define DBGV(xxx) +#define DBG(xxx) + + + + + + |