summaryrefslogtreecommitdiff
path: root/gamut/gamut.c
diff options
context:
space:
mode:
Diffstat (limited to 'gamut/gamut.c')
-rw-r--r--gamut/gamut.c682
1 files changed, 332 insertions, 350 deletions
diff --git a/gamut/gamut.c b/gamut/gamut.c
index 3fc1c97..dd4f43b 100644
--- a/gamut/gamut.c
+++ b/gamut/gamut.c
@@ -56,7 +56,7 @@
#define TRIANG_TOL 1e-10 /* [1e-10] Triangulation tollerance */
#define NORM_LOG_POW 0.25 /* [0.25] Normal, colorspace lopow value */
-#define RAST_LOG_POW 0.10 /* [0.10] Raster lopow value (is 0.05 too extreme ??) */
+#define RAST_LOG_POW 0.05 /* [0.05] Raster lopow value */
#undef TEST_CONVEX_HULL /* Use pure convex hull, not log hull */
@@ -156,7 +156,6 @@ static int getssvert(gamut *s, double *rad, double pos[3], double norm[3], int i
static void startnexttri(gamut *s);
static int getnexttri(gamut *s, int v[3]);
static double volume(gamut *s);
-static int write_to_vrml(gamut *s, vrml *wrl, double trans, int docusps);
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);
@@ -175,9 +174,8 @@ static int compute_vector_isect(gamut *s, double *p1, double *p2, double *min, d
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 nexpintersect(gamut *s, gamut *s1, gamut *s2);
-static int expdstbysrcmdst(gamut *s, gamut *s1, gamut *s2, gamut *s3,
- void (*cvect)(void *cntx, double *p2, double *p1), void *cntx);
+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);
static void compgawb(gamut *s);
@@ -658,8 +656,7 @@ int isRast /* Flag indicating Raster rather than colorspace */
s->getvert = getvert;
s->volume = volume;
s->intersect = intersect;
- s->nexpintersect = nexpintersect;
- s->expdstbysrcmdst = expdstbysrcmdst;
+ s->compdstgamut = compdstgamut;
s->radial = radial;
s->nradial = nradial;
s->nearest = nearest;
@@ -670,7 +667,6 @@ int isRast /* Flag indicating Raster rather than colorspace */
s->getwb = getwb;
s->setcusps = setcusps;
s->getcusps = getcusps;
- s->write_to_vrml = write_to_vrml;
s->write_vrml = write_vrml;
s->write_trans_vrml = write_trans_vrml;
s->write_gam = write_gam;
@@ -1066,13 +1062,44 @@ double pp[3] /* rectangular coordinate of point */
}
/* ------------------------------------ */
-
-/* intersect implementation */
-/* Assumes s has been initialised. */
-static void intersect_imp(gamut *s, gamut *sa, gamut *sb) {
+/* 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++) {
@@ -1092,7 +1119,6 @@ static void intersect_imp(gamut *s, gamut *sa, gamut *sb) {
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 */
@@ -1135,142 +1161,6 @@ static void intersect_imp(gamut *s, gamut *sa, gamut *sb) {
} END_FOR_ALL_ITEMS(tp1);
}
-}
-
-/* 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) {
- gamut *ss;
- int j;
-
- 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->sres = sa->sres > sb->sres ? sa->sres : sb->sres;
-
- s->isJab = sa->isJab;
-
- /* Clear some flags */
- s->cswbset = 0;
- s->cswbset = 0;
- s->dcuspixs = 0;
-
- /* If either is a raster gamut, make it a raster gamut */
- if (sa->isRast || sb->isRast)
- 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 */
- }
-
- for (j = 0; j < 3; j++)
- s->cent[j] = sa->cent[j];
-
- /* Grab white & black from whichever source has it */
- if (sb->cswbset)
- ss = sb;
- else if (sb->cswbset)
- ss = sa;
-
- if (ss->cswbset) {
- for (j = 0; j < 3; j++) {
- s->cs_wp[j] = ss->cs_wp[j];
- s->cs_bp[j] = ss->cs_bp[j];
- s->cs_kp[j] = ss->cs_kp[j];
- }
- s->cswbset = ss->cswbset;
- }
-
- /* Don't filter the points (gives a more accurate result ?) */
- s->nofilter = 1;
-
- intersect_imp(s, sa, sb);
-
- if (sa->gawbset) {
- compgawb(s);
- }
-
- s->nofilter = 0;
-
- return 0;
-}
-
-/* ------------------------------------ */
-
-/* Initialise this gamut with neutral axis points from sa, */
-/* and then intersected with sb. */
-/* Return NZ on error. */
-/* Expand sb with neutral points, and then intersect */
-/* with sa. */
-
-/* Return 1 if gamuts are not compatible */
-/* (We assume that the this gamut is currently empty) */
-static int nexpintersect(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->sres = sa->sres > sb->sres ? sa->sres : sb->sres;
-
- s->isJab = sa->isJab;
-
- /* If either is a raster gamut, make it a raster gamut */
- if (sa->isRast || sb->isRast)
- 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 */
- }
-
- for (j = 0; j < 3; j++)
- s->cent[j] = sa->cent[j];
-
- /* Clear some flags */
- s->cswbset = 0;
- s->cswbset = 0;
- s->dcuspixs = 0;
-
- /* Don't filter the points (gives a more accurate result ?) */
- s->nofilter = 1;
-
- /* Number of points to generate */
- k = 10;
- if (sa->nv <= 10)
- k = 5;
-
- /* Generate points from black to white */
- for (i = 0; i < k; i++) {
- double pp[3];
- double bf = i/(k-1.0);
-
- icmBlend3(pp, sa->cs_bp, sa->cs_wp, bf);
- expand_gamut(s, pp);
- }
-
- /* let intersect_imp to the hard work */
- intersect_imp(s, sa, sb);
s->nofilter = 0;
@@ -1279,36 +1169,49 @@ static int nexpintersect(gamut *s, gamut *sa, gamut *sb) {
/* ------------------------------------ */
/*
- Initialise this gamut with the image/destination gamut
- expanded by the amount that dest colorspace is outside
- the source colorspace gamut.
+ 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 as p1 -> p2 if it is supplied, or radially
- inwards if it is not. p2 is a "center" point to compute depth to.
+ 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)
+ (We assume that the this gamut is currently empty)
*/
#define MXNIS 40 /* Maximum raw intersections handled */
-static int expdstbysrcmdst(
-gamut *s, /* Gamut to be expanded */
-gamut *s1, /* (Image) destination gamut to be expanded */
+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 /* which returns p2 which is in desired direction from given p1 */
+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 */
- gispnt lp1[MXNIS], lp2[MXNIS], lp3[MXNIS]; /* Lists of intersections */
int ii, jj, kk; /* List indexes */
-//printf("\n~1 expdstbysrcmdst() called\n");
-
if (s1->compatible(s1, s2) == 0
|| s1->compatible(s2, s3) == 0)
return 1;
@@ -1323,14 +1226,6 @@ void *cntx /* which returns p2 which is in desired direction from given p1 */
s->isJab = s1->isJab;
s->isRast = s1->isRast;
- 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 */
- }
-
for (j = 0; j < 3; j++)
s->cent[j] = s1->cent[j];
@@ -1339,6 +1234,12 @@ void *cntx /* which returns p2 which is in desired direction from given p1 */
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;
@@ -1346,12 +1247,26 @@ void *cntx /* which returns p2 which is in desired direction from given p1 */
ss[1] = s2;
ss[2] = s3;
- /* Use all the triangle vertices from the two/three gamuts */
+//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 (k = 0; k < 3; k++) { /* For img, src & dst */
- /* For each vertex */
for (i = 0; i < ss[k]->nv; i++) {
double pp[3], ppv, p2[3];
double rr, r4;
@@ -1361,39 +1276,22 @@ void *cntx /* which returns p2 which is in desired direction from given p1 */
icmCpy3(pp, ss[k]->verts[i]->p); /* Point in question */
- if (k == 0) { /* Seed with all points from gamut */
- expand_gamut(s, pp); /* to ensure result can't be less. */
- }
-
//printf("\n~1 k %d, point %d: %f %f %f\n", k,i,pp[0],pp[1],pp[2]);
-
- /* Get the mapping vector */
if (cvect != NULL)
- cvect(cntx, p2, pp); /* Get mapping direction to center */
+ cvect(cntx, p2, pp); /* Get mapping direction */
else
- icmCpy3(p2, ss[k]->cent); /* Radial vector to center */
- icmNormalize33(pp, pp, p2, 1.0); /* Make p2->pp length 1.0 */
-
-//printf("~1 k %d, center %d: %f %f %f\n", k,i,p2[0],p2[1],p2[2]);
-
- /* Locate the intersecting segments for each gamut. */
- /* The returned parameter value will be >= 1.0 at and beyond the center point */
- /* and < 1.0 on the pp side. */
- /* We need intersections for all three for this to be a potential expansion point. */
- if ((ll1 = s1->vector_isectns(s1, pp, p2, lp1, MXNIS)) == 0) { /* Dest gamut */
-//printf("~1 no dst intersection\n");
- continue;
- }
- if ((ll2 = s2->vector_isectns(s2, pp, p2, lp2, MXNIS)) == 0) { /* Src space */
-//printf("~1 no sc intersection\n");
- continue;
- }
-
- if ((ll3 = s3->vector_isectns(s3, pp, p2, lp3, MXNIS)) == 0) { /* Dst space */
-//printf("~1 no dc intersection\n");
- continue;
- }
-
+ 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++)
@@ -1405,120 +1303,264 @@ void *cntx /* which returns p2 which is in desired direction from given p1 */
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) {
- /* We're only interested in the most outside intersection of each surface */
- /* on the side of the point in question. (we're ignoring and complex */
- /* topology) */
- if (lp1[0].pv > (1.0 - 1e-8)) {
-//printf("~1 dst point is on other side\n");
- continue;
- }
- if (lp2[0].pv > (1.0 - 1e-8)) {
-//printf("~1 sc point is on other side\n");
- continue;
- }
- if (lp3[0].pv > (1.0 - 1e-8)) {
-//printf("~1 sc point is on other side\n");
- continue;
- }
-
- /* Make sure that sc is inside dc, */
- /* and sc is on or above dst */
- if (lp2[0].pv > (lp3[0].pv - 1e-8)
- && lp2[0].pv <= (lp1[0].pv + 1e-8)) {
- double ex[3], sf;
-
- /* Expand point by up to dst - src */
- icmSub3(ex, lp3[0].ip, lp2[0].ip);
-
- /* Make expansion proportional to how far dst */
- /* is to sc */
- sf = (1.0 - lp1[0].pv)/(1.0 - lp2[0].pv);
- icmScale3(ex, ex, sf);
- icmAdd3(pp, lp1[0].ip, ex);
-
-//printf("~1 expanding point by sf %f = %f %f %f\nb",sf,ex[0],ex[1],ex[2]);
-//printf("~1 to %f %f %f\nb",pp[0],pp[1],pp[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);
- } else {
-//printf("~1 lp2 %f <= lp3 %f || lp2 %f > lp2 %f\n",lp2[0].pv,lp3[0].pv - 1e-8, lp2[0].pv, lp1[0].pv + 1e-8);
-//printf("~1 dc is inside sc or sc is inside dst\n");
}
}
}
- /* Generate points along the intersection of sc and dc, map them */
- /* to the image/dest gamut, and add them as well. This is to properly define */
- /* the edges of the expansion zone. */
+#ifdef STATS
+ tedges = edgestested = testhits = 0;
+#endif
- /* For sc on dc and then dc on sc */
- for (k = 0; k < 2; k++) {
- gamut *ss1, *ss2;
- gtri *tp1, *tp2; /* Triangle pointers */
+ /* 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) {
- ss1 = s2;
- ss2 = s3;
- } else {
- ss1 = s3;
- ss2 = s2;
+ 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 */
}
- /* Find the edges that intersect the other gamut */
- tp1 = ss1->tris;
- FOR_ALL_ITEMS(gtri, tp1) { /* For all ss1 triangles */
-
- for (j = 0; j < 3; j++) { /* For all edges in ss1 triangle */
- /* If edge passes through the other gamut */
- if ((tp1->e[j]->v[0]->f ^ tp1->e[j]->v[1]->f) & GVERT_ISOS) {
+ /* Now find the edges that intersect the other gamut */
+ tpa = sa->tris;
+ FOR_ALL_ITEMS(gtri, tpa) {
- /* Exhaustive search of other triangles in ss2, */
- /* to find the one that the edge intersects with. */
- tp2 = ss2->tris;
- FOR_ALL_ITEMS(gtri, tp2) {
+ 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 (tp2->mix[1][i] < tp1->mix[0][i]
- || tp2->mix[0][i] > tp1->mix[1][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(ss1, &pv, pp, tp1->e[j]->v[0]->p, tp1->e[j]->v[1]->p, tp2)
+ 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)) {
- double p2[3];
-
- /* Got intersection point pp. */
-
- /* Get the mapping vector */
+ /* 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 to center */
+ cvect(cntx, p2, pp); /* Get mapping direction */
else
- icmCpy3(p2, ss[k]->cent); /* Radial vector to center */
- icmNormalize33(pp, pp, p2, 1.0); /* Make p2->pp length 1.0 */
-
- /* Locate the intersecting segments for the img/dest gamut. */
- /* The returned parameter value will be >= 1.0 at and beyond */
- /* the center point and < 1.0 on the pp side. */
- /* We're only going to bother with the most outside point */
- if ((ll1 = s1->vector_isectns(s1, pp, p2, lp1, MXNIS)) == 0
- || lp1[0].pv > (1.0 - 1e-8)) {
- continue;
+ 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);
}
-
- expand_gamut(s, pp);
}
- } END_FOR_ALL_ITEMS(tp2);
+ } END_FOR_ALL_ITEMS(tpb);
}
}
- } END_FOR_ALL_ITEMS(tp1);
+ } 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;
@@ -3907,7 +3949,7 @@ double *nin /* Normalised center relative point */
if (fabs(denom) < 1e-9) {
/* Hmm. The ray is paralell to the triangle ? */
- error("radial_point: failed to intersect radial triangle, num %e, denom %e\n",num,denom);
+ error("radial_point: failed to intersect radial triangle\n");
}
rv = num/denom;
@@ -5582,7 +5624,7 @@ int ll /* Size of list. */
return 0; /* Too few to be useful */
}
- /* Now we need to turn the raw intersections into sanitized segment pairs. */
+ /* 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)
@@ -5759,67 +5801,7 @@ int ll /* Size of list. */
#endif /* INTERSECT_DEBUG */
/* ===================================================== */
-
-/* Append gamut to an open VRML/X3d file */
-/* Return non-zero on error */
-static int write_to_vrml(
- gamut *s,
- vrml *wrl,
- double trans, /* Transparency of gamut */
- int docusps /* Add cusps to vrml */
-) {
- int i;
- gtri *tp; /* Triangle pointer */
-
- if IS_LIST_EMPTY(s->tris)
- triangulate(s);
-
- 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++)
- wrl->add_marker(wrl, s->cusps[i], ccolors[i], 2.0);
- }
-
- wrl->start_line_set(wrl, 0);
-
- /* Spit out the vertex values, in order. */
- for (i = 0; i < s->nv; i++) {
- double out[3];
-
- if (!(s->verts[i]->f & GVERT_TRI))
- continue;
-
- /* 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];
- wrl->add_vertex(wrl, 0, out);
- }
-
- /* Add the surface triangles */
- tp = s->tris;
- FOR_ALL_ITEMS(gtri, tp) {
- int ix[3];
- ix[0] = tp->v[0]->tn;
- ix[1] = tp->v[1]->tn;
- ix[2] = tp->v[2]->tn;
- wrl->add_triangle(wrl, 0, ix);
- } END_FOR_ALL_ITEMS(tp);
-
- wrl->make_triangles_vc(wrl, 0, trans);
-
- return 0;
-}
-
-/* Write gamut to a VRML/X3d file */
+/* Write to a VRML/X3d file */
/* Return non-zero on error */
static int write_vrml(
gamut *s,
@@ -5830,7 +5812,7 @@ int docusps /* Non-zero if cusp points are to be marked */
return write_trans_vrml(s, filename, doaxes, docusps, NULL, NULL);
}
-/* Write gamut to a VRML/X3d file */
+/* Write to a VRML/X3d file */
/* Return non-zero on error */
static int write_trans_vrml(
gamut *s,
@@ -5944,7 +5926,7 @@ void *cntx
ix[2] = tp->v[2]->tn;
wrl->add_triangle(wrl, 0, ix);
} END_FOR_ALL_ITEMS(tp);
-#endif /* !SHOW_BUCKETS */
+#endif /* SHOW_BUCKETS */
{
double rgb[3];
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-07 10:06:21 +0000