diff options
Diffstat (limited to 'gamut/gamut.c')
| -rw-r--r--[-rwxr-xr-x] | gamut/gamut.c | 682 |
1 files changed, 332 insertions, 350 deletions
diff --git a/gamut/gamut.c b/gamut/gamut.c index 3fc1c97..dd4f43b 100755..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]; |