diff options
Diffstat (limited to 'gamut/gamut.c')
| -rw-r--r-- | gamut/gamut.c | 7136 |
1 files changed, 7136 insertions, 0 deletions
diff --git a/gamut/gamut.c b/gamut/gamut.c new file mode 100644 index 0000000..052b8d8 --- /dev/null +++ b/gamut/gamut.c @@ -0,0 +1,7136 @@ + +/* + * gamut + * + * Gamut support routines. + * + * Author: Graeme W. Gill + * Date: 9/3/2000 + * Version: 1.00 + * + * Copyright 2000 - 2006 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. + */ + +/* + The gamut boundary is comuted using a variation of + Jan Morovic's Segment Maximum approach. The variations + are: + + The segments are filtered with an adaptive depth structure, + so that approximately the same detail is kept on the gamut + surface. Multiple direction vectors at each point are retained. + The resultant points are used to create the overal convex + jull, but in an adaptive, non-linearly scaled radial space, + that allows for convexity in the PCS result. +*/ + +/* TTBD: + +*/ + +#include <stdio.h> +#include <stdlib.h> +#include <stdarg.h> +#include <fcntl.h> +#include <string.h> +#include <math.h> +#include <time.h> +#include "icc.h" +#include "numlib.h" +#include "vrml.h" +#include "gamut.h" +#include "cgats.h" +#include "sort.h" /* ../h sort macro */ +#include "counters.h" /* ../h counter macros */ +#include "xlist.h" /* ../h expandable list macros */ + +#define COLORED_VRML + +#define DO_TWOPASS /* Second pass with adjustment based on first pass */ + +#define FAKE_SEED_SIZE 0.1 /* Usually 0.1 */ +#define TRIANG_TOL 1e-10 /* Triangulation tollerance, usually 1e-10 */ + +#define NORM_LOG_POW 0.25 /* Normal, colorspace lopow value */ +#define RAST_LOG_POW 0.05 /* Raster lopow value */ + +#undef TEST_CONVEX_HULL /* Use pure convex hull, not log hull */ + +#undef DEBUG_TRIANG /* Enable detailed triangulation debugging */ +#undef DEBUG_TRIANG_VRML /* Create debug.wrl for each step of triangulation */ + /* (Only on second pass if #define DO_TWOPASS) */ +#undef DEBUG_TRIANG_VRML_STEP /* Wait for return after each step */ + +#undef DEBUG_SPLIT_VRML /* Create debug.wrl for each step of triangle plane split */ + +#undef TEST_LOOKUP +#undef TEST_NEAREST +#undef TEST_NOSORT /* Turn off sorted insersion of verticies */ + +#undef SHOW_BUCKETS /* Show vertex buckets as surface */ +#undef SHOW_SPHERE /* Show surface on sphere */ +#undef SHOW_HULL_PNTS /* Show log() length convex hull points */ + +#undef ASSERTS /* Do internal checking */ + +#undef INTERSECT_DEBUG /* Turn on compute_vector_isect debugging, inc isect.wrl plot */ +#undef INTERSECT_VERIFY /* Verify compute_vector_isect against brute force search */ + +/* These routines support: + + representing the 3D gamut boundary of a device or image as + radial surface height, described by a triangular poligon hull. + + Interogate the surface to find the point lying on the hull + in the same radial direction as the query point. + + Interogate the surface to find the point lying on the hull + that is the closest to the query point. + + Save the gamut as a vrml format viewable file. + Save the gamut as a CGATS format .gam file. + + */ + +/* TTBD: + * + * Would be nice to take the exact colorspace specification + * (ie. Lab vs. Jab + viewing conditions), and store them + * in the .gam file, so that a warning can be issued if + * the gamut colorspace is a mismatch in icclink, or to be able + * to translate the verticies into the correct colorespace. + * + * Add interface to transform all the nodes, while + * keeping structure (For use within gamut map, or + * transforming existing gamut from Lab to Jab etc.) + * + * Add inteface to fetch the triangle information ? + * + * Need to cleanup error handling. We just exit() at the moment. + * + * Replace BSP tree optmisation with ball tree, to speedup + * radial, nearest, and vector search ? + * + * The log surface stuff is a compromise, that ends up with + * some dings and nicks, and a not fully detailed/smooth surface. + * The fundamental limitation is the use of the Delaunay triangulation + * criteria, and the triangulation algorithm dependence on it for consistency. + * Want to switch to triangulation algorithm that doesn't + * depend on this, and can triangulate concave objects, + * so that something like alpha-shapes criteria can be + * used to filter out non surface points. Inserting verticies + * from largest radius to smallest seems to do the right thing + * with cusp ridges, and this property needs to be retained. + * + */ + + +#ifndef M_PI +#define M_PI (3.1415926535897932384626433832795) +#endif + +static void triangulate(gamut *s); +static void del_gamut(gamut *s); +static gvert *expand_gamut(gamut *s, double in[3]); +static double getsres(gamut *s); +static int getisjab(gamut *s); +static int getisrast(gamut *s); +static void setnofilt(gamut *s); +static void getcent(gamut *s, double *cent); +static void getrange(gamut *s, double *min, double *max); +static int compatible(gamut *s, gamut *t); +static int nrawverts(gamut *s); +static int getrawvert(gamut *s, double pos[3], int ix); +static int nraw0verts(gamut *s); +static int getraw0vert(gamut *s, double pos[3], int ix); +static int nverts(gamut *s); +static int getvert(gamut *s, double *rad, double pos[3], int ix); +static int nssverts(gamut *s, double xvra); +static int getssvert(gamut *s, double *rad, double pos[3], double norm[3], int ix); +static void startnexttri(gamut *s); +static int getnexttri(gamut *s, int v[3]); +static double volume(gamut *s); +static int write_trans_vrml(gamut *s, char *filename, int doaxes, int docusps, + void (*transform)(void *cntx, double out[3], double in[3]), void *cntx); +static int write_vrml(gamut *s, char *filename, int doaxes, int docusps); +static int write_gam(gamut *s, char *filename); +static int read_gam(gamut *s, char *filename); +static double radial(gamut *s, double out[3], double in[3]); +static double nradial(gamut *s, double out[3], double in[3]); +static void nearest(gamut *s, double out[3], double in[3]); +static void nearest_tri(gamut *s, double out[3], double in[3], gtri **ctri); +static void setwb(gamut *s, double *wp, double *bp, double *kp); +static int getwb(gamut *s, double *cswp, double *csbp, double *cskp, double *gawp, double *gabp, double *gakp); +static void setcusps(gamut *s, int flag, double in[3]); +static int getcusps(gamut *s, double cusps[6][3]); +static int compute_vector_isect(gamut *s, double *p1, double *p2, double *min, double *max, + double *mint, double *maxt, gtri **mntri, gtri **mxtri); +static int compute_vector_isectns(gamut *s, double *p1, double *p2, gispnt *lp, int ll); +static double log_scale(gamut *s, double ss); +static int intersect(gamut *s, gamut *s1, gamut *s2); +static int compdstgamut(gamut *s, gamut *img, gamut *src, gamut *dst, int docomp, int doexp, + gamut *nedst, void (*cvect)(void *cntx, double *vec, double *pos), void *cntx); +static int vect_intersect(gamut *s, double *rvp, double *ip, double *p1, double *p2, gtri *t); + +/* in isecvol.c: */ +extern double isect_volume(gamut *s1, gamut *s2); + +/* ------------------------------------ */ + +/* Generic hue directions in degrees for Lab and Jab */ +/* Must be in increasing order */ +double gam_hues[2][7] = { + { + /* Lab */ + 36.0, /* Red */ + 101.0, /* Yellow */ + 149.0, /* Green */ + 225.0, /* Cyan */ + 300.0, /* Blue */ + 337.0, /* Magenta */ + 36.0 + 360.0 /* Red */ + }, + { + /* Jab */ + 28.0, /* Red */ + 101.0, /* Yellow */ + 148.0, /* Green */ + 211.0, /* Cyan */ +// 269.0, /* Blue */ + 250.0, /* Blue */ + 346.0, /* Magenta */ + 28.0 + 360.0 /* Red */ + } +}; + + +/* ------------------------------------ */ +static +gvert *new_gvert( +gamut *s, +gquad *p, /* Parent quad (may be NULL) */ +int i, /* Intended node in quad */ +int f, /* Flag value to be OR'ed */ +double pp[3], /* Point in xyz rectangular coordinates, absolute */ +double rr[3], /* Radial coordinates */ +double lrr0, /* log scaled rr[0] */ +double sp[3], /* Point mapped to surface of unit sphere, relative to center */ +double ch[3] /* Point mapped for convex hull testing, relative to center */ +) { + gvert *v; + + if (s->doingfake == 0 && s->ul != NULL) { /* There is an unused one available */ + v = s->ul; + s->ul = v->ul; + + } else { /* Allocate a new one */ + + if (s->nv >= s->na) { /* We need a new slot in the list */ + if (s->na == 0) { + s->na = 5; + if ((s->verts = (gvert **)malloc(s->na * sizeof(gvert *))) == NULL) { + fprintf(stderr,"gamut: malloc failed on %d gvert pointer\n",s->na); + exit (-1); + } + } else { + s->na *= 2; + if ((s->verts = (gvert **)realloc(s->verts, s->na * sizeof(gvert *))) == NULL) { + fprintf(stderr,"gamut: realloc failed on %d gvert pointer\n",s->na); + exit (-1); + } + } + } + + if ((v = (gvert *)calloc(1, sizeof(gvert))) == NULL) { + fprintf(stderr,"gamut: malloc failed on gvert object\n"); + exit (-1); + } + s->verts[s->nv] = v; + v->n = s->nv++; + } + v->tag = 1; + + if (p != NULL) { + v->w = 0.5 * p->w; + v->h = 0.5 * p->h; + + v->hc = p->hc; + if (i & 1) + v->hc += 0.5 * v->w; + else + v->hc -= 0.5 * v->w; + + v->vc = p->vc; + if (i & 2) + v->vc += 0.5 * v->h; + else + v->vc -= 0.5 * v->h; + } else { + v->w = 0.0; + v->h = 0.0; + v->hc = 0.0; + v->vc = 0.0; + } + + v->f = GVERT_NONE | f; + v->ul = NULL; + v->rc = 1; + + v->p[0] = pp[0]; + v->p[1] = pp[1]; + v->p[2] = pp[2]; + v->r[0] = rr[0]; + v->r[1] = rr[1]; + v->r[2] = rr[2]; + v->lr0 = lrr0; + v->sp[0] = sp[0]; + v->sp[1] = sp[1]; + v->sp[2] = sp[2]; + v->ch[0] = ch[0]; + v->ch[1] = ch[1]; + v->ch[2] = ch[2]; + + return v; +} + +/* Set the size of gvert angular segment */ +static +void set_gvert( +gvert *v, /* gvert to set */ +gquad *p, /* Parent quad */ +int i /* Intended node in quad */ +) { + v->w = 0.5 * p->w; + v->h = 0.5 * p->h; + + v->hc = p->hc; + if (i & 1) + v->hc += 0.5 * v->w; + else + v->hc -= 0.5 * v->w; + + v->vc = p->vc; + if (i & 2) + v->vc += 0.5 * v->h; + else + v->vc -= 0.5 * v->h; +} + +/* Increment the reference count on a gvert */ +static gvert *inc_gvert(gamut *s, gvert *v) { + if (v == NULL) + return v; + v->rc++; +//printf("~1 incremented count on 0x%x to %d\n",v,v->rc); + return v; +} + +/* Decrement the reference count on a gvert */ +/* Copes with NULL gvert. */ +/* If the reference count goes to 0, place the */ +/* vert on the unused list. */ +static void dec_gvert(gamut *s, gvert *v) { + if (v == NULL) + return; + + v->rc--; +//printf("~1 decremended count on 0x%x to %d\n",v,v->rc); + +#ifdef ASSERTS + if (v->tag != 1) + error("Assert: doing decremented on gquad node"); + + if (v->rc < 0) + error("Assert: decremented gvert ref count too far"); +#endif + + if (v->rc <= 0) { /* Add it to the unused list */ + memset((void *)v, 0, sizeof(gvert)); + v->ul = s->ul; + s->ul = v; + } +} + +/* Delete all the gverts */ +static void del_gverts(gamut *s) { + int i; + for (i = 0; i < s->nv; i++) { + free(s->verts[i]); + } + if (s->verts != NULL) { + free(s->verts); + s->na = 0; + s->nv = 0; + } +} + +/* ------------------------------------ */ + +static +gquad *new_gquad( +gquad *p, /* Parent quad */ +int i /* Intended node in quad */ +) { + gquad *q; + if ((q = (gquad *)calloc(1, sizeof(gquad))) == NULL) { + fprintf(stderr,"gamut: calloc failed on gquad object\n"); + exit (-1); + } + q->tag = 2; + + q->w = 0.5 * p->w; + q->h = 0.5 * p->h; + + q->hc = p->hc; + if (i & 1) + q->hc += 0.5 * q->w; + else + q->hc -= 0.5 * q->w; + + q->vc = p->vc; + if (i & 2) + q->vc += 0.5 * q->h; + else + q->vc -= 0.5 * q->h; + + return q; +} + +/* Same as above, but create with explicit size */ +static +gquad *new_gquad2( +double l, /* Left border */ +double r, /* Right border */ +double b, /* Top border */ +double t /* Bottom border */ +) { + gquad *q; + if ((q = (gquad *)calloc(1, sizeof(gquad))) == NULL) { + fprintf(stderr,"gamut: calloc failed on gquad object\n"); + exit (-1); + } + q->tag = 2; + q->w = r - l; + q->h = t - b; + q->hc = (l + r) * 0.5; + q->vc = (t + b) * 0.5; + return q; +} + +static void +del_gquad(gquad *q) { + int i; + + if (q == NULL) + return; + for (i = 0; i < 4; i++) { + gnode *n = (gnode *)q->qt[i][0]; + if (n != NULL && n->tag == 2) + del_gquad((gquad *)n); /* Recurse */ + } + free(q); +} + +/* Helper functions */ + +/* Given a gquad and a location, decide which quandrant we're in */ +static int gquad_quadrant(gquad *q, double p[3]) { + int i; + +#ifdef ASSERTS + if (p[1] < (q->hc - q->w * 0.5 - 1e-10) + || p[1] > (q->hc + q->w * 0.5 + 1e-10) + || p[2] < (q->vc - q->h * 0.5 - 1e-10) + || p[2] > (q->vc + q->h * 0.5 + 1e-10)) { + fprintf(stderr,"error! point doesn't fall into bucket chosen for it!!!\n"); + fprintf(stderr,"point x: %f < %f\n", p[1], (q->hc - q->w * 0.5)); + fprintf(stderr,"point x: %f > %f\n", p[1], (q->hc + q->w * 0.5)); + fprintf(stderr,"point y: %f < %f\n", p[2], (q->vc - q->h * 0.5)); + fprintf(stderr,"point y: %f > %f\n", p[2], (q->vc + q->h * 0.5)); + exit(-1); + } +#endif + + i = 0; + if (p[1] >= q->hc) + i |= 1; + if (p[2] >= q->vc) + i |= 2; + + return i; +} + +/* ------------------------------------ */ +/* Allocate a BSP decision node structure */ +gbspn *new_gbspn(void) { + gbspn *t; + static int n = 0; /* Serial number */ + if ((t = (gbspn *) calloc(1, sizeof(gbspn))) == NULL) { + fprintf(stderr,"gamut: malloc failed - bspn node\n"); + exit(-1); + } + t->tag = 1; /* bspn decision node */ + t->n = n++; + + return t; +} + +/* Delete a BSP decision node struture */ +void del_gbspn(gbspn *t) { + free(t); +} + +/* ------------------------------------ */ +/* Allocate a BSP tree triangle list node structure */ +gbspl *new_gbspl( +int nt, /* Number of triangles in the list */ +gtri **t /* List of triangles to copy into structure */ +) { + gbspl *l; + int i; + static int n = 0; /* Serial number */ + if ((l = (gbspl *) calloc(1, sizeof(gbspl) + nt * sizeof(gtri *))) == NULL) { + fprintf(stderr,"gamut: malloc failed - bspl triangle tree node\n"); + exit(-1); + } + l->tag = 3; /* bspl triangle list node */ + l->n = n++; + l->nt = nt; + for (i = 0; i < nt; i++) + l->t[i] = t[i]; + + return l; +} + +/* Delete a BSP tree triangle list structure */ +void del_gbspl(gbspl *l) { + free(l); +} + +/* ------------------------------------ */ +/* Allocate a triangle structure */ +gtri *new_gtri(void) { + gtri *t; + static int n = 0; /* Serial number */ + if ((t = (gtri *) calloc(1, sizeof(gtri))) == NULL) { + fprintf(stderr,"gamut: malloc failed - gamut surface triangle\n"); + exit(-1); + } + t->tag = 2; /* Triangle */ + t->n = n++; + + return t; +} + +/* Delete a triangle struture */ +void del_gtri(gtri *t) { + free(t); +} + +/* ------------------------------------ */ +/* Allocate an edge structure */ +gedge *new_gedge(void) { + gedge *t; + static int n = 0; /* Serial number */ + if ((t = (gedge *) calloc(1, sizeof(gedge))) == NULL) { + fprintf(stderr,"gamut: malloc failed - triangle edge\n"); + exit(-1); + } + t->n = n++; + return t; +} + +/* Delete an edge struture */ +void del_gedge(gedge *t) { + free(t); +} + +/* ------------------------------------ */ + +/* Create a standard gamut map */ +gamut *new_gamut( +double sres, /* Resolution (in rect coord units) of surface triangles */ + /* 0.0 = default */ +int isJab, /* Flag indicating Jab space */ +int isRast /* Flag indicating Raster rather than colorspace */ +) { + gamut *s; + +#ifdef ASSERTS + fprintf(stderr,">>>>>>> ASSERTS ARE COMPILED INTO GAMUT.C <<<<<<<\n"); +#endif /* ASSERTS */ +#ifdef TEST_LOOKUP + fprintf(stderr,">>>>>>> TEST_LOOKUP IS COMPILED INTO GAMUT.C <<<<<<<\n"); +#endif +#ifdef TEST_NEAREST + fprintf(stderr,">>>>>>> TEST_NEAREST IS COMPILED INTO GAMUT.C <<<<<<<\n"); +#endif +#ifdef TEST_NOSORT /* Turn off sorted insersion of verticies */ + fprintf(stderr,">>>>>>> TEST_NOSORT IS COMPILED INTO GAMUT.C <<<<<<<\n"); +#endif +#ifdef INTERSECT_VERIFY + fprintf(stderr,">>>>>>> INTERSECT_VERIFY IS COMPILED INTO GAMUT.C <<<<<<<\n"); +#endif + + if ((s = (gamut *)calloc(1, sizeof(gamut))) == NULL) { + fprintf(stderr,"gamut: calloc failed on gamut object\n"); + exit (-1); + } + + if (sres <= 0.0) + sres = 10.0; /* default */ + if (sres > 15.0) /* Anything less is very poor */ + sres = 15.0; + s->sres = sres; + + if (isJab != 0) { + s->isJab = 1; + } + + if (isRast != 0) { + s->isRast = 1; + } + + if (s->isRast) { + s->logpow = RAST_LOG_POW; /* Wrap the surface more closely */ + s->no2pass = 1; /* Only do one pass */ + } else { + s->logpow = NORM_LOG_POW; /* Convex hull compression power */ + s->no2pass = 0; /* Do two passes */ + } + + /* Center point for radial values, surface creation etc. */ + /* To compare two gamuts using radial values, their cent must */ + /* be the same. */ + s->cent[0] = 50.0; + s->cent[1] = 0.0; + s->cent[2] = 0.0; + + s->mx[0] = -1e38; + s->mx[1] = -1e38; + s->mx[2] = -1e38; + s->mn[0] = 1e38; + s->mn[1] = 1e38; + s->mn[2] = 1e38; + + /* Create top level quadtree nodes */ + s->tl = new_gquad2(-M_PI, 0.0, -M_PI/2.0, M_PI/2.0); /* Left one */ + s->tr = new_gquad2(0.0, M_PI, -M_PI/2.0, M_PI/2.0); /* Right one */ + + INIT_LIST(s->tris); /* Init triangle list */ + INIT_LIST(s->edges); /* Init edge list (?) */ + s->read_inited = 0; + s->lu_inited = 0; + s->ne_inited = 0; + s->cswbset = 0; + s->gawbset = 0; + + /* Setup methods */ + s->del = del_gamut; + s->expand = expand_gamut; + s->getsres = getsres; + s->getisjab = getisjab; + s->getisrast = getisrast; + s->setnofilt = setnofilt; + s->getcent = getcent; + s->getrange = getrange; + s->compatible = compatible; + s->nrawverts = nrawverts; + s->getrawvert = getrawvert; + s->nraw0verts = nraw0verts; + s->getraw0vert = getraw0vert; + s->nssverts = nssverts; + s->getssvert = getssvert; + s->nverts = nverts; + s->getvert = getvert; + s->startnexttri = startnexttri; + s->getnexttri = getnexttri; + s->getvert = getvert; + s->volume = volume; + s->intersect = intersect; + s->compdstgamut = compdstgamut; + s->radial = radial; + s->nradial = nradial; + s->nearest = nearest; + s->nearest_tri = nearest_tri; + s->vector_isect = compute_vector_isect; + s->vector_isectns = compute_vector_isectns; + s->setwb = setwb; + s->getwb = getwb; + s->setcusps = setcusps; + s->getcusps = getcusps; + s->write_vrml = write_vrml; + s->write_trans_vrml = write_trans_vrml; + s->write_gam = write_gam; + s->read_gam = read_gam; + + return s; +} + +static void del_gnn(gnn *p); +static void del_gbsp(gbsp *n); + +/* Free and clear the triangulation structures, */ +/* and clear the triangulation vertex flags. */ +static void del_triang(gamut *s) { + int i; + gtri *tp; /* Triangle pointer */ + gedge *ep; + + /* Recursively free radial lookup acceleration structures */ + /* Do this before we delete triangles, because there may */ + /* be triangles in the tree. */ + if (s->lutree != NULL) { + del_gbsp(s->lutree); + s->lutree = NULL; + } + + if (s->tris != NULL) { + tp = s->tris; /* Delete all the triangles */ + FOR_ALL_ITEMS(gtri, tp) { + DEL_LINK(s->tris, tp); + del_gtri(tp); + } END_FOR_ALL_ITEMS(tp); + + INIT_LIST(s->tris); /* Init triangle list */ + } + + if (s->edges != NULL) { + ep = s->edges; /* Delete all the edges */ + FOR_ALL_ITEMS(gedge, ep) { + DEL_LINK(s->edges, ep); + del_gedge(ep); + } END_FOR_ALL_ITEMS(ep); + INIT_LIST(s->edges); /* Init edge list */ + } + + s->lu_inited = 0; + + if (s->nns != NULL) { + del_gnn(s->nns); + s->nns = NULL; + } + s->ne_inited = 0; + + /* Reset the vertex flags triangulation changes */ + for (i = 0; i < s->nv; i++) { + s->verts[i]->f &= ~GVERT_TRI; + s->verts[i]->f &= ~GVERT_INSIDE; + } +} + +static void +del_gamut(gamut *s) { + del_gquad(s->tl); + del_gquad(s->tr); + + del_triang(s); + del_gverts(s); + + if (s->ss != NULL) + s->ss->del(s->ss); + + free(s); +} + +/* ===================================================== */ +/* Segmented maxima filter code */ +/* ===================================================== */ + +/* + This implementation has two twists on the usual + segmented maxima filtering: + + * Rather than using a uniform radial grid, it + uses an adaptive, quadtree structure, so that + rather than having a constant angular resolution + for the segments, the sements are chosen so as to + aproximately correspond to a constant surface detail + level. [ A subtle problem with this approach is that + some points will be discarded early on, that wouldn't + be discarded later, when the quadtree is finer. A hack + is to run the points throught twice.] + + * Rather than keep the single sample with the + largest radius in each radial segment, + four samples are kept, each being the largest + in a different direction. This helps avoid + "knicks" in edges, where a non edge sample + displaces an edge sample within a segment. + + */ + +/* Helper function that returns nz if v1 should replace v2 */ +static int smreplace( +gamut *s, +int d, /* Slot number */ +gvert *v1, /* Candidate vertex */ +gvert *v2 /* Existing vertex */ +) { + double xx, w[3], c[3]; + int j; + if (v2 == NULL) + return 1; + + /* Filter out any points that are almost identical */ + /* This can cause numerical problems in the triangle BSP tree creation. */ + for (xx = 0.0, j = 0; j < 3; j++) { + double tt = v1->p[j] - v2->p[j]; + xx += tt * tt; + } + if (xx < (1e-4 * 1e-4)) + return 0; + + c[0] = s->cent[0]; + c[1] = s->cent[1]; + c[2] = s->cent[2]; + + /* Set L, a & b weighting depending on the slot */ + switch(d) { + case 1: + w[0] = 0.5; + w[1] = 1.0; + w[2] = 1.0; + break; + case 2: + w[0] = 2.0; + w[1] = 1.0; + w[2] = 1.0; + break; + case 3: + w[0] = 4.0; + w[1] = 1.0; + w[2] = 1.0; + break; + case 4: + w[0] = 0.0; + w[1] = 1.0; + w[2] = 0.1; + break; + case 5: + w[0] = 0.0; + w[1] = 0.1; + w[2] = 1.0; + break; + default: + w[0] = 1.0; + w[1] = 1.0; + w[2] = 1.0; + break; + } + w[0] *= w[0]; /* Because we're weighting the squares */ + w[1] *= w[1]; + w[2] *= w[2]; + return ( w[0] * (v1->p[0] - c[0]) * (v1->p[0] - c[0]) + + w[1] * (v1->p[1] - c[1]) * (v1->p[1] - c[1]) + + w[2] * (v1->p[2] - c[2]) * (v1->p[2] - c[2])) + > ( w[0] * (v2->p[0] - c[0]) * (v2->p[0] - c[0]) + + w[1] * (v2->p[1] - c[1]) * (v2->p[1] - c[1]) + + w[2] * (v2->p[2] - c[2]) * (v2->p[2] - c[2])); +} + + +/* Expand the gamut by adding a point. */ +/* If nofilter is set, return NULL if the point */ +/* is discarded, or the address of the point representing */ +/* the point added. If nofiler is not set, return NULL */ +static gvert *expand_gamut( +gamut *s, +double pp[3] /* rectangular coordinate of point */ +) { + gnode *n; /* Current node */ + gvert *nv, *ov; /* new vertex, old vertex */ + gquad *q; /* Parent quad */ + int i; /* Sub element within quad */ + int k; /* Index of direction slot */ + double rr[3]; /* Radial coordinate version of pp[] */ + double sp[3]; /* Unit shere mapped version of pp[] relative to center */ + double ch[3]; /* Convex hull testing mapped version of pp[] relative to center */ + double lrr0; /* log scaled rr[0] */ + double hang, vang; /* Critical angles for this points depth */ + double aa; + int j; + + if (s->tris != NULL || s->read_inited || s->lu_inited || s->ne_inited) { + fprintf(stderr,"Can't add points to gamut now!\n"); + exit(-1); + } + + if (s->doingfake == 0) + s->cu_inited = 0; /* Invalidate cust info */ + + /* Tracl bounding range */ + for (j = 0; j < 3; j++) { + if (pp[j] > s->mx[j]) + s->mx[j] = pp[j]; + if (pp[j] < s->mn[j]) + s->mn[j] = pp[j]; + } + + /* Convert to radial coords */ + gamut_rect2radial(s, rr, pp); + + if (rr[0] < 1e-6) /* Ignore a point right at the center */ + return NULL; + + /* Figure log scaled radius */ + lrr0 = log_scale(s, rr[0]); + + /* Compute unit shere mapped location */ + aa = 1.0/rr[0]; /* Adjustment to put in on unit sphere */ + for (j = 0; j < 3; j++) + sp[j] = (pp[j] - s->cent[j]) * aa; + + /* Compute hull testing mapped version */ + for (j = 0; j < 3; j++) + ch[j] = sp[j] * lrr0; + + /* compute twice angle resolution required (will compare to parent size) */ + hang = pow(rr[0], 1.01) * fabs(cos(rr[2])); + if (hang < 1e-9) + hang = 1e-9; + hang = 4.0 * s->sres/hang; + vang = 4.0 * s->sres/pow(rr[0], 1.01); + +//printf("~1 Point at %f %f %f, radial %f %f %f\n", pp[0], pp[1], pp[2], rr[0], rr[1], rr[2]); +//printf("~1 shere at %f %f %f, log %f %f %f, vhang %f %f\n", sp[0], sp[1], sp[2], ch[0], ch[1], ch[2], vang, hang); + + /* If nofilter flag is set, add all points as verticies */ + if (s->nofilter) { + + /* Filter out any points that are almost identical. */ + /* This can cause numerical problems in the triangle BSP tree creation. */ + for (i = 0; i < s->nv; i++) { + double xx; + + for (xx = 0.0, j = 0; j < 3; j++) { + double tt = pp[j] - s->verts[i]->p[j]; + xx += tt * tt; + } + if (xx < (1e-4 * 1e-4)) { + if (s->doingfake) + s->verts[i]->f |= GVERT_ESTP; + + return s->verts[i]; /* Existing point becomes added point */ + } + } + + /* Create a vertex for the point we're possibly adding */ + nv = new_gvert(s, NULL, 0, GVERT_SET | (s->doingfake ? (GVERT_FAKE | GVERT_ESTP) : 0), + pp, rr, lrr0, sp, ch); + + return nv; + } + /* else filter using adaptive segmented maxima */ + + /* Start by looking at the top level quads */ + if (rr[1] >= 0.0) { + q = s->tr; + } else { + q = s->tl; + } + n = (gnode *)q; +//printf("~1 Starting with quad 0x%x, width %f, height %f\n",q, q->w, q->h); + + /* Now recurse until we have a virtex at the right location and depth */ + for (;;) { + /* Recurse into quad node n */ + q = (gquad *)n; /* Parent node */ + i = gquad_quadrant(q, rr); /* Quadrand of parent */ + n = q->qt[i][0]; /* Child node in quadrant */ + +//printf("~1 Current quad 0x%x, width %f, height %f\n",q, q->w, q->h); +//printf("~1 Current child in quadrant %d, node 0x%x, type %d\n", i, n, n != NULL ? n->tag : 0); + + /* If we're at the right depth to create a vertex, break out of decent loop. */ + + if (n == NULL) { /* Create new node */ + if (q->w <= hang && q->h <= vang) { +//printf("~1 We're at the right depth to add vertex\n"); + break; + } + /* Else create a new quad */ + n = (gnode *)new_gquad(q, i); + q->qt[i][0] = n; +//printf("~1 Empty child node not deep enough, creating new quad node 0x%x\n",n); + + /* If we've found verticies at this node */ + } else if (n->tag == 1) { + int j; + gquad *qq; /* New child quad */ + gvert *vv[NSLOTS]; /* Existing verticies at this level */ + + if (q->w <= hang && q->h <= vang) { +//printf("~1 We're at the right depth to replace vertex\n"); + break; + } +//printf("~1 deepening verticies\n"); + + for (k = 0; k < NSLOTS; k++) + vv[k] = (gvert *)q->qt[i][k]; /* Save pointers to current verticies */ + +//printf("~1 existing verticies are 0x%x, 0x%x, 0x%x, 0x%x\n", vv[0], vv[1], vv[2], vv[3]); + + /* make a quad to replace the current verticies */ + qq = new_gquad(q, i); + n = (gnode *)qq; + q->qt[i][0] = n; + for (k = 1; k < NSLOTS; k++) + q->qt[i][k] = NULL; + +//printf("~1 added quad 0x%x to quadrant %d\n",i,q); + + /* Distribute verticies that were here, into new quad */ + for (j = 0; j < NSLOTS; j++) { /* For all existing verticies */ + + if (vv[j] == NULL) + continue; + +//printf("~1 re-distributing verticy 0x%x\n",vv[j]); + i = gquad_quadrant(qq, vv[j]->r); /* Quadrant for existing vertex */ + + set_gvert(vv[j], qq, i); /* Update vertex node location */ + + nv = vv[j]; + for (k = 0; nv != NULL && k < NSLOTS; k++) { /* For direction slot */ + ov = (gvert *)qq->qt[i][k]; + if (smreplace(s, k, nv, ov)) { + if (k == 0) { /* Track points that are in k == 0 direction */ + if (ov != NULL && ov->k0 > 0) + ov->k0--; + nv->k0++; + } +#ifndef NEVER + qq->qt[i][k] = (gnode *)inc_gvert(s, nv); + dec_gvert(s, ov); +#else + /* Use slots for best, 2nd best, etc */ + qq->qt[i][k] = (gnode *)inc_gvert(s, nv); + dec_gvert(s, nv); + nv = ov; +#endif +//printf("Node 0x%x rc %d at %f %f %f is replacing\n",nv, nv->rc, nv->p[0], nv->p[1], nv->p[2]); +//if (ov != NULL) printf(" replacing node 0x%x rc %d at %f %f %f\n",ov, ov->rc, ov->p[0], ov->p[1], ov->p[2]); +//else printf(" NULL\n"); + } + } + dec_gvert(s, nv); + } + } + /* Else it's a quad, and we will decend into it */ + + } /* keep decending until we find right depth */ + +//printf("~1 Got parent quad 0x%x, quadrant %d, vertex 0x%x\n", q, i, n); + + /* Create a vertex for the point we're possibly adding */ + nv = new_gvert(s, q, i, GVERT_SET, pp, rr, lrr0, sp, ch); + + /* Replace any existing gverts with this one */ + for (k = 0; k < NSLOTS; k++) { /* For direction slot */ + ov = (gvert *)q->qt[i][k]; + if (smreplace(s, k, nv, ov)) { + if (k == 0) { /* Track points that are in k == 0 direction */ + if (ov != NULL && ov->k0 > 0) + ov->k0--; + nv->k0++; + } +#ifndef NEVER + q->qt[i][k] = (gnode *)inc_gvert(s, nv); + dec_gvert(s, ov); +#else + /* Use slots for best, 2nd best, etc */ + q->qt[i][k] = (gnode *)inc_gvert(s, nv); + dec_gvert(s, nv); + nv = ov; +#endif + } + } + dec_gvert(s, nv); /* Make sure it's reclaimed if wasn't used */ + +//printf("~1 Point is done\n\n"); + + return NULL; +} + +/* ------------------------------------ */ +/* Initialise this gamut with the intersection of the */ +/* the two given gamuts. Return NZ on error. */ +/* Return 1 if gamuts are not compatible */ +/* (We assume that the this gamut is currently empty) */ +static int intersect(gamut *s, gamut *sa, gamut *sb) { + int i, j, k; + gamut *s1, *s2; + + if (sa->compatible(sa, sb) == 0) + return 1; + + if IS_LIST_EMPTY(sa->tris) + triangulate(sa); + if IS_LIST_EMPTY(sb->tris) + triangulate(sb); + + s->isJab = sa->isJab; + + if (sa->isRast || sb->isRast) + s->isRast = 1; + + for (j = 0; j < 3; j++) + s->cent[j] = sa->cent[j]; + + /* Clear some flags */ + s->cswbset = 0; + s->cswbset = 0; + s->dcuspixs = 0; + + /* If either is a raster gamut, make it a raster gamut */ + if (s->isRast) + s->logpow = RAST_LOG_POW; /* Wrap the surface more closely */ + else + s->logpow = NORM_LOG_POW; + + /* Don't filter the points (gives a more accurate result ?) */ + s->nofilter = 1; + + /* Add each source gamuts verticies that lie within */ + /* the other gamut */ + for (k = 0; k < 2; k++) { + gtri *tp1, *tp2; /* Triangle pointers */ + + if (k == 0) { + s1 = sa; + s2 = sb; + } else { + s1 = sb; + s2 = sa; + } + for (i = 0; i < s1->nv; i++) { + double pl; + + if (!(s1->verts[i]->f & GVERT_TRI)) + continue; + + pl = s2->nradial(s2, NULL, s1->verts[i]->p); + if (pl <= (1.0 + 1e-9)) { + expand_gamut(s, s1->verts[i]->p); + s1->verts[i]->f &= ~GVERT_ISOS; /* s1 vert is not outside s2 */ + } else { + s1->verts[i]->f |= GVERT_ISOS; /* s1 vert is outside s2 */ + } + } + + /* Now find the edges that intersect the other gamut */ + tp1 = s1->tris; + FOR_ALL_ITEMS(gtri, tp1) { /* For all s1 triangles */ + + for (j = 0; j < 3; j++) { /* For all edges in s1 triangle */ + /* If edge passes through the other gamut */ + if ((tp1->e[j]->v[0]->f ^ tp1->e[j]->v[1]->f) & GVERT_ISOS) { + + /* Exhaustive search of other triangles in s2, */ + /* to find the one that the edge intersects with. */ + tp2 = s2->tris; + FOR_ALL_ITEMS(gtri, tp2) { + double pv; + double tt[3]; + + /* Do a min/max intersection elimination test */ + for (i = 0; i < 3; i++) { + if (tp2->mix[1][i] < tp1->mix[0][i] + || tp2->mix[0][i] > tp1->mix[1][i]) + break; /* min/max don't overlap */ + } + if (i < 3) + continue; /* Skip this triangle, it can't intersect */ + + if (vect_intersect(s1, &pv, tt, tp1->e[j]->v[0]->p, tp1->e[j]->v[1]->p, tp2) + && pv >= (0.0 - 1e-10) && pv <= (1.0 + 1e-10)) { + expand_gamut(s, tt); + } + } END_FOR_ALL_ITEMS(tp2); + } + } + + } END_FOR_ALL_ITEMS(tp1); + } + + s->nofilter = 0; + + return 0; +} + +/* ------------------------------------ */ +/* + Initialise this gamut with a destination mapping gamut. + s1 is the image gamut (a possible sub-set of the src gamut) + s2 is the source gamut + s3 is the destination gamut + + if docomp: + this gamut will be set to the smaller of the img & dst gamuts + else + this gamut will be the img gamut + + if doexp + this gamut will be expanded by the amount dst is outside the src gamut. + + The vector direction of "inwards" is that returned by the + callback function if it is supplied, or radially inwards + if it is not. + + Return 1 if gamuts are not compatible. + (We assume that the this gamut is currently empty) + + */ +#define MXNIS 40 /* Maximum raw intersections handled */ + +static int compdstgamut( +gamut *s, /* This */ +gamut *s1, /* Image gamut, assumed to be a subset of source gamut */ +gamut *s2, /* The source space gamut */ +gamut *s3, /* The destination space gamut */ +int docomp, /* Flag, nz if compression is enabled */ +int doexp, /* Flag, nz if expansion is enabled. */ +gamut *nedst, /* Optional - if doexp, then expand nedst with non-expanded dst gamut */ +void (*cvect)(void *cntx, double *p2, double *p1), /* Compression direction callback */ +void *cntx /* Returns p2 which is in desired direction from given p1 */ +) { +#ifdef STATS + int tedges, edgestested, testhits; +#endif + int i, j, k; + gamut *ss[3]; + gispnt lp1[MXNIS], lp2[MXNIS], lp3[MXNIS]; /* Intersection lists */ + int ll1, ll2, ll3; /* Returned list length */ + int ii, jj, kk; /* List indexes */ + + if (s1->compatible(s1, s2) == 0 + || s1->compatible(s2, s3) == 0) + return 1; + + if IS_LIST_EMPTY(s1->tris) + triangulate(s1); + if IS_LIST_EMPTY(s2->tris) + triangulate(s2); + if IS_LIST_EMPTY(s3->tris) + triangulate(s3); + + s->isJab = s1->isJab; + s->isRast = s1->isRast; + + for (j = 0; j < 3; j++) + s->cent[j] = s1->cent[j]; + + /* Clear some flags */ + s->cswbset = 0; + s->cswbset = 0; + s->dcuspixs = 0; + + /* If s1 is a raster gamut, make output a raster gamut */ + if (s->isRast) + s->logpow = RAST_LOG_POW; /* Wrap the surface more closely */ + else + s->logpow = NORM_LOG_POW; + + /* Don't filter the points (gives a more accurate result ?) */ + s->nofilter = 1; + + ss[0] = s1; + ss[1] = s2; + ss[2] = s3; + +//printf("~1 compdstgamut docomp %d, doexp %d\n",docomp,doexp); + /* Reset the above/below flags */ + /* 1 = img above dst */ + /* 2 = img below dst */ + /* 4 = src above dst */ + /* 8 = src below dst */ + for (k = 0; k < 3; k++) { /* For img, src & dst */ + for (i = 0; i < ss[k]->nv; i++) { + if (!(ss[k]->verts[i]->f & GVERT_TRI)) + continue; + + ss[k]->verts[i]->as = 0; + } + } + + /* Use all the triangle vertices from the three gamuts */ + /* as candidate points, because any of them might */ + /* determine a surface feature. */ + for (k = 0; k < 3; k++) { /* For img, src & dst */ + + for (i = 0; i < ss[k]->nv; i++) { + double pp[3], ppv, p2[3]; + double rr, r4; + + if (!(ss[k]->verts[i]->f & GVERT_TRI)) + continue; + + icmCpy3(pp, ss[k]->verts[i]->p); /* Point in question */ + +//printf("\n~1 k %d, point %d: %f %f %f\n", k,i,pp[0],pp[1],pp[2]); + if (cvect != NULL) + cvect(cntx, p2, pp); /* Get mapping direction */ + else + icmCpy3(p2, ss[k]->cent); /* Radial vector */ + icmNormalize33(p2, p2, pp, 1.0); + + /* Locate the intersecting segments for each gamut */ + ll1 = ll2 = ll3 = 0; + ll1 = s1->vector_isectns(s1, pp, p2, lp1, MXNIS); /* Image gamut */ + if (doexp) /* For dst - src expansion */ + ll2 = s2->vector_isectns(s2, pp, p2, lp2, MXNIS); /* Src gamut */ + if (docomp || doexp) /* For img ^ dst or dst - src */ + ll3 = s3->vector_isectns(s3, pp, p2, lp3, MXNIS); /* Dst gamut */ + +//printf("~1 ll1 %d, ll2 %d, ll3 %d\n",ll1,ll2,ll3); +#ifdef NEVER + printf("img segments:\n"); + for (ii = 0; ii < ll1; ii++) + printf("Isect %d: pv %f, dir %d, edge %d, tri %d\n",ii,lp1[ii].pv,lp1[ii].dir,lp1[ii].edge,lp1[ii].tri->n); + printf("src segments:\n"); + for (ii = 0; ii < ll2; ii++) + printf("Isect %d: pv %f, dir %d, edge %d, tri %d\n",ii,lp2[ii].pv,lp2[ii].dir,lp2[ii].edge,lp2[ii].tri->n); + printf("dst segments:\n"); + for (ii = 0; ii < ll3; ii++) + printf("Isect %d: pv %f, dir %d, edge %d, tri %d\n",ii,lp3[ii].pv,lp3[ii].dir,lp3[ii].edge,lp3[ii].tri->n); +#endif + /* Now go through each image segment */ + for (ii = 0; ii < ll1; ii += 2) { + + icmCpy3(pp, lp1[ii].ip); + ppv = lp1[ii].pv; + +//printf("~1 img pnt at %f\n",lp1[ii].pv); + if (docomp) { + + /* Locate a dst segment around or below img point */ + for (kk = 0; kk < ll3; kk += 2) { + if ((lp1[kk+1].pv + 1e-8) >= ppv) + break; + } + + if (kk >= ll3) { /* No dst segments or none below */ + ss[k]->verts[i]->as |= 1; /* img above dst */ +//printf("~1 img pnt is outside dst\n"); + continue; + } else { +//printf("~1 ing %f - %f, dst %f - %f\n", lp1[ii].pv,lp1[ii+1].pv,lp3[kk].pv,lp3[kk+1].pv); + /* Use the dst point if it is smaller */ + if (ppv < lp3[kk].pv) { + icmCpy3(pp, lp3[kk].ip); + ppv = lp3[kk].pv; +//printf("~1 dst is smaller %f\n",ppv); + ss[k]->verts[i]->as |= 1; /* img above dst */ + } else { +//printf("~1 ing is smaller %f\n",ppv); + ss[k]->verts[i]->as |= 2; /* img below dst */ + } + } + } + + if (nedst != NULL) + expand_gamut(nedst, pp); + + while (doexp) { + /* Locate a src segment that ing point lies in. */ + for (jj = 0; jj < ll2; jj += 2) { + if (lp1[ii].pv >= (lp2[jj].pv - 1e-8) + && lp1[ii].pv <= (lp2[jj+1].pv + 1e-8)) + break; + } + if (jj >= ll2) { + ss[k]->verts[i]->as |= 4; /* src above dst ?? */ + break; /* No overlapping segment */ + } + + /* Locate a dst segment that src point lies in */ + for (kk = 0; kk < ll3; kk += 2) { + if (lp2[jj].pv >= (lp3[kk].pv - 1e-8) + && lp2[jj].pv <= (lp3[kk+1].pv + 1e-8)) + break; + } + if (kk >= ll2) { + ss[k]->verts[i]->as |= 4; /* src above dst ?? */ + break; /* No overlapping segment */ + } + + /* if src is outside dst, do nothing */ + if (lp3[kk].pv >= lp2[jj].pv) { + ss[k]->verts[i]->as |= 4; /* src above dst */ + break; + } + + /* Expand point by dst - src */ + icmAdd3(pp, pp, lp3[kk].ip); + icmSub3(pp, pp, lp2[jj].ip); + ss[k]->verts[i]->as |= 8; /* src below dst */ +//printf("~1 expanding point by %f - %f\nb",lp3[kk].pv,lp2[jj].pv); + break; + } + expand_gamut(s, pp); + } + } + } + +#ifdef STATS + tedges = edgestested = testhits = 0; +#endif + + /* Add any intersection points between img/dst, and src/dst */ + for (k = 0; k < 4; k++) { + int mask; + gamut *sa, *sb; + gtri *tpa, *tpb; /* Triangle pointers */ + + if (k == 0) { + if (!docomp) + continue; + mask = 3; + sa = s1; /* img */ + sb = s3; /* dst */ + } else if (k == 1) { + if (!docomp) + continue; + mask = 3; + sa = s3; /* dst */ + sb = s1; /* img */ + } else if (k == 2) { + if (!doexp) + continue; + mask = 12; + sa = s2; /* src */ + sb = s3; /* dst */ + } else if (k == 3) { + if (!doexp) + continue; + mask = 12; + sa = s3; /* dst */ + sb = s2; /* src */ + } + + /* Now find the edges that intersect the other gamut */ + tpa = sa->tris; + FOR_ALL_ITEMS(gtri, tpa) { + + for (j = 0; j < 3; j++) { /* For each edge */ +#ifdef STATS + tedges++; +#endif + /* If edge disposition flags aren't the same */ + if ((tpa->e[j]->v[0]->as ^ tpa->e[j]->v[1]->as) & mask + || (tpa->e[j]->v[0]->as & mask) == 3 || (tpa->e[j]->v[0]->as & mask) == 12 + || (tpa->e[j]->v[1]->as & mask) == 3 || (tpa->e[j]->v[1]->as & mask) == 12) { +#ifdef STATS + edgestested++; +#endif + /* Exhaustive search of other triangles */ + tpb = sb->tris; + FOR_ALL_ITEMS(gtri, tpb) { + double pv; + double pp[3]; + + /* Do a min/max intersection elimination test */ + for (i = 0; i < 3; i++) { + if (tpb->mix[1][i] < tpa->mix[0][i] + || tpb->mix[0][i] > tpa->mix[1][i]) + break; /* min/max don't overlap */ + } + if (i < 3) + continue; /* Skip this triangle, it can't intersect */ + + if (vect_intersect(sa, &pv, pp, tpa->e[j]->v[0]->p, tpa->e[j]->v[1]->p, tpb) + && pv >= (0.0 - 1e-10) && pv <= (1.0 + 1e-10)) { + /* Process intersection point pp the same as the first loop */ + double ppv, p2[3]; + double rr, r4; +#ifdef STATS + testhits++; +#endif +//printf("\n~1 tri isxn point %f %f %f\n", pp[0],pp[1],pp[2]); + if (cvect != NULL) + cvect(cntx, p2, pp); /* Get mapping direction */ + else + icmCpy3(p2, sa->cent); /* Radial vector */ + icmNormalize33(p2, p2, pp, 1.0); + + /* Locate the intersecting segments for each gamut */ + ll1 = ll2 = ll3 = 0; + ll1 = s1->vector_isectns(s1, pp, p2, lp1, MXNIS); /* Image gamut */ + if (doexp) /* For dst - src expansion */ + ll2 = s2->vector_isectns(s2, pp, p2, lp2, MXNIS); /* Src gamut */ + if (docomp || doexp) /* For img ^ dst or dst - src */ + ll3 = s3->vector_isectns(s3, pp, p2, lp3, MXNIS); /* Dst gamut */ + +//printf("~1 ll1 %d, ll2 %d, ll3 %d\n",ll1,ll2,ll3); +#ifdef NEVER + printf("img segments:\n"); + for (ii = 0; ii < ll1; ii++) + printf("Isect %d: pv %f, dir %d, edge %d, tri %d\n",ii,lp1[ii].pv,lp1[ii].dir,lp1[ii].edge,lp1[ii].tri->n); + printf("src segments:\n"); + for (ii = 0; ii < ll2; ii++) + printf("Isect %d: pv %f, dir %d, edge %d, tri %d\n",ii,lp2[ii].pv,lp2[ii].dir,lp2[ii].edge,lp2[ii].tri->n); + printf("dst segments:\n"); + for (ii = 0; ii < ll3; ii++) + printf("Isect %d: pv %f, dir %d, edge %d, tri %d\n",ii,lp3[ii].pv,lp3[ii].dir,lp3[ii].edge,lp3[ii].tri->n); +#endif + /* Now go through each image segment */ + for (ii = 0; ii < ll1; ii += 2) { + + icmCpy3(pp, lp1[ii].ip); + ppv = lp1[ii].pv; + +//printf("~1 img pnt at %f\n",lp1[ii].pv); + if (docomp) { + + /* Locate a dst segment around or below img point */ + for (kk = 0; kk < ll3; kk += 2) { + if ((lp1[kk+1].pv + 1e-8) >= ppv) + break; + } + + if (kk >= ll3) { /* No dst segments or none below */ +//printf("~1 img pnt is outside dst\n"); + continue; + } else { +//printf("~1 ing %f - %f, dst %f - %f\n", lp1[ii].pv,lp1[ii+1].pv,lp3[kk].pv,lp3[kk+1].pv); + /* Use the dst point if it is smaller */ + if (ppv < lp3[kk].pv) { + icmCpy3(pp, lp3[kk].ip); + ppv = lp3[kk].pv; +//printf("~1 dst is smaller %f\n",ppv); + } else { +//printf("~1 ing is smaller %f\n",ppv); + } + } + } + + if (nedst != NULL) + expand_gamut(nedst, pp); + + while (doexp) { + /* Locate a src segment that ing point lies in */ + for (jj = 0; jj < ll2; jj += 2) { + if (lp1[ii].pv >= (lp2[jj].pv - 1e-8) + && lp1[ii].pv <= (lp2[jj+1].pv + 1e-8)) + break; + } + if (jj >= ll2) { + break; /* No overlapping segment */ + } + + /* Locate a dst segment that src point lies in */ + for (kk = 0; kk < ll3; kk += 2) { + if (lp2[jj].pv >= (lp3[kk].pv - 1e-8) + && lp2[jj].pv <= (lp3[kk+1].pv + 1e-8)) + break; + } + if (kk >= ll2) { + break; /* No overlapping segment */ + } + + /* if src is outside dst, do nothing */ + if (lp3[kk].pv >= lp2[jj].pv) { + break; + } + + /* Expand point by dst - src */ + icmAdd3(pp, pp, lp3[kk].ip); + icmSub3(pp, pp, lp2[jj].ip); +//printf("~1 expanding point by %f - %f\nb",lp3[kk].pv,lp2[jj].pv); + break; + } + expand_gamut(s, pp); + } + } + } END_FOR_ALL_ITEMS(tpb); + } + } + + } END_FOR_ALL_ITEMS(tpa); + } + +#ifdef STATS + printf("Total edges %d, edges tested %d, edge hits %d\n", tedges, edgestested, testhits); +#endif + s->nofilter = 0; + + return 0; +} +#undef MXNIS + +/* ------------------------------------ */ +/* Locate the vertices most likely to correspond to the */ +/* primary and secondary colors (cusps) */ +/* + * Notes: + * + * To better support gamuts of devices with more than 4 colorants, + * it may be necessary to add another flag type that expands + * the cusps to lie on the actual gamut surface, as for + * some devices this lies outside the pure colorant combinations. + * + * Flinging a grid of values at this doesn't always + * return sensible results. Sometimes two "cusps" might be + * unreasonable close to each other (ie. one isn't a real cusp). + * This can cause gammut mapping to fail ... + * + * How could this be made more robust ? + * + */ + +static void setcusps(gamut *s, int flag, double in[3]) { + int i, j; + + /* - - - - - - - - - - - - - - - - - - - - - - */ + if (flag == 0) { /* Reset */ + for (j = 0; j < 6; j++) { + s->cusps[j][0] = 0.0; /* Marker values */ + s->cusps[j][1] = 0.0; + s->cusps[j][2] = 0.0; + } + s->dcuspixs = 0; + s->cu_inited = 0; + return; + + /* - - - - - - - - - - - - - - - - - - - - - - */ + } else if (flag == 2) { /* Finalize */ + + if (s->dcuspixs > 0) { + double JCh[3]; + double hues[6]; + int r, br = 0.0; + double berr; + + /* Figure out where to put the ones we got */ + for (j = 0; j < 6; j++) { /* Compute the hues */ + icmLab2LCh(JCh, s->dcusps[j]); + hues[j] = JCh[2]; + } + + /* Sort them into hue order */ + for (j = 0; j < 5; j++) { + for (i = j+1; i < 6; i++) { + if (hues[j] > hues[i]) { + double tt; + tt = hues[j]; hues[j] = hues[i]; hues[i] = tt; + tt = s->dcusps[j][0]; s->dcusps[j][0] = s->dcusps[i][0]; s->dcusps[i][0] = tt; + tt = s->dcusps[j][1]; s->dcusps[j][1] = s->dcusps[i][1]; s->dcusps[i][1] = tt; + tt = s->dcusps[j][2]; s->dcusps[j][2] = s->dcusps[i][2]; s->dcusps[i][2] = tt; + } + } + } + + /* Figure out which is the best match by rotation */ + berr = 1e6; + for (r = 0; r < 6; r++) { + double terr = 0.0; + for (j = 0; j < 6; j++) { + double tt; + + tt = fabs(gam_hues[s->isJab][j] - hues[(j + r) % 6]); + if (tt > 180.0) + tt = 360.0 - tt; + terr += tt; + } + if (terr < berr) { + br = r; + berr = terr; + } + } + /* Place them at that rotation */ + for (j = 0; j < 6; j++) { /* Compute the hues */ +//printf("~1 placing hue %f ix %d into hue %f ix %d\n", hues[(j + br) % 6],(j + br) % 6, gam_hues[s->isJab][j] ,j); + + s->cusps[j][0] = s->dcusps[(j + br) % 6][0]; + s->cusps[j][1] = s->dcusps[(j + br) % 6][1]; + s->cusps[j][2] = s->dcusps[(j + br) % 6][2]; + } + } + + /* Check we've got a cusp */ + for (j = 0; j < 6; j++) { + if (s->cusps[j][0] == 0.0 + && s->cusps[j][1] == 0.0 + && s->cusps[j][2] == 0.0) { + s->cu_inited = 0; + return; /* Not all have been set */ + } + } + + { + double JCh[3]; + double hues[6]; + + /* Check how far appart the cusps are in hue angle */ + // ~~999 + for (j = 0; j < 6; j++) { + icmLab2LCh(JCh, s->cusps[j]); + hues[j] = JCh[2]; +//printf("~1 cusp %d = hue %f\n",j,hues[j]); + } + for (j = 0; j < 6; j++) { + int k = j < 5 ? j + 1 : 0; + double rh, h; + rh = gam_hues[s->isJab][k] - gam_hues[s->isJab][j]; + if (rh < 0.0) + rh = 360 + rh; + h = hues[k] - hues[j]; + if (h < 0.0) + h = 360 + h; +//printf("~1 cusp %d - %d = ref dh %f, dh %f\n",j,k,rh,h); + + /* if our delta is less than half reference, */ + /* assume the cusps are bad. */ + if ((2.0 * h) < rh) { + + s->cu_inited = 0; /* Not trustworthy */ +//printf("~1 cusps are not trustworthy\n"); + return; + } + } + } + + s->cu_inited = 1; + return; + + /* - - - - - - - - - - - - - - - - - - - - - - */ + } else if (flag == 3) { /* Definite 1/6 cusp */ + + if (s->dcuspixs >= 6) { +//printf("~1 too many cusp values added\n"); + return; /* Error - extra cusp ignored */ + } + s->dcusps[s->dcuspixs][0] = in[0]; + s->dcusps[s->dcuspixs][1] = in[1]; + s->dcusps[s->dcuspixs++][2] = in[2]; + + } else { /* Consider another point as a cusp point */ + double JCh[3]; + int bj = 0, sbj = 0; + double bh = 1e6, sbh = 1e6; + double ns, es; + + icmLab2LCh(JCh, in); + +//printf("~1 cusp at %f %f %f\n",JCh[0],JCh[1],JCh[2]); + /* See which hue it is closest and 2nd closet to cusp hue. */ + for (j = 0; j < 6; j++) { + double tt; + + tt = fabs(gam_hues[s->isJab][j] - JCh[2]); + if (tt > 180.0) + tt = 360.0 - tt; + + if (tt < bh) { + if (bh < sbh) { + sbh = bh; + sbj = bj; + } + bh = tt; + bj = j; + } else if (tt < sbh) { + sbh = tt; + sbj = j; + } + } + + /* Compute distance of existing and new */ + es = s->cusps[bj][1] * s->cusps[bj][1] + s->cusps[bj][2] * s->cusps[bj][2]; + ns = in[1] * in[1] + in[2] * in[2]; +//printf("~1 chroma dist of existing %f, new %f\n",es,ns); + if (ns > es) { +//printf("~1 New closest\n"); + s->cusps[bj][0] = in[0]; + s->cusps[bj][1] = in[1]; + s->cusps[bj][2] = in[2]; + + } else { /* If 2nd best has no entry, use this to fill it */ + if (s->cusps[sbj][0] == 0.0 + && s->cusps[sbj][1] == 0.0 + && s->cusps[sbj][2] == 0.0) { +//printf("~1 Fill with 2nd closest\n"); + s->cusps[sbj][0] = in[0]; + s->cusps[sbj][1] = in[1]; + s->cusps[sbj][2] = in[2]; + } + } + } +} + + +/* Get the cusp values for red, yellow, green, cyan, blue & magenta */ +/* Return nz if there are no cusps available */ +static int getcusps( +gamut *s, +double cusps[6][3] +) { + int i, j; + + if (s->cu_inited == 0) { + return 1; + } + + for (i = 0; i < 6; i++) + for (j = 0; j < 3; j++) + cusps[i][j] = s->cusps[i][j]; + + return 0; +} + +/* ===================================================== */ +/* Triangulation code */ +/* ===================================================== */ + +static double ne_point_on_tri(gamut *s, gtri *t, double *out, double *in); + +/* Given three points, compute the normalised plane equation */ +/* of a surface through them. */ +/* Return non-zero on error */ +static int plane_equation( +double *eq, /* Return equation parameters */ +double *p0, /* The three points */ +double *p1, +double *p2 +) { + double ll, v1[3], v2[3]; + + /* Compute vectors along edges */ + v1[0] = p1[0] - p0[0]; + v1[1] = p1[1] - p0[1]; + v1[2] = p1[2] - p0[2]; + + v2[0] = p2[0] - p0[0]; + v2[1] = p2[1] - p0[1]; + v2[2] = p2[2] - p0[2]; + + /* Compute cross products v1 x v2, which will be the normal */ + eq[0] = v1[1] * v2[2] - v1[2] * v2[1]; + eq[1] = v1[2] * v2[0] - v1[0] * v2[2]; + eq[2] = v1[0] * v2[1] - v1[1] * v2[0]; + + /* Normalise the equation */ + ll = sqrt(eq[0] * eq[0] + eq[1] * eq[1] + eq[2] * eq[2]); + if (ll < 1e-10) { + return 1; + } + eq[0] /= ll; + eq[1] /= ll; + eq[2] /= ll; + + /* Compute the plane equation constant */ + eq[3] = - (eq[0] * p0[0]) + - (eq[1] * p0[1]) + - (eq[2] * p0[2]); + +#ifdef NEVER + /* Veritify the plane equation */ + { + double c; + c = eq[0] * p0[0] + + eq[1] * p0[1] + + eq[2] * p0[2] + + eq[3]; + if (fabs(c) > 1e-10) { + printf("Plane equation check 0 failed by %f\n",c); + } + c = eq[0] * p1[0] + + eq[1] * p1[1] + + eq[2] * p1[2] + + eq[3]; + if (fabs(c) > 1e-10) { + printf("Plane equation check 1 failed by %f\n",c); + } + c = eq[0] * p2[0] + + eq[1] * p2[1] + + eq[2] * p2[2] + + eq[3]; + if (fabs(c) > 1e-10) { + printf("Plane equation check 2 failed by %f\n",c); + } + } +#endif /* NEVER */ + return 0; +} + +/* Compute the log surface plane equation for the triangle */ +/* and other triangle attributes. (Doesn't depend on edge info.) */ +void +comptriattr( +gamut *s, +gtri *t +) { + int j; + static double v0[3] = {0.0, 0.0, 0.0}; + double cp[3]; /* Closest point - not used */ + + /* Compute the plane equation for the absolute triangle. */ + /* This is used for testing if a point is inside the gamut hull. */ + plane_equation(t->pe, t->v[0]->p, t->v[1]->p, t->v[2]->p); + + /* Compute the plane equation for the triangle */ + /* based on the log compressed convex hull verticies. */ + /* This is used for convex hull construction. */ + plane_equation(t->che, t->v[0]->ch, t->v[1]->ch, t->v[2]->ch); + + /* Compute the plane equation for the triangle */ + /* mapped to the surface of the sphere */ + /* This can be used for point in triangle testing ?? */ + plane_equation(t->spe, t->v[0]->sp, t->v[1]->sp, t->v[2]->sp); + + /* Compute the plane equations of the spherical mapped vertex */ + /* values with regard to the center of the sphere, so that */ + /* a point in triangle test can be performed, and baricentric, */ + /* coordinates can be computed. */ + plane_equation(t->ee[0], v0, t->v[1]->sp, t->v[2]->sp); + plane_equation(t->ee[1], v0, t->v[2]->sp, t->v[0]->sp); + plane_equation(t->ee[2], v0, t->v[0]->sp, t->v[1]->sp); + + /* Compute the radius range of the triangle to the center */ + /* Compute the maximum from the vertexes */ + t->rs1 = -1.0; + for (j = 0; j < 3; j++) { + int k; + double rs, tt; + for (rs = 0.0, k = 0;k < 3; k++) { + tt = t->v[j]->p[k] - s->cent[k]; + rs += tt * tt; + } + if (rs > t->rs1) + t->rs1 = rs; + } + /* The minimum may be on the plane, an edge or a vertex, */ + /* so use closest point in triangle function. */ + t->rs0 = ne_point_on_tri(s, t, cp, s->cent); + + /* Allow a tollerance around the radius squareds */ + t->rs0 -= 1e-4; + t->rs1 += 1e-4; + +#ifdef NEVER // ??? +#ifdef ASSERTS + { + double tt[3]; /* Triangle test point */ + double ds; + for (j = 0; j < 3; j++) { + tt[j] = (t->v[0]->p[j] + t->v[1]->p[j] + t->v[2]->p[j])/3.0; + tt[j] -= s->cent[j]; /* Make it center relative */ + } + for (j = 0; j < 3; j++) { + ds = t->ee[j][0] * tt[0] /* Point we know is inside */ + + t->ee[j][1] * tt[1] + + t->ee[j][2] * tt[2] + + t->ee[j][3]; + if (ds > 1e-8) + break; /* Not within triangle */ + } + if (j < 3) { + fprintf(stderr,"Assert: point expected to be within triangle %d (vx %d %d %d) is not\n", + t->n, t->v[0]->n, t->v[1]->n, t->v[2]->n); + fprintf(stderr,"Known point is %f, expect -ve\n",ds); + exit(-1); + } + } +#endif /* ASSERTS */ +#endif /* NEVER */ + +} + +/* By using the pow() or log() of the radial distance, */ +/* blended with a sphere surface, we try and strike a compromise */ +/* between a pure convex hull surface, and a pure Delaunay triangulation */ +/* the latter which would show dings and nicks from points */ +/* that fall withing the "real" gamut. */ +static double log_scale(gamut *s, double rr) { + double aa; + +#ifdef TEST_CONVEX_HULL + return rr; +#else +#ifdef NEVER /* (Not using this version, doesn't work reliably) */ + aa = (2.0 + rr)/3.0; /* Blend with sphere */ + aa = log(aa); /* Allow for concave slope */ + if (aa < 0.0) /* but constrain to be +ve */ + aa = 0.0; +#else /* (Using simpler version) */ + aa = 20.0 * pow(rr, s->logpow); /* Default 0.25 */ +#endif +#endif /* TEST_CONVEX_HULL */ + + return aa; +} + +/* Comput r[], lr0, sp[] and ch[] from p[] for all verticies */ +/* (Note that lr0 will be the first cut, log_scale() value) */ +static void +compute_vertex_coords( +gamut *s +) { + int i, j; + + for (i = 0; i < s->nv; i++) { + gamut_rect2radial(s, s->verts[i]->r, s->verts[i]->p); + + if (s->verts[i]->r[0] < 1e-6) { /* Ignore a point right at the center */ + s->verts[i]->lr0 = 0.0; + for (j = 0; j < 3; j++) { + s->verts[i]->sp[j] = 0.0; + s->verts[i]->ch[j] = 0.0; + } + } else { + double aa; + + /* Figure log scaled radius */ + s->verts[i]->lr0 = log_scale(s, s->verts[i]->r[0]); + + /* Compute unit shere mapped location */ + aa = 1.0/s->verts[i]->r[0]; /* Adjustment to put in on unit sphere */ + for (j = 0; j < 3; j++) + s->verts[i]->sp[j] = (s->verts[i]->p[j] - s->cent[j]) * aa; + + /* Compute hull testing mapped version */ + for (j = 0; j < 3; j++) + s->verts[i]->ch[j] = s->verts[i]->p[j] * s->verts[i]->lr0; + } + } +} + +/* Sort the verticies from maximum radius */ +static void sort_verticies( +gamut *s +) { + int i; + +#ifndef TEST_NOSORT + + /* Sort them */ +//#define HEAP_COMPARE(A,B) (A->r[0] > B->r[0]) +#define HEAP_COMPARE(A,B) (A->lr0 > B->lr0) + HEAPSORT(gvert *, s->verts, s->nv) +#undef HEAP_COMPARE + +#endif /* !TEST_NOSORT */ + + /* Renumber them */ + for (i = 0; i < s->nv; i++) { + s->verts[i]->n = i; + } +} + +/* Number just the verticies that have been set, */ +/* and those that have been used in the convex hull */ +static void renumber_verticies( +gamut *s +) { + int i, j; + + for (j = i = 0; i < s->nv; i++) { + if (!(s->verts[i]->f & GVERT_SET)) + continue; + + s->verts[i]->sn = j; + j++; + } + s->nsv = j; + + for (j = i = 0; i < s->nv; i++) { + if (!(s->verts[i]->f & GVERT_TRI)) + continue; + + s->verts[i]->tn = j; + j++; + } + s->ntv = j; +} + +#ifdef ASSERTS + +/* Diagnpostic aid */ +/* Check that the triangulation adjacenty info is OK */ +static void check_triangulation(gamut *s, int final) { + int i, j; + gtri *tp; /* Triangle pointer */ + gedge *ep; /* Edge pointer */ + int failed = 0; + + tp = s->tris; + FOR_ALL_ITEMS(gtri, tp) { + + /* Check verticies for duplication */ + for (i = 0; i < 3; i++) { + for (j = i+1; j < 3; j++) { + if (tp->v[i] == tp->v[j]) { + failed = 1; + printf("Validation failed - duplicate verticies:\n"); + printf("Triangle %d, has verticies %d %d %d\n", tp->n, tp->v[0]->n, tp->v[1]->n, tp->v[2]->n); + fflush(stdout); + } + } + } + + /* Check edges for duplication */ + for (i = 0; i < 3; i++) { + for (j = i+1; j < 3; j++) { + if (tp->e[i] == tp->e[j]) { + failed = 1; + printf("Validation failed - duplicate connectivity:\n"); + printf("Triangle %d, has verticies %d %d %d\n", tp->n, tp->v[0]->n, tp->v[1]->n, tp->v[2]->n); + printf("Triangle %d, has edges %d %d %d\n", tp->n, tp->e[0]->n, tp->e[1]->n, tp->e[2]->n); + fflush(stdout); + } + } + } + + /* Check connectivity */ + for (i = 0; i < 3; i++) { + gtri *t1, *t2; + gedge *e; + int ei1, ei2; + int tei; /* Edges index for this triangle [0..1] */ + + e = tp->e[i]; /* The edge in question */ + tei = tp->ei[i]; + ei1 = e->ti[tei]; /* Edges record of edge index withing this triangle */ + + /* Check that the edges reconing of what index edge it is */ + /* for this triangle is correct */ + if (ei1 != i) { + failed = 1; + printf("Validation failed - triangle edge index doesn't match record withing edge:\n"); + printf("Triangle %d, edge index %d edge %d has record %d\n", tp->n, i, e->n, ei1); + fflush(stdout); + } + + /* Check that the edges pointer to the triangle is this triangle */ + if (tp != e->t[tei]) { + failed = 1; + printf("Validation failed - edge doesn't point back to triangle:\n"); + printf("Triangle %d, edge index %d is edge %d\n",tp->n, i, e->n); + printf("Edge %d, triangle index %d is triangle %d\n", e->n, tei, e->t[tei]->n); + printf("Edge %d, triangle index %d is triangle %d\n", e->n, tei^1, e->t[tei^1]->n); + fflush(stdout); + } + + /* Check the verticies for this edge match edge record */ + if ((e->v[0] != tp->v[i] || e->v[1] != tp->v[(i+1) % 3]) + && (e->v[1] != tp->v[i] || e->v[0] != tp->v[(i+1) % 3])) { + failed = 1; + printf("Validation failed - edge doesn't have same verticies as triangle expects:\n"); + printf("Triangle %d, has verticies %d %d\n", tp->n, tp->v[i]->n, tp->v[(i+1) % 3]->n); + printf("Edge %d, has verticies %d %d\n", e->n, e->v[0]->n, e->v[1]->n); + fflush(stdout); + } + + t2 = e->t[tei ^ 1]; /* The other triangle */ + ei2 = e->ti[tei ^ 1]; /* Edges index number withing triangle t2 */ + + if (t2 == tp) { + failed = 1; + printf("Validation failed - connects to itself:\n"); + printf("Triangle %d, has edges %d %d %d\n", tp->n, tp->e[0]->n, tp->e[1]->n, tp->e[2]->n); + fflush(stdout); + } + + /* Check that the connection is reflective */ + if (e != t2->e[ei2]) { + failed = 1; + printf("Validation failed - connectivity not reflected:\n"); + printf("Triangle %d, edge index %d points to edge %d\n",tp->n, i, e->n); + printf("Triangle %d, edge index %d points to edge %d\n",t2->n, ei2, t2->e[ei2]->n); + fflush(stdout); + } + } + } END_FOR_ALL_ITEMS(tp); + if (failed) { + exit(-1); + } + + /* - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + /* Check that every point is part of a triangle and edge */ + + for (i = 0; i < s->nv; i++) { /* Reset the assert flag */ + gvert *v = s->verts[i]; + + v->as = 0; + + /* Check out the flags */ + if (!(v->f & GVERT_SET)) { + if ((v->f & GVERT_TRI) + || (v->f & GVERT_INSIDE)) { + printf("Validation failed - vertex %d has strange flags 0x%x\n",i, v->f); + fflush(stdout); + failed = 1; + } + } else { + if ((v->f & GVERT_TRI) && (v->f & GVERT_INSIDE)) { + printf("Validation failed - vertex %d has strange flags 0x%x\n",i, v->f); + fflush(stdout); + failed = 1; + } + } + } + + FOR_ALL_ITEMS(gtri, tp) { + for (i = 0; i < 3; i++) + tp->v[i]->as |= 1; /* Vertex is in a triangle */ + } END_FOR_ALL_ITEMS(tp); + + ep = s->edges; + FOR_ALL_ITEMS(gedge, ep) { + ep->v[0]->as |= 2; /* Vertex is in an edge */ + ep->v[1]->as |= 2; + ep->as = 0; /* Reset the assert flag */ + } END_FOR_ALL_ITEMS(ep); + + for (i = 0; i < s->nv; i++) { + if (s->verts[i]->f & GVERT_TRI) { + if ((s->verts[i]->as & 1) == 0) { + printf("Validation failed - vertex %d is not in any triangles\n",i); + fflush(stdout); + failed = 1; + } + if ((s->verts[i]->as & 2) == 0) { + printf("Validation failed - vertex %d is not in any edge\n",i); + fflush(stdout); + failed = 1; + } + } + } + + + /* - - - - - - - - - - - - - - - - - - - - - - */ + /* Check that every edge is part of a triangle */ + + /* as flag in triangle was reset above */ + FOR_ALL_ITEMS(gtri, tp) { + for (i = 0; i < 3; i++) + tp->e[i]->as |= 1; /* Mark edge used in triangle */ + } END_FOR_ALL_ITEMS(tp); + + ep = s->edges; + FOR_ALL_ITEMS(gedge, ep) { + if (ep->as != 1) { + printf("Validation failed - edge %d is not in any triangle\n",ep->n); + fflush(stdout); + failed = 1; + } + } END_FOR_ALL_ITEMS(ep); + + if (failed) { + exit(-1); + } +} + +#endif /* ASSERTS */ + +/* -------------------------------------- */ +/* Add a face to the hit list, if it is not a duplicate. */ +static void add_to_hit_list( +gamut *s, +gtri **hlp, /* Hit list */ +gtri *cf /* Face to be added (triangle verts 0, 1) */ +) { + gtri *tp; /* Triangle pointer */ + gvert *c0 = cf->v[0]; + gvert *c1 = cf->v[1]; + +//printf("Adding face to hit list %d: %d %d\n", +//cf->n, cf->v[0]->n, cf->v[1]->n); + + tp = *hlp; + /* Search current faces in hit list */ + FOR_ALL_ITEMS(gtri, tp) { + gvert *v0 = tp->v[0]; + gvert *v1 = tp->v[1]; + if ((c0 == v0 && c1 == v1) /* Same face from other side */ + || (c0 == v1 && c1 == v0)) { + /* Duplicate found */ +//printf("Duplicate found %d: %d %d\n", +//tp->n, tp->v[0]->n, tp->v[1]->n); + DEL_LINK(*hlp, tp); /* Delete from the hit list */ + + /* Check face is common */ + if (cf->e[0] != tp->e[0]) { + fprintf(stderr,"gamut: internal error - face match inconsistency\n"); + exit(-1); + } + /* Delete edge */ + DEL_LINK(s->edges, cf->e[0]); + del_gedge(cf->e[0]); + + /* Delete the two faces */ + del_gtri(tp); + del_gtri(cf); + return; + } + } END_FOR_ALL_ITEMS(tp); + + /* Safe to add it to face hit list */ + /* This removes triangle from triangles list ? */ + ADD_ITEM_TO_BOT(*hlp, cf); +//printf("Face added\n"); +} + +/* Add a triangles faces to the hit list. */ +static void add_tri_to_hit_list( +gamut *s, +gtri **hlp, /* Hit list */ +gtri *tp /* Triangle faces to be added */ +) { + int j; + gtri *t1, *t2; + + /* In case some verticies disapear below the log surface, */ + /* and don't remain part of the triangulation, we mark them off. */ + for (j = 0; j < 3 ; j++) { + tp->v[j]->f &= ~GVERT_TRI; + tp->v[j]->f |= GVERT_INSIDE; + } + + /* Decompose the triangle into three faces, each face being stored */ + /* into a triangle created on the hit list, using verticices 0, 1. */ + /* The edges adjacency info remains valid for the three faces, */ + /* as does the edge plane equation. */ + DEL_LINK(s->tris, tp); /* Delete it from the triangulation list */ + t1 = new_gtri(); + t1->v[0] = tp->v[1]; /* Duplicate with rotated faces */ + t1->v[1] = tp->v[2]; + t1->e[0] = tp->e[1]; /* Edge adjacency for this edge */ + t1->ei[0] = tp->ei[1]; /* Edge index of this triangle */ + t1->e[0]->t[t1->ei[0]] = t1; /* Fixup reverse adjacency for valid edge */ + t1->e[0]->ti[t1->ei[0]] = 0; /* Rotated index of new triangles edge */ + t1->e[1] = t1->e[2] = NULL; /* be safe */ + for (j = 0; j < 4; j++) /* Copy edge plane equation */ + t1->ee[2][j] = tp->ee[0][j]; + + t2 = new_gtri(); + t2->v[0] = tp->v[2]; /* Duplicate with rotated faces */ + t2->v[1] = tp->v[0]; + t2->e[0] = tp->e[2]; /* Edge adjacency for this edge */ + t2->ei[0] = tp->ei[2]; /* Edge index of this triangle */ + t2->e[0]->t[t2->ei[0]] = t2; /* Fixup reverse adjacency for valid edge */ + t2->e[0]->ti[t2->ei[0]] = 0; /* Rotated index of new triangles edge */ + t2->e[1] = t2->e[2] = NULL; /* be safe */ + for (j = 0; j < 4; j++) /* Copy edge plane equation */ + t2->ee[2][j] = tp->ee[1][j]; + + tp->e[1] = tp->e[2] = NULL; /* be safe */ + add_to_hit_list(s, hlp, tp); /* Add edge 0 to hit list as is */ + add_to_hit_list(s, hlp, t1); /* Add edge 1 to hit list */ + add_to_hit_list(s, hlp, t2); /* Add edge 2 to hit list */ +} + +#ifdef DEBUG_TRIANG + typedef struct { + int tix[3]; /* Triangle indexes */ + int type; /* 0 = hit, 1 = added */ + } tidxs; +#endif + +/* Insert a vertex into the triangulation */ +static void insert_vertex( +gamut *s, +gvert *v /* Vertex to insert */ +) { + gtri *tp, *tp2; /* Triangle pointers */ + gtri *hl; /* Triangle face hit list (polygon faces) */ + double tol = TRIANG_TOL; + int hit = 0; /* Vertex expands hull flag */ +#ifdef DEBUG_TRIANG + int intri = 0; /* Vertex landed in a triangle */ + XLIST(tidxs, hittris) + tidxs xxs; + + XLIST_INIT(tidxs, &hittris); +#endif + +#ifdef DEBUG_TRIANG + printf("Adding vertex %d: %f %f %f to triangles\n", v->n, v->p[0], v->p[1], v->p[2]); +#endif + + /* First we search the current triangles, and convert */ + /* any trianges that are visible from the new point, */ + /* into a list of faces stored on the face */ + /* hit list. */ + /* We are using a brute force search, which will make the */ + /* algorithm speed proportional to n^2. For better performance */ + /* with a large number of verticies, an acceleration structure */ + /* should be used to speed circumradius hit detection. */ + v->f &= ~GVERT_INSIDE; /* Reset flags */ + v->f &= ~GVERT_TRI; + INIT_LIST(hl); + hit = 0; + tp = s->tris; + FOR_ALL_ITEMS(gtri, tp) { + double c; + + /* Check the depth out compared to this triangle log plane equation */ + c = tp->che[0] * v->ch[0] + + tp->che[1] * v->ch[1] + + tp->che[2] * v->ch[2] + + tp->che[3]; + + /* If vertex is above the log hull surface, add triangle to the hit list. */ + if (c < -tol) { +#ifdef DEBUG_TRIANG + int j; + double bds = -1e10; +#endif + hit = 1; + +#ifdef DEBUG_TRIANG + printf("Got a hit on triangle %d: %d %d %d by %f\n", + tp->n, tp->v[0]->n, tp->v[1]->n, tp->v[2]->n,c); +#endif + +#ifdef DEBUG_TRIANG + for (j = 0; j < 3; j++) { + double ds; + ds = tp->ee[j][0] * v->ch[0] + + tp->ee[j][1] * v->ch[1] + + tp->ee[j][2] * v->ch[2] + + tp->ee[j][3]; + if (ds > tol) { + printf("Vertex is not in triangle by %e\n",ds); + break; + } + if (ds > bds) + bds = ds; + } + if (j >= 3) { + printf("Vertex is in triangle by %e\n",bds); + intri = 1; /* Landed in this triangle */ + } + + xxs.tix[0] = tp->v[0]->n, xxs.tix[1] = tp->v[1]->n, xxs.tix[2] = tp->v[2]->n; + xxs.type = 0; + XLIST_ADD(&hittris, xxs) +#endif + add_tri_to_hit_list(s, &hl, tp); + } + } END_FOR_ALL_ITEMS(tp); + + if (hit == 0) { + +//printf("No hits - must be inside the log hull\n"); + v->f |= GVERT_INSIDE; /* This point is inside the log hull */ + v->f &= ~GVERT_TRI; + } else { + int changed = 1; + +#ifdef DEBUG_TRIANG + /* Point doesn't lie radially within any of the triangles it is */ + /* above the plane of. This is a geometric conundrum. (?) */ +if (!intri) printf("~1 ###### vertex didn't land in any triangle! ########\n"); +#endif + +//printf("Checking out hit polygon edges:\n"); + /* Now we must make a pass though the hit list, checking that each */ + /* hit list face will make a correctly oriented, non-sliver triangle */ + /* when joined to the vertex. */ + /* Do this check until there are no changes */ + for (;changed != 0 ;) { + tp = hl; + changed = 0; + FOR_ALL_ITEMS(gtri, tp) { + /* Check which side of the edge our vertex is */ + double ds; + ds = tp->ee[2][0] * v->ch[0] + + tp->ee[2][1] * v->ch[1] + + tp->ee[2][2] * v->ch[2] + + tp->ee[2][3]; +//printf("Vertex margin to edge = %e\n",ds); + /* If vertex is not to the right of this edge by tol */ + /* add associated triangle to the hit list. */ + if (ds > -tol) { + gtri *xtp; +//printf("~1 ###### vertex on wrong side by %e - expand hit list ######\n",ds); + if (tp->e[0]->t[0] != tp) + xtp = tp->e[0]->t[0]; + else + xtp = tp->e[0]->t[1]; +//printf("Got a hit on triangle %d: %d %d %d\n", xtp->n, xtp->v[0]->n, xtp->v[1]->n, xtp->v[2]->n); + +#ifdef DEBUG_TRIANG + xxs.tix[0] = xtp->v[0]->n, xxs.tix[1] = xtp->v[1]->n, xxs.tix[2] = xtp->v[2]->n; + xxs.type = 1; + XLIST_ADD(&hittris, xxs) +#endif + + add_tri_to_hit_list(s, &hl, xtp); + changed = 1; + break; + } + } END_FOR_ALL_ITEMS(tp); + } + +#ifdef DEBUG_TRIANG_VRML + write_diag_vrml(s, v->ch, hittris.no, hittris.list, hl); +#endif + +//printf("About to turn polygon faces into triangles\n"); + /* Turn all the faces that made it to the */ + /* hit list, into triangles using the new vertex. */ + tp = hl; + FOR_ALL_ITEMS(gtri, tp) { + tp->v[2] = v; /* Add third vertex to face to make triangle */ + comptriattr(s, tp); /* Compute triangle attributes */ + + /* Find the new adjacent triangles from the triangles being formed, */ + /* to maintain edge adjacency information. */ + /* Do only one edge at a time, since each adjacency */ + /* will be visited twice. */ + tp2 = hl; + FOR_ALL_ITEMS(gtri, tp2) { + if (tp2->v[0] == tp->v[1]) { /* Found 1/2 tp/tp2 edge adjacency */ + gedge *e; + e = new_gedge(); + ADD_ITEM_TO_BOT(s->edges, e); /* Append to edge list */ + tp->e[1] = e; /* Point to edge */ + tp->ei[1] = 0; /* edges 0th triangle */ + e->t[0] = tp; /* triangles 1st edge */ + e->ti[0] = 1; /* triangles 1st edge */ + tp2->e[2] = e; /* Point to edge */ + tp2->ei[2] = 1; /* edges 1st triangle */ + e->t[1] = tp2; /* Triangles 2nd edge */ + e->ti[1] = 2; /* Triangles 2nd edge */ + e->v[0] = v; /* Add the two verticies */ + e->v[1] = tp->v[1]; + } + } END_FOR_ALL_ITEMS(tp2); + +//printf("~1 Creating new triangle %d: %d %d %d\n", tp->n, tp->v[0]->n, tp->v[1]->n, tp->v[2]->n); + } END_FOR_ALL_ITEMS(tp); + +#ifdef DEBUG_TRIANG_VRML + tp = hl; + hittris.no = 0; + FOR_ALL_ITEMS(gtri, tp) { + xxs.tix[0] = tp->v[0]->n, xxs.tix[1] = tp->v[1]->n, xxs.tix[2] = tp->v[2]->n; + xxs.type = 2; + XLIST_ADD(&hittris, xxs) + } END_FOR_ALL_ITEMS(tp); + write_diag_vrml(s, v->ch, hittris.no, hittris.list, NULL); +#ifdef DEBUG_TRIANG_VRML_STEP +#ifdef DO_TWOPASS + if (s->pass > 0) +#endif /* DO_TWOPASS */ + getchar(); +#endif +#endif + + /* Move them to the triangulation. */ + tp = hl; + FOR_ALL_ITEMS(gtri, tp) { + int j; + DEL_LINK(hl, tp); /* Gone from the hit list */ + ADD_ITEM_TO_BOT(s->tris, tp); /* Append to triangulation list */ + for (j = 0; j < 3 ; j++) { /* Verticies weren't dropped from triangulation */ + tp->v[j]->f |= GVERT_TRI; + tp->v[j]->f &= ~GVERT_INSIDE; + } + } END_FOR_ALL_ITEMS(tp); + + v->f |= GVERT_TRI; /* This vertex has been added to triangulation */ + v->f &= ~GVERT_INSIDE; /* and it's not inside */ + } + +#ifdef DEBUG_TRIANG + XLIST_FREE(&hittris); +#endif +} + +/* - - - - - - - - - - - - - - - - */ + +/* Create the convex hull surface triangulation */ +static void triangulate_ch( +gamut *s +) { + /* Establish the base triangulation */ + { + int i, j; + gvert *tvs[4]; /* Initial verticies */ + gtri *tr[4]; /* Initial triangles */ + gedge *ed[6]; /* Initial edges */ + double fsz = FAKE_SEED_SIZE; /* Initial tetra size */ + double ff[3]; + int onf; + static double foffs[4][3] = { /* tetrahedral offsets */ + { 0.0, 0.0, 1.0 }, + { 0.0, 0.80254, -0.5 }, + { 0.75, -0.330127, -0.5 }, + { -0.75, -0.330127, -0.5 } + }; + + /* Delete any current fake verticies */ + for (j = i = 0; i < s->nv; i++) { + s->verts[i]->f &= ~GVERT_ESTP; /* Unmark non-fake establishment points */ + if (!(s->verts[i]->f & GVERT_FAKE)) + s->verts[j++] = s->verts[i]; + else + free(s->verts[i]); + } + s->nv = j; + + /* Re-create fake points on each pass */ + onf = s->nofilter; + s->nofilter = 1; /* Turn off filtering */ + s->doingfake = 1; /* Adding fake points */ + + for (j = i = 0; i < 4; i++) { + ff[0] = fsz * foffs[i][2] + s->cent[0]; + ff[1] = fsz * foffs[i][0] + s->cent[1]; + ff[2] = fsz * foffs[i][1] + s->cent[2]; + if ((tvs[j++] = expand_gamut(s, ff)) == NULL) { + fprintf(stderr,"gamut: internal error - failed to register a fake initial verticies!\n"); + exit (-1); + } + } + + s->nofilter = onf; + s->doingfake = 0; + +#ifdef NEVER + printf("Initial verticies:\n"); + for (i = 0; i < 4; i++) { + printf(" %d: %f %f %f\n",tvs[i]->n, tvs[i]->p[0], tvs[i]->p[1], tvs[i]->p[2]); + } +#endif + /* Setup the initial triangulation */ + for (i = 0; i < 4; i++) { + tr[i] = new_gtri(); + } + + for (i = 0; i < 6; i++) { + ed[i] = new_gedge(); + ADD_ITEM_TO_BOT(s->edges, ed[i]); + } + + /* Enter the edge verticies */ + ed[0]->v[0] = tvs[0]; + ed[0]->v[1] = tvs[1]; + ed[1]->v[0] = tvs[1]; + ed[1]->v[1] = tvs[2]; + ed[2]->v[0] = tvs[0]; + ed[2]->v[1] = tvs[2]; + ed[3]->v[0] = tvs[0]; + ed[3]->v[1] = tvs[3]; + ed[4]->v[0] = tvs[1]; + ed[4]->v[1] = tvs[3]; + ed[5]->v[0] = tvs[2]; + ed[5]->v[1] = tvs[3]; + + /* Triangle facing in the +x, +y +z direction */ + tr[0]->v[0] = tvs[0]; + tr[0]->v[1] = tvs[1]; + tr[0]->v[2] = tvs[2]; + + tr[0]->e[0] = ed[0]; /* Should make edge joining a function ? */ + tr[0]->ei[0] = 0; + ed[0]->t[0] = tr[0]; + ed[0]->ti[0] = 0; + + tr[0]->e[1] = ed[1]; + tr[0]->ei[1] = 0; + ed[1]->t[0] = tr[0]; + ed[1]->ti[0] = 1; + + tr[0]->e[2] = ed[2]; + tr[0]->ei[2] = 0; + ed[2]->t[0] = tr[0]; + ed[2]->ti[0] = 2; + + comptriattr(s, tr[0]); /* Compute triangle attributes */ + ADD_ITEM_TO_BOT(s->tris, tr[0]); /* Append to list */ + + /* Triangle facing in the -x, +y +z direction */ + tr[1]->v[0] = tvs[0]; + tr[1]->v[1] = tvs[3]; + tr[1]->v[2] = tvs[1]; + + tr[1]->e[0] = ed[3]; + tr[1]->ei[0] = 0; + ed[3]->t[0] = tr[1]; + ed[3]->ti[0] = 0; + + tr[1]->e[1] = ed[4]; + tr[1]->ei[1] = 0; + ed[4]->t[0] = tr[1]; + ed[4]->ti[0] = 1; + + tr[1]->e[2] = ed[0]; + tr[1]->ei[2] = 1; + ed[0]->t[1] = tr[1]; + ed[0]->ti[1] = 2; + + comptriattr(s, tr[1]); /* Compute triangle attributes */ + ADD_ITEM_TO_BOT(s->tris, tr[1]); /* Append to list */ + + /* Triangle facing in the -y +z direction */ + tr[2]->v[0] = tvs[0]; + tr[2]->v[1] = tvs[2]; + tr[2]->v[2] = tvs[3]; + + tr[2]->e[0] = ed[2]; + tr[2]->ei[0] = 1; + ed[2]->t[1] = tr[2]; + ed[2]->ti[1] = 0; + + tr[2]->e[1] = ed[5]; + tr[2]->ei[1] = 0; + ed[5]->t[0] = tr[2]; + ed[5]->ti[0] = 1; + + tr[2]->e[2] = ed[3]; + tr[2]->ei[2] = 1; + ed[3]->t[1] = tr[2]; + ed[3]->ti[1] = 2; + + comptriattr(s, tr[2]); /* Compute triangle attributes */ + ADD_ITEM_TO_BOT(s->tris, tr[2]); /* Append to list */ + + /* Triangle facing in the -z direction */ + tr[3]->v[0] = tvs[1]; + tr[3]->v[1] = tvs[3]; + tr[3]->v[2] = tvs[2]; + + tr[3]->e[0] = ed[4]; + tr[3]->ei[0] = 1; + ed[4]->t[1] = tr[3]; + ed[4]->ti[1] = 0; + + tr[3]->e[1] = ed[5]; + tr[3]->ei[1] = 1; + ed[5]->t[1] = tr[3]; + ed[5]->ti[1] = 1; + + tr[3]->e[2] = ed[1]; + tr[3]->ei[2] = 1; + ed[1]->t[1] = tr[3]; + ed[1]->ti[1] = 2; + + comptriattr(s, tr[3]); /* Compute triangle attributes */ + ADD_ITEM_TO_BOT(s->tris, tr[3]); /* Append to list */ + + /* The four used verticies are now part of the triangulation */ + for (i = 0; i < 4; i++) { + tvs[i]->f |= GVERT_TRI; +//printf("Base triangle %d: %d %d %d (Verticies 0x%x, 0x%x, 0x%x, 0x%x)\n", tr[i]->n, tr[i]->v[0]->n, tr[i]->v[1]->n, tr[i]->v[2]->n, tr[i]->v[0], tr[i]->v[1], tr[i]->v[2]); + } +#ifdef ASSERTS + check_triangulation(s, 0); +#endif + } + + /* Sort the verticies from maximum radius, */ + /* to make our log convex hull logic work */ + sort_verticies(s); + + { + int i; + /* Complete the triangulation by adding all the remaining verticies */ + /* in order of decreasing radius, so that those below the log */ + /* convex hull get discarded. */ + for (i = 0; i < s->nv; i++) { + if (!(s->verts[i]->f & GVERT_SET) + || (s->verts[i]->f & GVERT_TRI) + || (s->verts[i]->f & GVERT_INSIDE)) { + continue; + } + + insert_vertex(s, s->verts[i]); +#ifdef ASSERTS + check_triangulation(s, 0); +#endif + } + } + + /* Number the used verticies */ + renumber_verticies(s); + +#ifdef ASSERTS + check_triangulation(s, 1); +#endif +} + +/* - - - - - - - - - - - - - - - - - - - - - - - - - */ +/* Compute a new convex hull mapping radius for the sample points, */ +/* on the basis of the initial mapping. */ +static void compute_smchrad( +gamut *s +) { + int i, j; + double zz[3] = { 0.0, 0.0, 1.0 }; + double rot[3][3]; + double ssr = 0.5 * s->sres; /* Sample size radius in delta E */ +// double ssr = 10.0; + int res = 4; /* resolution of sample grid */ + +//printf("~1 computing smoothed chrads\n"); + + /* Compute the new log surface value */ + for (i = 0; i < s->nv; i++) { + double pr; /* Current surface radius at this vertex */ + double rad, rw; /* Smoothed radius, weight */ + double out[3]; /* Point on surface for this vertex */ + int x, y; + + if (!(s->verts[i]->f & GVERT_SET)) { + continue; + } + +//printf("~1 vertex %d, %f %f %f\n",i, s->verts[i]->p[0], s->verts[i]->p[1], s->verts[i]->p[2]); + + /* Find the average surface level near this point */ + pr = s->radial(s, out, s->verts[i]->p); + +//printf("~1 surface radius %f, location %f %f %f\n",pr, out[0],out[1],out[2]); + + /* Compute a rotation that lines Z up with the radial direction */ + out[0] -= s->cent[0]; /* Radial vector through point to surface */ + out[1] -= s->cent[1]; + out[2] -= s->cent[2]; + zz[2] = pr; /* Z vector of same length */ + icmRotMat(rot, zz, out); /* Compute vector from Z to radial */ + out[0] += s->cent[0]; + out[1] += s->cent[1]; + out[2] += s->cent[2]; + rad = rw = 0.0; + + /* Sample a rectangular array orthogonal to radial vector, */ + /* and weight samples appropriately */ + for (x = 0; x < res; x++) { + for (y = 0; y < res; y++) { + double tt, off[3], rv, in[3]; + + off[0] = 2.0 * (x/(res-1.0) - 0.5); /* -1.0 to 1.0 */ + off[1] = 2.0 * (y/(res-1.0) - 0.5); /* -1.0 to 1.0 */ + off[2] = 0.0; + + rv = off[0] * off[0] + off[1] * off[1]; + if (rv > 1.0) + continue; /* Outside circle */ +#ifdef NEVER + rv = 1.0 - sqrt(rv); /* Radius from center */ + rv = rv * rv * (3.0 - 2.0 * rv); /* Spline weight it */ +#else + rv = 1.0; /* Moving average weighting */ +#endif + + off[0] *= ssr; /* Scale offset to sampling radius */ + off[1] *= ssr; + + /* Rotate offset to be orthogonal to radial vector */ + icmMulBy3x3(off, rot, off); + + /* Add offset to surface point over current vertex */ + in[0] = out[0] + off[0]; + in[1] = out[1] + off[1]; + in[2] = out[2] + off[2]; + +//printf("~1 grid %d %d, weight %f, offset %f %f %f\n",x,y,rv,off[0],off[1],off[2]); + + /* Sum weighted radius at sample point */ + tt = s->radial(s, NULL, in); + tt = log_scale(s, tt); + rad += rv * tt; + rw += rv; + + } + } + /* Compute sample filtered radius at the sample point */ + rad /= rw; +//printf("~1 sampled radius = %f\n\n",rad); + + /* Now compute new hull mapping radius for this point, */ + /* based on dividing out the sampled radius */ + +// s->verts[i]->lr0 = 40.0 + s->verts[i]->lr0 - rad; + s->verts[i]->lr0 = 40.0 + log_scale(s, s->verts[i]->r[0]) - rad; + /* Prevent silliness */ + if (s->verts[i]->lr0 < (2.0 * FAKE_SEED_SIZE)) + s->verts[i]->lr0 = (2.0 * FAKE_SEED_SIZE); + +//printf("~1 new lr0 = %f\n\n",s->verts[i]->lr0); + + /* recompute ch[] for new lr0 */ + for (j = 0; j < 3; j++) + s->verts[i]->ch[j] = s->verts[i]->sp[j] * s->verts[i]->lr0; + } +} + +/* ===================================================== */ +/* Overall triangulation */ +static void triangulate( +gamut *s +) { + + /* Create the convex hull */ + triangulate_ch(s); + +#if defined(DO_TWOPASS) && !defined(TEST_CONVEX_HULL) + if (s->no2pass == 0) { +#ifdef DEBUG_TRIANG + printf("############ Starting second pass ###################\n"); +#endif + compute_smchrad(s); + del_triang(s); + s->pass++; + triangulate_ch(s); + + /* Three passes is typically slightly better, but slower... */ +// compute_smchrad(s); +// del_triang(s); +// s->pass++; +// triangulate_ch(s); + } +#endif /* DO_TWOPASS && !TEST_CONVEX_HULL */ +} + +/* ===================================================== */ +/* Code that makes use of the triangulation */ +/* ===================================================== */ + +/* return the current surface resolution */ +static double getsres( +gamut *s +) { + return s->sres; +} + +/* return the isJab flag value */ +static int getisjab( +gamut *s +) { + return s->isJab; +} + +/* return the isRast flag value */ +static int getisrast( +gamut *s +) { + return s->isRast; +} + +/* Disable segmented maxima filtering */ +static void setnofilt(gamut *s) { + s->nofilter = 1; +} + +/* return the gamut center value */ +static void getcent(gamut *s, double *cent) { + cent[0] = s->cent[0]; + cent[1] = s->cent[1]; + cent[2] = s->cent[2]; +} + +/* Return the gamut min/max range */ +static void getrange(gamut *s, double *min, double *max) { + + if (min != NULL) { + min[0] = s->mn[0]; + min[1] = s->mn[1]; + min[2] = s->mn[2]; + } + if (max != NULL) { + max[0] = s->mx[0]; + max[1] = s->mx[1]; + max[2] = s->mx[2]; + } +} + +/* return nz if the two gamut are compatible */ +static int compatible( +gamut *s, struct _gamut *t) { + int j; + + /* The same colorspace ? */ + if ((s->isJab && !t->isJab) + || (!s->isJab && t->isJab)) { + return 0; + } + + /* The same gamut center ? */ + for (j = 0; j < 3; j++) { + if (fabs(s->cent[j] - t->cent[j]) > 1e-9) { + return 0; + } + } + return 1; +} + + +/* Return the number of raw verticies used to construct surface */ +static int nrawverts( +gamut *s +) { + int i, nrv = 0; + + /* Sort them so that triangulate doesn't mess indexing up */ + sort_verticies(s); + + /* Count them */ + for (i = 0; i < s->nv; i++) { + if (s->verts[i]->f & GVERT_SET) + nrv++; + } + + return nrv; +} + +/* Return the raw (triangle and non-triangle surface) verticies */ +/* location given its index. */ +/* return the next (sparse) index, or -1 if beyond last */ +static int getrawvert( +gamut *s, +double pos[3], /* Return absolute position */ +int ix /* Input index */ +) { + if (ix < 0) + return -1; + + /* Find then next used in the triangulation */ + for (; ix < s->nv; ix++) { + if (!(s->verts[ix]->f & GVERT_SET)) + continue; + break; + } + + if (ix >= s->nv) + return -1; + + pos[0] = s->verts[ix]->p[0]; + pos[1] = s->verts[ix]->p[1]; + pos[2] = s->verts[ix]->p[2]; + + return ix+1; +} + +/* Return the number of raw direction 0 verticies used */ +/* to construct surface. (Direction 0 is radial direction maxima) */ +static int nraw0verts( +gamut *s +) { + int i, nrv = 0; + + /* Sort them so that triangulate doesn't mess indexing up */ + sort_verticies(s); + + /* Count them */ + for (i = 0; i < s->nv; i++) { + if ((s->verts[i]->f & GVERT_SET) + && (s->verts[i]->k0 > 0)) + nrv++; + } + + return nrv; +} + +/* Return the raw (triangle and non-triangle surface) direction 0 */ +/* verticies location given its index. (Direction 0 is radial direction maxima) */ +/* return the next (sparse) index, or -1 if beyond last */ +static int getraw0vert( +gamut *s, +double pos[3], /* Return absolute position */ +int ix /* Input index */ +) { + if (ix < 0) + return -1; + + /* Find then next used in the triangulation and direction 0 */ + for (; ix < s->nv; ix++) { + if (!(s->verts[ix]->f & GVERT_SET) + || !(s->verts[ix]->k0 > 0)) + continue; + break; + } + + if (ix >= s->nv) + return -1; + + pos[0] = s->verts[ix]->p[0]; + pos[1] = s->verts[ix]->p[1]; + pos[2] = s->verts[ix]->p[2]; + + return ix+1; +} + +/* Return the number of stratified sampling surface verticies, */ +/* for the given verticies per unit area parameter. */ +static int nssverts( +gamut *s, +double xvra /* Extra vertex ratio */ +) { + + if IS_LIST_EMPTY(s->tris) + triangulate(s); + +//printf("~1 nssverts called with xvra = %f\n",xvra); + if (s->xvra != xvra) { + int i, j; + gtri *tp; /* Triangle pointer */ + double tarea; /* Total area */ + double tnverts; /* Target number of verticies */ + int anverts; /* Actual number of verticies */ + + /* Calculate the total surface area of the triangulation */ + tarea = 0.0; + tp = s->tris; + FOR_ALL_ITEMS(gtri, tp) { + double sp, ss[3]; /* Triangle side lengths */ + double dp; /* Dot product of point in triangle and normal */ + + for (i = 0; i < 3; i++) { /* For each edge */ + for (ss[i] = 0.0, j = 0; j < 3; j++) { + double dd = tp->e[i]->v[1]->p[j] - tp->e[i]->v[0]->p[j]; + ss[i] += dd * dd; + } + ss[i] = sqrt(ss[i]); + } + + /* semi-perimeter */ + sp = 0.5 * (ss[0] + ss[1] + ss[2]); + + /* Area of triangle */ + tp->area = sqrt(sp * (sp - ss[0]) * (sp - ss[1]) * (sp - ss[2])); + + tarea += tp->area; + } END_FOR_ALL_ITEMS(tp); + + /* target number of vectors */ + tnverts = xvra * s->ntv; +//printf("~1 total area = %f, tnverts = %f\n",tarea, tnverts); + + /* Number that need to be added using stratified sampling */ + tnverts -= (double)s->ntv; + anverts = 0; + + /* Compute number of extra verticies for each triangle */ + if (tnverts > 0.0) { + double exvpua; /* Extra verticies per unit area to create */ + + exvpua = tnverts/tarea; +//printf("~1 extra verts = %f\n",exvpua); + tp = s->tris; + FOR_ALL_ITEMS(gtri, tp) { + tp->ssverts = (int)(exvpua * tp->area + 0.5); + anverts += tp->ssverts; + } END_FOR_ALL_ITEMS(tp); + } + anverts += s->ntv; + s->xvra = xvra; + s->ssnverts = anverts; + } + +//printf("~1 returning total verts %d\n",s->ssnverts); + return s->ssnverts; +} + +/* Return the stratified sampling surface verticies */ +/* location and radius. nssverts() sets vpua */ +static int getssvert( +gamut *s, +double *rad, /* Return radial radius */ +double pos[3], /* Return absolute position */ +double norm[3], /* Return normal of triangle it orginates from */ +int ix /* Input index */ +) { + int sskip = 0; /* Number of points to skip after each reset of pseudo rand */ + +//printf("getssvert called\n"); + + if (ix < 0) + return -1; + + if (ix < s->nv) { + + /* Find then next used vertex in the triangulation */ + for (; ix < s->nv; ix++) { + if (!(s->verts[ix]->f & GVERT_TRI)) + continue; + break; + } + } + + if (ix < s->nv) { /* returning verticies */ + + if (rad != NULL) + *rad = s->verts[ix]->r[0]; + if (pos != NULL) { + pos[0] = s->verts[ix]->p[0]; + pos[1] = s->verts[ix]->p[1]; + pos[2] = s->verts[ix]->p[2]; + } + if (norm != NULL) { + gvert *vp = s->verts[ix]; + gtri *tp; + int i, j, nt = 0; + for (j = 0; j < 3; j++) + norm[j] = 0.0; + + /* Slow, but search all triangles for this vertex. */ + /* Return the average normal of all the triangles it is part of */ + tp = s->tris; + FOR_ALL_ITEMS(gtri, tp) { + for (i = 0; i < 3; i++) { + if (tp->v[i] == vp) { + for (j = 0; j < 3; j++) + norm[j] += tp->pe[j]; + nt++; + break; + } + } + } END_FOR_ALL_ITEMS(tp); + if (nt == 0) + error("gamut::getssvert() vertex doesn't have a triangle"); + for (j = 0; j < 3; j++) + norm[j] /= (double)nt; + } +//printf("~1 returning tri vertex %f %f %f\n", pos[0],pos[1],pos[2]); + + } else { /* We're generating ss points for each triangle */ + int i, j; + double tt; + double uv[2]; + double tr[3]; /* Baricentric weighting */ + double vv[3]; + + if (s->ss == NULL) { + if ((s->ss = new_sobol(2)) == NULL) + error("gamut::getssvert() new_sobol() failed"); + for (i = 0; i < sskip; i++) + s->ss->next(s->ss, uv); + } + if (ix == s->nv) { /* Start of generating verticies in triangles */ + +//printf("~1 setting up for scan through triangles\n"); + s->nexttri = s->tris; + if (s->nexttri == NULL) + return -1; + s->ssvertn = 0; + s->ss->reset(s->ss); + } + if (s->ssvertn >= s->nexttri->ssverts) { + do { +//printf("~1 skipping to next triangle\n"); + s->nexttri = NEXT_FWD(s->nexttri); + if (s->nexttri == s->tris) + return -1; + } while(s->nexttri->ssverts <= 0); + s->ssvertn = 0; + s->ss->reset(s->ss); + for (i = 0; i < sskip; i++) + s->ss->next(s->ss, uv); + } +//printf("~1 generating ss vert %d out of %d\n",s->ssvertn+1,s->nexttri->ssverts); + s->ss->next(s->ss, uv); + + tt = sqrt(uv[0]); + tr[0] = 1 - tt; + tr[1] = uv[1] * tt; + tr[2] = 1.0 - tr[0] - tr[1]; + + vv[0] = vv[1] = vv[2] = 0.0; + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) + vv[j] += s->nexttri->v[i]->p[j] * tr[i]; + } + + if (rad != NULL) + *rad = icmNorm33(vv, s->cent); + if (pos != NULL) { + pos[0] = vv[0]; + pos[1] = vv[1]; + pos[2] = vv[2]; + } + if (norm != NULL) { + norm[0] = s->nexttri->pe[0]; + norm[1] = s->nexttri->pe[1]; + norm[2] = s->nexttri->pe[2]; + } + s->ssvertn++; +//printf("~1 returning ss vertex %f %f %f\n", pos[0],pos[1],pos[2]); + } + + return ix+1; +} + +/* Return the number of verticies in the triangulated surface */ +static int nverts( +gamut *s +) { + if IS_LIST_EMPTY(s->tris) + triangulate(s); + + return s->ntv; +} + +/* Return the verticies location and radius given its index. */ +/* return the next (sparse) index, or -1 if beyond last */ +static int getvert( +gamut *s, +double *rad, /* Return radial radius */ +double pos[3], /* Return absolute position */ +int ix /* Input index */ +) { + if (ix >= s->nv) + return -1; + + /* Find then next used in the triangulation */ + for (; ix < s->nv; ix++) { + if (!(s->verts[ix]->f & GVERT_TRI)) + continue; + break; + } + if (ix >= s->nv) + return -1; + + if (rad != NULL) + *rad = s->verts[ix]->r[0]; + if (pos != NULL) { + pos[0] = s->verts[ix]->p[0]; + pos[1] = s->verts[ix]->p[1]; + pos[2] = s->verts[ix]->p[2]; + } + + return ix+1; +} + + +/* Reset indexing through triangles for getnexttri() */ +static void startnexttri(gamut *s) { + if IS_LIST_EMPTY(s->tris) + triangulate(s); + + s->nexttri = NULL; +} + +/* Return the next surface triange, nz on no more */ +static int getnexttri( +gamut *s, +int v[3] /* Return indexes for same order as getvert() */ +) { + if IS_LIST_EMPTY(s->tris) + triangulate(s); + + if (s->nexttri == NULL) { + s->nexttri = s->tris; + if (s->nexttri == NULL) + return 1; + } else { + s->nexttri = NEXT_FWD(s->nexttri); + if (s->nexttri == s->tris) + return 1; + } + + v[0] = s->nexttri->v[0]->tn; + v[1] = s->nexttri->v[1]->tn; + v[2] = s->nexttri->v[2]->tn; + return 0; +} + +/* ===================================================== */ + +/* Return the total volume of the gamut */ +/* [ We use the formula from "Area of planar polygons and */ +/* volume of polyhedra" by Ronald N. Goldman, */ +/* Graphics Gems II, pp 170 ] */ +static double volume( +gamut *s +) { + int i, j; + gtri *tp; /* Triangle pointer */ + double vol; /* Gamut volume */ + + if IS_LIST_EMPTY(s->tris) + triangulate(s); + + /* Compute the area of each triangle in the list, */ + /* and accumulate the gamut volume. */ + tp = s->tris; + vol = 0.0; + FOR_ALL_ITEMS(gtri, tp) { + double sp, ss[3]; /* Triangle side lengths */ + double area; /* Area of this triangle */ + double dp; /* Dot product of point in triangle and normal */ + + for (i = 0; i < 3; i++) { /* For each edge */ + for (ss[i] = 0.0, j = 0; j < 3; j++) { + double dd = tp->e[i]->v[1]->p[j] - tp->e[i]->v[0]->p[j]; + ss[i] += dd * dd; + } + ss[i] = sqrt(ss[i]); + } + + /* semi-perimeter */ + sp = 0.5 * (ss[0] + ss[1] + ss[2]); + + /* Area of triangle */ + area = sqrt(sp * (sp - ss[0]) * (sp - ss[1]) * (sp - ss[2])); + + /* Dot product between first vertex in triangle and the unit normal vector */ + dp = tp->v[0]->p[0] * tp->pe[0] + + tp->v[0]->p[1] * tp->pe[1] + + tp->v[0]->p[2] * tp->pe[2]; + + /* Accumulate gamut volume */ + vol += dp * area; + + } END_FOR_ALL_ITEMS(tp); + + vol = fabs(vol)/3.0; + + return vol; +} + +/* ===================================================== */ +/* ===================================================== */ +/* Given a point, */ +/* return the distance to the gamut surface. */ + +static void init_lu(gamut *s); +static gtri *radial_point_triang(gamut *s, gbsp *np, double in[3]); +static double radial_point(gamut *s, gbsp *np, double in[3]); + +/* Given a point, return the point in that direction */ +/* that lies on the gamut surface. Return the radial */ +/* radius to the surface point */ +/* Brute force search version. */ +static double +radial_bf( +gamut *s, +double *out, /* result point (absolute)*/ +double *in /* input point (absolute)*/ +) { + gtri *tp; + int j; + double ss, rv; + double nin[3]; /* Normalised input vector */ + +//printf("~1 radial called with %f %f %f\n", in[0], in[1], in[2]); + if IS_LIST_EMPTY(s->tris) + triangulate(s); + + /* Compute vector length to center point */ + for (ss = 0.0, j = 0; j < 3; j++) + ss += (in[j] - s->cent[j]) * (in[j] - s->cent[j]); + ss = 1.0/sqrt(ss); /* Normalising factor */ + for (ss = 0.0, j = 0; j < 3; j++) + nin[j] = s->cent[j] + (in[j] - s->cent[j]) * ss; + + tp = s->tris; + FOR_ALL_ITEMS(gtri, tp) { + if (vect_intersect(s, &rv, out, s->cent, nin, tp)) { + if (rv > 0.0) /* Expect only one intersection */ + break; + } + } END_FOR_ALL_ITEMS(tp); + +//printf("~1 result = %f %f %f\n",out[0], out[1], out[2]); + + return rv; +} + +/* Implementation for following two functions: */ +/* Given a point, return the point in that direction */ +/* that lies on the gamut surface. Use the BSP accellerated search. */ +/* Return the radial length of the input and radial length of result */ +static void +_radial( +gamut *s, +double *ir, /* return input radius (may be NULL) */ +double *or, /* return output radius (may be NULL) */ +double *out, /* result point (absolute) (may be NULL) */ +double *in /* input point (absolute)*/ +) { + int j; + double ss, rv; + double nin[3]; /* Normalised input vector */ + + if IS_LIST_EMPTY(s->tris) + triangulate(s); + + /* We have to find out which triangle the point is in */ + if (s->lu_inited == 0) { + init_lu(s); /* Init BSP search tree */ + } +//if (trace) printf("~1 radial called with %f %f %f\n", in[0], in[1], in[2]); + + for (j = 0; j < 3; j++) + nin[j] = in[j] - s->cent[j]; /* relative to gamut center */ + + for (ss = 0.0, j = 0; j < 3; j++) + ss += nin[j] * nin[j]; + ss = sqrt(ss); + if (ss > 1e-9) { /* Normalise to 1.0 */ + for (j = 0; j < 3; j++) + nin[j] /= ss; + } else { + nin[0] = 1.0; + nin[1] = nin[2] = 0.0; + } + +//if (trace) printf("~1 Normalised in = %f %f %f\n", nin[0], nin[1], nin[2]); + rv = radial_point(s, s->lutree, nin); + + if (rv < 0.0) { + error("gamut: radial internal error - failed to find triangle\n"); + } + + if (out != NULL) { + for (j = 0; j < 3; j++) + out[j] = nin[j] * rv + s->cent[j]; /* Scale out to surface length, absolute */ +//if (trace) printf("~1 result = %f %f %f\n",out[0], out[1], out[2]); + } + + if (ir != NULL) { +//if (trace) printf("~1 input radius res = %f\n",ss); + *ir = ss; + } + + if (or != NULL) { +//if (trace) printf("~1 output radius res = %f\n",rv); + *or = rv; + } +} + +/* Given a point, return the point in that direction */ +/* that lies on the gamut surface */ +/* Return the normalised radial radius to the surface point */ +static double +nradial( +gamut *s, +double *out, /* result point (absolute) (May be NULL) */ +double *in /* input point (absolute)*/ +) { + double ss, rv; + + _radial(s, &ss, &rv, out, in); + return ss/rv; +} + +/* Given a point, return the point in that direction */ +/* that lies on the gamut surface */ +/* Return the radial radius to the surface point */ +static double +radial( +gamut *s, +double *out, /* result point (absolute) (May be NULL) */ +double *in /* input point (absolute)*/ +) { + double ss, rv; + + _radial(s, &ss, &rv, out, in); + return rv; +} + +void lu_split(gamut *s, gbsp **np, int rdepth, gtri **list, int llen); + +/* Setup the radial lookup function acceleration structure */ +static void +init_lu( +gamut *s +) { + static double v0[3] = {0.0, 0.0, 0.0}; + static + gedge *ep; /* Edge pointer */ + gtri *tp; /* Triangle pointer */ + gtri **tlist; + int ntris; + +//printf("~1 init_lu called\n"); + + /* Create mean angle dividing plane equations */ + ep = s->edges; + FOR_ALL_ITEMS(gedge, ep) { + plane_equation(ep->re, v0, ep->v[0]->sp, ep->v[1]->sp); + } END_FOR_ALL_ITEMS(ep); + + /* Create the initial triangle list */ + /* First count them */ + ntris = 0; + tp = s->tris; + FOR_ALL_ITEMS(gtri, tp) { + ntris++; + } END_FOR_ALL_ITEMS(tp); + + /* Allocate a list */ + if ((tlist = (gtri **) malloc(ntris * sizeof(gtri *))) == NULL) { + fprintf(stderr,"gamut: malloc failed - top level triangle list (%d entries)\n",ntris); + exit(-1); + } + + /* Then add them to the list */ + ntris = 0; + tp = s->tris; + FOR_ALL_ITEMS(gtri, tp) { + tlist[ntris] = tp; + ntris++; + } END_FOR_ALL_ITEMS(tp); + + /* Recursively split them, and add objects to leaves */ + lu_split(s, &s->lutree, 0, tlist, ntris); + + free(tlist); + +//printf("~1 init_lu done\n"); + s->lu_inited = 1; +} + +/* + * BSP accellerator: + * This is setup specifically to accellerate finding + * the radial point on the gamut surface. To do this, all + * the BSP plains pass through the gamut center, creating + * wedge shaped sub-division regions. + * + * For accellerating the vector intersect code, this isn't + * so fabulous, and a general unconstrained BSP tree would + * be better. To address this, an orhogonal element to the + * radial BSP's is provided in the radius squared range + * of each set of elements below a BSP node. + */ + +/* Recursive routine to choose a partition plane, */ +/* and then split the triangle list between the */ +/* +ve and -ve sides, or add triangles as leaves. */ +void +lu_split( +gamut *s, +gbsp **np, /* Address of node pointer to be set */ +int rdepth, /* Current recursion depth */ +gtri **list, /* Current triangle list */ +int llen /* Number of triangles in the list */ +) { + double rs0, rs1; /* Radius squared range of elements */ + int ii, jj; /* Progress through edges */ + int pcount; /* Current best scored try */ + int ncount; + int bcount; + int mcount; + double peqs[4] = { 0.0, 0.0, 0.0, 0.0 }; + gtri **plist, **nlist; /* New sub-lists */ + int pix, nix; /* pos/ned sublist indexes */ + gbspn *bspn; /* BSP decision node */ + +//printf("~1\nlu_split called at depth %d with %d triangles\n",rdepth, llen); +#ifdef DEBUG + if (llen <= 3) { + int i; + for (i = 0; i < llen; i++) { + printf("Triang index %d = %d\n",i, list[i]->n); + printf("Triang verts %d %d %d\n", + list[i]->v[0]->tn, list[i]->v[1]->tn, list[i]->v[2]->tn); + printf("Vert 0 at %.18f %.18f %.18f\n",list[i]->v[0]->sp[0], list[i]->v[0]->sp[1], list[i]->v[0]->sp[2]); + printf("Vert 1 at %.18f %.18f %.18f\n",list[i]->v[1]->sp[0], list[i]->v[1]->sp[1], list[i]->v[1]->sp[2]); + printf("Vert 2 at %.18f %.18f %.18f\n",list[i]->v[2]->sp[0], list[i]->v[2]->sp[1], list[i]->v[2]->sp[2]); + } + } +#endif /* DEBUG */ + + if ((rdepth+1) >= BSPDEPTH) { /* Oops */ + printf("gamut internal error: ran out of recursion depth in BSP\n"); + exit (-1); + } + + /* Scan our list or triangles and figure out radius squared range */ + { + int i, j, e; + double rs; + + rs0 = 1e120; + rs1 = -1.0; + for (i = 0; i < llen; i++) { + if (list[i]->rs0 < rs0) + rs0 = list[i]->rs0; + if (list[i]->rs1 > rs1) + rs1 = list[i]->rs1; + } +//printf("~1 no triangs %d, rs range %f - %f\n",llen,rs0,rs1); + } + + pcount = ncount = bcount = -1; + mcount = 0; + /* test every edge in turn */ + for (ii = jj = 0;ii < llen;) { + double eqs[4]; + int i; + gedge *ep; /* Edge pointer */ + int pc, nc, bc; /* Score a try, postive count, negative count, both count */ + int mc; /* Minumum count */ + + ep = list[ii]->e[jj]; + eqs[0] = ep->re[0]; /* Use this edge */ + eqs[1] = ep->re[1]; + eqs[2] = ep->re[2]; + eqs[3] = ep->re[3]; + if (++jj > 2) { + jj = 0; + ii++; + } + + /* Do the trial split */ + pc = nc = bc = 0; + for (i = 0; i < llen; i++) { + int j; + int po, ne; + + /* Compute distance from plane of all verticies in triangle */ + po = ne = 0; + for (j = 0; j < 3; j++) { /* For triangle verticies */ + double ds; + /* Compute distance to dividing plane of this vertex */ + ds = eqs[0] * list[i]->v[j]->sp[0] + + eqs[1] * list[i]->v[j]->sp[1] + + eqs[2] * list[i]->v[j]->sp[2] + + eqs[3]; + /* Figure if the verticies are clearly to one side of the plane */ + if (ds > 1e-10) { + po++; + } else if (ds < -1e-10) { + ne++; + } + } + /* Score this split */ + if (po) { + pc++; + if (ne) { + nc++; + bc++; + list[i]->sort = 3; /* Both */ + } else { + list[i]->sort = 1; /* +ve */ + } + } else if (ne) { + nc++; + list[i]->sort = 2; /* -ve */ + } else { /* Hmm. Neither */ + bc++; + list[i]->sort = 3; /* Assume both */ + } + } + mc = pc < nc ? pc : nc; /* Size of smallest group */ + mc -= bc; +//printf("~1 lu_split trial %d, mc %d, pc %d, nc %d, bc %d\n",ii * 3 + jj, mc, pc, nc, bc); + if (mc > mcount) { /* New largest small group */ + mcount = mc; + pcount = pc; + ncount = nc; + bcount = bc; + peqs[0] = eqs[0]; + peqs[1] = eqs[1]; + peqs[2] = eqs[2]; + peqs[3] = eqs[3]; +//printf("~1 new best - plane mc = %d, %f %f %f %f\n",mc, peqs[0], peqs[1], peqs[2], peqs[3]); + for (i = 0; i < llen; i++) { + list[i]->bsort = list[i]->sort; + } + } + } + +#ifdef DEBUG_SPLIT_VRML + write_split_diag_vrml(s, list, llen); + getchar(); +#endif /* DEBUG_SPLIT_VRML */ + + if (ii >= llen && bcount < 0) { /* We failed to find a split plane. */ + /* This is usually a result of the list being 2 or more triangles */ + /* that do not share any edges (disconected from each other), and */ + /* lying so that any split plane formed from an edge of one, */ + /* intersects one of the others. */ + /* In theory we could solve this by picking some */ + /* other radial split plane ? */ + + /* Instead leave our list of triangles as the leaf node, */ + /* and let the search algorithms deal with this. */ + + *np = (gbsp *)new_gbspl(llen, list); + (*np)->rs0 = rs0; /* Radius squared range */ + (*np)->rs1 = rs1; +//printf("~1 lu_split returning with a non split list of %d triangles\n",llen); + return; + } + + /* Divide the triangles into two lists */ + bspn = new_gbspn(); /* Next node */ + *np = (gbsp *)bspn; /* Put it in place */ + bspn->rs0 = rs0; /* Radius squared range */ + bspn->rs1 = rs1; + bspn->pe[0] = peqs[0]; /* Plane equation */ + bspn->pe[1] = peqs[1]; + bspn->pe[2] = peqs[2]; + bspn->pe[3] = peqs[3]; + + /* Allocate the sub lists */ + if ((plist = (gtri **) malloc(pcount * sizeof(gtri *))) == NULL) { + fprintf(stderr,"gamut: malloc failed - pos sub-list\n"); + exit(-1); + } + if ((nlist = (gtri **) malloc(ncount * sizeof(gtri *))) == NULL) { + fprintf(stderr,"gamut: malloc failed - neg sub-list\n"); + exit(-1); + } + + /* Fill them in */ + for (pix = nix = ii = 0; ii < llen; ii++) { + if (list[ii]->bsort & 1) { /* Positive */ + plist[pix] = list[ii]; + pix++; + } + if (list[ii]->bsort & 2) { /* Negative */ + nlist[nix] = list[ii]; + nix++; + } + } + + /* Recurse if there are more triangles to split */ + if (pix == 1) { + bspn->po = (gbsp *)plist[0]; /* leaf node */ +//printf("~1 pos leaf with triangle %d\n",plist[0]->n); + } else if (pix > 1) { +//printf("~1 About to recurse on positive with list of %d\n",pix); + lu_split(s, &bspn->po, rdepth+1, plist, pix); + } + + if (nix == 1) { +//printf("~1 neg leaf with triangle %d\n",nlist[0]->n); + bspn->ne = (gbsp *)nlist[0]; /* leaf node */ + } else if (nix > 1) { +//printf("~1 About to recurse on negative with list of %d\n",nix); + lu_split(s, &bspn->ne, rdepth+1, nlist, nix); + } + + free(plist); + free(nlist); +//printf("~1 lu_split returning\n"); +} + +/* Given a point and a node in the BSP tree, recurse down */ +/* the tree, or return the triangle it lies in. */ +/* Return NULL if it wasn't in any triangle (shouldn't happen with a closed gamut ?). */ +static gtri *radial_point_triang( +gamut *s, +gbsp *np, /* BSP node pointer we're at */ +double *nin /* Normalised center relative point */ +) { + gtri *rv; +//if (trace) printf("~1 rad_pnt_tri: BSP 0x%x tag = %d, point %f %f %f\n", np,np->tag,nin[0],nin[1],nin[2]); + if (np->tag == 1) { /* It's a BSP node */ + gbspn *n = (gbspn *)np; + double ds; + + ds = n->pe[0] * nin[0] + + n->pe[1] * nin[1] + + n->pe[2] * nin[2] + + n->pe[3]; + +//if (trace) printf("~1 checking against BSP plane, ds = %e\n",ds); + /* Recurse down both sides it might be in */ + if (ds > -1e-12) { + if ((rv = radial_point_triang(s, n->po, nin)) != NULL) + return rv; + } + if (ds < 1e-12) { + if ((rv = radial_point_triang(s, n->ne, nin)) != NULL) + return rv; + } + return NULL; /* Hmm */ + + } else { /* It's a triangle or list of triangles */ + int nt; /* Number of triangles in list */ + gtri **tpp; /* Pointer to list of triangles */ + int i, j; + + if (np->tag == 2) { /* It's a triangle */ + tpp = (gtri **)&np; + nt = 1; + } else if (np->tag == 3) { /* It's a triangle list */ + gbspl *n = (gbspl *)np; + tpp = n->t; + nt = n->nt; + } + + /* Go through the list and stop at the first triangle */ + /* that the node lies in. */ + for (i = 0; i < nt; i++, tpp++) { + gtri *t = *tpp; + + /* Check if the point is within this triangle */ + for (j = 0; j < 3; j++) { + double ds; + ds = t->ee[j][0] * nin[0] + + t->ee[j][1] * nin[1] + + t->ee[j][2] * nin[2] + + t->ee[j][3]; + if (ds > 1e-10) + break; /* Not within triangle */ + } + if (j >= 3) { +//if (trace) printf("~1 located triangle from list that we're in %d\n",n->t[i]->n); + return t; + } + } + /* Hmm. */ + } + +//if (trace) printf("~1 failed to find a triangle\n"); + return NULL; +} + +/* Return the location on the surface of the triangle */ +/* that is intersected by the radial direction */ +/* of the given relative point. Return the distance to */ +/* the gamut surface. */ +static double radial_point( +gamut *s, +gbsp *np, /* BSP node pointer we're at */ +double *nin /* Normalised center relative point */ +) { + gtri *t; + double rv; + +//if (trace) printf("~1 radial_point: BSP 0x%x tag = %d, point %f %f %f\n", np,np->tag,nin[0],nin[1],nin[2]); + + t = radial_point_triang(s, np, nin); + + /* If we failed to find a triangle, or the result was incorrect, do a */ + /* brute force search to be sure of the result. */ + if (t == NULL) { + error("rspl.radial: failed to find radial triangle\n"); + } + + /* Compute the intersection of the input vector with the triangle plane */ + /* (Since nin[] is already relative, we don't need to subtract cent[] from it) */ + rv = -(t->pe[0] * s->cent[0] + t->pe[1] * s->cent[1] + t->pe[2] * s->cent[2] + t->pe[3])/ + (t->pe[0] * nin[0] + t->pe[1] * nin[1] + t->pe[2] * nin[2]); + +#ifdef ASSERTS + /* check the result */ + { + double tt[3]; + double ds; + int j; + for (j = 0; j < 3; j++) /* Compute result, absolute */ + tt[j] = nin[j] * rv + s->cent[j]; + + ds = t->pe[0] * tt[0] + + t->pe[1] * tt[1] + + t->pe[2] * tt[2] + + t->pe[3]; + + if (fabs(ds) > 1e-6) { + fprintf(stderr,"radial: distance to plane not zero! %e\n",ds); + exit(-1); + } + + /* Check if the closest point is within this triangle */ + for (j = 0; j < 3; j++) { + double ds; + ds = t->ee[j][0] * (tt[0] - s->cent[0]) + + t->ee[j][1] * (tt[1] - s->cent[1]) + + t->ee[j][2] * (tt[2] - s->cent[2]) + + t->ee[j][3]; + if (ds > 1e-8) { + fprintf(stderr,"radial: lookup point wasn't within its triangle (%f) !!\n",ds); + exit(-1); + } + } + } +#endif /* ASSERTS */ + +//if (trace) printf("~1 radial_point: rv = %f\n",rv); + return rv; +} + +/* Recursively free a gbsp node and all its children */ +static void del_gbsp(gbsp *n) { + int tag = n->tag; + + if (tag == 1) { /* Another decision node */ + gbspn *dn = (gbspn *)n; + del_gbsp(dn->po); /* Delete children */ + del_gbsp(dn->ne); + del_gbspn(dn); /* And itself */ + + } else if (tag == 3) { /* If a triangle list */ + gbspl *dl = (gbspl *)n; + del_gbspl(dl); /* Delete itself */ + } + + /* Don't delete triangles (tag == 2) since they */ + /* have their own linked list, and may have already been deleted. */ + /* Note we need to be called _before_ triangles are deleted though, */ + /* since we access them to get the tag. */ +} + +/* =================================== */ +/* Given a point, */ +/* return the nearest point on the gamut surface. */ + +#define GNN_INF 1e307 +static void init_ne(gamut *s); + +/* Given an absolute point, return the point on the gamut */ +/* surface that is closest to it. */ +/* Use a brute force search */ +static void +nearest_bf( +gamut *s, +double *out, /* result point (absolute) */ +double *q /* Target point (absolute) */ +) { + gtri *tp; + double bdist = 1e308; /* Best possible distance to an object outside the window */ + + + if IS_LIST_EMPTY(s->tris) + triangulate(s); + + tp = s->tris; + FOR_ALL_ITEMS(gtri, tp) { + double r[3]; /* Possible solution point */ + double tdist; + + /* Compute distance from query point to this object */ + tdist = ne_point_on_tri(s, tp, r, q); + + if (tdist < bdist) { /* New best point */ + bdist = tdist; + out[0] = r[0]; + out[1] = r[1]; + out[2] = r[2]; + } + } END_FOR_ALL_ITEMS(tp); +} + +/* Using nearest neighbourhood accelleration structure: */ + +/* Given an absolute point, return the point on the gamut */ +/* surface that is closest to it. */ +static void +nearest_tri( +gamut *s, +double *rout, /* result point (absolute) */ +double *q, /* Target point (absolute) */ +gtri **ctri /* If not NULL, return pointer to nearest triangle */ +) { + gnn *p; /* Pointer to nearest neighbor structure */ + int e, i; + double r[3] = {0.0, 0.0, 0.0 }; /* Possible solution point */ + double out[3] = {0.0, 0.0, 0.0}; /* Current best output value */ + int wex[3 * 2]; /* Current window edge indexes */ + double wed[3 * 2]; /* Current window edge distances */ + /* Indexes are axis * 2 +0 for lower edge, */ + /* +1 for upper edge of search box. */ + /* We are comparing lower edge of search box */ + /* with upper edge of bounding box etc. */ + +//printf("~1 nearest called\n"); + if IS_LIST_EMPTY(s->tris) + triangulate(s); + + /* We have to find out which triangle the point will be nearest */ + if (s->ne_inited == 0) { + init_ne(s); /* Init nn structure */ + } + p = s->nns; + + if ((p->tbase + 3) < p->tbase) { /* Overflow of touch count */ + for (i = 0; i < p->n; i++) + p->sax[0][i]->touch = 0; /* reset it in all the objects */ + p->tbase = 0; + } + p->ttarget = p->tbase + 3; /* Target touch value */ + +//printf("\n"); +//printf("Query point is %f %f %f\n",q[0], q[1], q[2]); + + /* Find starting indexes within axis arrays */ + for (e = 0; e < (2 * 3); e++) { /* For all axes min & max */ + int f = e/2; /* Axis */ + int mm = (e ^ 1) & 1; /* Min/Max index used for edges */ + int i0, i1, i2; + double v0, v1, v2; + double qf, ww; + + /* Binary search this edge */ + qf = q[f]; /* strength reduced q[f] */ + +//printf("\n"); +//printf("isearching axis %d %s for %f\n",f, e & 1 ? "max" : "min", qf); + i0 = 0; + i2 = p->n - 1; + v0 = p->sax[e][i0]->mix[mm][f]; + v2 = p->sax[e][i2]->mix[mm][f]; +//printf("start points %d - %d, bound %f - %f\n",i0, i2, v0, v2); + + if (qf <= v0) { + i2 = i0; + v2 = v0; + } else if (qf >= v2) { + i0 = i2; + v0 = v2; + } else { + do { + i1 = (i2 + i0)/2; /* Trial point */ + v1 = p->sax[e][i1]->mix[mm][f]; /* Value at trial */ + if (v1 < qf) { + i0 = i1; /* Take top half */ + v0 = v1; + } else { + i2 = i1; /* Take bottom half */ + v2 = v1; + } +//printf("current point %d - %d, bound %f - %f\n",i0, i2, v0, v2); + } while ((i2 - i0) > 1); + } + + if (e & 1) { /* Max side of window */ + int tc; /* total object count */ + + ww = v2 - qf; + wed[e] = fabs(ww) * ww; + wex[e] = i2; + + /* Check that min and max together will cover at least p->n objects */ + tc = p->n - i2 + wex[e ^ 1] + 1; +//printf("got %d, expected %d\n",tc, p->n); + + /* (I don't really understand why this works!) */ + if (tc < p->n) { /* We haven't accounted for all the objects */ + int el = e ^ 1; /* Low side sax */ + int ti0, ti2; + double tv0, tv2; + + ti0 = wex[el]; + ti2 = i2; +//printf("We have straddling objects, initial indexes are %d - %d\n",ti0, ti2); + + /* While straddling objects remain undiscovered: */ + while (tc < p->n) { + tv0 = GNN_INF; /* Guard values */ + tv2 = -GNN_INF; + + /* Increment low side until we find a straddler */ + while (ti0 < (p->n-1)) { + ww = p->sax[el][++ti0]->mix[0][f]; /* Position of the other end */ + if (ww < qf) { +//printf("found low object %d at index %d that straddles\n",p->sax[el][ti0]->n,ti0); + tv0 = qf - p->sax[el][ti0]->mix[1][f]; + break; + } + } + + /* Decrement high side until we find a straddler */ + while (ti2 > 0) { + ww = p->sax[e][--ti2]->mix[1][f]; /* Position of the other end */ + if (ww > qf) { +//printf("found high object %d at index %d that straddles\n",p->sax[e][ti2]->n,ti2); + tv2 = p->sax[e][ti2]->mix[0][f] - qf; + break; + } + } + /* Choose the closest */ + if (tv0 > tv2) { + wed[el] = fabs(tv0) * tv0; + wex[el] = ti0; + tc++; + } else { + wed[e] = fabs(tv2) * tv2; + wex[e] = ti2; + tc++; + } + } +//printf("After correction we have %d - %d\n",wex[e^1], wex[e]); + } + } else { /* Min side of window */ + ww = q[f] - v0; + wed[e] = fabs(ww) * ww; + wex[e] = i0; + } + } + + /* Expand a 3 dimenstional cube centered on the target point, */ + /* jumping to the next nearest point on any axis, discovering */ + /* any bounding boxes that are within the expanding window */ + /* by checking their touch count. */ + + /* The first point found establishes the initial best distance. */ + /* When the window expands beyond the point where it can improve */ + /* the best distance, stop */ + + { + double bw = 0.0; /* Current window distance */ + double bdist = 1e308; /* Best possible distance to an object outside the window */ + gtri *bobj = NULL; + int ptested = 0; /* Stats */ + int pcalced = 0; /* Stats */ + + /* Until we're done */ + for (;;ptested++) { + int ee; /* Axis & expanding box edge */ + int ff; /* Axis */ + int ii; /* Index of chosen point */ + gtri *ob; /* Current object */ + unsigned int ctv; /* Current touch value */ +//printf("\n"); +//printf("wwidth = %f, bdist = %f, window = %d-%d, %d-%d, %d-%d\n", +//bw, bobj == NULL ? 0.0 : bdist, wex[0], wex[1], wex[2], wex[3], wex[4], wex[5]); +//printf("window edge distances are = %f-%f, %f-%f, %f-%f\n", +//wed[0], wed[1], wed[2], wed[3], wed[4], wed[5]); + + /* find next (smallest) window increment axis and direction */ + ee = 0; + ii = wex[ee]; + bw = wed[ee]; + for (e = 1; e < (2 * 3); e++) { + if (wed[e] < bw) { + ee = e; + ii = wex[e]; + bw = wed[e]; + } + } +//printf("Next best is axisdir %d, object %d, axis index %d, best possible dist %f\n", +//ee, p->sax[ee][ii]->n, ii, bw); + + if (bw == GNN_INF || bw > bdist) { + break; /* Can't got any further, or further points will be worse */ + } + +#ifdef ASSERTS + if (ii < 0 || ii >= p->n) { + printf("Assert: went out of bounds of sorted axis array\n"); + exit(0); + } +#endif + /* Chosen point on ee axis/direction, index ii */ + ff = ee / 2; /* Axis only */ + + ob = p->sax[ee][ii]; + + /* Touch value of current object */ + ctv = ob->touch; + + if (ctv < p->ttarget) { /* Not been dealt with before */ + + /* Touch this new window boundary point */ + ob->touch = ctv = ((ctv < p->tbase) ? p->tbase : ctv) + 1; + +//printf("New touch count on %d is %d, target %d\n", ob->n, p->sax[ee][ii]->touch, p->ttarget); + + /* Check the point out */ + if (ctv == (p->tbase + 3)) { /* Is within window on all axes */ + double tdist; + + pcalced++; /* Stats */ + + /* Compute distance from query point to this object */ + tdist = ne_point_on_tri(s, ob, r, q); + +//printf("Got new best point %d, dist %f\n",i,tdist); + if (tdist < bdist) { /* New best point */ + bobj = ob; + bdist = tdist; + out[0] = r[0]; + out[1] = r[1]; + out[2] = r[2]; + } + } + } + + /* Increment next window edge candidate, and figure new edge distance */ + if (ee & 1) { /* Top */ + if (++wex[ee] >= p->n) { + wed[ee] = GNN_INF; + wex[ee]--; + } else { + double ww = p->sax[ee][wex[ee]]->mix[0][ff] - q[ff]; + wed[ee] = fabs(ww) * ww; + } + } else { + if (--wex[ee] < 0) { + wed[ee] = GNN_INF; + wex[ee]++; + } else { + double ww = q[ff] - p->sax[ee][wex[ee]]->mix[1][ff]; + wed[ee] = fabs(ww) * ww; + } + } + } + +//printf("Searched %d points out of %d = %f%%\n",ptested, p->n, 100.0 * ptested/p->n); + + p->tbase += 3; /* Next touch */ + + if (rout != NULL) { + rout[0] = out[0]; /* Copy results to output */ + rout[1] = out[1]; + rout[2] = out[2]; + } + + if (ctri != NULL) + *ctri = bobj; + + return; + } +} + +/* Given an absolute point, return the point on the gamut */ +/* surface that is closest to it. */ +static void +nearest( +gamut *s, +double *rout, /* result point (absolute) */ +double *q /* Target point (absolute) */ +) { + nearest_tri(s, rout, q, NULL); +} + +/* Perturb the containment points to avoid */ +/* numerical co-incidence */ +double perturb[21] = { + 8.9919295344233395e-283, 1.1639766020018968e+224, 1.2554893023590904e+232, + 2.3898157055642966e+190, 1.5697612415774029e-076, 6.6912978722191457e+281, + 1.2369092402930559e+277, 1.4430907501246712e-153, 3.0017439193018232e+238, + 1.2978311824382444e+161, 5.5068703318775818e-311, 7.7791723264448801e-260, + 4.4296571592384350e+281, 8.9481529920968425e+165, 1.2845894914769635e-153, + 2.0835868791190880e-076, 5.4310198502711138e+241, 4.8689849775675438e+275, + 9.2709981544886391e+122, 3.7958270103353899e-153, 7.1366083837501666e-154 +}; + +/* Setup the nearest function acceleration structure */ +static void +init_ne( +gamut *s +) { + gnn *p; + int i, k; + gtri *tp; /* Triangle pointer */ + int ntris; + double psf; + +//printf("~1 init_ne called\n"); + + /* Allocate the nearest neighbor acceleration structure */ + if ((s->nns = p = (gnn *) calloc(1, sizeof(gnn))) == NULL) { + fprintf(stderr,"gamut: calloc failed - gnn structure\n"); + exit(-1); + } + + /* Count triangles */ + ntris = 0; + tp = s->tris; + FOR_ALL_ITEMS(gtri, tp) { + ntris++; + } END_FOR_ALL_ITEMS(tp); + + p->n = ntris; + p->tbase = 0; /* Initialse touch flag */ + + /* Allocate the arrays spaces */ + for (k = 0; k < (3 * 2); k++) { + if ((p->sax[k] = (gtri **)malloc(sizeof(gtri *) * ntris)) == NULL) + error("Failed to allocate sorted index array"); + } + + /* Compute pertbation factor */ + for (psf = 0.0, i = 1; i < 21; i++) + psf += perturb[i]; + psf *= perturb[0]; + + /* For each triangle, create the triangle bounding box values, */ + /* and add them to the axis lists. */ + tp = s->tris; + i = 0; + FOR_ALL_ITEMS(gtri, tp) { + int j; + for (j = 0; j < 3; j++) { /* Init */ + tp->mix[0][j] = 1e38; + tp->mix[1][j] = -1e38; + } + for (k = 0; k < 3; k++) { + for (j = 0; j < 3; j++) { + if (tp->v[k]->p[j] < tp->mix[0][j]) /* New min */ + tp->mix[0][j] = psf * tp->v[k]->p[j]; + if (tp->v[k]->p[j] > tp->mix[1][j]) /* New max */ + tp->mix[1][j] = psf * tp->v[k]->p[j]; + } + p->sax[k * 2 + 0][i] = tp; + p->sax[k * 2 + 1][i] = tp; + } +//printf("~1 tri %d has bb %f - %f, %f - %f, %f - %f\n", i, tp->mix[0][0], tp->mix[1][0], tp->mix[0][1], tp->mix[1][1], tp->mix[0][2], tp->mix[1][2]); + i++; + } END_FOR_ALL_ITEMS(tp); + + + /* Sort the axis arrays */ + for (k = 0; k < 3; k++) { + + /* Sort upper edges of bounding box */ +#define HEAP_COMPARE(A,B) (A->mix[1][k] < B->mix[1][k]) + HEAPSORT(gtri *, &p->sax[k * 2 + 0][0], ntris) +#undef HEAP_COMPARE + + /* Sort lower edges of bounding box */ +#define HEAP_COMPARE(A,B) (A->mix[0][k] < B->mix[0][k]) + HEAPSORT(gtri *, &p->sax[k * 2 + 1][0], ntris) +#undef HEAP_COMPARE + } + s->ne_inited = 1; + +//printf("~1 init_ne done\n"); +} + +/* Free everything */ +static void del_gnn(gnn *p) { + int k; + + for (k = 0; k < (3 * 2); k++) { + free (p->sax[k]); + } + + free(p); +} + +/* ===================================================== */ +/* Define the colorspaces white and black point. May be NULL if unknown. */ +/* Note that as in all of the gamut library, we assume that we are in */ +/* an L*a*b* or Jab type color space. */ +static void setwb( +gamut *s, +double *wp, +double *bp, +double *kp +) { + if (wp != NULL) { + s->cs_wp[0] = wp[0]; + s->cs_wp[1] = wp[1]; + s->cs_wp[2] = wp[2]; + } else { + s->cs_wp[0] = 100.0; + s->cs_wp[1] = 0.0; + s->cs_wp[2] = 0.0; + } + + if (bp != NULL) { + s->cs_bp[0] = bp[0]; + s->cs_bp[1] = bp[1]; + s->cs_bp[2] = bp[2]; + } else { + s->cs_bp[0] = 0.0; + s->cs_bp[1] = 0.0; + s->cs_bp[2] = 0.0; + } + + if (kp != NULL) { + s->cs_kp[0] = kp[0]; + s->cs_kp[1] = kp[1]; + s->cs_kp[2] = kp[2]; + } else { + s->cs_kp[0] = s->cs_bp[0]; + s->cs_kp[1] = s->cs_bp[1]; + s->cs_kp[2] = s->cs_bp[2]; + } + + s->cswbset = 1; +} + + +/* Compute the gamut white/black points, assuming */ +/* that the colorspace white/black points have been set. */ +/* The gamut white/black are the points on the colorspace */ +/* white/black axis that have the same L values as the */ +/* extremes within the gamut. */ +static void compgawb(gamut *s) { + int i; + double ff, Lmax, Lmin, LKmin; + + if (s->cswbset == 0 || s->gawbset != 0) + return; /* Nothing to do */ + + Lmax = -1000.0; + Lmin = 1000.0; + + /* Discover min and max L values */ + for (i = 0; i < s->nv; i++) { + if ((s->verts[i]->f & GVERT_SET) == 0 ) + continue; + + if (s->verts[i]->p[0] > Lmax) + Lmax = s->verts[i]->p[0]; + if (s->verts[i]->p[0] < Lmin) + Lmin = s->verts[i]->p[0]; + } + + LKmin = Lmin; + + if (Lmax > s->cs_wp[0]) /* Slightly Strange */ + Lmax = s->cs_wp[0]; + if (Lmin < s->cs_bp[0]) /* Also Slightly strange */ + Lmin = s->cs_bp[0]; + if (LKmin < s->cs_kp[0]) /* Expected */ + LKmin = s->cs_kp[0]; + + /* Locate points along colorspace grey axis */ + /* that correspond to the L extremes */ + ff = (Lmax - s->cs_bp[0])/(s->cs_wp[0] - s->cs_bp[0]); + s->ga_wp[0] = Lmax; + s->ga_wp[1] = ff * (s->cs_wp[1] - s->cs_bp[1]) + s->cs_bp[1]; + s->ga_wp[2] = ff * (s->cs_wp[2] - s->cs_bp[2]) + s->cs_bp[2]; + + ff = (Lmin - s->cs_bp[0])/(s->cs_wp[0] - s->cs_bp[0]); + s->ga_bp[0] = Lmin; + s->ga_bp[1] = ff * (s->cs_wp[1] - s->cs_bp[1]) + s->cs_bp[1]; + s->ga_bp[2] = ff * (s->cs_wp[2] - s->cs_bp[2]) + s->cs_bp[2]; + + ff = (LKmin - s->cs_kp[0])/(s->cs_wp[0] - s->cs_kp[0]); + s->ga_kp[0] = LKmin; + s->ga_kp[1] = ff * (s->cs_wp[1] - s->cs_kp[1]) + s->cs_kp[1]; + s->ga_kp[2] = ff * (s->cs_wp[2] - s->cs_kp[2]) + s->cs_kp[2]; + + s->gawbset = 1; +} + +/* Get the colorspace and gamut white & black points. */ +/* Return pointers may be NULL */ +/* Return non-zero if not possible. */ +static int getwb( +gamut *s, +double *cswp, /* Color space */ +double *csbp, +double *cskp, +double *gawp, /* Gamut */ +double *gabp, +double *gakp +) { +//printf("~1 getwb() called\n"); + if (s->cswbset == 0) { + return 1; + } + + if (cswp != NULL) { + cswp[0] = s->cs_wp[0]; + cswp[1] = s->cs_wp[1]; + cswp[2] = s->cs_wp[2]; + } + + if (csbp != NULL) { + csbp[0] = s->cs_bp[0]; + csbp[1] = s->cs_bp[1]; + csbp[2] = s->cs_bp[2]; + } + + if (cskp != NULL) { + cskp[0] = s->cs_kp[0]; + cskp[1] = s->cs_kp[1]; + cskp[2] = s->cs_kp[2]; + } +//printf("~1 colorspace white %f %f %f, black %f %f %f, kblack %f %f %f\n", s->cs_wp[0], s->cs_wp[1], s->cs_wp[2], s->cs_bp[0], s->cs_bp[1], s->cs_bp[2], s->cs_kp[0],s->cs_kp[1],s->cs_kp[2]); + + if (gawp != NULL || gabp != NULL) { +//printf("~1 computing gamut white & black\n"); + compgawb(s); /* make sure we have gamut white/black available */ + } + + if (gawp != NULL) { + gawp[0] = s->ga_wp[0]; + gawp[1] = s->ga_wp[1]; + gawp[2] = s->ga_wp[2]; + } + + if (gabp != NULL) { + gabp[0] = s->ga_bp[0]; + gabp[1] = s->ga_bp[1]; + gabp[2] = s->ga_bp[2]; + } + + if (gakp != NULL) { + gakp[0] = s->ga_kp[0]; + gakp[1] = s->ga_kp[1]; + gakp[2] = s->ga_kp[2]; + } +//printf("~1 gamut white %f %f %f, black %f %f %f, kblack %f %f %f\n", s->ga_wp[0], s->ga_wp[1], s->ga_wp[2], s->ga_bp[0], s->ga_bp[1], s->ga_bp[2], s->ga_kp[0],s->ga_kp[1],s->ga_kp[2]); + + return 0; +} + + +/* ---------------------------------------------------- */ +/* Per-triangle primitives used to compute brute force */ +/* radial & vector intersection and nearest point, */ +/* as well as gamut surface intersections. */ +/* See if the given triangle intersect the given vector. */ +/* Return 1 if it does, 0 if it doesn't */ +static int vect_intersect( +gamut *s, +double *rvp, /* parameter, 0.0 = p1, 1.0 = p2 */ +double *ip, /* return intersection point */ +double *p1, /* First point of vector (ie black) */ +double *p2, /* Second point of vector (ie white) */ +gtri *t /* Triangle in question */ +) { + double ti; /* Axis parameter value */ + double vv[3]; /* vector vector */ + double ival[3]; /* Intersection value */ + double den; + int j; + + vv[0] = p2[0] - p1[0]; + vv[1] = p2[1] - p1[1]; + vv[2] = p2[2] - p1[2]; + + den = t->pe[0] * vv[0] + t->pe[1] * vv[1] + t->pe[2] * vv[2]; + if (fabs(den) < 1e-10) { + return 0; + } + + /* Compute the intersection of the grey axis vector with the triangle plane */ + ti = -(t->pe[0] * p1[0] + t->pe[1] * p1[1] + t->pe[2] * p1[2] + t->pe[3])/den; + + /* Compute the actual intersection point */ + ival[0] = p1[0] + ti * vv[0]; + ival[1] = p1[1] + ti * vv[1]; + ival[2] = p1[2] + ti * vv[2]; + + /* Check if the intersection point is within the triangle */ + for (j = 0; j < 3; j++) { + double ds; + ds = t->ee[j][0] * (ival[0] - s->cent[0]) /* Convert to relative for edge check */ + + t->ee[j][1] * (ival[1] - s->cent[1]) + + t->ee[j][2] * (ival[2] - s->cent[2]) + + t->ee[j][3]; + if (ds > 1e-8) { + return 0; /* Not within triangle */ + } + } +//printf("~1 vect_intersect got intersection with tri %d at %f\n",t->n,ti); + + /* Got an intersection point */ + ip[0] = ival[0]; + ip[1] = ival[1]; + ip[2] = ival[2]; + + *rvp = ti; + + return 1; +} + +/* Given a point and a triangle, return the closest point on */ +/* the triangle closest to the given point. Also return the distance squared */ +/* (Doesn't depend on triangle edge info) */ +static double ne_point_on_tri( +gamut *s, +gtri *t, /* Triangle to use */ +double *out, /* Absolute output point */ +double *in /* Absolute input point */ +) { + int j; + double rv; + double bdist; + + /* Compute the point on the triangles plane, that is orthogonal */ + /* (closest) to the target point. */ + rv = (t->pe[0] * in[0] + t->pe[1] * in[1] + t->pe[2] * in[2] + t->pe[3])/ + (t->pe[0] * t->pe[0] + t->pe[1] * t->pe[1] + t->pe[2] * t->pe[2]); + + out[0] = in[0] - rv * t->pe[0]; + out[1] = in[1] - rv * t->pe[1]; + out[2] = in[2] - rv * t->pe[2]; + + /* Check if the closest point is within this triangle */ + for (j = 0; j < 3; j++) { + double ds; + ds = t->ee[j][0] * (out[0] - s->cent[0]) /* Convert to relative for edge check */ + + t->ee[j][1] * (out[1] - s->cent[1]) + + t->ee[j][2] * (out[2] - s->cent[2]) + + t->ee[j][3]; + if (ds > 1e-8) { + break; /* Not within triangle */ + } + } + if (j >= 3) { /* It's OK */ + return rv * rv; /* rv is distance since pe length is 1.0 */ + } + + /* Not in triangle, so find closest point along any edge, */ + /* or at the verticies. (don't use edge info, it may not be set up) */ + bdist = 1e38; + for (j = 0; j < 3; j++) { /* For each edge */ + gvert *v0 = t->v[j], *v1 = t->v[j >= 2 ? 0 : j+1]; + int k; + double nu, de, ds; + for (de = 0.0, k = 0; k < 3; k++) { + double tt = v1->p[k] - v0->p[k]; + de += tt * tt; + } + for (nu = 0.0, k = 0; k < 3; k++) + nu += (v1->p[k] - v0->p[k]) * (in[k] - v0->p[k]); + + ds = nu/de; + + if (ds >= 0.0 && ds <= 1.0) { /* Valid edge */ + double tout[3], ss; + for (ss = 0.0, k = 0; k < 3; k++) { + tout[k] = v0->p[k] + ds * (v1->p[k] - v0->p[k]); + ss += (in[k] - tout[k]) * (in[k] - tout[k]); + } + if (ss < bdist) { + bdist = ss; + out[0] = tout[0]; + out[1] = tout[1]; + out[2] = tout[2]; + } + } + } + + for (j = 0; j < 3; j++) { /* For each vertex */ + int k; + double ss; + for (ss = 0.0, k = 0; k < 3; k++) { + double tt; + tt = in[k] - t->v[j]->p[k]; + ss += tt * tt; + } + + if (ss < bdist) { + bdist = ss; + out[0] = t->v[j]->p[0]; + out[1] = t->v[j]->p[1]; + out[2] = t->v[j]->p[2]; + } + } + + return bdist; +} + +/* ----------------------------------------------------- */ +/* Arbitrary vector intersect */ + +/* Given a vector, find the two extreme intersection with */ +/* the gamut surface using a brute force search. */ +/* Return 0 if there is no intersection */ +static int compute_vector_isect_bf( +gamut *s, +double *p1, /* First point (ie black) */ +double *p2, /* Second point (ie white) */ +double *omin, /* Return gamut surface points, min = closest to p1 */ +double *omax, /* max = farthest from p1 */ +double *omnt, /* Return parameter values for p1 and p2, 0 being at p1, */ +double *omxt, /* and 1 being at p2 */ +gtri **omntri, /* Return the intersection triangles */ +gtri **omxtri +) { + gtri *tp, *t0, *t1; + double ip[3], min[3], max[3]; + double mint, maxt; + int j; + + if IS_LIST_EMPTY(s->tris) + triangulate(s); + + maxt = -1e68; /* Setup to find min/max */ + mint = 1e68; + + tp = s->tris; + FOR_ALL_ITEMS(gtri, tp) { + double rv; + if (vect_intersect(s, &rv, ip, p1, p2, tp)) { + if (rv < mint) { + min[0] = ip[0]; + min[1] = ip[1]; + min[2] = ip[2]; + mint = rv; + t0 = tp; + } + if (rv > maxt) { + max[0] = ip[0]; + max[1] = ip[1]; + max[2] = ip[2]; + maxt = rv; + t1 = tp; + } + } + } END_FOR_ALL_ITEMS(tp); + + if (((omin != NULL || omnt != NULL || omntri != NULL) && mint == 1e68) + || ((omax != NULL || omxt != NULL || omxtri != NULL) && maxt == -1e68)) { + return 0; + } + + if (omin != NULL) + for (j = 0; j < 3; j++) + omin[j] = min[j]; + + if (omax != NULL) + for (j = 0; j < 3; j++) + omax[j] = max[j]; + + if (omnt != NULL) + *omnt = mint; + + if (omxt != NULL) + *omxt = maxt; + + if (omntri != NULL) + *omntri = t0; + + if (omxtri != NULL) + *omxtri = t1; + + return 1; +} + + +#ifdef INTERSECT_DEBUG + +#define ISDBG(xxx) if (deb_insect) printf xxx + +int deb_insect = 0; /* Do vrml plot */ + +/* Debug - given a BSP node, add all the triangles vertexes indexes */ +/* below this node to the diagnosti wrl */ +static void debug_bsp_triangl_wrl( +gamut *s, +gbsp *np, /* BSP node pointer we're at */ +vrml *wrl /* Diagnostic plot */ +) { + if (np->tag == 1) { /* It's a BSP node */ + gbspn *n = (gbspn *)np; + + debug_bsp_triangl_wrl(s, n->po, wrl); + debug_bsp_triangl_wrl(s, n->ne, wrl); + + } else { + int nt; /* Number of triangles in list */ + gtri **tpp; /* Pointer to list of triangles */ + int i, j; + + if (np->tag == 2) { /* It's a triangle */ + tpp = &((gtri *)np); + nt = 1; + } else if (np->tag == 3) { /* It's a triangle list */ + gbspl *n = (gbspl *)np; + tpp = n->t; + nt = n->nt; + } + /* Go through the list of triangles and intersect with vector */ + for (i = 0; i < nt; i++, tpp++) { + gtri *t = *tpp; + int ix[3]; + + ix[0] = t->v[0]->tn; + ix[1] = t->v[1]->tn; + ix[2] = t->v[2]->tn; + + wrl->add_triangle(wrl, 0, ix); + } + } +} +#else /* !INTERSECT_DEBUG */ +# define ISDBG(xxx) +#endif /* !INTERSECT_DEBUG */ + +/* Recursive vector intersect using BSP accelleration. */ +static void vector_isect_rec( +gamut *s, +gbsp *np, /* BSP node pointer we're at */ +double *vb, /* Center relative base point of vector */ +double *vv, /* Vector direction from base */ +double t0, /* Start parameter value of line */ +double rs0, /* line start point radius squared */ +double t1, /* End parameter value of line */ +double rs1, /* line end point radius squared */ +double tc, /* Parameter value of closest point on line to center */ +double rsc, /* Radius squared of closest point on line to center */ +double rse0, /* Effective radius squared minimum */ +double rse1, /* Effective radius squared maximum */ +gispnt *lp, /* List to set intersections in. */ +int ll, /* Size of list. 0 == 2, min and max only */ +int *lu /* Number used in list */ +) { + double den; /* Intersection denominator */ + double ti; /* Intersection parameter value */ + double rsi; /* Radius squared of intersection point */ + double ip[3]; /* Intersection point */ + +#ifdef INTERSECT_DEBUG + if (ll == 0) + printf("\nvector_isect_rec got seg %f - %f (%e), best %f, %f\n",t0,t1,t1-t0,lp[0].pv,lp[1].pv); + else + printf("\nvector_isect_rec got seg %f - %f (%e), no isects %d\n",t0,t1,t1-t0,*lu); + printf(" rs %f, %f, tc %f, rsc %f, rse0 %f, rse1 %f, BSP rs %f - %f\n",rs0, rs1, tc, rsc, rse0, rse1, np->rs0, np->rs1); +#endif +#ifdef INTERSECT_DEBUG + if (deb_insect) { + vrml *wrl = NULL; + double cc[3] = { 1.0, 1.0, 0.0 }; + double red[3] = { 1.0, 0.0, 0.0 }; + double green[3] = { 0.0, 1.0, 0.0 }; + double blue[3] = { 0.0, 0.0, 1.0 }; + double p1[3], p2[3]; + int i; + + unlink("isect2.wrl"); + rename("isect.wrl", "isect2.wrl"); + + if ((wrl = new_vrml("isect.wrl", 0)) == NULL) + error("New vrml failed"); + + /* The triangles below the BSP */ + for (i = 0; i < s->nv; i++) { + if (!(s->verts[i]->f & GVERT_TRI)) + continue; + wrl->add_vertex(wrl, 0, s->verts[i]->p); + } + debug_bsp_triangl_wrl(s, np, wrl); + wrl->make_triangles(wrl, 0, 0.0, cc); + + /* The segment. The vrml browser can go crazy if the line */ + /* is too long, so limit it to +/- 100 */ + { + double vbl = icmNorm3(vb); + double tt0 = t0, tt1 = t1; + if (tt0 > 0.0) { + if ((tt0 * vbl) > 100.0) + tt0 = 100.0/vbl; + } else { + if ((tt0 * vbl) < -100.0) + tt0 = -100.0/vbl; + } + if (tt1 > 0.0) { + if ((tt1 * vbl) > 100.0) + tt1 = 100.0/vbl; + } else { + if ((tt1 * vbl) < -100.0) + tt1 = -100.0/vbl; + } + for (i = 0; i < 3; i++) { + p1[i] = s->cent[i] + vb[i] + tt0 * vv[i]; + p2[i] = s->cent[i] + vb[i] + tt1 * vv[i]; + } + } + wrl->add_col_vertex(wrl, 1, p1, red); + wrl->add_col_vertex(wrl, 1, p2, red); + wrl->make_lines(wrl, 1, 2); + + /* Add two initial points */ + for (i = 0; i < 3; i++) { + p1[i] = s->cent[i] + vb[i] + 0.0 * vv[i]; + p2[i] = s->cent[i] + vb[i] + 1.0 * vv[i]; + } + wrl->add_marker(wrl, p1, green, 0.5); + wrl->add_marker(wrl, p2, blue, 0.5); + + wrl->del(wrl); + printf("Waiting for input after writing 'isect.wrl':\n"); + getchar(); + + } +#endif + + if (np->tag == 1) { /* It's a BSP node */ + int j; + gbspn *n = (gbspn *)np; + double ds; + gbsp *n1, *n2; /* next bsp to recurse to */ + double ti1, ti2; /* Intersection parameters for recursion */ + double rse1_0, rse1_1, rse2_0, rse2_1; + + ISDBG(("vector_isect_rec at bsp node %d\n")); + + /* Try and compute intersection with BSP */ + + den = n->pe[0] * vv[0] + n->pe[1] * vv[1] + n->pe[2] * vv[2]; + if (fabs(den) > 1e-12) { + /* Compute the intersection point */ + ti = -(n->pe[0] * vb[0] + n->pe[1] * vb[1] + n->pe[2] * vb[2] + n->pe[3])/den; + ISDBG(("intersects BSP plane at ti %f\n",ti)); + } + if (fabs(den) < 1e-12 || ti < (t0 - 1e-6) || ti > (t1 + 1e-6)) { /* Doesn't intersect */ + double ds; /* or intersects outside segment */ + + ISDBG(("doesn't intersect BSP plane within segment\n")); + /* Figure which side of the BSP segment is on */ + ds = n->pe[0] * (vb[0] + 0.5 * (t0 + t1) * vv[0]) + + n->pe[1] * (vb[1] + 0.5 * (t0 + t1) * vv[1]) + + n->pe[2] * (vb[2] + 0.5 * (t0 + t1) * vv[2]) + + n->pe[3]; + + ISDBG(("recursing down side that segment is on\n")); + /* And recurse approproately if it can be improved */ + /* ???? Shouldn't we recurse down both if ds ~= 0.0 ? */ + if (ds >= 0.0) { + if (rse0 <= n->po->rs1 + && rse1 >= n->po->rs0 + && (ll > 0 || t0 < lp[0].pv || t1 > lp[1].pv)) + vector_isect_rec(s, n->po, vb, vv, t0, rs0, t1, rs1, tc, rsc, rse0, rse1, + lp, ll, lu); + } else { + if (rse0 <= n->ne->rs1 + && rse1 >= n->ne->rs0 + && (ll > 0 || t0 < lp[0].pv || t1 > lp[1].pv)) + vector_isect_rec(s, n->ne, vb, vv, t0, rs0, t1, rs1, tc, rsc, rse0, rse1, + lp, ll, lu); + } + ISDBG(("vector_isect_rec returning\n")); + return; + } + + /* Compute radius squared to center point at split point */ + for (rsi = 0.0, j = 0; j < 3; j++) { + den = vb[j] + ti * vv[j]; + rsi += den * den; + } + + /* Compute the effective radius squared range for each segment */ + rse1_0 = rs0; + rse1_1 = rs0; + if (rsi < rse1_0) + rse1_0 = rsi; + if (rsi > rse1_1) + rse1_1 = rsi; + if (tc >= t0 && tc <= ti) { /* Closest point is within segment */ + if (rsc < rse1_0) + rse1_0 = rsc; + if (rsc > rse1_1) + rse1_1 = rsc; + } + + rse2_0 = rsi; + rse2_1 = rsi; + if (rs1 < rse2_0) + rse2_0 = rs1; + if (rs1 > rse2_1) + rse2_1 = rs1; + if (tc >= ti && tc <= t1) { /* Closest point is within segment */ + if (rsc < rse2_0) + rse2_0 = rsc; + if (rsc > rse2_1) + rse2_1 = rsc; + } + + /* Test t0-1.0 to see what side of the BSP t0..ti is. */ + ip[0] = vb[0] + (t0-1.0) * vv[0]; + ip[1] = vb[1] + (t0-1.0) * vv[1]; + ip[2] = vb[2] + (t0-1.0) * vv[2]; + + ds = n->pe[0] * ip[0] + + n->pe[1] * ip[1] + + n->pe[2] * ip[2] + + n->pe[3]; + + /* Because we're intersecting a line segment, we don't */ + /* have to recurse down both sides of the BSP tree ? */ + if (ds >= 0.0) { + n1 = n->po; + n2 = n->ne; + } else { + n1 = n->ne; + n2 = n->po; + } + + /* Make sure that touching segments get properly tested */ + ti1 = ti + 1e-7; + ti2 = ti - 1e-7; + + /* Split the line into two segments and recurse for each. */ + /* Don't recurse if the line can't improve either min or max. */ + if (rse1_0 <= n1->rs1 + && rse1_1 >= n1->rs0 + && (ll > 0 || t0 < lp[0].pv || ti1 > lp[1].pv)) { + ISDBG(("recursing segment 1/2 %f .. %f\n",t0,ti1)); + vector_isect_rec(s, n1, vb, vv, t0, rs0, ti1, rsi, tc, rsc, rse1_0, rse1_1, + lp, ll, lu); + } +#ifdef INTERSECT_DEBUG + else if (deb_insect) { + printf("Skipped seg 1/2 because rse1_0 %f > n1->rs1 %f ? || rse1_1 %f < n1->rs0 %f\n",rse1_0,n1->rs1,rse1_1, n1->rs0); + if (ll == 0) + printf("|| ti1 %f <= t0 %f ? || ti1 %f <= omxt %f && t0 %f <= omnt %f ?\n",ti1,t0,ti1,lp[1].pv,t0,lp[0].pv); + else + printf("|| ti1 %f <= t0 %f ?\n",ti1,t0); + } +#endif + if (rse2_0 <= n2->rs1 + && rse2_1 >= n2->rs0 + && (ll > 0 || ti2 < lp[0].pv || t1 > lp[1].pv)) { + ISDBG(("recursing segment 2/2 %f .. %f\n",ti2,t1)); + vector_isect_rec(s, n2, vb, vv, ti2, rsi, t1, rs1, tc, rsc, rse2_0, rse2_1, + lp, ll, lu); + } +#ifdef INTERSECT_DEBUG + else if (deb_insect) { + printf("Skipped seg 2/2 because rse2_0 %f > n2->rs1 %f ? || rse2_1 %f < n2->rs0 %f\n",rse2_0,n2->rs1,rse2_1, n2->rs0); + if (ll == 0) + printf("|| t1 %f <= ti2 %f ? || t1 %f <= omxt %f && ti2 %f <= omnt %f ?\n",t1,ti2,t1,lp[1].pv,ti2,lp[0].pv); + else + printf("|| t1 %f <= ti2 %f ?\n",t1,ti2); + } +#endif + + ISDBG(("vector_isect_rec returning\n")); + return; + + /* It's a list of triangles */ + } else { + int nt; /* Number of triangles in list */ + gtri **tpp; /* Pointer to list of triangles */ + int i, j; + + if (np->tag == 2) { /* It's a triangle */ + tpp = (gtri **)&np; + nt = 1; + } else if (np->tag == 3) { /* It's a triangle list */ + gbspl *n = (gbspl *)np; + tpp = n->t; + nt = n->nt; + } + ISDBG(("vector_isect_rec at triangle(s) %d\n",nt)); + /* Go through the list of triangles and intersect with vector */ + for (i = 0; i < nt; i++, tpp++) { + double bds; + gtri *t = *tpp; + + ISDBG(("triangle no %d\n",t->n)); + + den = t->pe[0] * vv[0] + t->pe[1] * vv[1] + t->pe[2] * vv[2]; + if (fabs(den) < 1e-12) { + ISDBG(("segment is tangent to triangle\n")); + continue; + } + + /* Compute the intersection of vector with the BSP plane */ + ti = -(t->pe[0] * (vb[0] + s->cent[0]) + + t->pe[1] * (vb[1] + s->cent[1]) + + t->pe[2] * (vb[2] + s->cent[2]) + + t->pe[3])/den; + ISDBG(("segment intersects at %f\n",ti)); + + /* Compute the actual (center relative) intersection point */ + ip[0] = vb[0] + ti * vv[0]; + ip[1] = vb[1] + ti * vv[1]; + ip[2] = vb[2] + ti * vv[2]; + ISDBG(("triangle intersection point %f %f %f\n",ip[0]+s->cent[0],ip[1]+s->cent[1],ip[2]+s->cent[2])); + + /* Check if the intersection point is within the triangle */ + bds = -1e6; + for (j = 0; j < 3; j++) { + double ds; + ds = t->ee[j][0] * ip[0] + + t->ee[j][1] * ip[1] + + t->ee[j][2] * ip[2] + + t->ee[j][3]; + if (ds > 1e-8) + break; /* Not within triangle */ + if (ds > bds) + bds = ds; + } + if (j < 3) { + ISDBG(("intersection not within triangle\n")); + continue; /* Not within triangle, so ignore */ + } + + /* Add intersection to list */ + if (ll > 0) { /* List of all */ + if (*lu < ll) { + lp[*lu].pv = ti; + icmAdd3(lp[*lu].ip,ip,s->cent); /* Abs. intersection point */ + lp[*lu].dir = den > 0.0 ? 1 : 0; + lp[*lu].edge = bds > 0.0 ? 1 : 0; + lp[*lu].tri = t; + ISDBG(("new isect %d: pv %f, dir %d, edge %d\n",*lu,ti,lp[*lu].dir,lp[*lu].edge)); + (*lu)++; + } else { + ISDBG(("new isect %d: List Too Short %d!!!\n",*lu,ll)); + } + } else { /* Bigest/smallest list of 2 */ + if (ti < lp[0].pv) { + ISDBG(("new min %f\n",ti)); + lp[0].pv = ti; + icmAdd3(lp[0].ip,ip,s->cent); /* Abs. intersection point */ + lp[0].dir = den > 0.0 ? 1 : 0; + lp[0].edge = bds > 0.0 ? 1 : 0; + lp[0].tri = t; + } + if (ti > lp[1].pv) { + ISDBG(("new max %f\n",ti)); + lp[1].pv = ti; + icmAdd3(lp[1].ip,ip,s->cent); /* Abs. intersection point */ + lp[1].dir = den > 0.0 ? 1 : 0; + lp[1].edge = bds > 0.0 ? 1 : 0; + lp[1].tri = t; + } + } + } + ISDBG(("vector_isect_rec returning\n")); + return; + } +} + +/* Re-evaluate the intersections with an offset */ +static void reevaluate_isectns( +gamut *s, +gispnt *lp, /* List to set intersections in. */ +int n, /* Number to re-evaluate */ +double offset, /* Amount to offset vector */ +double *ivb, /* Center relative base point of vector */ +double *vv /* Vector direction from base */ +) { + double vb[3]; /* Offset vb */ + double sv; + int sj, j, i; + + /* Decide which way to offset base */ + sv = -1e20; + for (j = 0; j < 3; j++) { + if (fabs(vv[j]) > sv) { /* Locate largest direction */ + sv = fabs(vv[j]); + sj = j; + } + } + + /* Apply the offset to other than largest */ + for (j = 0; j < 3; j++) { + vb[j] = ivb[j]; + if (j == sj) + continue; + vb[j] += offset; + } + + for (i = 0; i < n; i++) { + double den, ti, ip[3], bds; + gtri *t = lp[i].tri; + + lp[i].dir = 0; + lp[i].edge = 2; /* Assume no intersect */ + + den = t->pe[0] * vv[0] + t->pe[1] * vv[1] + t->pe[2] * vv[2]; + if (fabs(den) < 1e-12) { + continue; + } + + /* Compute the intersection of vector with the BSP plane */ + ti = -(t->pe[0] * (vb[0] + s->cent[0]) + + t->pe[1] * (vb[1] + s->cent[1]) + + t->pe[2] * (vb[2] + s->cent[2]) + + t->pe[3])/den; + + /* Compute the actual intersection point */ + ip[0] = vb[0] + ti * vv[0]; + ip[1] = vb[1] + ti * vv[1]; + ip[2] = vb[2] + ti * vv[2]; + + /* Check if the intersection point is within the triangle */ + bds = -1e6; + for (j = 0; j < 3; j++) { + double ds; + ds = t->ee[j][0] * ip[0] + + t->ee[j][1] * ip[1] + + t->ee[j][2] * ip[2] + + t->ee[j][3]; + if (ds > 1e-8) + break; /* Not within triangle */ + if (ds > bds) + bds = ds; + } + if (j < 3) { + continue; /* Not within triangle, so ignore */ + } + + /* Update intersection info */ + lp[i].dir = den > 0.0 ? 1 : 0; + lp[i].edge = bds > 0.0 ? 1 : 0; + } +} + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ +#ifdef INTERSECT_VERIFY +#define compute_vector_isect _compute_vector_isect +#endif + +/* Given a vector, find the two extreme intersection with */ +/* the gamut surface. */ +/* BSP accellerated version */ +/* Return 0 if there is no intersection */ +static int compute_vector_isect( +gamut *s, +double *p1, /* First point (ie param value 0.0) */ +double *p2, /* Second point (ie param value 1.0) */ +double *omin, /* Return gamut surface points, min = closest to p1 */ +double *omax, /* max = farthest from p1 */ +double *omnt, /* Return parameter values for p1 and p2, 0 being at p1, */ +double *omxt, /* and 1 being at p2 */ +gtri **omntri, /* Return the intersection triangles */ +gtri **omxtri +) { + gtri *tp; + double vb[3], vv[3]; /* Center relative base of vector, vector of vector */ + double tt, t0, rs0, t1, rs1, tc, rsc; + double rse0, rse1; /* Effective radius squared min & max */ + gispnt islist[2]; /* min and max result */ + int lu = 0, j; + int rv = 0; + + if IS_LIST_EMPTY(s->tris) + triangulate(s); + + if (s->lu_inited == 0) + init_lu(s); /* Init BSP search tree */ + + /* Convert twp points to center relative base + vector direction */ + for (tt = 0.0, j = 0; j < 3; j++) { + vv[j] = p2[j] - p1[j]; + tt += vv[j] * vv[j]; + vb[j] = p1[j] - s->cent[j]; /* relative to gamut center */ + } + /* If vector is too small to have a valid direction */ + if (tt < 1e-12) { + return 0; + } + + islist[0].pv = 1e68; + islist[1].pv = -1e68; /* Setup to find min/max */ + t0 = -1e6; + t1 = 1e6; + + /* Compute radius range of segment */ + for (rs0 = rs1 = 0.0, j = 0; j < 3; j++) { + tt = vb[j] + t0 * vv[j]; + rs0 += tt * tt; + tt = vb[j] + t1 * vv[j]; + rs1 += tt * tt; + } + + /* Compute the point closest to the center */ + tc = -(vv[0] * vb[0] + vv[1] * vb[1] + vv[2] * vb[2]) + / (vv[0] * vv[0] + vv[1] * vv[1] + vv[2] * vv[2]); + +#ifdef INTERSECT_VERIFY + /* Check that this is correct */ + for (tt = 0.0, j = 0; j < 3; j++) { + double pp; + pp = vb[j] + tc * vv[j]; + tt += pp * vv[j]; + } + if (fabs(tt) > 1e-5) + error("Failed to locate closest point on vector"); +#endif /* INTERSECT_VERIFY */ + + for (rsc = 0.0, j = 0; j < 3; j++) { + tt = vb[j] + tc * vv[j]; + rsc += tt * tt; + } + + /* Compute the effective min/max radius squared */ + rse0 = rs0; + rse1 = rs0; + if (rs1 < rse0) + rse0 = rs1; + if (rs1 > rse1) + rse1 = rs1; + if (tc >= t0 && tc <= t1) { /* Closest point is within segment */ + if (rsc < rse0) + rse0 = rsc; + if (rsc > rse1) + rse1 = rsc; + } + + vector_isect_rec(s, s->lutree, vb, vv, t0, rs0, t1, rs1, tc, rsc, rse0, rse1, + islist, 0, &lu); + + /* If we failed to locate a requested intersection */ + if (((omin != NULL || omnt != NULL || omntri != NULL) && islist[0].pv == 1e68) + || ((omax != NULL || omxt != NULL || omxtri != NULL) && islist[1].pv == -1e68)) { + rv = 0; + + } else { + + if (omin != NULL) { + for (j = 0; j < 3; j++) + icmCpy3(omin,islist[0].ip); + ISDBG(("Fast min = %f %f %f\n", omin[0], omin[1], omin[2])); + } + + if (omax != NULL) { + icmCpy3(omax,islist[1].ip); + ISDBG(("Fast max = %f %f %f\n", omax[0], omax[1], omax[2])); + } + + if (omnt != NULL) + *omnt = islist[0].pv; + + if (omxt != NULL) + *omxt = islist[1].pv; + + if (omntri != NULL) + *omntri = islist[0].tri; + + if (omxtri != NULL) + *omxtri = islist[1].tri; + + rv = 1; + } + + return rv; +} + +#ifdef INTERSECT_VERIFY +#undef compute_vector_isect + +/* Verifying version of above */ +static int compute_vector_isect( +gamut *s, +double *p1, /* First point (ie param value 0.0) */ +double *p2, /* Second point (ie param value 1.0) */ +double *omin, /* Return gamut surface points, min = closest to p1 */ +double *omax, /* max = farthest from p1 */ +double *omnt, /* Return parameter values for p1 and p2, 0 being at p1, */ +double *omxt, /* and 1 being at p2 */ +gtri **omintri, /* Return the intersection triangles */ +gtri **omaxtri +) { + int rv, _rv; + double _omin[3]; + double _omax[3]; + double _omnt; + double _omxt; + gtri *_omintri; + gtri *_omaxtri; + int fail = 0; + + ISDBG(("\n\n###########################################\n")); + + /* Call the routine we're checking */ + rv = _compute_vector_isect(s, p1, p2, omin, omax, omnt, omxt, omintri, omaxtri); + + _rv = compute_vector_isect_bf(s, p1, p2, _omin, _omax, &_omnt, &_omxt, &_omintri, &_omaxtri); + + if (rv != _rv) { + warning("compute_vector_isect verify: rv %d != _rv %d\n",rv,_rv); + fail = 1; + } + + if (rv == 1) { + int j; + if (omnt != NULL) + if (fabs (*omnt - _omnt) > 1e-4) { + warning("compute_vector_isect verify:\n omnt %f != _omnt %f\n",*omnt,_omnt); + fail = 1; + } + if (omxt != NULL) + if (fabs (*omxt - _omxt) > 1e-4) { + warning("compute_vector_isect verify:\n omxt %f != _omxt %f\n",*omxt,_omxt); + fail = 1; + } + if (omin != NULL) { + ISDBG(("bf min = %f %f %f\n", _omin[0], _omin[1], _omin[2])); + for (j = 0; j < 3; j++) { + if (fabs (omin[j] - _omin[j]) > 1e-4) + break; + } + if (j < 3) { + warning("compute_vector_isect verify:\n omin %f %f %f != _omin %f %f %f\n", omin[0], omin[1], omin[2], _omin[0], _omin[1], _omin[2]); + fail = 1; + } + } + if (omax != NULL) { + ISDBG(("bf max = %f %f %f\n", _omax[0], _omax[1], _omax[2])); + for (j = 0; j < 3; j++) { + if (fabs (omax[j] - _omax[j]) > 1e-4) + break; + } + if (j < 3) { + warning("compute_vector_isect verify:\n omax %f %f %f != _omax %f %f %f\n", omax[0], omax[1], omax[2], _omax[0], _omax[1], _omax[2]); + fail = 1; + } + } +#ifdef NEVER + if (omintri != NULL) + if (*omintri != _omintri) { + warning("compute_vector_isect verify:\n omintri %d != _omintri %d\n",(*omintri)->n, _omintri->n); + } + if (omaxtri != NULL) + if (*omaxtri != _omaxtri) { + warning("compute_vector_isect verify:\n omaxtri %d != _omaxtri %d\n",(*omaxtri)->n, _omaxtri->n); + } +#endif /* NEVER */ + } + if (fail) { +#ifdef INTERSECT_DEBUG + printf("Re-running intersect with debug trace on\n"); + deb_insect = 1; + _compute_vector_isect(s, p1, p2, _omin, _omax, &_omnt, &_omxt, &_omintri, &_omaxtri); +#endif /* INTERSECT_DEBUG */ + error("Verify failed"); + } + return rv; +} + +#endif /* INTERSECT_VERIFY */ + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ +/* Compute all the intersection pairs of the vector p1->p2 with */ +/* the gamut surface. lp points to an array of ll gispnt to be */ +/* filled in. If the list is too small, intersections will be */ +/* arbitrarily ignored. */ +/* Return the number of intersections set in list. This will always be even, */ +/* and there will be pairs of in/out intersections in the direction p1->p2. */ +static int compute_vector_isectns( +gamut *s, +double *p1, /* First point (ie param value 0.0) */ +double *p2, /* Second point (ie param value 1.0) */ +gispnt *lp, /* List to return in/out intersection pairs */ +int ll /* Size of list. */ +) { + gtri *tp; + double vb[3], vv[3]; /* Center relative base of vector, vector of vector */ + double tt, t0, rs0, t1, rs1, tc, rsc, vscale; + double rse0, rse1; /* Effective radius squared min & max */ + int lu = 0, i, j, k, m, pdir; + int rv = 0; + +#ifdef INTERSECT_DEBUG + printf("compute_vector_isectns p1 %f %f %f, p2 %f %f %f\n", p1[0], p1[1], p1[2], p2[0], p2[1], p2[2]); +#endif + if IS_LIST_EMPTY(s->tris) + triangulate(s); + + if (s->lu_inited == 0) + init_lu(s); /* Init BSP search tree */ + + /* Convert twp points to relative base + vector direction */ + for (tt = 0.0, j = 0; j < 3; j++) { + vv[j] = p2[j] - p1[j]; + tt += vv[j] * vv[j]; + vb[j] = p1[j] - s->cent[j]; /* relative to gamut center */ + } + /* If vector is too small to have a valid direction, give up */ + if (tt < 1e-12) { +#ifdef INTERSECT_DEBUG + printf("p2 too close to p1\n"); +#endif + return 0; + } + + /* Scale factor to make parameter delta = gamut space unit */ + vscale = 1.0/sqrt(tt); + + t0 = -1e6 * vscale; /* Set the parameter search space */ + t1 = 1e6 * vscale; + + /* Compute radius range of segment */ + for (rs0 = rs1 = 0.0, j = 0; j < 3; j++) { + tt = vb[j] + t0 * vv[j]; + rs0 += tt * tt; + tt = vb[j] + t1 * vv[j]; + rs1 += tt * tt; + } + + /* Compute the point closest to the center */ + tc = -(vv[0] * vb[0] + vv[1] * vb[1] + vv[2] * vb[2]) + / (vv[0] * vv[0] + vv[1] * vv[1] + vv[2] * vv[2]); + +#ifdef INTERSECT_DEBUG + /* Check that this is correct */ + for (tt = 0.0, j = 0; j < 3; j++) { + double pp; + pp = vb[j] + tc * vv[j]; + tt += pp * vv[j]; + } + if (fabs(tt) > 1e-5) + error("Failed to locate closest point on vector"); +#endif /* INTERSECT_DEBUG */ + + for (rsc = 0.0, j = 0; j < 3; j++) { + tt = vb[j] + tc * vv[j]; + rsc += tt * tt; + } + + /* Compute the effective min/max radius squared */ + rse0 = rs0; + rse1 = rs0; + if (rs1 < rse0) + rse0 = rs1; + if (rs1 > rse1) + rse1 = rs1; + if (tc >= t0 && tc <= t1) { /* Closest point is within segment */ + if (rsc < rse0) + rse0 = rsc; + if (rsc > rse1) + rse1 = rsc; + } + + /* Recursively locate all the triangle intersections using BSP */ + vector_isect_rec(s, s->lutree, vb, vv, t0, rs0, t1, rs1, tc, rsc, rse0, rse1, + lp, ll, &lu); + + if (lu <= 1) { +#ifdef INTERSECT_DEBUG + printf("%d intersections found\n",lu); +#endif + return 0; /* Too few to be useful */ + } + + /* Now we need to turn th raw intersections into sanitized segment pairs. */ + + /* Sort the intersections by parameter value */ +#define HEAP_COMPARE(A,B) (A.pv < B.pv) + HEAPSORT(gispnt, lp, lu) +#undef HEAP_COMPARE + +#ifdef INTERSECT_DEBUG + printf("Before sanitizing %d\n",lu); + for (i = 0; i < lu; i++) + printf("Isect %d: pv %f, dir %d, edge %d, tri %d\n",i,lp[i].pv,lp[i].dir,lp[i].edge,lp[i].tri->n); +#endif + + /* Remove any duplicate intersections (triangles) */ + for (j = i = 0; i < lu; i++) { + + for (k = i+1; k < lu; k++) { + if (lp[k].tri == lp[i].tri) { + lp[k].edge &= lp[i].edge; /* Keep non-edge status */ + break; + } + } + if (k < lu) + continue; /* Skip this one */ + + /* Accept this intersection */ + memmove(&lp[j], &lp[i], sizeof(gispnt)); + j++; + } + lu = j; + + if (lu <= 1) { +#ifdef INTERSECT_DEBUG + printf("%d intersections after removing duplicates\n",lu); +#endif + return 0; /* Too few to be useful */ + } + +#ifdef INTERSECT_DEBUG + printf("After removing duplicates %d\n",lu); + for (i = 0; i < lu; i++) { + printf("Isect %d: pv %f, dir %d, edge %d, tri %d\n",i,lp[i].pv,lp[i].dir,lp[i].edge,lp[i].tri->n); + } +#endif + + /* Sanitize the intersections. */ + /* We must end up with in/out segment pairs. */ + /* j = output index, i = current index */ + pdir = 0; /* Previous isection direction = "out" */ + for (j = i = 0; i < lu;) { + int nin, nout; /* Number fully in/out */ + int inx, outx; /* Indexes of representative in/out */ + int npin, npout; /* Number partially in/out */ + int pinx, poutx; /* Indexes of representative in/out */ + +//printf("~1 at %d out of %d, %d saved\n",i,lu,j); + /* Two tries, re-evaluate before second try */ + for (m = 0; m < 2; m++) { + + /* See how many we have at the next pv, and */ + /* decide if they're in, or out or both. */ + nin = nout = npin = npout = 0; + for (k = i; k < lu; k++) { + if (i != k && fabs((lp[i].pv - lp[k].pv) * vscale) >= 0.0001) + break; + if (lp[k].dir) { /* In */ + if (lp[k].edge == 0) { + nin++; + inx = k; + } else if (lp[k].edge == 1) { + npin++; + pinx = k; + } + } else { /* Out */ + if (lp[k].edge == 0) { + nout++; + outx = k; + } else if (lp[k].edge == 1) { + npout++; + poutx = k; + } + } + } + + if (m == 1 /* We've already re-evaluated */ + || (k - i) <= 2 /* Not worth re-evaluating */ + || (npin == 0 && npout == 0)) /* All definite */ + break; + + /* Re-evaluate the intersections with an offset. */ + /* (We need this because often the point of interest */ + /* is a vertex, and there are lots of intersections */ + /* with the edges of the triangles that share that vertex) */ + reevaluate_isectns(s, lp+i, k-i, 1e-5, vb, vv); + +#ifdef INTERSECT_DEBUG + { + int ii; + printf("After re-evaluating intersections\n"); + for (ii = i; ii < k; ii++) + printf("Isect %d: pv %f, dir %d, edge %d, tri %d\n",ii,lp[ii].pv,lp[ii].dir,lp[ii].edge,lp[ii].tri->n); + } +#endif + } + +//printf("~1 nin %d, nout %d, npin %d, npout %d\n",nin,nout,npin,npout); + /* Create a zero length segment */ + if ((k-i) > 1 + && ((nin >= 1 && nout >= 1) /* both are definite */ + || (nin == 0 && nout == 0 && npin >= 1 && npout >= 1) /* Both glancing */ + || (nin == 0 && nout == 0 && npin == 0 && npout == 0)) /* No intersections now */ + ) { + if (pdir != 0) { /* Not correct for segment */ + i = k; +//printf("~1 neither or both or uncertain\n"); + continue; /* Discard them all */ + } +//printf("~1 creating zero length segment\n"); + /* Hmm. For reasonable triangles we should really */ + /* grab in/out from original evaluation... */ + memmove(&lp[j], &lp[i], sizeof(gispnt)); + lp[j].dir = 1; + lp[j].edge = 1; + j++; + memmove(&lp[j], &lp[i+1], sizeof(gispnt)); + lp[j].dir = 0; + lp[j].edge = 1; + j++; + i = k; + continue; + } + + /* We expect one conclusion */ + if (nin >= 1) + i = inx; + else if (nout >= 1) + i = outx; + else if (npin >= 1) + i = pinx; + else /* npout >= 1 */ + i = poutx; +//printf("~1 using %d\n",i); + + if ((lp[i].dir ^ pdir) == 0) { /* Not opposite to previous */ +//printf("~1 not opposite, discard it\n"); + /* This shouldn't happen. */ + i = k; + continue; /* Discard it */ + } +//printf("~1 save %d\n",i); + /* Accept this intersection */ + memmove(&lp[j], &lp[i], sizeof(gispnt)); + pdir = lp[j].dir; + j++; + i = k; + } + if (j & 1) /* Hmm. We ended up odd. This shouldn't happen. */ + j--; + rv = j; + +#ifdef INTERSECT_DEBUG + if (rv == 0) + printf("No intersections left\n"); + else { + printf("After sanitizing %d\n",rv); + for (i = 0; i < rv; i++) + printf("Isect %d: pv %f, dir %d, edge %d, tri %d\n",i,lp[i].pv,lp[i].dir,lp[i].edge,lp[i].tri->n); + } +#endif + return rv; +} + +#ifdef INTERSECT_DEBUG +#undef ISDBG +#endif /* INTERSECT_DEBUG */ + +/* ===================================================== */ +/* Write to a VRML .wrl file */ +/* Return non-zero on error */ +static int write_vrml( +gamut *s, +char *filename, +int doaxes, /* Non-zero if axes are to be written */ +int docusps /* Non-zero if cusp points are to be marked */ +) { + return write_trans_vrml(s, filename, doaxes, docusps, NULL, NULL); +} + +/* Write to a VRML .wrl file */ +/* Return non-zero on error */ +static int write_trans_vrml( +gamut *s, +char *filename, +int doaxes, /* Non-zero if axes are to be written */ +int docusps, /* Non-zero if cusp points are to be marked */ +void (*transform)(void *cntx, double out[3], double in[3]), /* Optional transformation callback */ +void *cntx +) { + int i; + gtri *tp; /* Triangle pointer */ + FILE *wrl; + struct { + double x, y, z; + double wx, wy, wz; + double r, g, b; + } axes[5] = { + { 0 - s->cent[1], 0 - s->cent[2], 50 - s->cent[0], 2, 2, 100, .7, .7, .7 }, + /* L axis */ + { 50 - s->cent[1], 0 - s->cent[2], 0 - s->cent[0], 100, 2, 2, 1, 0, 0 }, + /* +a (red) axis */ + { 0 - s->cent[1], -50 - s->cent[2], 0 - s->cent[0], 2, 100, 2, 0, 0, 1 }, + /* -b (blue) axis */ + { -50 - s->cent[1], 0 - s->cent[2], 0 - s->cent[0], 100, 2, 2, 0, 1, 0 }, + /* -a (green) axis */ + { 0 - s->cent[1], 50 - s->cent[2], 0 - s->cent[0], 2, 100, 2, 1, 1, 0 }, + /* +b (yellow) axis */ + }; + + /* Define the labels */ + struct { + double x, y, z; + double size; + char *string; + double r, g, b; + } labels[6] = { + { -2 - s->cent[1], 2 - s->cent[2], - s->cent[0] + 100 + 10, 10, "+L*", .7, .7, .7 }, + /* Top of L axis */ + { -2 - s->cent[1], 2 - s->cent[2], - s->cent[0] - 10, 10, "0", .7, .7, .7 }, + /* Bottom of L axis */ + { 100 + 5 - s->cent[1], -3 - s->cent[2], 0 - s->cent[0], 10, "+a*", 1, 0, 0 }, + /* +a (red) axis */ + { -5 - s->cent[1], -100 - 10 - s->cent[2], 0 - s->cent[0], 10, "-b*", 0, 0, 1 }, + /* -b (blue) axis */ + { -100 - 15 - s->cent[1], -3 - s->cent[2], 0 - s->cent[0], 10, "-a*", 0, 0, 1 }, + /* -a (green) axis */ + { -5 - s->cent[1], 100 + 5 - s->cent[2], 0 - s->cent[0], 10, "+b*", 1, 1, 0 }, + /* +b (yellow) axis */ + }; + + if IS_LIST_EMPTY(s->tris) + triangulate(s); + + if ((wrl = fopen(filename,"w")) == NULL) { + fprintf(stderr,"Error opening output file '%s'\n",filename); + return 2; + } + + /* Spit out a VRML 2 Object surface of gamut */ + + fprintf(wrl,"#VRML V2.0 utf8\n"); + fprintf(wrl,"\n"); + fprintf(wrl,"# Created by the Argyll CMS\n"); + fprintf(wrl,"Transform {\n"); + fprintf(wrl,"children [\n"); + fprintf(wrl," NavigationInfo {\n"); + fprintf(wrl," type \"EXAMINE\" # It's an object we examine\n"); + fprintf(wrl," } # We'll add our own light\n"); + fprintf(wrl,"\n"); + fprintf(wrl," DirectionalLight {\n"); + fprintf(wrl," intensity 0.2\n"); + fprintf(wrl," ambientIntensity 0.1\n"); + fprintf(wrl," direction -1 -1 -1\n"); + fprintf(wrl," }\n"); + fprintf(wrl," DirectionalLight {\n"); + fprintf(wrl," intensity 0.6\n"); + fprintf(wrl," ambientIntensity 0.2\n"); + fprintf(wrl," direction 1 1 1\n"); + fprintf(wrl," }\n"); + fprintf(wrl,"\n"); + fprintf(wrl," Viewpoint {\n"); + fprintf(wrl," position 0 0 250 # Position we view from\n"); + fprintf(wrl," }\n"); + fprintf(wrl,"\n"); + if (doaxes != 0) { + fprintf(wrl,"# Lab axes as boxes:\n"); + for (i = 0; i < 5; i++) { + fprintf(wrl,"Transform { translation %f %f %f\n", axes[i].x, axes[i].y, axes[i].z); + fprintf(wrl,"\tchildren [\n"); + fprintf(wrl,"\t\tShape {\n"); + fprintf(wrl,"\t\t\tgeometry Box { size %f %f %f }\n", + axes[i].wx, axes[i].wy, axes[i].wz); + fprintf(wrl,"\t\t\tappearance Appearance { material Material "); + fprintf(wrl,"{ diffuseColor %f %f %f} }\n", axes[i].r, axes[i].g, axes[i].b); + fprintf(wrl,"\t\t}\n"); + fprintf(wrl,"\t]\n"); + fprintf(wrl,"}\n"); + } + fprintf(wrl,"# Axes identification:\n"); + for (i = 0; i < 6; i++) { + fprintf(wrl,"Transform { translation %f %f %f\n", labels[i].x, labels[i].y, labels[i].z); + fprintf(wrl,"\tchildren [\n"); + fprintf(wrl,"\t\tShape {\n"); + fprintf(wrl,"\t\t\tgeometry Text { string [\"%s\"]\n",labels[i].string); + fprintf(wrl,"\t\t\t\tfontStyle FontStyle { family \"SANS\" style \"BOLD\" size %f }\n", + labels[i].size); + fprintf(wrl,"\t\t\t\t}\n"); + fprintf(wrl,"\t\t\tappearance Appearance { material Material "); + fprintf(wrl,"{ diffuseColor %f %f %f} }\n", labels[i].r, labels[i].g, labels[i].b); + fprintf(wrl,"\t\t}\n"); + fprintf(wrl,"\t]\n"); + fprintf(wrl,"}\n"); + } + fprintf(wrl,"\n"); + } + fprintf(wrl," Transform {\n"); + fprintf(wrl," translation 0 0 0\n"); + fprintf(wrl," children [\n"); + fprintf(wrl," Shape { \n"); + fprintf(wrl," geometry IndexedFaceSet {\n"); + fprintf(wrl," ccw FALSE\n"); + fprintf(wrl," convex TRUE\n"); + fprintf(wrl,"\n"); + fprintf(wrl," coord Coordinate { \n"); + fprintf(wrl," point [ # Verticy coordinates\n"); + + /* Spit out the point values, in order. */ + /* Note that a->x, b->y, L->z */ + for (i = 0; i < s->nv; i++) { + double out[3]; + +#ifdef SHOW_BUCKETS /* Show vertex buckets as surface */ + if (!(s->verts[i]->f & GVERT_SET)) +#else + if (!(s->verts[i]->f & GVERT_TRI)) +#endif + continue; + +#ifdef SHOW_BUCKETS /* Show vertex buckets as surface */ + { + double cc[3], rr[3]; +# ifdef SHOW_SPHERE /* Show surface on sphere */ + rr[0] = 50.0; /* Sphere radius */ +# else + rr[0] = s->verts[i]->r[0], +# endif /* SHOW_SPHERE */ + + rr[1] = s->verts[i]->hc - 0.5 * s->verts[i]->w; + rr[2] = s->verts[i]->vc - 0.5 * s->verts[i]->h; + gamut_radial2rect(s, cc, rr); + fprintf(wrl,"%f %f %f,\n",cc[1], cc[2], cc[0]); + + rr[1] = s->verts[i]->hc - 0.5 * s->verts[i]->w; + rr[2] = s->verts[i]->vc + 0.5 * s->verts[i]->h; + gamut_radial2rect(s, cc, rr); + fprintf(wrl,"%f %f %f,\n",cc[1], cc[2], cc[0]); + + rr[1] = s->verts[i]->hc + 0.5 * s->verts[i]->w; + rr[2] = s->verts[i]->vc + 0.5 * s->verts[i]->h; + gamut_radial2rect(s, cc, rr); + fprintf(wrl,"%f %f %f,\n",cc[1], cc[2], cc[0]); + + rr[1] = s->verts[i]->hc + 0.5 * s->verts[i]->w; + rr[2] = s->verts[i]->vc - 0.5 * s->verts[i]->h; + gamut_radial2rect(s, cc, rr); + fprintf(wrl,"%f %f %f,\n",cc[1], cc[2], cc[0]); + } + +#else /* Show point data */ + +# ifdef SHOW_SPHERE /* Show surface on sphere */ + fprintf(wrl,"%f %f %f,\n",s->verts[i]->sp[1], s->verts[i]->sp[2], + s->verts[i]->sp[0]); +# else +# ifdef SHOW_HULL_PNTS + fprintf(wrl,"%f %f %f,\n",s->verts[i]->ch[1], s->verts[i]->ch[2], + s->verts[i]->ch[0]); +# else + /* Show normal gamut surface */ + out[0] = s->verts[i]->p[0]; + out[1] = s->verts[i]->p[1]; + out[2] = s->verts[i]->p[2]; + + if (transform) + transform(cntx, out, out); /* Do transform */ + + fprintf(wrl,"%f %f %f,\n",out[1]-s->cent[1], out[2]-s->cent[2], out[0]-s->cent[0]); + +# endif /* SHOW_HULL_PNTS */ +# endif /* SHOW_SPHERE */ + +#endif /* SHOW_BUCKETS */ + + } + fprintf(wrl," ]\n"); + fprintf(wrl," }\n"); + fprintf(wrl,"\n"); + fprintf(wrl," coordIndex [ # Indexes of poligon Verticies \n"); + +#ifdef SHOW_BUCKETS /* Show vertex buckets as surface */ + for (i = 0; i < s->nv; i++) { + int j = s->verts[i]->sn; + if (!(s->verts[i]->f & GVERT_SET)) + continue; + fprintf(wrl,"%d, %d, %d, %d, -1\n", j * 4, j * 4 + 1, j * 4 + 2, j * 4 + 3); + } +#else /* Show gamut triangular surface */ + tp = s->tris; + FOR_ALL_ITEMS(gtri, tp) { + fprintf(wrl,"%d, %d, %d, -1\n", tp->v[0]->tn, tp->v[1]->tn, tp->v[2]->tn); + } END_FOR_ALL_ITEMS(tp); +#endif /* SHOW_BUCKETS */ + + fprintf(wrl," ]\n"); + fprintf(wrl,"\n"); + fprintf(wrl," colorPerVertex TRUE\n"); + fprintf(wrl," color Color {\n"); + fprintf(wrl," color [ # RGB colors of each vertex\n"); + + /* Spit out the colors for each vertex */ + for (i = 0; i < s->nv; i++) { + double rgb[3]; +#ifdef SHOW_BUCKETS /* Show vertex buckets as surface */ + if (!(s->verts[i]->f & GVERT_SET)) +#else + if (!(s->verts[i]->f & GVERT_TRI)) +#endif + continue; + +#ifdef COLORED_VRML + gamut_Lab2RGB(rgb, s->verts[i]->p); +#else + rgb[0] = rgb[1] = rgb[2] = 1.0; +#endif + fprintf(wrl,"%f %f %f,\n", rgb[0], rgb[1], rgb[2]); +#ifdef SHOW_BUCKETS /* Show vertex buckets as surface */ + fprintf(wrl,"%f %f %f,\n", rgb[0], rgb[1], rgb[2]); + fprintf(wrl,"%f %f %f,\n", rgb[0], rgb[1], rgb[2]); + fprintf(wrl,"%f %f %f,\n", rgb[0], rgb[1], rgb[2]); +#endif /* SHOW_BUCKETS */ + } + fprintf(wrl," ] \n"); + fprintf(wrl," }\n"); + fprintf(wrl," }\n"); + fprintf(wrl," appearance Appearance { \n"); + fprintf(wrl," material Material {\n"); + fprintf(wrl," transparency 0.0\n"); + fprintf(wrl," ambientIntensity 0.3\n"); + fprintf(wrl," shininess 0.5\n"); + fprintf(wrl," }\n"); + fprintf(wrl," }\n"); + fprintf(wrl," } # end Shape\n"); + fprintf(wrl," ]\n"); + fprintf(wrl," }\n"); + + + if (s->gawbset && doaxes) { + + /* Show the gamut white and black points */ + fprintf(wrl,"\n"); + fprintf(wrl," Transform {\n"); + fprintf(wrl," translation %f %f %f\n",s->ga_wp[1]-s->cent[1], s->ga_wp[2]-s->cent[2], s->ga_wp[0]-s->cent[0]); + fprintf(wrl," children [\n"); + fprintf(wrl," Shape { \n"); + fprintf(wrl," geometry Sphere { radius 2.0 }\n"); + fprintf(wrl," appearance Appearance { material Material { diffuseColor 0.9 0.9 0.9 } }\n"); + fprintf(wrl," } \n"); + fprintf(wrl," ]\n"); + fprintf(wrl," }\n"); + fprintf(wrl,"\n"); + fprintf(wrl," Transform {\n"); + fprintf(wrl," translation %f %f %f\n",s->ga_bp[1]-s->cent[1], s->ga_bp[2]-s->cent[2], s->ga_bp[0]-s->cent[0]); + fprintf(wrl," children [\n"); + fprintf(wrl," Shape { \n"); + fprintf(wrl," geometry Sphere { radius 2.0 }\n"); + fprintf(wrl," appearance Appearance { material Material { diffuseColor 0.9 0.9 0.9 } }\n"); + fprintf(wrl," } \n"); + fprintf(wrl," ]\n"); + fprintf(wrl," }\n"); + } + + if (docusps && s->cu_inited != 0) { + double ccolors[6][3] = { + { 1.0, 0.1, 0.1 }, /* Red */ + { 1.0, 1.0, 0.1 }, /* Yellow */ + { 0.1, 1.0, 0.1 }, /* Green */ + { 0.1, 1.0, 1.0 }, /* Cyan */ + { 0.1, 0.1, 1.0 }, /* Blue */ + { 1.0, 0.1, 1.0 } /* Magenta */ + }; + + for (i = 0; i < 6; i++) { + fprintf(wrl,"\n"); + fprintf(wrl," Transform {\n"); + fprintf(wrl," translation %f %f %f\n",s->cusps[i][1]-s->cent[1], s->cusps[i][2]-s->cent[2], s->cusps[i][0]-s->cent[0]); + fprintf(wrl," children [\n"); + fprintf(wrl," Shape { \n"); + fprintf(wrl," geometry Sphere { radius 2.0 }\n"); + fprintf(wrl," appearance Appearance { material Material { diffuseColor %f %f %f } }\n", ccolors[i][0],ccolors[i][1],ccolors[i][2]); + fprintf(wrl," } \n"); + fprintf(wrl," ]\n"); + fprintf(wrl," }\n"); + } + } + +#ifdef TEST_LOOKUP + { + int i, j; + double in[3], out[3]; + + fprintf(wrl,"\n"); + fprintf(wrl,"Shape {\n"); + fprintf(wrl," geometry PointSet { \n"); + fprintf(wrl," coord Coordinate { \n"); + fprintf(wrl," point [\n"); + + for (i = 0; i < 10; i++) { + double ss; + /* Create random vector relative to center, absolute */ + in[0] = (rand() / (double)RAND_MAX) - 0.5 + s->cent[0]; + in[1] = (rand() / (double)RAND_MAX) - 0.5 + s->cent[1]; + in[2] = (rand() / (double)RAND_MAX) - 0.5 + s->cent[2]; + + s->radial(s, out, in); /* Lookup point on gamut surface */ + + out[0] = (out[0] - s->cent[0]) * 1.01 + s->cent[0]; + out[1] = (out[1] - s->cent[1]) * 1.01 + s->cent[1]; + out[2] = (out[2] - s->cent[2]) * 1.01 + s->cent[2]; + fprintf(wrl,"%f %f %f,\n",out[1]-s->cent[1], out[2]-s->cent[2], out[0]-s->cent[0]); + } + fprintf(wrl," ]\n"); + fprintf(wrl," }\n"); + fprintf(wrl," }\n"); + fprintf(wrl,"} # end shape\n"); + + } +#endif /* TEST_LOOKUP */ + +#ifdef TEST_NEAREST + { +#define NTPTS 500 + int i, j; + double in[3], out[3]; + + fprintf(wrl,"\n"); + fprintf(wrl,"Shape {\n"); + fprintf(wrl," geometry IndexedLineSet { \n"); + fprintf(wrl," coord Coordinate { \n"); + fprintf(wrl," point [\n"); + + for (i = 0; i < NTPTS; i++) { + double ss; + /* Create random vector relative to center */ + in[0] = (rand() / (double)RAND_MAX) - 0.5; + in[1] = (rand() / (double)RAND_MAX) - 0.5; + in[2] = (rand() / (double)RAND_MAX) - 0.5; + +#ifndef NEVER /* Make points just above surface */ + in[0] += s->cent[0]; /* Make absolute */ + in[1] += s->cent[1]; + in[2] += s->cent[2]; + s->radial(s, in, in); /* Lookup point on gamut surface */ + in[0] = (in[0] - s->cent[0]) * 1.20 + s->cent[0]; /* Extend by 10% */ + in[1] = (in[1] - s->cent[1]) * 1.20 + s->cent[1]; + in[2] = (in[2] - s->cent[2]) * 1.20 + s->cent[2]; +#else + /* Make distance 150 */ + ss = sqrt(in[0] * in[0] + in[1] * in[1] + in[2] * in[2]); + in[0] = 60.0/ss * in[0] + s->cent[0]; + in[1] = 60.0/ss * in[1] + s->cent[1]; + in[2] = 60.0/ss * in[2] + s->cent[2]; +#endif + +// s->radial(s, out, in); /* Lookup point on gamut surface */ + + s->nearest(s, out, in); /* Nearest point on gamut surface */ + + fprintf(wrl,"%f %f %f,\n",in[1]-s->cent[1], in[2]-s->cent[2], in[0]-s->cent[0]); + fprintf(wrl,"%f %f %f,\n",out[1]-s->cent[1], out[2]-s->cent[2], out[0]-s->cent[0]); + } + + fprintf(wrl," ]\n"); + fprintf(wrl," }\n"); + fprintf(wrl," coordIndex [\n"); + + for (i = 0; i < NTPTS; i++) { + fprintf(wrl,"%d, %d, -1,\n", i * 2, i * 2 + 1); + } + fprintf(wrl," ]\n"); + fprintf(wrl," }\n"); + fprintf(wrl,"} # end shape\n"); + + } +#endif /* TEST_NEAREST */ + + fprintf(wrl,"\n"); + fprintf(wrl," ] # end of children for world\n"); + fprintf(wrl,"}\n"); + + if (fclose(wrl) != 0) { + fprintf(stderr,"Error closing output file '%s'\n",filename); + return 2; + } + + return 0; +} + + +/* ----------------------------------- */ +/* Write to a CGATS .gam file */ +/* Return non-zero on error */ +static int write_gam( +gamut *s, +char *filename +) { + time_t clk = time(0); + struct tm *tsp = localtime(&clk); + char *atm = asctime(tsp); /* Ascii time */ + int i; + gtri *tp; /* Triangle pointer */ + cgats *gam; + char buf[100]; + + if IS_LIST_EMPTY(s->tris) + triangulate(s); + + gam = new_cgats(); /* Create a CGATS structure */ + gam->add_other(gam, "GAMUT"); + + gam->add_table(gam, tt_other, 0); /* Start the first table as type "GAMUT" */ + + gam->add_kword(gam, 0, "DESCRIPTOR", "Argyll Gamut surface poligon data", NULL); + gam->add_kword(gam, 0, "ORIGINATOR", "Argyll CMS gamut library", NULL); + atm[strlen(atm)-1] = '\000'; /* Remove \n from end */ + gam->add_kword(gam, 0, "CREATED",atm, NULL); + +#ifdef NEVER + /* would be nice to add extra info like description, source (ie icc filename) etc. */ + gam->add_kword(gam, 0, "DEVICE_CLASS","INPUT", NULL); /* What sort of device this is */ +#endif + if (s->isJab) + gam->add_kword(gam, 0, "COLOR_REP","JAB", NULL); + else + gam->add_kword(gam, 0, "COLOR_REP","LAB", NULL); + + if (s->isRast) + gam->add_kword(gam, 0, "SURF_TYPE","RASTER", NULL); + + sprintf(buf,"%f %f %f", s->cent[0], s->cent[1], s->cent[2]); + gam->add_kword(gam, 0, "GAMUT_CENTER",buf, NULL); + + /* If the white and black points are known, put them in the file */ + if (s->cswbset) { + + compgawb(s); /* make sure we have gamut white/black available */ + + sprintf(buf,"%f %f %f", s->cs_wp[0], s->cs_wp[1], s->cs_wp[2]); + gam->add_kword(gam, 0, "CSPACE_WHITE",buf, NULL); + + sprintf(buf,"%f %f %f", s->ga_wp[0], s->ga_wp[1], s->ga_wp[2]); + gam->add_kword(gam, 0, "GAMUT_WHITE",buf, NULL); + + sprintf(buf,"%f %f %f", s->cs_bp[0], s->cs_bp[1], s->cs_bp[2]); + gam->add_kword(gam, 0, "CSPACE_BLACK",buf, NULL); + + sprintf(buf,"%f %f %f", s->ga_bp[0], s->ga_bp[1], s->ga_bp[2]); + gam->add_kword(gam, 0, "GAMUT_BLACK",buf, NULL); + } + + /* If cusp values are known, put them in the file */ + if (s->cu_inited != 0) { + char buf1[50], buf2[100]; + char *cnames[6] = { "RED", "YELLOW", "GREEN", "CYAN", "BLUE", "MAGENTA" }; + + for (i = 0; i < 6; i++) { + sprintf(buf1,"CUSP_%s", cnames[i]); + sprintf(buf2,"%f %f %f", s->cusps[i][0], s->cusps[i][1], s->cusps[i][2]); + gam->add_kword(gam, 0, buf1, buf2, NULL); + } + } + + gam->add_kword(gam, 0, NULL, NULL, "First come the triangle verticy location"); + + gam->add_field(gam, 0, "VERTEX_NO", i_t); + gam->add_field(gam, 0, "LAB_L", r_t); + gam->add_field(gam, 0, "LAB_A", r_t); + gam->add_field(gam, 0, "LAB_B", r_t); + + /* Spit out the vertex values, in order. */ + for (i = 0; i < s->nv; i++) { + if (!(s->verts[i]->f & GVERT_TRI)) + continue; + gam->add_set(gam, 0, s->verts[i]->tn, + s->verts[i]->p[0], s->verts[i]->p[1], s->verts[i]->p[2]); + } + + gam->add_table(gam, tt_other, 0); /* Start the second table */ + gam->set_table_flags(gam, 1, 1, 1, 0); /* Suppress id & kwords */ + gam->add_kword(gam, 1, NULL, NULL, "And then come the triangles"); + + gam->add_field(gam, 1, "VERTEX_0", i_t); + gam->add_field(gam, 1, "VERTEX_1", i_t); + gam->add_field(gam, 1, "VERTEX_2", i_t); + + tp = s->tris; + FOR_ALL_ITEMS(gtri, tp) { + gam->add_set(gam, 1, tp->v[0]->tn, tp->v[1]->tn, tp->v[2]->tn); + } END_FOR_ALL_ITEMS(tp); + + + if (gam->write_name(gam, filename)) { + fprintf(stderr,"Error writing to file '%s' : '%s'\n",filename, gam->err); + return 2; + } + + gam->del(gam); /* Clean up */ + return 0; +} + +/* ----------------------------------- */ +/* Read from a CGATS .gam file */ +/* Return non-zero on error */ +static int read_gam( +gamut *s, +char *filename +) { + int i; + cgats *gam; + gtri *tp; + int nverts; + int ntris; + int Lf, af, bf; /* Fields holding L, a & b data */ + int v0f, v1f, v2f; /* Fields holding verticies 0, 1 & 2 */ + int cw, cb; /* Colorspace white, black keyword indexes */ + int gw, gb; /* Gamut white, black keyword indexes */ + + if (s->tris != NULL || s->read_inited || s->lu_inited || s->ne_inited) { + fprintf(stderr,"Can't add read into gamut after it is initialised!\n"); + return 1; + } + + gam = new_cgats(); /* Create a CGATS structure */ + + gam->add_other(gam, "GAMUT"); /* Setup to cope with a gamut file */ + + if (gam->read_name(gam, filename)) { + fprintf(stderr,"Input file '%s' error : %s",filename, gam->err); + return 1; + } + + if (gam->t[0].tt != tt_other || gam->t[0].oi != 0) { + fprintf(stderr,"Input file isn't a GAMUT format file"); + return 1; + } + if (gam->ntables != 2) { + fprintf(stderr,"Input file doesn't contain exactly two tables"); + return 1; + } + + /* Figure the basic colorspace information */ + s->isJab = 0; + if ((cw = gam->find_kword(gam, 0, "COLOR_REP")) >= 0) { + if (strcmp(gam->t[0].kdata[cw], "JAB") == 0) + s->isJab = 1; + } + + /* Figure the surface type */ + s->isRast = 0; + if ((cw = gam->find_kword(gam, 0, "SURF_TYPE")) >= 0) { + if (strcmp(gam->t[0].kdata[cw], "RASTER") == 0) + s->isRast = 1; + } + if (s->isRast) { + s->logpow = RAST_LOG_POW; /* Wrap the surface more closely */ + s->no2pass = 1; /* Only do one pass */ + } else { + s->logpow = NORM_LOG_POW; /* Convex hull compression power */ + s->no2pass = 0; /* Do two passes */ + } + + /* If we can find the the colorspace white and black points, add them to the gamut */ + cw = gam->find_kword(gam, 0, "CSPACE_WHITE"); + cb = gam->find_kword(gam, 0, "CSPACE_BLACK"); + if (cw >= 0 && cb >= 0) { + int ok = 1; + if (sscanf(gam->t[0].kdata[cw], "%lf %lf %lf", + &s->cs_wp[0], &s->cs_wp[1], &s->cs_wp[2]) != 3) { + ok = 0; + } + + if (sscanf(gam->t[0].kdata[cb], "%lf %lf %lf", + &s->cs_bp[0], &s->cs_bp[1], &s->cs_bp[2]) != 3) { + ok = 0; + } + + if (ok) { + s->cswbset = 1; + } + } + + /* If we can find the the gamut white and black points, add them to the gamut */ + gw = gam->find_kword(gam, 0, "GAMUT_WHITE"); + gb = gam->find_kword(gam, 0, "GAMUT_BLACK"); + if (gw >= 0 && gb >= 0) { + int ok = 1; + if (sscanf(gam->t[0].kdata[gw], "%lf %lf %lf", + &s->ga_wp[0], &s->ga_wp[1], &s->ga_wp[2]) != 3) { + ok = 0; + } + + if (sscanf(gam->t[0].kdata[gb], "%lf %lf %lf", + &s->ga_bp[0], &s->ga_bp[1], &s->ga_bp[2]) != 3) { + ok = 0; + } + + if (ok) { + s->gawbset = 1; + } + } + + /* See if there are cusp values */ + { + int kk; + char buf1[50]; + char *cnames[6] = { "RED", "YELLOW", "GREEN", "CYAN", "BLUE", "MAGENTA" }; + + for (i = 0; i < 6; i++) { + sprintf(buf1,"CUSP_%s", cnames[i]); + if ((kk = gam->find_kword(gam, 0, buf1)) < 0) + break; + + if (sscanf(gam->t[0].kdata[kk], "%lf %lf %lf", + &s->cusps[i][0], &s->cusps[i][1], &s->cusps[i][2]) != 3) { + break; + } + } + if (i >= 6) + s->cu_inited = 1; + } + + + if ((nverts = gam->t[0].nsets) <= 0) { + fprintf(stderr,"No verticies"); + return 1; + } + if ((ntris = gam->t[1].nsets) <= 0) { + fprintf(stderr,"No triangles"); + return 1; + } + + /* Get ready to read the verticy data */ + if ((Lf = gam->find_field(gam, 0, "LAB_L")) < 0) { + fprintf(stderr,"Input file doesn't contain field LAB_L"); + return 1; + } + if (gam->t[0].ftype[Lf] != r_t) { + fprintf(stderr,"Field LAB_L is wrong type"); + return 1; + } + if ((af = gam->find_field(gam, 0, "LAB_A")) < 0) { + fprintf(stderr,"Input file doesn't contain field LAB_A"); + return 1; + } + if (gam->t[0].ftype[af] != r_t) { + fprintf(stderr,"Field LAB_A is wrong type"); + return 1; + } + if ((bf = gam->find_field(gam, 0, "LAB_B")) < 0) { + fprintf(stderr,"Input file doesn't contain field LAB_B"); + return 1; + } + if (gam->t[0].ftype[bf] != r_t) { + fprintf(stderr,"Field LAB_B is wrong type"); + return 1; + } + + /* Allocate an array to point at the verts */ + if ((s->verts = (gvert **)malloc(nverts * sizeof(gvert *))) == NULL) { + fprintf(stderr,"gamut: malloc failed on gvert pointer\n"); + return 2; + } + s->nv = s->na = nverts; + + for (i = 0; i < nverts; i++) { + gvert *v; + + /* Allocate and fill in each verticies basic information */ + if ((v = (gvert *)calloc(1, sizeof(gvert))) == NULL) { + fprintf(stderr,"gamut: malloc failed on gvert object\n"); + return 2; + } + s->verts[i] = v; + v->tag = 1; + v->tn = v->n = i; + v->f = GVERT_SET | GVERT_TRI; /* Will be part of the triangulation */ + + v->p[0] = *((double *)gam->t[0].fdata[i][Lf]); + v->p[1] = *((double *)gam->t[0].fdata[i][af]); + v->p[2] = *((double *)gam->t[0].fdata[i][bf]); + + gamut_rect2radial(s, v->r, v->p); + } + s->ntv = i; + + /* Compute the other vertex values */ + compute_vertex_coords(s); + + /* Get ready to read the triangle data */ + if ((v0f = gam->find_field(gam, 1, "VERTEX_0")) < 0) { + fprintf(stderr,"Input file doesn't contain field VERTEX_0"); + return 1; + } + if (gam->t[1].ftype[v0f] != i_t) { + fprintf(stderr,"Field VERTEX_0 is wrong type"); + return 1; + } + if ((v1f = gam->find_field(gam, 1, "VERTEX_1")) < 0) { + fprintf(stderr,"Input file doesn't contain field VERTEX_1"); + return 1; + } + if (gam->t[1].ftype[v1f] != i_t) { + fprintf(stderr,"Field VERTEX_1 is wrong type"); + return 1; + } + if ((v2f = gam->find_field(gam, 1, "VERTEX_2")) < 0) { + fprintf(stderr,"Input file doesn't contain field VERTEX_2"); + return 1; + } + if (gam->t[1].ftype[v2f] != i_t) { + fprintf(stderr,"Field VERTEX_2 is wrong type"); + return 1; + } + + /* Create all the triangles */ + for (i = 0; i < ntris; i++) { + gtri *t; + int v0, v1, v2; + + t = new_gtri(); + ADD_ITEM_TO_BOT(s->tris, t); /* Append to triangulation list */ + + v0 = *((int *)gam->t[1].fdata[i][v0f]); + v1 = *((int *)gam->t[1].fdata[i][v1f]); + v2 = *((int *)gam->t[1].fdata[i][v2f]); + + t->v[0] = s->verts[v0]; + t->v[1] = s->verts[v1]; + t->v[2] = s->verts[v2]; + + comptriattr(s, t); /* Compute triangle attributes */ + } + + /* Connect edge information */ + tp = s->tris; + FOR_ALL_ITEMS(gtri, tp) { + int en; + + for (en = 0; en < 3; en++) { /* For each edge */ + gedge *e; + gvert *v0, *v1; /* The two verticies of the edge */ + gtri *tp2; /* The other triangle */ + int em; /* The other edge */ + gvert *w0, *w1; /* The other verticies */ + + v0 = tp->v[en]; + v1 = tp->v[en < 2 ? en+1 : 0]; + + if (v0->n > v1->n) + continue; /* Skip every other edge */ + + /* Find the corresponding edge of the other triangle */ + w0 = w1 = NULL; + tp2 = s->tris; + FOR_ALL_ITEMS(gtri, tp2) { + for (em = 0; em < 3; em++) { /* For each edge */ + w0 = tp2->v[em]; + w1 = tp2->v[em < 2 ? em+1 : 0]; + if (v0 == w1 && v1 == w0) /* Found it */ + break; + } + if (em < 3) + break; + } END_FOR_ALL_ITEMS(tp2); + if (w0 == NULL) { + /* Should clean up ? */ + fprintf(stderr,".gam file triangle data is not consistent\n"); + return 1; + } + + if (tp->e[en] != NULL + || tp2->e[em] != NULL) { + fprintf(stderr,".gam file triangle data is not consistent\n"); + fprintf(stderr,"tp1->e[%d] = 0x%p, tp2->e[%d]= 0x%p\n",en, + (void *)tp->e[en],em,(void *)tp2->e[em]); + return 1; + } + + /* Creat the edge structure */ + e = new_gedge(); + ADD_ITEM_TO_BOT(s->edges, e); /* Append to edge list */ + tp->e[en] = e; /* This edge */ + tp->ei[en] = 0; /* 0th triangle in edge */ + e->t[0] = tp; /* 0th triangle is tp */ + e->ti[0] = en; /* 0th triangles en edge */ + tp2->e[em] = e; /* This edge */ + tp2->ei[em] = 1; /* 1st triangle in edge */ + e->t[1] = tp2; /* 1st triangle is tp2 */ + e->ti[1] = em; /* 1st triangles em edge */ + e->v[0] = v0; /* The two verticies */ + e->v[1] = v1; + } + } END_FOR_ALL_ITEMS(tp); + + gam->del(gam); /* Clean up */ + + s->read_inited = 1; /* It's now valid */ + +#ifdef ASSERTS + check_triangulation(s, 1); /* Check out our work */ +#endif + + return 0; +} + +/* ===================================================== */ +/* ===================================================== */ + +/* Convert from rectangular to radial coordinates */ +void +gamut_rect2radial( +gamut *s, +double out[3], /* Radius, longitude, lattitude out */ +double in[3] /* Lab in */ +) { + double L, a, b; /* Lab values */ + double R, g, t; /* Radial value */ + double c; /* Chromatic length */ + + + L = in[0] - s->cent[0]; /* Offset value */ + a = in[1] - s->cent[1]; + b = in[2] - s->cent[2]; + c = a * a + b * b; + R = c + L * L; + c = sqrt(c); /* Saturation */ + R = sqrt(R); /* Vector length */ + + if (R < 1e-6) { /* Hmm, a point at the center */ + g = t = 0.0; + } else { + + /* Figure out the longitude, -pi to +pi */ + if (c < 1e-6) { + g = 0.0; + } else { + g = asin(b/c); + if (a < 0.0) { + if (b >= 0.0) + g = M_PI - g; + else + g = -g - M_PI; + } + } + + /* Figure out the lattitude, -pi/2 to +pi/2 */ + t = asin(L/R); + } + out[0] = R; + out[1] = g; + out[2] = t; +} + +/* Convert from radial to rectangular coordinates */ +void +gamut_radial2rect( +gamut *s, +double out[3], /* Lab out */ +double in[3] /* Radius, longitude, lattitude in */ +) { + double R, g, t; /* Radial value */ + double L, a, b; /* Lab values */ + double c; /* Chromatic length */ + + R = in[0]; + g = in[1]; + t = in[2]; + + L = R * sin(t); + c = R * cos(t); + + a = c * cos(g); + b = c * sin(g); + + out[0] = L + s->cent[0]; + out[1] = a + s->cent[1]; + out[2] = b + s->cent[2]; +} + + +/* -------------------------------------------------- */ + +/* Convert a gamut Lab value to an RGB value for display purposes */ +void +gamut_Lab2RGB(double *out, double *in) { + double L = in[0], a = in[1], b = in[2]; + double x,y,z,fx,fy,fz; + double R, G, B; + + /* Scale so that black is visible */ + L = L * (100 - 40.0)/100.0 + 40.0; + + /* First convert to XYZ using D50 white point */ + if (L > 8.0) { + fy = (L + 16.0)/116.0; + y = pow(fy,3.0); + } else { + y = L/903.2963058; + fy = 7.787036979 * y + 16.0/116.0; + } + + fx = a/500.0 + fy; + if (fx > 24.0/116.0) + x = pow(fx,3.0); + else + x = (fx - 16.0/116.0)/7.787036979; + + fz = fy - b/200.0; + if (fz > 24.0/116.0) + z = pow(fz,3.0); + else + z = (fz - 16.0/116.0)/7.787036979; + + x *= 0.9642; /* Multiply by white point, D50 */ + y *= 1.0; + z *= 0.8249; + + /* Now convert to sRGB values */ + R = x * 3.2410 + y * -1.5374 + z * -0.4986; + G = x * -0.9692 + y * 1.8760 + z * 0.0416; + B = x * 0.0556 + y * -0.2040 + z * 1.0570; + + if (R < 0.0) + R = 0.0; + else if (R > 1.0) + R = 1.0; + + if (G < 0.0) + G = 0.0; + else if (G > 1.0) + G = 1.0; + + if (B < 0.0) + B = 0.0; + else if (B > 1.0) + B = 1.0; + + R = pow(R, 1.0/2.2); + G = pow(G, 1.0/2.2); + B = pow(B, 1.0/2.2); + + out[0] = R; + out[1] = G; + out[2] = B; +} + + +/* -------------------------------------------------- */ + +#ifdef DEBUG_TRIANG + +/* Write a surface contrsuction diagnostic VRML .wrl file */ +static int write_diag_vrml( +gamut *s, +double vv[3], /* Vertex being added */ +int nh, /* Number of hit triangles */ +tidxs *hixs, /* verticy indexes of hit triangles */ +gtri *hl /* Edge hit list (may be NULL) */ +) { + char *filename; + int i, j; + gtri *tp; /* Triangle pointer */ + FILE *wrl; + + if (hl) + filename = "diag1.wrl"; /* Triangles hit */ + else + filename = "diag2.wrl"; /* Triangles formed */ + + if ((wrl = fopen(filename,"w")) == NULL) { + fprintf(stderr,"Error opening output file '%s'\n",filename); + return 2; + } + + /* Spit out a VRML 2 Object surface of gamut */ + + fprintf(wrl,"#VRML V2.0 utf8\n"); + fprintf(wrl,"\n"); + fprintf(wrl,"# Created by the Argyll CMS\n"); + fprintf(wrl,"Transform {\n"); + fprintf(wrl,"children [\n"); + fprintf(wrl," NavigationInfo {\n"); + fprintf(wrl," type \"EXAMINE\" # It's an object we examine\n"); + fprintf(wrl," } # We'll add our own light\n"); + fprintf(wrl,"\n"); + fprintf(wrl," DirectionalLight {\n"); + fprintf(wrl," direction 0 0 -1 # Light illuminating the scene\n"); + fprintf(wrl," direction 0 -1 0 # Light illuminating the scene\n"); + fprintf(wrl," }\n"); + fprintf(wrl,"\n"); + fprintf(wrl," Viewpoint {\n"); + fprintf(wrl," position 0 0 200 # Position we view from\n"); + fprintf(wrl," }\n"); + fprintf(wrl,"\n"); + + fprintf(wrl," Transform {\n"); + fprintf(wrl," translation 0 0 0\n"); + fprintf(wrl," children [\n"); + fprintf(wrl," Shape { \n"); + fprintf(wrl," geometry IndexedFaceSet {\n"); + fprintf(wrl," ccw FALSE\n"); + fprintf(wrl," convex TRUE\n"); + fprintf(wrl,"\n"); + fprintf(wrl," coord Coordinate { \n"); + fprintf(wrl," point [ # Verticy coordinates\n"); + + /* Spit out the vertex values, in order. */ + /* Note that a->x, b->y, L->z */ + for (i = 0; i < s->nv; i++) { + double out[3]; + + /* Show normal gamut surface */ + out[0] = s->verts[i]->ch[0]; + out[1] = s->verts[i]->ch[1]; + out[2] = s->verts[i]->ch[2]; + + fprintf(wrl,"%f %f %f,\n",out[1], out[2], out[0]); + } + fprintf(wrl," ]\n"); + fprintf(wrl," }\n"); + fprintf(wrl,"\n"); + fprintf(wrl," coordIndex [ # Indexes of poligon Verticies \n"); + + tp = s->tris; + FOR_ALL_ITEMS(gtri, tp) { + fprintf(wrl,"%d, %d, %d, -1\n", tp->v[0]->n, tp->v[1]->n, tp->v[2]->n); + } END_FOR_ALL_ITEMS(tp); + + for (i = 0; i < nh; i++) { + fprintf(wrl,"%d, %d, %d, -1\n", hixs[i].tix[0], hixs[i].tix[1], hixs[i].tix[2]); + } + + fprintf(wrl," ]\n"); + fprintf(wrl,"\n"); + fprintf(wrl," colorPerVertex FALSE\n"); + fprintf(wrl," color Color {\n"); + fprintf(wrl," color [ # RGB colors of each vertex\n"); + + /* Spit out the colors for each face */ + tp = s->tris; + FOR_ALL_ITEMS(gtri, tp) { + fprintf(wrl,"%f %f %f,\n", 0.7, 0.7, 0.7); + } END_FOR_ALL_ITEMS(tp); + for (i = 0; i < nh; i++) { + if (hixs[i].type == 0) + fprintf(wrl,"%f %f %f,\n", 0.4, 1.0, 0.4); /* Green for hit */ + else if (hixs[i].type == 1) + fprintf(wrl,"%f %f %f,\n", 0.4, 0.4, 1.0); /* Blue for extra */ + else + fprintf(wrl,"%f %f %f,\n", 0.8, 0.8, 0.2); /* Yellow for new */ + } + fprintf(wrl," ] \n"); + fprintf(wrl," }\n"); + fprintf(wrl," } # end IndexedFaceSet\n"); + + fprintf(wrl," appearance Appearance { \n"); + fprintf(wrl," material Material {\n"); + fprintf(wrl," transparency 0.0\n"); + fprintf(wrl," ambientIntensity 0.3\n"); + fprintf(wrl," shininess 0.5\n"); + fprintf(wrl," }\n"); + fprintf(wrl," }\n"); + fprintf(wrl," } # end Shape\n"); + fprintf(wrl," ]\n"); + fprintf(wrl," } # end of transform\n"); + + /* center of gamut */ + fprintf(wrl,"\n"); + fprintf(wrl," Transform {\n"); + fprintf(wrl," translation %f %f %f\n",0.0, 0.0, 0.0); + fprintf(wrl," children [\n"); + fprintf(wrl," Shape { \n"); + fprintf(wrl," geometry Sphere { radius 1.5 }\n"); + fprintf(wrl," appearance Appearance { material Material { diffuseColor %f %f %f } }\n", 1.0, 1.0, 0.0); + fprintf(wrl," } \n"); + fprintf(wrl," ]\n"); + fprintf(wrl," }\n"); + + /* vertex being added */ + fprintf(wrl,"\n"); + fprintf(wrl," Transform {\n"); + fprintf(wrl," translation %f %f %f\n",vv[1], vv[2], vv[0]); + fprintf(wrl," children [\n"); + fprintf(wrl," Shape { \n"); + fprintf(wrl," geometry Sphere { radius 1.5 }\n"); + fprintf(wrl," appearance Appearance { material Material { diffuseColor %f %f %f } }\n", 1.0, 0.0, 0.0); + fprintf(wrl," } \n"); + fprintf(wrl," ]\n"); + fprintf(wrl," }\n"); + + + /* Verticies for Polygon edges, marked by directional cones */ + if (hl != NULL) { + double base[3] = { 0.0, 0.0, 1.0 }; /* Default orientation of cone is b axis */ + tp = hl; + FOR_ALL_ITEMS(gtri, tp) { + double len; + double loc[3]; + double vec[3]; + double axis[3]; /* Axis to rotate around */ + double rot; /* In radians */ + +//printf("~1 edge vert %d to %d\n",tp->v[0]->n, tp->v[1]->n); +//printf("~1 edge %f %f %f to %f %f %f\n", +//tp->v[0]->ch[0], tp->v[0]->ch[1], tp->v[0]->ch[2], +//tp->v[1]->ch[0], tp->v[1]->ch[1], tp->v[1]->ch[2]); + + icmAdd3(loc, tp->v[1]->ch, tp->v[0]->ch); + icmScale3(loc, loc, 0.5); + icmSub3(vec, tp->v[1]->ch, tp->v[0]->ch); + len = icmNorm3(vec); + + if (len < 1.0) + len = 1.0; + + icmNormalize3(base, base, 1.0); + icmNormalize3(vec, vec, 1.0); + icmCross3(axis, base, vec); + rot = icmDot3(base, vec); +//printf("~1 Axis = %f %f %f\n",axis[0],axis[1],axis[2]); + if (icmNorm3sq(axis) < 1e-10) { /* 0 or 180 degrees */ + double base2[3]; + int mxi = 0; + base2[0] = vec[1]; /* Comute vector in a different direction */ + base2[1] = vec[2]; + base2[2] = vec[0]; + for (j = 1; j < 3; j++) { + if (fabs(base2[j]) > fabs(base2[mxi])) + mxi = j; + } + base2[mxi] = -base2[mxi]; + + icmCross3(axis, base2, vec); + if (icmNorm3sq(axis) < 1e-10) { /* 0 or 180 degrees */ + error("VRML rotate axis still too small"); + } + if (rot < 0.0) + rot = 3.1415926; + else + rot = 0.0; + } else { + rot = acos(rot); +//printf("~1 rotation %f\n",rot); + } + + fprintf(wrl,"\n"); + fprintf(wrl," Transform {\n"); + fprintf(wrl," rotation %f %f %f %f\n",axis[1], axis[2], axis[0], rot); + fprintf(wrl," translation %f %f %f\n",loc[1], loc[2], loc[0]); + fprintf(wrl," children [\n"); + fprintf(wrl," Shape { \n"); + fprintf(wrl," geometry Cone { bottomRadius 0.5 height %f }\n",len); + fprintf(wrl," appearance Appearance { material Material { diffuseColor 0.7 0.0 1.0 } }\n"); + fprintf(wrl," } \n"); + fprintf(wrl," ]\n"); + fprintf(wrl," }\n"); + } END_FOR_ALL_ITEMS(tp); + } + + fprintf(wrl,"\n"); + fprintf(wrl," ] # end of children for world\n"); + fprintf(wrl,"}\n"); + + if (fclose(wrl) != 0) { + fprintf(stderr,"Error closing output file '%s'\n",filename); + return 2; + } + + return 0; +} + +#endif /* DEBUG_TRIANG */ + + +#ifdef DEBUG_SPLIT_VRML + +/* Write a triangle split diagnostic VRML .wrl file */ +static int write_split_diag_vrml( +gamut *s, +gtri **list, /* Triangle list */ +int llen /* Number of triangles in the list */ +) { + char *filename; + int i, j; + FILE *wrl; + + filename = "diag3.wrl"; /* Triangles split */ + + if ((wrl = fopen(filename,"w")) == NULL) { + fprintf(stderr,"Error opening output file '%s'\n",filename); + return 2; + } + + /* Spit out a VRML 2 Object surface of gamut */ + + fprintf(wrl,"#VRML V2.0 utf8\n"); + fprintf(wrl,"\n"); + fprintf(wrl,"# Created by the Argyll CMS\n"); + fprintf(wrl,"Transform {\n"); + fprintf(wrl,"children [\n"); + fprintf(wrl," NavigationInfo {\n"); + fprintf(wrl," type \"EXAMINE\" # It's an object we examine\n"); + fprintf(wrl," } # We'll add our own light\n"); + fprintf(wrl,"\n"); + fprintf(wrl," DirectionalLight {\n"); + fprintf(wrl," ambientIntensity 0.3 # Ambient light illuminating the scene\n"); + fprintf(wrl," direction 0 0 -1 # Light illuminating the scene\n"); + fprintf(wrl," direction 0 -1 0 # Light illuminating the scene\n"); + fprintf(wrl," }\n"); + fprintf(wrl,"\n"); + fprintf(wrl," Viewpoint {\n"); + fprintf(wrl," position 0 0 5 # Position we view from\n"); + fprintf(wrl," }\n"); + fprintf(wrl,"\n"); + + fprintf(wrl," Transform {\n"); + fprintf(wrl," translation 0 0 0\n"); + fprintf(wrl," children [\n"); + fprintf(wrl," Shape { \n"); + fprintf(wrl," geometry IndexedFaceSet {\n"); + fprintf(wrl," ccw FALSE\n"); + fprintf(wrl," convex TRUE\n"); + fprintf(wrl," solid FALSE\n"); + fprintf(wrl,"\n"); + fprintf(wrl," coord Coordinate { \n"); + fprintf(wrl," point [ # Verticy coordinates\n"); + + /* Spit out the vertex values, in order. */ + for (i = 0; i < llen; i++) { + fprintf(wrl,"%f %f %f,\n",list[i]->v[0]->sp[0], list[i]->v[0]->sp[1], list[i]->v[0]->sp[2]); + fprintf(wrl,"%f %f %f,\n",list[i]->v[1]->sp[0], list[i]->v[1]->sp[1], list[i]->v[1]->sp[2]); + fprintf(wrl,"%f %f %f,\n",list[i]->v[2]->sp[0], list[i]->v[2]->sp[1], list[i]->v[2]->sp[2]); + } + fprintf(wrl," ]\n"); + fprintf(wrl," }\n"); + fprintf(wrl,"\n"); + fprintf(wrl," coordIndex [ # Indexes of poligon Verticies \n"); + + for (i = 0; i < llen; i++) { + fprintf(wrl,"%d, %d, %d, -1\n", i * 3 + 0, i * 3 + 1, i * 3 + 2); + } + fprintf(wrl," ]\n"); + fprintf(wrl,"\n"); + fprintf(wrl," colorPerVertex FALSE\n"); + fprintf(wrl," color Color {\n"); + fprintf(wrl," color [ # RGB colors of each vertex\n"); + + /* Spit out the colors for each face */ + for (i = 0; i < llen; i++) { + if (list[i]->bsort == 1) { /* Positive */ + fprintf(wrl,"%f %f %f,\n", 1.0, 0.3, 0.3); /* Red */ + } else if (list[i]->bsort == 2) { /* Negative */ + fprintf(wrl,"%f %f %f,\n", 0.3, 1.0, 0.3); /* Green */ + } else if (list[i]->bsort == 3) { /* Both */ + fprintf(wrl,"%f %f %f,\n", 1.0, 1.0, 0.3); /* Yellow */ + } else { /* Neither */ + fprintf(wrl,"%f %f %f,\n", 0.3, 0.3, 1.0); /* Blue */ + } + } + fprintf(wrl," ] \n"); + fprintf(wrl," }\n"); + fprintf(wrl," } # end IndexedFaceSet\n"); + + fprintf(wrl," appearance Appearance { \n"); + fprintf(wrl," material Material {\n"); + fprintf(wrl," transparency 0.0\n"); + fprintf(wrl," ambientIntensity 0.3\n"); + fprintf(wrl," shininess 0.5\n"); + fprintf(wrl," }\n"); + fprintf(wrl," }\n"); + fprintf(wrl," } # end Shape\n"); + fprintf(wrl," ]\n"); + fprintf(wrl," } # end of transform\n"); + + /* center of gamut */ + fprintf(wrl,"\n"); + fprintf(wrl," Transform {\n"); + fprintf(wrl," translation %f %f %f\n",0.0, 0.0, 0.0); + fprintf(wrl," children [\n"); + fprintf(wrl," Shape { \n"); + fprintf(wrl," geometry Sphere { radius 0.05 }\n"); + fprintf(wrl," appearance Appearance { material Material { diffuseColor %f %f %f } }\n", 1.0, 1.0, 0.0); + fprintf(wrl," } \n"); + fprintf(wrl," ]\n"); + fprintf(wrl," }\n"); + + fprintf(wrl,"\n"); + fprintf(wrl," ] # end of children for world\n"); + fprintf(wrl,"}\n"); + + if (fclose(wrl) != 0) { + fprintf(stderr,"Error closing output file '%s'\n",filename); + return 2; + } + + return 0; +} + +#endif /* DEBUG_SPLIT_VRML */ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + |