thirdparty/liblcms2/src/cmssm.c

   1 //---------------------------------------------------------------------------------
   2 //
   3 //  Little Color Management System
   4 //  Copyright (c) 1998-2016 Marti Maria Saguer
   5 //
   6 // Permission is hereby granted, free of charge, to any person obtaining
   7 // a copy of this software and associated documentation files (the "Software"),
   8 // to deal in the Software without restriction, including without limitation
   9 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  10 // and/or sell copies of the Software, and to permit persons to whom the Software
  11 // is furnished to do so, subject to the following conditions:
  12 //
  13 // The above copyright notice and this permission notice shall be included in
  14 // all copies or substantial portions of the Software.
  15 //
  16 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  17 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO
  18 // THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  19 // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  20 // LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  21 // OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  22 // WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  23 //
  24 //---------------------------------------------------------------------------------
  25 //
  26
  27 #include "lcms2_internal.h"
  28
  29
  30 // ------------------------------------------------------------------------
  31
  32 // Gamut boundary description by using Jan Morovic's Segment maxima method
  33 // Many thanks to Jan for allowing me to use his algorithm.
  34
  35 // r = C*
  36 // alpha = Hab
  37 // theta = L*
  38
  39 #define SECTORS 16      // number of divisions in alpha and theta
  40
  41 // Spherical coordinates
  42 typedef struct {
  43
  44     cmsFloat64Number r;
  45     cmsFloat64Number alpha;
  46     cmsFloat64Number theta;
  47
  48 } cmsSpherical;
  49
  50 typedef  enum {
  51         GP_EMPTY,
  52         GP_SPECIFIED,
  53         GP_MODELED
  54
  55     } GDBPointType;
  56
  57
  58 typedef struct {
  59
  60     GDBPointType Type;
  61     cmsSpherical p;         // Keep also alpha & theta of maximum
  62
  63 } cmsGDBPoint;
  64
  65
  66 typedef struct {
  67
  68     cmsContext ContextID;
  69     cmsGDBPoint Gamut[SECTORS][SECTORS];
  70
  71 } cmsGDB;
  72
  73
  74 // A line using the parametric form
  75 // P = a + t*u
  76 typedef struct {
  77
  78     cmsVEC3 a;
  79     cmsVEC3 u;
  80
  81 } cmsLine;
  82
  83
  84 // A plane using the parametric form
  85 // Q = b + r*v + s*w
  86 typedef struct {
  87
  88     cmsVEC3 b;
  89     cmsVEC3 v;
  90     cmsVEC3 w;
  91
  92 } cmsPlane;
  93
  94
  95
  96 // --------------------------------------------------------------------------------------------
  97
  98 // ATAN2() which always returns degree positive numbers
  99
 100 static
 101 cmsFloat64Number _cmsAtan2(cmsFloat64Number y, cmsFloat64Number x)
 102 {
 103     cmsFloat64Number a;
 104
 105     // Deal with undefined case
 106     if (x == 0.0 && y == 0.0) return 0;
 107
 108     a = (atan2(y, x) * 180.0) / M_PI;
 109
 110     while (a < 0) {
 111         a += 360;
 112     }
 113
 114     return a;
 115 }
 116
 117 // Convert to spherical coordinates
 118 static
 119 void ToSpherical(cmsSpherical* sp, const cmsVEC3* v)
 120 {
 121
 122     cmsFloat64Number L, a, b;
 123
 124     L = v ->n[VX];
 125     a = v ->n[VY];
 126     b = v ->n[VZ];
 127
 128     sp ->r = sqrt( L*L + a*a + b*b );
 129
 130    if (sp ->r == 0) {
 131         sp ->alpha = sp ->theta = 0;
 132         return;
 133     }
 134
 135     sp ->alpha = _cmsAtan2(a, b);
 136     sp ->theta = _cmsAtan2(sqrt(a*a + b*b), L);
 137 }
 138
 139
 140 // Convert to cartesian from spherical
 141 static
 142 void ToCartesian(cmsVEC3* v, const cmsSpherical* sp)
 143 {
 144     cmsFloat64Number sin_alpha;
 145     cmsFloat64Number cos_alpha;
 146     cmsFloat64Number sin_theta;
 147     cmsFloat64Number cos_theta;
 148     cmsFloat64Number L, a, b;
 149
 150     sin_alpha = sin((M_PI * sp ->alpha) / 180.0);
 151     cos_alpha = cos((M_PI * sp ->alpha) / 180.0);
 152     sin_theta = sin((M_PI * sp ->theta) / 180.0);
 153     cos_theta = cos((M_PI * sp ->theta) / 180.0);
 154
 155     a = sp ->r * sin_theta * sin_alpha;
 156     b = sp ->r * sin_theta * cos_alpha;
 157     L = sp ->r * cos_theta;
 158
 159     v ->n[VX] = L;
 160     v ->n[VY] = a;
 161     v ->n[VZ] = b;
 162 }
 163
 164
 165 // Quantize sector of a spherical coordinate. Saturate 360, 180 to last sector
 166 // The limits are the centers of each sector, so
 167 static
 168 void QuantizeToSector(const cmsSpherical* sp, int* alpha, int* theta)
 169 {
 170     *alpha = (int) floor(((sp->alpha * (SECTORS)) / 360.0) );
 171     *theta = (int) floor(((sp->theta * (SECTORS)) / 180.0) );
 172
 173     if (*alpha >= SECTORS)
 174         *alpha = SECTORS-1;
 175     if (*theta >= SECTORS)
 176         *theta = SECTORS-1;
 177 }
 178
 179
 180 // Line determined by 2 points
 181 static
 182 void LineOf2Points(cmsLine* line, cmsVEC3* a, cmsVEC3* b)
 183 {
 184
 185     _cmsVEC3init(&line ->a, a ->n[VX], a ->n[VY], a ->n[VZ]);
 186     _cmsVEC3init(&line ->u, b ->n[VX] - a ->n[VX],
 187                             b ->n[VY] - a ->n[VY],
 188                             b ->n[VZ] - a ->n[VZ]);
 189 }
 190
 191
 192 // Evaluate parametric line
 193 static
 194 void GetPointOfLine(cmsVEC3* p, const cmsLine* line, cmsFloat64Number t)
 195 {
 196     p ->n[VX] = line ->a.n[VX] + t * line->u.n[VX];
 197     p ->n[VY] = line ->a.n[VY] + t * line->u.n[VY];
 198     p ->n[VZ] = line ->a.n[VZ] + t * line->u.n[VZ];
 199 }
 200
 201
 202
 203 /*
 204     Closest point in sector line1 to sector line2 (both are defined as 0 <=t <= 1)
 205     http://softsurfer.com/Archive/algorithm_0106/algorithm_0106.htm
 206
 207     Copyright 2001, softSurfer (www.softsurfer.com)
 208     This code may be freely used and modified for any purpose
 209     providing that this copyright notice is included with it.
 210     SoftSurfer makes no warranty for this code, and cannot be held
 211     liable for any real or imagined damage resulting from its use.
 212     Users of this code must verify correctness for their application.
 213
 214 */
 215
 216 static
 217 cmsBool ClosestLineToLine(cmsVEC3* r, const cmsLine* line1, const cmsLine* line2)
 218 {
 219     cmsFloat64Number a, b, c, d, e, D;
 220     cmsFloat64Number sc, sN, sD;
 221     //cmsFloat64Number tc; // left for future use
 222     cmsFloat64Number tN, tD;
 223     cmsVEC3 w0;
 224
 225     _cmsVEC3minus(&w0, &line1 ->a, &line2 ->a);
 226
 227     a  = _cmsVEC3dot(&line1 ->u, &line1 ->u);
 228     b  = _cmsVEC3dot(&line1 ->u, &line2 ->u);
 229     c  = _cmsVEC3dot(&line2 ->u, &line2 ->u);
 230     d  = _cmsVEC3dot(&line1 ->u, &w0);
 231     e  = _cmsVEC3dot(&line2 ->u, &w0);
 232
 233     D  = a*c - b * b;      // Denominator
 234     sD = tD = D;           // default sD = D >= 0
 235
 236     if (D <  MATRIX_DET_TOLERANCE) {   // the lines are almost parallel
 237
 238         sN = 0.0;        // force using point P0 on segment S1
 239         sD = 1.0;        // to prevent possible division by 0.0 later
 240         tN = e;
 241         tD = c;
 242     }
 243     else {                // get the closest points on the infinite lines
 244
 245         sN = (b*e - c*d);
 246         tN = (a*e - b*d);
 247
 248         if (sN < 0.0) {       // sc < 0 => the s=0 edge is visible
 249
 250             sN = 0.0;
 251             tN = e;
 252             tD = c;
 253         }
 254         else if (sN > sD) {   // sc > 1 => the s=1 edge is visible
 255             sN = sD;
 256             tN = e + b;
 257             tD = c;
 258         }
 259     }
 260
 261     if (tN < 0.0) {           // tc < 0 => the t=0 edge is visible
 262
 263         tN = 0.0;
 264         // recompute sc for this edge
 265         if (-d < 0.0)
 266             sN = 0.0;
 267         else if (-d > a)
 268             sN = sD;
 269         else {
 270             sN = -d;
 271             sD = a;
 272         }
 273     }
 274     else if (tN > tD) {      // tc > 1 => the t=1 edge is visible
 275
 276         tN = tD;
 277
 278         // recompute sc for this edge
 279         if ((-d + b) < 0.0)
 280             sN = 0;
 281         else if ((-d + b) > a)
 282             sN = sD;
 283         else {
 284             sN = (-d + b);
 285             sD = a;
 286         }
 287     }
 288     // finally do the division to get sc and tc
 289     sc = (fabs(sN) < MATRIX_DET_TOLERANCE ? 0.0 : sN / sD);
 290     //tc = (fabs(tN) < MATRIX_DET_TOLERANCE ? 0.0 : tN / tD); // left for future use.
 291
 292     GetPointOfLine(r, line1, sc);
 293     return TRUE;
 294 }
 295
 296
 297
 298 // ------------------------------------------------------------------ Wrapper
 299
 300
 301 // Allocate & free structure
 302 cmsHANDLE  CMSEXPORT cmsGBDAlloc(cmsContext ContextID)
 303 {
 304     cmsGDB* gbd = (cmsGDB*) _cmsMallocZero(ContextID, sizeof(cmsGDB));
 305     if (gbd == NULL) return NULL;
 306
 307     gbd -> ContextID = ContextID;
 308
 309     return (cmsHANDLE) gbd;
 310 }
 311
 312
 313 void CMSEXPORT cmsGBDFree(cmsHANDLE hGBD)
 314 {
 315     cmsGDB* gbd = (cmsGDB*) hGBD;
 316     if (hGBD != NULL)
 317         _cmsFree(gbd->ContextID, (void*) gbd);
 318 }
 319
 320
 321 // Auxiliar to retrieve a pointer to the segmentr containing the Lab value
 322 static
 323 cmsGDBPoint* GetPoint(cmsGDB* gbd, const cmsCIELab* Lab, cmsSpherical* sp)
 324 {
 325     cmsVEC3 v;
 326     int alpha, theta;
 327
 328     // Housekeeping
 329     _cmsAssert(gbd != NULL);
 330     _cmsAssert(Lab != NULL);
 331     _cmsAssert(sp != NULL);
 332
 333     // Center L* by substracting half of its domain, that's 50
 334     _cmsVEC3init(&v, Lab ->L - 50.0, Lab ->a, Lab ->b);
 335
 336     // Convert to spherical coordinates
 337     ToSpherical(sp, &v);
 338
 339     if (sp ->r < 0 || sp ->alpha < 0 || sp->theta < 0) {
 340          cmsSignalError(gbd ->ContextID, cmsERROR_RANGE, "spherical value out of range");
 341          return NULL;
 342     }
 343
 344     // On which sector it falls?
 345     QuantizeToSector(sp, &alpha, &theta);
 346
 347     if (alpha < 0 || theta < 0 || alpha >= SECTORS || theta >= SECTORS) {
 348          cmsSignalError(gbd ->ContextID, cmsERROR_RANGE, " quadrant out of range");
 349          return NULL;
 350     }
 351
 352     // Get pointer to the sector
 353     return &gbd ->Gamut[theta][alpha];
 354 }
 355
 356 // Add a point to gamut descriptor. Point to add is in Lab color space.
 357 // GBD is centered on a=b=0 and L*=50
 358 cmsBool CMSEXPORT cmsGDBAddPoint(cmsHANDLE hGBD, const cmsCIELab* Lab)
 359 {
 360     cmsGDB* gbd = (cmsGDB*) hGBD;
 361     cmsGDBPoint* ptr;
 362     cmsSpherical sp;
 363
 364
 365     // Get pointer to the sector
 366     ptr = GetPoint(gbd, Lab, &sp);
 367     if (ptr == NULL) return FALSE;
 368
 369     // If no samples at this sector, add it
 370     if (ptr ->Type == GP_EMPTY) {
 371
 372         ptr -> Type = GP_SPECIFIED;
 373         ptr -> p    = sp;
 374     }
 375     else {
 376
 377
 378         // Substitute only if radius is greater
 379         if (sp.r > ptr -> p.r) {
 380
 381                 ptr -> Type = GP_SPECIFIED;
 382                 ptr -> p    = sp;
 383         }
 384     }
 385
 386     return TRUE;
 387 }
 388
 389 // Check if a given point falls inside gamut
 390 cmsBool CMSEXPORT cmsGDBCheckPoint(cmsHANDLE hGBD, const cmsCIELab* Lab)
 391 {
 392     cmsGDB* gbd = (cmsGDB*) hGBD;
 393     cmsGDBPoint* ptr;
 394     cmsSpherical sp;
 395
 396     // Get pointer to the sector
 397     ptr = GetPoint(gbd, Lab, &sp);
 398     if (ptr == NULL) return FALSE;
 399
 400     // If no samples at this sector, return no data
 401     if (ptr ->Type == GP_EMPTY) return FALSE;
 402
 403     // In gamut only if radius is greater
 404
 405     return (sp.r <= ptr -> p.r);
 406 }
 407
 408 // -----------------------------------------------------------------------------------------------------------------------
 409
 410 // Find near sectors. The list of sectors found is returned on Close[].
 411 // The function returns the number of sectors as well.
 412
 413 // 24   9  10  11  12
 414 // 23   8   1   2  13
 415 // 22   7   *   3  14
 416 // 21   6   5   4  15
 417 // 20  19  18  17  16
 418 //
 419 // Those are the relative movements
 420 // {-2,-2}, {-1, -2}, {0, -2}, {+1, -2}, {+2,  -2},
 421 // {-2,-1}, {-1, -1}, {0, -1}, {+1, -1}, {+2,  -1},
 422 // {-2, 0}, {-1,  0}, {0,  0}, {+1,  0}, {+2,   0},
 423 // {-2,+1}, {-1, +1}, {0, +1}, {+1,  +1}, {+2,  +1},
 424 // {-2,+2}, {-1, +2}, {0, +2}, {+1,  +2}, {+2,  +2}};
 425
 426
 427 static
 428 const struct _spiral {
 429
 430     int AdvX, AdvY;
 431
 432     } Spiral[] = { {0,  -1}, {+1, -1}, {+1,  0}, {+1, +1}, {0,  +1}, {-1, +1},
 433                    {-1,  0}, {-1, -1}, {-1, -2}, {0,  -2}, {+1, -2}, {+2, -2},
 434                    {+2, -1}, {+2,  0}, {+2, +1}, {+2, +2}, {+1, +2}, {0,  +2},
 435                    {-1, +2}, {-2, +2}, {-2, +1}, {-2, 0},  {-2, -1}, {-2, -2} };
 436
 437 #define NSTEPS (sizeof(Spiral) / sizeof(struct _spiral))
 438
 439 static
 440 int FindNearSectors(cmsGDB* gbd, int alpha, int theta, cmsGDBPoint* Close[])
 441 {
 442     int nSectors = 0;
 443     int a, t;
 444     cmsUInt32Number i;
 445     cmsGDBPoint* pt;
 446
 447     for (i=0; i < NSTEPS; i++) {
 448
 449         a = alpha + Spiral[i].AdvX;
 450         t = theta + Spiral[i].AdvY;
 451
 452         // Cycle at the end
 453         a %= SECTORS;
 454         t %= SECTORS;
 455
 456         // Cycle at the begin
 457         if (a < 0) a = SECTORS + a;
 458         if (t < 0) t = SECTORS + t;
 459
 460         pt = &gbd ->Gamut[t][a];
 461
 462         if (pt -> Type != GP_EMPTY) {
 463
 464             Close[nSectors++] = pt;
 465         }
 466     }
 467
 468     return nSectors;
 469 }
 470
 471
 472 // Interpolate a missing sector. Method identifies whatever this is top, bottom or mid
 473 static
 474 cmsBool InterpolateMissingSector(cmsGDB* gbd, int alpha, int theta)
 475 {
 476     cmsSpherical sp;
 477     cmsVEC3 Lab;
 478     cmsVEC3 Centre;
 479     cmsLine ray;
 480     int nCloseSectors;
 481     cmsGDBPoint* Close[NSTEPS + 1];
 482     cmsSpherical closel, templ;
 483     cmsLine edge;
 484     int k, m;
 485
 486     // Is that point already specified?
 487     if (gbd ->Gamut[theta][alpha].Type != GP_EMPTY) return TRUE;
 488
 489     // Fill close points
 490     nCloseSectors = FindNearSectors(gbd, alpha, theta, Close);
 491
 492
 493     // Find a central point on the sector
 494     sp.alpha = (cmsFloat64Number) ((alpha + 0.5) * 360.0) / (SECTORS);
 495     sp.theta = (cmsFloat64Number) ((theta + 0.5) * 180.0) / (SECTORS);
 496     sp.r     = 50.0;
 497
 498     // Convert to Cartesian
 499     ToCartesian(&Lab, &sp);
 500
 501     // Create a ray line from centre to this point
 502     _cmsVEC3init(&Centre, 50.0, 0, 0);
 503     LineOf2Points(&ray, &Lab, &Centre);
 504
 505     // For all close sectors
 506     closel.r = 0.0;
 507     closel.alpha = 0;
 508     closel.theta = 0;
 509
 510     for (k=0; k < nCloseSectors; k++) {
 511
 512         for(m = k+1; m < nCloseSectors; m++) {
 513
 514             cmsVEC3 temp, a1, a2;
 515
 516             // A line from sector to sector
 517             ToCartesian(&a1, &Close[k]->p);
 518             ToCartesian(&a2, &Close[m]->p);
 519
 520             LineOf2Points(&edge, &a1, &a2);
 521
 522             // Find a line
 523             ClosestLineToLine(&temp, &ray, &edge);
 524
 525             // Convert to spherical
 526             ToSpherical(&templ, &temp);
 527
 528
 529             if ( templ.r > closel.r &&
 530                  templ.theta >= (theta*180.0/SECTORS) &&
 531                  templ.theta <= ((theta+1)*180.0/SECTORS) &&
 532                  templ.alpha >= (alpha*360.0/SECTORS) &&
 533                  templ.alpha <= ((alpha+1)*360.0/SECTORS)) {
 534
 535                 closel = templ;
 536             }
 537         }
 538     }
 539
 540     gbd ->Gamut[theta][alpha].p = closel;
 541     gbd ->Gamut[theta][alpha].Type = GP_MODELED;
 542
 543     return TRUE;
 544
 545 }
 546
 547
 548 // Interpolate missing parts. The algorithm fist computes slices at
 549 // theta=0 and theta=Max.
 550 cmsBool CMSEXPORT cmsGDBCompute(cmsHANDLE hGBD, cmsUInt32Number dwFlags)
 551 {
 552     int alpha, theta;
 553     cmsGDB* gbd = (cmsGDB*) hGBD;
 554
 555     _cmsAssert(hGBD != NULL);
 556
 557     // Interpolate black
 558     for (alpha = 0; alpha < SECTORS; alpha++) {
 559
 560         if (!InterpolateMissingSector(gbd, alpha, 0)) return FALSE;
 561     }
 562
 563     // Interpolate white
 564     for (alpha = 0; alpha < SECTORS; alpha++) {
 565
 566         if (!InterpolateMissingSector(gbd, alpha, SECTORS-1)) return FALSE;
 567     }
 568
 569
 570     // Interpolate Mid
 571     for (theta = 1; theta < SECTORS; theta++) {
 572         for (alpha = 0; alpha < SECTORS; alpha++) {
 573
 574             if (!InterpolateMissingSector(gbd, alpha, theta)) return FALSE;
 575         }
 576     }
 577
 578     // Done
 579     return TRUE;
 580
 581     cmsUNUSED_PARAMETER(dwFlags);
 582 }
 583
 584
 585
 586
 587 // --------------------------------------------------------------------------------------------------------
 588
 589 // Great for debug, but not suitable for real use
 590
 591 #if 0
 592 cmsBool cmsGBDdumpVRML(cmsHANDLE hGBD, const char* fname)
 593 {
 594     FILE* fp;
 595     int   i, j;
 596     cmsGDB* gbd = (cmsGDB*) hGBD;
 597     cmsGDBPoint* pt;
 598
 599     fp = fopen (fname, "wt");
 600     if (fp == NULL)
 601         return FALSE;
 602
 603     fprintf (fp, "#VRML V2.0 utf8\n");
 604
 605     // set the viewing orientation and distance
 606     fprintf (fp, "DEF CamTest Group {\n");
 607     fprintf (fp, "\tchildren [\n");
 608     fprintf (fp, "\t\tDEF Cameras Group {\n");
 609     fprintf (fp, "\t\t\tchildren [\n");
 610     fprintf (fp, "\t\t\t\tDEF DefaultView Viewpoint {\n");
 611     fprintf (fp, "\t\t\t\t\tposition 0 0 340\n");
 612     fprintf (fp, "\t\t\t\t\torientation 0 0 1 0\n");
 613     fprintf (fp, "\t\t\t\t\tdescription \"default view\"\n");
 614     fprintf (fp, "\t\t\t\t}\n");
 615     fprintf (fp, "\t\t\t]\n");
 616     fprintf (fp, "\t\t},\n");
 617     fprintf (fp, "\t]\n");
 618     fprintf (fp, "}\n");
 619
 620     // Output the background stuff
 621     fprintf (fp, "Background {\n");
 622     fprintf (fp, "\tskyColor [\n");
 623     fprintf (fp, "\t\t.5 .5 .5\n");
 624     fprintf (fp, "\t]\n");
 625     fprintf (fp, "}\n");
 626
 627     // Output the shape stuff
 628     fprintf (fp, "Transform {\n");
 629     fprintf (fp, "\tscale .3 .3 .3\n");
 630     fprintf (fp, "\tchildren [\n");
 631
 632     // Draw the axes as a shape:
 633     fprintf (fp, "\t\tShape {\n");
 634     fprintf (fp, "\t\t\tappearance Appearance {\n");
 635     fprintf (fp, "\t\t\t\tmaterial Material {\n");
 636     fprintf (fp, "\t\t\t\t\tdiffuseColor 0 0.8 0\n");
 637     fprintf (fp, "\t\t\t\t\temissiveColor 1.0 1.0 1.0\n");
 638     fprintf (fp, "\t\t\t\t\tshininess 0.8\n");
 639     fprintf (fp, "\t\t\t\t}\n");
 640     fprintf (fp, "\t\t\t}\n");
 641     fprintf (fp, "\t\t\tgeometry IndexedLineSet {\n");
 642     fprintf (fp, "\t\t\t\tcoord Coordinate {\n");
 643     fprintf (fp, "\t\t\t\t\tpoint [\n");
 644     fprintf (fp, "\t\t\t\t\t0.0 0.0 0.0,\n");
 645     fprintf (fp, "\t\t\t\t\t%f 0.0 0.0,\n",  255.0);
 646     fprintf (fp, "\t\t\t\t\t0.0 %f 0.0,\n",  255.0);
 647     fprintf (fp, "\t\t\t\t\t0.0 0.0 %f]\n",  255.0);
 648     fprintf (fp, "\t\t\t\t}\n");
 649     fprintf (fp, "\t\t\t\tcoordIndex [\n");
 650     fprintf (fp, "\t\t\t\t\t0, 1, -1\n");
 651     fprintf (fp, "\t\t\t\t\t0, 2, -1\n");
 652     fprintf (fp, "\t\t\t\t\t0, 3, -1]\n");
 653     fprintf (fp, "\t\t\t}\n");
 654     fprintf (fp, "\t\t}\n");
 655
 656
 657     fprintf (fp, "\t\tShape {\n");
 658     fprintf (fp, "\t\t\tappearance Appearance {\n");
 659     fprintf (fp, "\t\t\t\tmaterial Material {\n");
 660     fprintf (fp, "\t\t\t\t\tdiffuseColor 0 0.8 0\n");
 661     fprintf (fp, "\t\t\t\t\temissiveColor 1 1 1\n");
 662     fprintf (fp, "\t\t\t\t\tshininess 0.8\n");
 663     fprintf (fp, "\t\t\t\t}\n");
 664     fprintf (fp, "\t\t\t}\n");
 665     fprintf (fp, "\t\t\tgeometry PointSet {\n");
 666
 667     // fill in the points here
 668     fprintf (fp, "\t\t\t\tcoord Coordinate {\n");
 669     fprintf (fp, "\t\t\t\t\tpoint [\n");
 670
 671     // We need to transverse all gamut hull.
 672     for (i=0; i < SECTORS; i++)
 673         for (j=0; j < SECTORS; j++) {
 674
 675             cmsVEC3 v;
 676
 677             pt = &gbd ->Gamut[i][j];
 678             ToCartesian(&v, &pt ->p);
 679
 680             fprintf (fp, "\t\t\t\t\t%g %g %g", v.n[0]+50, v.n[1], v.n[2]);
 681
 682             if ((j == SECTORS - 1) && (i == SECTORS - 1))
 683                 fprintf (fp, "]\n");
 684             else
 685                 fprintf (fp, ",\n");
 686
 687         }
 688
 689         fprintf (fp, "\t\t\t\t}\n");
 690
 691
 692
 693     // fill in the face colors
 694     fprintf (fp, "\t\t\t\tcolor Color {\n");
 695     fprintf (fp, "\t\t\t\t\tcolor [\n");
 696
 697     for (i=0; i < SECTORS; i++)
 698         for (j=0; j < SECTORS; j++) {
 699
 700            cmsVEC3 v;
 701
 702             pt = &gbd ->Gamut[i][j];
 703
 704
 705             ToCartesian(&v, &pt ->p);
 706
 707
 708         if (pt ->Type == GP_EMPTY)
 709             fprintf (fp, "\t\t\t\t\t%g %g %g", 0.0, 0.0, 0.0);
 710         else
 711             if (pt ->Type == GP_MODELED)
 712                 fprintf (fp, "\t\t\t\t\t%g %g %g", 1.0, .5, .5);
 713             else {
 714                 fprintf (fp, "\t\t\t\t\t%g %g %g", 1.0, 1.0, 1.0);
 715
 716             }
 717
 718         if ((j == SECTORS - 1) && (i == SECTORS - 1))
 719                 fprintf (fp, "]\n");
 720             else
 721                 fprintf (fp, ",\n");
 722     }
 723     fprintf (fp, "\t\t\t}\n");
 724
 725
 726     fprintf (fp, "\t\t\t}\n");
 727     fprintf (fp, "\t\t}\n");
 728     fprintf (fp, "\t]\n");
 729     fprintf (fp, "}\n");
 730
 731     fclose (fp);
 732
 733     return TRUE;
 734 }
 735 #endif
 736