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