[gnome-color-manager] Add gcm_profile_generate_gamut_hull() to generate a 3D gamut hull of any 3 color profile

[Richard Hughes <rhughes src gnome org>]

commit bfd8017c356f4cc75cb19b67e69ff50f54f07afc
Author: Richard Hughes <richard hughsie com>
Date:   Sat Mar 26 18:21:58 2011 +0000

    Add gcm_profile_generate_gamut_hull() to generate a 3D gamut hull of any 3 color profile
    
    The original code was taken from icc_examin copyright 2004-2009 Kai-Uwe Behrmann <ku b gmx de>,
    although quite modified from how it was found. icc_examin is also GPLv2+ licenced.

 src/gcm-profile.c   |  428 +++++++++++++++++++++++++++++++++++++++++++++++++++
 src/gcm-profile.h   |    3 +
 src/gcm-self-test.c |   21 +++
 3 files changed, 452 insertions(+), 0 deletions(-)
---
diff --git a/src/gcm-profile.c b/src/gcm-profile.c
index 1e3edb2..6bfd8dc 100644
--- a/src/gcm-profile.c
+++ b/src/gcm-profile.c
@@ -38,6 +38,7 @@
 
 #include "gcm-profile.h"
 #include "gcm-color.h"
+#include "gcm-hull.h"
 
 static void     gcm_profile_finalize	(GObject     *object);
 
@@ -99,6 +100,8 @@ enum {
 
 G_DEFINE_TYPE (GcmProfile, gcm_profile, G_TYPE_OBJECT)
 
+#define HYP(a,b)		(sqrt((a)*(a) + (b)*(b)))
+
 static void gcm_profile_file_monitor_changed_cb (GFileMonitor *monitor, GFile *file, GFile *other_file, GFileMonitorEvent event_type, GcmProfile *profile);
 
 /**
@@ -1472,6 +1475,431 @@ out:
 }
 
 /**
+ * gcm_profile_create_lab_cube:
+ *
+ * The original code was taken from icc_examin,
+ *  Copyright 2004-2009 Kai-Uwe Behrmann <ku b gmx de>
+ **/
+static gdouble *
+gcm_profile_create_lab_cube (guint res)
+{
+	gdouble *lab = NULL;
+	gdouble max = 0.99;
+	gdouble min = 0.01;
+	gint area;
+	gint channels_n = 3;
+	gint pos;
+	gsize size;
+	guint x, y;
+
+	size = 4 * res * (res+1) + 2 * (res-1) * (res-1);
+	lab = g_new0 (gdouble, size * channels_n);
+	if (lab == NULL)
+		goto out;
+
+	g_debug ("created 2*%ix%i array", (guint)size, (guint)channels_n);
+
+	/* side squares */
+	for (y = 0; y <= res; ++y) {
+		for (x = 0; x < 4 * res; ++x) {
+			area = 0;
+			pos = (y * 4 * res + x) * channels_n;
+
+			lab[pos + 0] = pow(0.9999 - (gdouble)y / (gdouble)res, 2.0) + 0.0001;
+			if (area * res <= x && x < ++area * res) {
+				lab[pos + 1] = min + (x - (area - 1) * res) / (gdouble)res * (max-min);
+				lab[pos + 2] = min;
+			} else if (area * res <= x && x < ++area * res) {
+				lab[pos + 1] = max;
+				lab[pos + 2] = min + (x - (area - 1) * res) / (gdouble)res * (max-min);
+			} else if (area * res <= x && x < ++area * res) {
+				lab[pos + 1] = max - (x - (area - 1) * res) / (gdouble)res * (max-min);
+				lab[pos + 2] = max;
+			} else if (area * res <= x && x < ++area * res) {
+				lab[pos + 1] = min;
+				lab[pos + 2] = max - (x - (area - 1) * res) / (double)res * (max-min);
+			}
+		}
+	}
+
+	/* bottom and top square */
+	for (y = 0; y < (res - 1); ++y) {
+		for (x = 0; x < 2 * (res - 1); ++x) {
+			gint x_pos;
+			gint y_pos;
+			gdouble val;
+
+			pos = (4 * res * (res + 1) + y * 2 * (res - 1) + x) * channels_n;
+			area = 1;
+			x_pos = x + 1;
+			y_pos = y + 1;
+			val = (gdouble)y_pos/(gdouble)res * (max-min);
+
+			if (/*0 <= x &&*/ x < res - 1) {
+				lab[pos + 0] = 1.0;
+				lab[pos + 1] = min + (x_pos - (area - 1) * (res - 1)) / (gdouble)res * (max-min);
+				lab[pos + 2] = min + val;
+			} else if (res - 1 <= x && x < 2 * res - 2) {
+				++area;
+				lab[pos + 1] = min + (x_pos - (area - 1) * (res - 1)) / (gdouble)res * (max-min);
+				lab[pos + 2] = min + val;
+				lab[pos + 0] = HYP (lab[pos + 1] - 0.5, lab[pos + 2] - 0.5)/100.; /* 0.0 */
+			}
+		}
+	}
+out:
+	return lab;
+}
+
+/**
+ * gcm_profile_create_hull_for_data:
+ *
+ * The original code was taken from icc_examin,
+ *  Copyright 2004-2009 Kai-Uwe Behrmann <ku b gmx de>
+ **/
+static GcmHull *
+gcm_profile_create_hull_for_data (guint res, gdouble *lab, gdouble *rgb)
+{
+	GcmColorRGB color;
+	GcmColorXYZ xyz;
+	GcmHull *hull = NULL;
+	gint channels_n = 3;
+	gint off;
+	gsize i;
+	gsize size;
+	guint face[3];
+	guint x, y;
+
+	size = 4 * res * (res+1) + 2 * (res-1) * (res-1);
+
+	hull = gcm_hull_new ();
+
+	/* collect colour points */
+	for (i = 0; i < size; ++i) {
+		xyz.X = lab[i*channels_n+0];
+		xyz.Y = lab[i*channels_n+1];
+		xyz.Z = lab[i*channels_n+2];
+		color.R = rgb[i*channels_n+0];
+		color.G = rgb[i*channels_n+1];
+		color.B = rgb[i*channels_n+2];
+		gcm_hull_add_vertex (hull, &xyz, &color);
+	}
+
+	for (y = 0; y < res; ++y) {
+		for (x = 0; x < 4 * res; ++x) {
+			gint x_ = x;
+			if (x == 4 * res - 1)
+				x_ = -1;
+			face[0] = y * 4*res+x;
+			face[1] = y * 4*res+x_+1;
+			face[2] = (y+1)*4*res+x;
+			gcm_hull_add_face (hull, face, 3);
+
+			face[0] = y * 4*res+x_+1;
+			face[1] = (y+1)*4*res+x_+1;
+			face[2] = (y+1)*4*res+x;
+			gcm_hull_add_face (hull, face, 3);
+		}
+	}
+
+	off = 4 * res * (res + 1);
+
+	/* 1 0 0 (L res b) */
+	face[0] = 4*res-1;
+	face[1] = off;
+	face[2] = 0;
+	gcm_hull_add_face (hull, face, 3);
+
+	face[0] = off;
+	face[2] = 0;
+	face[1] = 1;
+	gcm_hull_add_face (hull, face, 3);
+
+	/* 0 0 0 */
+	face[1] = off-1;
+	face[0] = off+res-1;
+	face[2] = off-4*res;
+	gcm_hull_add_face (hull, face, 3);
+
+	face[1] = off+res-1;
+	face[2] = off-4*res;
+	face[0] = off - 4*res+1;
+	gcm_hull_add_face (hull, face, 3);
+
+	/* 0 0 1 */
+	face[2] = off-res;
+	face[1] = off-res-1;
+	face[0] = off+2*(res-1)*(res-1)-res+1;
+	gcm_hull_add_face (hull, face, 3);
+
+	face[0] = off-res;
+	face[1] = off-res+1;
+	face[2] = off+2*(res-1)*(res-1)-res+1;
+	gcm_hull_add_face (hull, face, 3);
+
+	/* 0 1 1 */
+	face[0] = off-2*res+1;
+	face[2] = off-2*res;
+	face[1] = off+2*(res-1)*(res-1)-1;
+	gcm_hull_add_face (hull, face, 3);
+
+	face[1] = off-2*res;
+	face[2] = off+2*(res-1)*(res-1)-1;
+	face[0] = off-2*res-1;
+	gcm_hull_add_face (hull, face, 3);
+
+	/* 1 1 1 */
+	face[0] = 2*res-1;
+	face[2] = 2*res;
+	face[1] = off+2*(res-1)*(res-1)-res;
+	gcm_hull_add_face (hull, face, 3);
+
+	face[1] = 2*res;
+	face[2] = off+2*(res-1)*(res-1)-res;
+	face[0] = 2*res+1;
+	gcm_hull_add_face (hull, face, 3);
+
+	/* 1 0 1 */
+	face[2] = 3*res;
+	face[0] = 3*res-1;
+	face[1] = off+2*(res-1)*(res-1)-2*res+2;
+	gcm_hull_add_face (hull, face, 3);
+
+	face[2] = 3*res;
+	face[1] = 3*res+1;
+	face[0] = off+2*(res-1)*(res-1)-2*res+2;
+	gcm_hull_add_face (hull, face, 3);
+
+	/* 1 1 0 */
+	face[0] = off+res-2;
+	face[1] = res + 1;
+	face[2] = res - 1;
+	gcm_hull_add_face (hull, face, 3);
+
+	face[0] = res + 1;
+	face[2] = res - 1;
+	face[1] = res;
+	gcm_hull_add_face (hull, face, 3);
+
+	/* 0 1 0 */
+	face[0] = off+2*(res-1)-1;
+	face[1] = off-3*res-1;
+	face[2] = off-3*res;
+	gcm_hull_add_face (hull, face, 3);
+
+	face[1] = off+2*(res-1)-1;
+	face[0] = off-3*res+1;
+	face[2] = off-3*res+0;
+	gcm_hull_add_face (hull, face, 3);
+
+	for (y = 0; y < res; ++y) {
+		if (0 < y && y < res - 1) {
+			/* 0 0 . */
+			face[2] = off-y;
+			face[0] = off+(y+1)*2*(res-1)-res+1;
+			face[1] = off-y-1;
+			gcm_hull_add_face (hull, face, 3);
+
+			face[0] = off+(y+0)*2*(res-1)-res+1;
+			face[2] = off-y;
+			face[1] = off+(y+1)*2*(res-1)-res+1;
+			gcm_hull_add_face (hull, face, 3);
+
+			/* 0 1 . */
+			face[1] = off+(y+1)*2*(res-1)-1;
+			face[0] = off-3*res+y+1;
+			face[2] = off+(y)*2*(res-1)-1;
+			gcm_hull_add_face (hull, face, 3);
+
+			face[1] = off-3*res+y+1;
+			face[2] = off+(y)*2*(res-1)-1;
+			face[0] = off-3*res+y;
+			gcm_hull_add_face (hull, face, 3);
+
+			/* 1 0 . */
+			face[0] = off+2*(res-1)*(res-1)-(y+1)*2*(res-1);
+			face[1] = 3*res+y+1;
+			face[2] = off+2*(res-1)*(res-1)-y*2*(res-1);
+			gcm_hull_add_face (hull, face, 3);
+
+			face[0] = 3*res+y+1;
+			face[2] = off+2*(res-1)*(res-1)-y*2*(res-1);
+			face[1] = 3*res+y;
+			gcm_hull_add_face (hull, face, 3);
+
+			/* 1 1 . */
+			face[0] = off+2*(res-1)*(res-1)-(y+1)*2*(res-1)+res-2;
+			face[1] = off+2*(res-1)*(res-1)-(y+0)*2*(res-1)+res-2;
+			face[2] = 2*res-y;
+			gcm_hull_add_face (hull, face, 3);
+
+			face[0] = 2*res-y-1;
+			face[1] = off+2*(res-1)*(res-1)-(y+1)*2*(res-1)+res-2;
+			face[2] = 2*res-y;
+			gcm_hull_add_face (hull, face, 3);
+		}
+
+		for (x = 0; x < 2 * res; ++x) {
+			gint x_ = x + off;
+
+			/* lower border */
+			if ( y == 0 ) {
+				if (x == 0) {
+				} else if (x == res - 1) {
+				} else if (x < res - 1) {
+					/* 1 . 0 */
+					face[0] = off + x - 1;
+					face[1] = off + x;
+					face[2] = x;
+					gcm_hull_add_face (hull, face, 3);
+
+					face[0] = off + x;
+					face[2] = x;
+					face[1] = x + 1;
+					gcm_hull_add_face (hull, face, 3);
+
+					/* 0 . 1 */
+					face[0] = off-res-x;
+					face[2] = off-res-x-1;
+					face[1] = off+2*(res-1)*(res-1)-res+x;
+					gcm_hull_add_face (hull, face, 3);
+
+					face[2] = off-res-x-1;
+					face[0] = off+2*(res-1)*(res-1)-res+x;
+					face[1] = off+2*(res-1)*(res-1)-res+x+1;
+					gcm_hull_add_face (hull, face, 3);
+
+					/* 1 . 1 */
+					face[0] = 3*res - x;
+					face[1] = 3*res - x-1;
+					face[2] = off+2*(res-1)*(res-1)-2*(res-1)+x-1;
+					gcm_hull_add_face (hull, face, 3);
+
+					face[0] = 3*res - x-1;
+					face[2] = off+2*(res-1)*(res-1)-2*(res-1)+x-1;
+					face[1] = off+2*(res-1)*(res-1)-2*(res-1)+x;
+					gcm_hull_add_face (hull, face, 3);
+
+				} else if (x > res + 1) {
+					/* 0 . 0 */
+					face[0] = off+x-3;
+					face[2] = off+x-3+1;
+					face[1] = 4*res*(res+1)-4*res + x-res-1;
+					gcm_hull_add_face (hull, face, 3);
+
+					face[1] = off+x-3+1;
+					face[2] = 4*res*(res+1)-4*res + x-res-1;
+					face[0] = 4*res*(res+1)-4*res + x-res;
+					gcm_hull_add_face (hull, face, 3);
+				}
+
+			/* upper border */
+			} else if ( y == res - 1 ) {
+				if (x == 0) {
+				}
+			} else if (/*0 <= x &&*/ x < res - 1 - 1) {
+
+				/* upper middle field (*L=0.0) */
+				face[0] = (y-1) * 2*(res-1)+x_;
+				face[2] = (y-1)*2*(res-1)+x_+1;
+				face[1] = (y+0)*2*(res-1)+x_;
+				gcm_hull_add_face (hull, face, 3);
+
+				face[2] = (y-1)*2*(res-1)+x_+1;
+				face[0] = (y+0)*2*(res-1)+x_;
+				face[1] = (y+0)*2*(res-1)+x_+1;
+				gcm_hull_add_face (hull, face, 3);
+
+			} else if (res - 1 <= x && x < 2 * res - 2 - 1) {
+
+				/* lower middle field (*L=1.0) */
+				face[0] = (y-1) * 2*(res-1)+x_;
+				face[1] = (y-1)*2*(res-1)+x_+1;
+				face[2] = (y+0)*2*(res-1)+x_;
+				gcm_hull_add_face (hull, face, 3);
+
+				face[0] = (y-1)*2*(res-1)+x_+1;
+				face[2] = (y+0)*2*(res-1)+x_;
+				face[1] = (y+0)*2*(res-1)+x_+1;
+				gcm_hull_add_face (hull, face, 3);
+			}
+		}
+	}
+
+	return hull;
+}
+
+/**
+ * gcm_profile_generate_gamut_hull:
+ * @profile: a #GcmProfile
+ * @res: The resolution. 10 is quick, 20 is more precise. 12 is a good default.
+ *
+ * A cube from six squares with the range of the Lab cube will be
+ * transformed to a profile colour space and then converted to a
+ * mesh.
+ *
+ * The original code was taken from icc_examin,
+ *  Copyright 2004-2009 Kai-Uwe Behrmann <ku b gmx de>
+ **/
+GcmHull *
+gcm_profile_generate_gamut_hull (GcmProfile *profile, guint res)
+{
+	cmsHPROFILE lab_profile = NULL;
+	cmsHPROFILE srgb_profile = NULL;
+	cmsHTRANSFORM lab_transform = NULL;
+	cmsHTRANSFORM srgb_transform = NULL;
+	GcmHull *hull = NULL;
+	GcmProfilePrivate *priv = profile->priv;
+	gdouble *lab = NULL;
+	gdouble *rgb = NULL;
+	gint channels_n = 3;
+	gsize size = 4 * res * (res+1) + 2 * (res-1) * (res-1);
+
+	/* create data array */
+	lab = gcm_profile_create_lab_cube (res);
+	rgb = g_new0 (gdouble, size * channels_n);
+	if (rgb == NULL)
+		goto out;
+
+	/* run the cube through the Lab profile */
+	lab_profile = cmsCreateLab4Profile (cmsD50_xyY ());
+	lab_transform = cmsCreateTransform (priv->lcms_profile, TYPE_RGB_DBL,
+					    lab_profile, TYPE_Lab_DBL,
+					    INTENT_ABSOLUTE_COLORIMETRIC, 0);
+	if (lab_transform == NULL) {
+		g_warning ("failed to create Lab transform");
+		goto out;
+	}
+	cmsDoTransform (lab_transform, lab, lab, size);
+
+	/* run the cube through the sRGB profile */
+	srgb_profile = cmsCreate_sRGBProfile ();
+	srgb_transform = cmsCreateTransform (lab_profile, TYPE_Lab_DBL,
+					     srgb_profile, TYPE_RGB_DBL,
+					     INTENT_PERCEPTUAL, 0);
+	if (srgb_transform == NULL) {
+		g_warning ("failed to create sRGB transform");
+		goto out;
+	}
+	cmsDoTransform (srgb_transform, lab, rgb, size);
+
+	/* create gamut hull */
+	hull = gcm_profile_create_hull_for_data (res, lab, rgb);
+out:
+	g_free (rgb);
+	g_free (lab);
+	if (lab_profile != NULL)
+		cmsCloseProfile (lab_profile);
+	if (srgb_profile != NULL)
+		cmsCloseProfile (srgb_profile);
+	if (lab_transform != NULL)
+		cmsDeleteTransform (lab_transform);
+	if (srgb_transform != NULL)
+		cmsDeleteTransform (srgb_transform);
+	return hull;
+}
+
+/**
  * gcm_profile_file_monitor_changed_cb:
  **/
 static void
diff --git a/src/gcm-profile.h b/src/gcm-profile.h
index 11e6a7f..d01073e 100644
--- a/src/gcm-profile.h
+++ b/src/gcm-profile.h
@@ -28,6 +28,7 @@
 
 #include "gcm-clut.h"
 #include "gcm-color.h"
+#include "gcm-hull.h"
 
 G_BEGIN_DECLS
 
@@ -79,6 +80,8 @@ GcmClut		*gcm_profile_generate_vcgt		(GcmProfile	*profile,
 							 guint		 size);
 GcmClut		*gcm_profile_generate_curve		(GcmProfile	*profile,
 							 guint		 size);
+GcmHull		*gcm_profile_generate_gamut_hull	(GcmProfile	*profile,
+							 guint		 res);
 gboolean	 gcm_profile_set_vcgt_from_data		(GcmProfile	*profile,
 							 guint16	*red,
 							 guint16	*green,
diff --git a/src/gcm-self-test.c b/src/gcm-self-test.c
index 7ff812f..bbf3577 100644
--- a/src/gcm-self-test.c
+++ b/src/gcm-self-test.c
@@ -372,6 +372,8 @@ gcm_test_profile_func (void)
 	GcmColorYxy green;
 	GcmColorYxy blue;
 	GcmColorYxy white;
+	GcmHull *hull;
+	gchar *data;
 
 	/* bluish test */
 	profile = gcm_profile_new ();
@@ -475,6 +477,25 @@ gcm_test_profile_func (void)
 
 	g_object_unref (file);
 	g_object_unref (profile);
+
+	/* get gamut hull */
+	profile = gcm_profile_new ();
+	file = g_file_new_for_path (TESTDATADIR "/ibm-t61.icc");
+	ret = gcm_profile_parse (profile, file, &error);
+	g_assert_no_error (error);
+	g_assert (ret);
+	hull = gcm_profile_generate_gamut_hull (profile, 12);
+	g_assert (hull != NULL);
+
+	/* save as PLY file */
+	data = gcm_hull_export_to_ply (hull);
+	ret = g_file_set_contents ("/tmp/gamut.ply", data, -1, NULL);
+	g_assert (ret);
+
+	g_free (data);
+	g_object_unref (hull);
+	g_object_unref (file);
+	g_object_unref (profile);
 }
 
 static void