1 //---------------------------------------------------------------------------------
2 //
3 //  Little Color Management System
4 //  Copyright (c) 1998-2011 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
_cmsAtan2(cmsFloat64Number y,cmsFloat64Number x)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
ToSpherical(cmsSpherical * sp,const cmsVEC3 * v)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
ToCartesian(cmsVEC3 * v,const cmsSpherical * sp)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
QuantizeToSector(const cmsSpherical * sp,int * alpha,int * theta)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
LineOf2Points(cmsLine * line,cmsVEC3 * a,cmsVEC3 * b)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
GetPointOfLine(cmsVEC3 * p,const cmsLine * line,cmsFloat64Number t)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
ClosestLineToLine(cmsVEC3 * r,const cmsLine * line1,const cmsLine * line2)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, tN, tD;
222     cmsVEC3 w0;
223 
224     _cmsVEC3minus(&w0, &line1 ->a, &line2 ->a);
225 
226     a  = _cmsVEC3dot(&line1 ->u, &line1 ->u);
227     b  = _cmsVEC3dot(&line1 ->u, &line2 ->u);
228     c  = _cmsVEC3dot(&line2 ->u, &line2 ->u);
229     d  = _cmsVEC3dot(&line1 ->u, &w0);
230     e  = _cmsVEC3dot(&line2 ->u, &w0);
231 
232     D  = a*c - b * b;      // Denominator
233     sD = tD = D;           // default sD = D >= 0
234 
235     if (D <  MATRIX_DET_TOLERANCE) {   // the lines are almost parallel
236 
237         sN = 0.0;        // force using point P0 on segment S1
238         sD = 1.0;        // to prevent possible division by 0.0 later
239         tN = e;
240         tD = c;
241     }
242     else {                // get the closest points on the infinite lines
243 
244         sN = (b*e - c*d);
245         tN = (a*e - b*d);
246 
247         if (sN < 0.0) {       // sc < 0 => the s=0 edge is visible
248 
249             sN = 0.0;
250             tN = e;
251             tD = c;
252         }
253         else if (sN > sD) {   // sc > 1 => the s=1 edge is visible
254             sN = sD;
255             tN = e + b;
256             tD = c;
257         }
258     }
259 
260     if (tN < 0.0) {           // tc < 0 => the t=0 edge is visible
261 
262         tN = 0.0;
263         // recompute sc for this edge
264         if (-d < 0.0)
265             sN = 0.0;
266         else if (-d > a)
267             sN = sD;
268         else {
269             sN = -d;
270             sD = a;
271         }
272     }
273     else if (tN > tD) {      // tc > 1 => the t=1 edge is visible
274 
275         tN = tD;
276 
277         // recompute sc for this edge
278         if ((-d + b) < 0.0)
279             sN = 0;
280         else if ((-d + b) > a)
281             sN = sD;
282         else {
283             sN = (-d + b);
284             sD = a;
285         }
286     }
287     // finally do the division to get sc and tc
288     sc = (fabs(sN) < MATRIX_DET_TOLERANCE ? 0.0 : sN / sD);
289     tc = (fabs(tN) < MATRIX_DET_TOLERANCE ? 0.0 : tN / tD);
290 
291     GetPointOfLine(r, line1, sc);
292     return TRUE;
293 }
294 
295 
296 
297 // ------------------------------------------------------------------ Wrapper
298 
299 
300 // Allocate & free structure
cmsGBDAlloc(cmsContext ContextID)301 cmsHANDLE  CMSEXPORT cmsGBDAlloc(cmsContext ContextID)
302 {
303     cmsGDB* gbd = (cmsGDB*) _cmsMallocZero(ContextID, sizeof(cmsGDB));
304     if (gbd == NULL) return NULL;
305 
306     gbd -> ContextID = ContextID;
307 
308     return (cmsHANDLE) gbd;
309 }
310 
311 
cmsGBDFree(cmsHANDLE hGBD)312 void CMSEXPORT cmsGBDFree(cmsHANDLE hGBD)
313 {
314     cmsGDB* gbd = (cmsGDB*) hGBD;
315     if (hGBD != NULL)
316         _cmsFree(gbd->ContextID, (void*) gbd);
317 }
318 
319 
320 // Auxiliar to retrieve a pointer to the segmentr containing the Lab value
321 static
GetPoint(cmsGDB * gbd,const cmsCIELab * Lab,cmsSpherical * sp)322 cmsGDBPoint* GetPoint(cmsGDB* gbd, const cmsCIELab* Lab, cmsSpherical* sp)
323 {
324     cmsVEC3 v;
325     int alpha, theta;
326 
327     // Housekeeping
328     _cmsAssert(gbd != NULL);
329     _cmsAssert(Lab != NULL);
330     _cmsAssert(sp != NULL);
331 
332     // Center L* by substracting half of its domain, that's 50
333     _cmsVEC3init(&v, Lab ->L - 50.0, Lab ->a, Lab ->b);
334 
335     // Convert to spherical coordinates
336     ToSpherical(sp, &v);
337 
338     if (sp ->r < 0 || sp ->alpha < 0 || sp->theta < 0) {
339          cmsSignalError(gbd ->ContextID, cmsERROR_RANGE, "spherical value out of range");
340          return NULL;
341     }
342 
343     // On which sector it falls?
344     QuantizeToSector(sp, &alpha, &theta);
345 
346     if (alpha < 0 || theta < 0 || alpha >= SECTORS || theta >= SECTORS) {
347          cmsSignalError(gbd ->ContextID, cmsERROR_RANGE, " quadrant out of range");
348          return NULL;
349     }
350 
351     // Get pointer to the sector
352     return &gbd ->Gamut[theta][alpha];
353 }
354 
355 // Add a point to gamut descriptor. Point to add is in Lab color space.
356 // GBD is centered on a=b=0 and L*=50
cmsGDBAddPoint(cmsHANDLE hGBD,const cmsCIELab * Lab)357 cmsBool CMSEXPORT cmsGDBAddPoint(cmsHANDLE hGBD, const cmsCIELab* Lab)
358 {
359     cmsGDB* gbd = (cmsGDB*) hGBD;
360     cmsGDBPoint* ptr;
361     cmsSpherical sp;
362 
363 
364     // Get pointer to the sector
365     ptr = GetPoint(gbd, Lab, &sp);
366     if (ptr == NULL) return FALSE;
367 
368     // If no samples at this sector, add it
369     if (ptr ->Type == GP_EMPTY) {
370 
371         ptr -> Type = GP_SPECIFIED;
372         ptr -> p    = sp;
373     }
374     else {
375 
376 
377         // Substitute only if radius is greater
378         if (sp.r > ptr -> p.r) {
379 
380                 ptr -> Type = GP_SPECIFIED;
381                 ptr -> p    = sp;
382         }
383     }
384 
385     return TRUE;
386 }
387 
388 // Check if a given point falls inside gamut
cmsGDBCheckPoint(cmsHANDLE hGBD,const cmsCIELab * Lab)389 cmsBool CMSEXPORT cmsGDBCheckPoint(cmsHANDLE hGBD, const cmsCIELab* Lab)
390 {
391     cmsGDB* gbd = (cmsGDB*) hGBD;
392     cmsGDBPoint* ptr;
393     cmsSpherical sp;
394 
395     // Get pointer to the sector
396     ptr = GetPoint(gbd, Lab, &sp);
397     if (ptr == NULL) return FALSE;
398 
399     // If no samples at this sector, return no data
400     if (ptr ->Type == GP_EMPTY) return FALSE;
401 
402     // In gamut only if radius is greater
403 
404     return (sp.r <= ptr -> p.r);
405 }
406 
407 // -----------------------------------------------------------------------------------------------------------------------
408 
409 // Find near sectors. The list of sectors found is returned on Close[].
410 // The function returns the number of sectors as well.
411 
412 // 24   9  10  11  12
413 // 23   8   1   2  13
414 // 22   7   *   3  14
415 // 21   6   5   4  15
416 // 20  19  18  17  16
417 //
418 // Those are the relative movements
419 // {-2,-2}, {-1, -2}, {0, -2}, {+1, -2}, {+2,  -2},
420 // {-2,-1}, {-1, -1}, {0, -1}, {+1, -1}, {+2,  -1},
421 // {-2, 0}, {-1,  0}, {0,  0}, {+1,  0}, {+2,   0},
422 // {-2,+1}, {-1, +1}, {0, +1}, {+1,  +1}, {+2,  +1},
423 // {-2,+2}, {-1, +2}, {0, +2}, {+1,  +2}, {+2,  +2}};
424 
425 
426 static
427 const struct _spiral {
428 
429     int AdvX, AdvY;
430 
431     } Spiral[] = { {0,  -1}, {+1, -1}, {+1,  0}, {+1, +1}, {0,  +1}, {-1, +1},
432                    {-1,  0}, {-1, -1}, {-1, -2}, {0,  -2}, {+1, -2}, {+2, -2},
433                    {+2, -1}, {+2,  0}, {+2, +1}, {+2, +2}, {+1, +2}, {0,  +2},
434                    {-1, +2}, {-2, +2}, {-2, +1}, {-2, 0},  {-2, -1}, {-2, -2} };
435 
436 #define NSTEPS (sizeof(Spiral) / sizeof(struct _spiral))
437 
438 static
FindNearSectors(cmsGDB * gbd,int alpha,int theta,cmsGDBPoint * Close[])439 int FindNearSectors(cmsGDB* gbd, int alpha, int theta, cmsGDBPoint* Close[])
440 {
441     int nSectors = 0;
442     int a, t;
443     cmsUInt32Number i;
444     cmsGDBPoint* pt;
445 
446     for (i=0; i < NSTEPS; i++) {
447 
448         a = alpha + Spiral[i].AdvX;
449         t = theta + Spiral[i].AdvY;
450 
451         // Cycle at the end
452         a %= SECTORS;
453         t %= SECTORS;
454 
455         // Cycle at the begin
456         if (a < 0) a = SECTORS + a;
457         if (t < 0) t = SECTORS + t;
458 
459         pt = &gbd ->Gamut[t][a];
460 
461         if (pt -> Type != GP_EMPTY) {
462 
463             Close[nSectors++] = pt;
464         }
465     }
466 
467     return nSectors;
468 }
469 
470 
471 // Interpolate a missing sector. Method identifies whatever this is top, bottom or mid
472 static
InterpolateMissingSector(cmsGDB * gbd,int alpha,int theta)473 cmsBool InterpolateMissingSector(cmsGDB* gbd, int alpha, int theta)
474 {
475     cmsSpherical sp;
476     cmsVEC3 Lab;
477     cmsVEC3 Centre;
478     cmsLine ray;
479     int nCloseSectors;
480     cmsGDBPoint* Close[NSTEPS + 1];
481     cmsSpherical closel, templ;
482     cmsLine edge;
483     int k, m;
484 
485     // Is that point already specified?
486     if (gbd ->Gamut[theta][alpha].Type != GP_EMPTY) return TRUE;
487 
488     // Fill close points
489     nCloseSectors = FindNearSectors(gbd, alpha, theta, Close);
490 
491 
492     // Find a central point on the sector
493     sp.alpha = (cmsFloat64Number) ((alpha + 0.5) * 360.0) / (SECTORS);
494     sp.theta = (cmsFloat64Number) ((theta + 0.5) * 180.0) / (SECTORS);
495     sp.r     = 50.0;
496 
497     // Convert to Cartesian
498     ToCartesian(&Lab, &sp);
499 
500     // Create a ray line from centre to this point
501     _cmsVEC3init(&Centre, 50.0, 0, 0);
502     LineOf2Points(&ray, &Lab, &Centre);
503 
504     // For all close sectors
505     closel.r = 0.0;
506     closel.alpha = 0;
507     closel.theta = 0;
508 
509     for (k=0; k < nCloseSectors; k++) {
510 
511         for(m = k+1; m < nCloseSectors; m++) {
512 
513             cmsVEC3 temp, a1, a2;
514 
515             // A line from sector to sector
516             ToCartesian(&a1, &Close[k]->p);
517             ToCartesian(&a2, &Close[m]->p);
518 
519             LineOf2Points(&edge, &a1, &a2);
520 
521             // Find a line
522             ClosestLineToLine(&temp, &ray, &edge);
523 
524             // Convert to spherical
525             ToSpherical(&templ, &temp);
526 
527 
528             if ( templ.r > closel.r &&
529                  templ.theta >= (theta*180.0/SECTORS) &&
530                  templ.theta <= ((theta+1)*180.0/SECTORS) &&
531                  templ.alpha >= (alpha*360.0/SECTORS) &&
532                  templ.alpha <= ((alpha+1)*360.0/SECTORS)) {
533 
534                 closel = templ;
535             }
536         }
537     }
538 
539     gbd ->Gamut[theta][alpha].p = closel;
540     gbd ->Gamut[theta][alpha].Type = GP_MODELED;
541 
542     return TRUE;
543 
544 }
545 
546 
547 // Interpolate missing parts. The algorithm fist computes slices at
548 // theta=0 and theta=Max.
cmsGDBCompute(cmsHANDLE hGBD,cmsUInt32Number dwFlags)549 cmsBool CMSEXPORT cmsGDBCompute(cmsHANDLE hGBD, cmsUInt32Number dwFlags)
550 {
551     int alpha, theta;
552     cmsGDB* gbd = (cmsGDB*) hGBD;
553 
554     _cmsAssert(hGBD != NULL);
555 
556     // Interpolate black
557     for (alpha = 0; alpha < SECTORS; alpha++) {
558 
559         if (!InterpolateMissingSector(gbd, alpha, 0)) return FALSE;
560     }
561 
562     // Interpolate white
563     for (alpha = 0; alpha < SECTORS; alpha++) {
564 
565         if (!InterpolateMissingSector(gbd, alpha, SECTORS-1)) return FALSE;
566     }
567 
568 
569     // Interpolate Mid
570     for (theta = 1; theta < SECTORS; theta++) {
571         for (alpha = 0; alpha < SECTORS; alpha++) {
572 
573             if (!InterpolateMissingSector(gbd, alpha, theta)) return FALSE;
574         }
575     }
576 
577     // Done
578     return TRUE;
579 
580     cmsUNUSED_PARAMETER(dwFlags);
581 }
582 
583 
584 
585 
586 // --------------------------------------------------------------------------------------------------------
587 
588 // Great for debug, but not suitable for real use
589 
590 #if 0
591 cmsBool cmsGBDdumpVRML(cmsHANDLE hGBD, const char* fname)
592 {
593     FILE* fp;
594     int   i, j;
595     cmsGDB* gbd = (cmsGDB*) hGBD;
596     cmsGDBPoint* pt;
597 
598     fp = fopen (fname, "wt");
599     if (fp == NULL)
600         return FALSE;
601 
602     fprintf (fp, "#VRML V2.0 utf8\n");
603 
604     // set the viewing orientation and distance
605     fprintf (fp, "DEF CamTest Group {\n");
606     fprintf (fp, "\tchildren [\n");
607     fprintf (fp, "\t\tDEF Cameras Group {\n");
608     fprintf (fp, "\t\t\tchildren [\n");
609     fprintf (fp, "\t\t\t\tDEF DefaultView Viewpoint {\n");
610     fprintf (fp, "\t\t\t\t\tposition 0 0 340\n");
611     fprintf (fp, "\t\t\t\t\torientation 0 0 1 0\n");
612     fprintf (fp, "\t\t\t\t\tdescription \"default view\"\n");
613     fprintf (fp, "\t\t\t\t}\n");
614     fprintf (fp, "\t\t\t]\n");
615     fprintf (fp, "\t\t},\n");
616     fprintf (fp, "\t]\n");
617     fprintf (fp, "}\n");
618 
619     // Output the background stuff
620     fprintf (fp, "Background {\n");
621     fprintf (fp, "\tskyColor [\n");
622     fprintf (fp, "\t\t.5 .5 .5\n");
623     fprintf (fp, "\t]\n");
624     fprintf (fp, "}\n");
625 
626     // Output the shape stuff
627     fprintf (fp, "Transform {\n");
628     fprintf (fp, "\tscale .3 .3 .3\n");
629     fprintf (fp, "\tchildren [\n");
630 
631     // Draw the axes as a shape:
632     fprintf (fp, "\t\tShape {\n");
633     fprintf (fp, "\t\t\tappearance Appearance {\n");
634     fprintf (fp, "\t\t\t\tmaterial Material {\n");
635     fprintf (fp, "\t\t\t\t\tdiffuseColor 0 0.8 0\n");
636     fprintf (fp, "\t\t\t\t\temissiveColor 1.0 1.0 1.0\n");
637     fprintf (fp, "\t\t\t\t\tshininess 0.8\n");
638     fprintf (fp, "\t\t\t\t}\n");
639     fprintf (fp, "\t\t\t}\n");
640     fprintf (fp, "\t\t\tgeometry IndexedLineSet {\n");
641     fprintf (fp, "\t\t\t\tcoord Coordinate {\n");
642     fprintf (fp, "\t\t\t\t\tpoint [\n");
643     fprintf (fp, "\t\t\t\t\t0.0 0.0 0.0,\n");
644     fprintf (fp, "\t\t\t\t\t%f 0.0 0.0,\n",  255.0);
645     fprintf (fp, "\t\t\t\t\t0.0 %f 0.0,\n",  255.0);
646     fprintf (fp, "\t\t\t\t\t0.0 0.0 %f]\n",  255.0);
647     fprintf (fp, "\t\t\t\t}\n");
648     fprintf (fp, "\t\t\t\tcoordIndex [\n");
649     fprintf (fp, "\t\t\t\t\t0, 1, -1\n");
650     fprintf (fp, "\t\t\t\t\t0, 2, -1\n");
651     fprintf (fp, "\t\t\t\t\t0, 3, -1]\n");
652     fprintf (fp, "\t\t\t}\n");
653     fprintf (fp, "\t\t}\n");
654 
655 
656     fprintf (fp, "\t\tShape {\n");
657     fprintf (fp, "\t\t\tappearance Appearance {\n");
658     fprintf (fp, "\t\t\t\tmaterial Material {\n");
659     fprintf (fp, "\t\t\t\t\tdiffuseColor 0 0.8 0\n");
660     fprintf (fp, "\t\t\t\t\temissiveColor 1 1 1\n");
661     fprintf (fp, "\t\t\t\t\tshininess 0.8\n");
662     fprintf (fp, "\t\t\t\t}\n");
663     fprintf (fp, "\t\t\t}\n");
664     fprintf (fp, "\t\t\tgeometry PointSet {\n");
665 
666     // fill in the points here
667     fprintf (fp, "\t\t\t\tcoord Coordinate {\n");
668     fprintf (fp, "\t\t\t\t\tpoint [\n");
669 
670     // We need to transverse all gamut hull.
671     for (i=0; i < SECTORS; i++)
672         for (j=0; j < SECTORS; j++) {
673 
674             cmsVEC3 v;
675 
676             pt = &gbd ->Gamut[i][j];
677             ToCartesian(&v, &pt ->p);
678 
679             fprintf (fp, "\t\t\t\t\t%g %g %g", v.n[0]+50, v.n[1], v.n[2]);
680 
681             if ((j == SECTORS - 1) && (i == SECTORS - 1))
682                 fprintf (fp, "]\n");
683             else
684                 fprintf (fp, ",\n");
685 
686         }
687 
688         fprintf (fp, "\t\t\t\t}\n");
689 
690 
691 
692     // fill in the face colors
693     fprintf (fp, "\t\t\t\tcolor Color {\n");
694     fprintf (fp, "\t\t\t\t\tcolor [\n");
695 
696     for (i=0; i < SECTORS; i++)
697         for (j=0; j < SECTORS; j++) {
698 
699            cmsVEC3 v;
700 
701             pt = &gbd ->Gamut[i][j];
702 
703 
704             ToCartesian(&v, &pt ->p);
705 
706 
707         if (pt ->Type == GP_EMPTY)
708             fprintf (fp, "\t\t\t\t\t%g %g %g", 0.0, 0.0, 0.0);
709         else
710             if (pt ->Type == GP_MODELED)
711                 fprintf (fp, "\t\t\t\t\t%g %g %g", 1.0, .5, .5);
712             else {
713                 fprintf (fp, "\t\t\t\t\t%g %g %g", 1.0, 1.0, 1.0);
714 
715             }
716 
717         if ((j == SECTORS - 1) && (i == SECTORS - 1))
718                 fprintf (fp, "]\n");
719             else
720                 fprintf (fp, ",\n");
721     }
722     fprintf (fp, "\t\t\t}\n");
723 
724 
725     fprintf (fp, "\t\t\t}\n");
726     fprintf (fp, "\t\t}\n");
727     fprintf (fp, "\t]\n");
728     fprintf (fp, "}\n");
729 
730     fclose (fp);
731 
732     return TRUE;
733 }
734 #endif
735