1// The colorful package provides all kinds of functions for working with colors.
2package colorful
3
4import (
5	"fmt"
6	"image/color"
7	"math"
8)
9
10// A color is stored internally using sRGB (standard RGB) values in the range 0-1
11type Color struct {
12	R, G, B float64
13}
14
15// Implement the Go color.Color interface.
16func (col Color) RGBA() (r, g, b, a uint32) {
17	r = uint32(col.R*65535.0 + 0.5)
18	g = uint32(col.G*65535.0 + 0.5)
19	b = uint32(col.B*65535.0 + 0.5)
20	a = 0xFFFF
21	return
22}
23
24// Constructs a colorful.Color from something implementing color.Color
25func MakeColor(col color.Color) (Color, bool) {
26	r, g, b, a := col.RGBA()
27	if a == 0 {
28		return Color{0, 0, 0}, false
29	}
30
31	// Since color.Color is alpha pre-multiplied, we need to divide the
32	// RGB values by alpha again in order to get back the original RGB.
33	r *= 0xffff
34	r /= a
35	g *= 0xffff
36	g /= a
37	b *= 0xffff
38	b /= a
39
40	return Color{float64(r) / 65535.0, float64(g) / 65535.0, float64(b) / 65535.0}, true
41}
42
43// Might come in handy sometimes to reduce boilerplate code.
44func (col Color) RGB255() (r, g, b uint8) {
45	r = uint8(col.R*255.0 + 0.5)
46	g = uint8(col.G*255.0 + 0.5)
47	b = uint8(col.B*255.0 + 0.5)
48	return
49}
50
51// This is the tolerance used when comparing colors using AlmostEqualRgb.
52const Delta = 1.0 / 255.0
53
54// This is the default reference white point.
55var D65 = [3]float64{0.95047, 1.00000, 1.08883}
56
57// And another one.
58var D50 = [3]float64{0.96422, 1.00000, 0.82521}
59
60// Checks whether the color exists in RGB space, i.e. all values are in [0..1]
61func (c Color) IsValid() bool {
62	return 0.0 <= c.R && c.R <= 1.0 &&
63		0.0 <= c.G && c.G <= 1.0 &&
64		0.0 <= c.B && c.B <= 1.0
65}
66
67func clamp01(v float64) float64 {
68	return math.Max(0.0, math.Min(v, 1.0))
69}
70
71// Returns Clamps the color into valid range, clamping each value to [0..1]
72// If the color is valid already, this is a no-op.
73func (c Color) Clamped() Color {
74	return Color{clamp01(c.R), clamp01(c.G), clamp01(c.B)}
75}
76
77func sq(v float64) float64 {
78	return v * v
79}
80
81func cub(v float64) float64 {
82	return v * v * v
83}
84
85// DistanceRgb computes the distance between two colors in RGB space.
86// This is not a good measure! Rather do it in Lab space.
87func (c1 Color) DistanceRgb(c2 Color) float64 {
88	return math.Sqrt(sq(c1.R-c2.R) + sq(c1.G-c2.G) + sq(c1.B-c2.B))
89}
90
91// Check for equality between colors within the tolerance Delta (1/255).
92func (c1 Color) AlmostEqualRgb(c2 Color) bool {
93	return math.Abs(c1.R-c2.R)+
94		math.Abs(c1.G-c2.G)+
95		math.Abs(c1.B-c2.B) < 3.0*Delta
96}
97
98// You don't really want to use this, do you? Go for BlendLab, BlendLuv or BlendHcl.
99func (c1 Color) BlendRgb(c2 Color, t float64) Color {
100	return Color{c1.R + t*(c2.R-c1.R),
101		c1.G + t*(c2.G-c1.G),
102		c1.B + t*(c2.B-c1.B)}
103}
104
105// Utility used by Hxx color-spaces for interpolating between two angles in [0,360].
106func interp_angle(a0, a1, t float64) float64 {
107	// Based on the answer here: http://stackoverflow.com/a/14498790/2366315
108	// With potential proof that it works here: http://math.stackexchange.com/a/2144499
109	delta := math.Mod(math.Mod(a1-a0, 360.0)+540, 360.0) - 180.0
110	return math.Mod(a0+t*delta+360.0, 360.0)
111}
112
113/// HSV ///
114///////////
115// From http://en.wikipedia.org/wiki/HSL_and_HSV
116// Note that h is in [0..360] and s,v in [0..1]
117
118// Hsv returns the Hue [0..360], Saturation and Value [0..1] of the color.
119func (col Color) Hsv() (h, s, v float64) {
120	min := math.Min(math.Min(col.R, col.G), col.B)
121	v = math.Max(math.Max(col.R, col.G), col.B)
122	C := v - min
123
124	s = 0.0
125	if v != 0.0 {
126		s = C / v
127	}
128
129	h = 0.0 // We use 0 instead of undefined as in wp.
130	if min != v {
131		if v == col.R {
132			h = math.Mod((col.G-col.B)/C, 6.0)
133		}
134		if v == col.G {
135			h = (col.B-col.R)/C + 2.0
136		}
137		if v == col.B {
138			h = (col.R-col.G)/C + 4.0
139		}
140		h *= 60.0
141		if h < 0.0 {
142			h += 360.0
143		}
144	}
145	return
146}
147
148// Hsv creates a new Color given a Hue in [0..360], a Saturation and a Value in [0..1]
149func Hsv(H, S, V float64) Color {
150	Hp := H / 60.0
151	C := V * S
152	X := C * (1.0 - math.Abs(math.Mod(Hp, 2.0)-1.0))
153
154	m := V - C
155	r, g, b := 0.0, 0.0, 0.0
156
157	switch {
158	case 0.0 <= Hp && Hp < 1.0:
159		r = C
160		g = X
161	case 1.0 <= Hp && Hp < 2.0:
162		r = X
163		g = C
164	case 2.0 <= Hp && Hp < 3.0:
165		g = C
166		b = X
167	case 3.0 <= Hp && Hp < 4.0:
168		g = X
169		b = C
170	case 4.0 <= Hp && Hp < 5.0:
171		r = X
172		b = C
173	case 5.0 <= Hp && Hp < 6.0:
174		r = C
175		b = X
176	}
177
178	return Color{m + r, m + g, m + b}
179}
180
181// You don't really want to use this, do you? Go for BlendLab, BlendLuv or BlendHcl.
182func (c1 Color) BlendHsv(c2 Color, t float64) Color {
183	h1, s1, v1 := c1.Hsv()
184	h2, s2, v2 := c2.Hsv()
185
186	// We know that h are both in [0..360]
187	return Hsv(interp_angle(h1, h2, t), s1+t*(s2-s1), v1+t*(v2-v1))
188}
189
190/// HSL ///
191///////////
192
193// Hsl returns the Hue [0..360], Saturation [0..1], and Luminance (lightness) [0..1] of the color.
194func (col Color) Hsl() (h, s, l float64) {
195	min := math.Min(math.Min(col.R, col.G), col.B)
196	max := math.Max(math.Max(col.R, col.G), col.B)
197
198	l = (max + min) / 2
199
200	if min == max {
201		s = 0
202		h = 0
203	} else {
204		if l < 0.5 {
205			s = (max - min) / (max + min)
206		} else {
207			s = (max - min) / (2.0 - max - min)
208		}
209
210		if max == col.R {
211			h = (col.G - col.B) / (max - min)
212		} else if max == col.G {
213			h = 2.0 + (col.B-col.R)/(max-min)
214		} else {
215			h = 4.0 + (col.R-col.G)/(max-min)
216		}
217
218		h *= 60
219
220		if h < 0 {
221			h += 360
222		}
223	}
224
225	return
226}
227
228// Hsl creates a new Color given a Hue in [0..360], a Saturation [0..1], and a Luminance (lightness) in [0..1]
229func Hsl(h, s, l float64) Color {
230	if s == 0 {
231		return Color{l, l, l}
232	}
233
234	var r, g, b float64
235	var t1 float64
236	var t2 float64
237	var tr float64
238	var tg float64
239	var tb float64
240
241	if l < 0.5 {
242		t1 = l * (1.0 + s)
243	} else {
244		t1 = l + s - l*s
245	}
246
247	t2 = 2*l - t1
248	h /= 360
249	tr = h + 1.0/3.0
250	tg = h
251	tb = h - 1.0/3.0
252
253	if tr < 0 {
254		tr++
255	}
256	if tr > 1 {
257		tr--
258	}
259	if tg < 0 {
260		tg++
261	}
262	if tg > 1 {
263		tg--
264	}
265	if tb < 0 {
266		tb++
267	}
268	if tb > 1 {
269		tb--
270	}
271
272	// Red
273	if 6*tr < 1 {
274		r = t2 + (t1-t2)*6*tr
275	} else if 2*tr < 1 {
276		r = t1
277	} else if 3*tr < 2 {
278		r = t2 + (t1-t2)*(2.0/3.0-tr)*6
279	} else {
280		r = t2
281	}
282
283	// Green
284	if 6*tg < 1 {
285		g = t2 + (t1-t2)*6*tg
286	} else if 2*tg < 1 {
287		g = t1
288	} else if 3*tg < 2 {
289		g = t2 + (t1-t2)*(2.0/3.0-tg)*6
290	} else {
291		g = t2
292	}
293
294	// Blue
295	if 6*tb < 1 {
296		b = t2 + (t1-t2)*6*tb
297	} else if 2*tb < 1 {
298		b = t1
299	} else if 3*tb < 2 {
300		b = t2 + (t1-t2)*(2.0/3.0-tb)*6
301	} else {
302		b = t2
303	}
304
305	return Color{r, g, b}
306}
307
308/// Hex ///
309///////////
310
311// Hex returns the hex "html" representation of the color, as in #ff0080.
312func (col Color) Hex() string {
313	// Add 0.5 for rounding
314	return fmt.Sprintf("#%02x%02x%02x", uint8(col.R*255.0+0.5), uint8(col.G*255.0+0.5), uint8(col.B*255.0+0.5))
315}
316
317// Hex parses a "html" hex color-string, either in the 3 "#f0c" or 6 "#ff1034" digits form.
318func Hex(scol string) (Color, error) {
319	format := "#%02x%02x%02x"
320	factor := 1.0 / 255.0
321	if len(scol) == 4 {
322		format = "#%1x%1x%1x"
323		factor = 1.0 / 15.0
324	}
325
326	var r, g, b uint8
327	n, err := fmt.Sscanf(scol, format, &r, &g, &b)
328	if err != nil {
329		return Color{}, err
330	}
331	if n != 3 {
332		return Color{}, fmt.Errorf("color: %v is not a hex-color", scol)
333	}
334
335	return Color{float64(r) * factor, float64(g) * factor, float64(b) * factor}, nil
336}
337
338/// Linear ///
339//////////////
340// http://www.sjbrown.co.uk/2004/05/14/gamma-correct-rendering/
341// http://www.brucelindbloom.com/Eqn_RGB_to_XYZ.html
342
343func linearize(v float64) float64 {
344	if v <= 0.04045 {
345		return v / 12.92
346	}
347	return math.Pow((v+0.055)/1.055, 2.4)
348}
349
350// LinearRgb converts the color into the linear RGB space (see http://www.sjbrown.co.uk/2004/05/14/gamma-correct-rendering/).
351func (col Color) LinearRgb() (r, g, b float64) {
352	r = linearize(col.R)
353	g = linearize(col.G)
354	b = linearize(col.B)
355	return
356}
357
358// A much faster and still quite precise linearization using a 6th-order Taylor approximation.
359// See the accompanying Jupyter notebook for derivation of the constants.
360func linearize_fast(v float64) float64 {
361	v1 := v - 0.5
362	v2 := v1 * v1
363	v3 := v2 * v1
364	v4 := v2 * v2
365	//v5 := v3*v2
366	return -0.248750514614486 + 0.925583310193438*v + 1.16740237321695*v2 + 0.280457026598666*v3 - 0.0757991963780179*v4 //+ 0.0437040411548932*v5
367}
368
369// FastLinearRgb is much faster than and almost as accurate as LinearRgb.
370// BUT it is important to NOTE that they only produce good results for valid colors r,g,b in [0,1].
371func (col Color) FastLinearRgb() (r, g, b float64) {
372	r = linearize_fast(col.R)
373	g = linearize_fast(col.G)
374	b = linearize_fast(col.B)
375	return
376}
377
378func delinearize(v float64) float64 {
379	if v <= 0.0031308 {
380		return 12.92 * v
381	}
382	return 1.055*math.Pow(v, 1.0/2.4) - 0.055
383}
384
385// LinearRgb creates an sRGB color out of the given linear RGB color (see http://www.sjbrown.co.uk/2004/05/14/gamma-correct-rendering/).
386func LinearRgb(r, g, b float64) Color {
387	return Color{delinearize(r), delinearize(g), delinearize(b)}
388}
389
390func delinearize_fast(v float64) float64 {
391	// This function (fractional root) is much harder to linearize, so we need to split.
392	if v > 0.2 {
393		v1 := v - 0.6
394		v2 := v1 * v1
395		v3 := v2 * v1
396		v4 := v2 * v2
397		v5 := v3 * v2
398		return 0.442430344268235 + 0.592178981271708*v - 0.287864782562636*v2 + 0.253214392068985*v3 - 0.272557158129811*v4 + 0.325554383321718*v5
399	} else if v > 0.03 {
400		v1 := v - 0.115
401		v2 := v1 * v1
402		v3 := v2 * v1
403		v4 := v2 * v2
404		v5 := v3 * v2
405		return 0.194915592891669 + 1.55227076330229*v - 3.93691860257828*v2 + 18.0679839248761*v3 - 101.468750302746*v4 + 632.341487393927*v5
406	} else {
407		v1 := v - 0.015
408		v2 := v1 * v1
409		v3 := v2 * v1
410		v4 := v2 * v2
411		v5 := v3 * v2
412		// You can clearly see from the involved constants that the low-end is highly nonlinear.
413		return 0.0519565234928877 + 5.09316778537561*v - 99.0338180489702*v2 + 3484.52322764895*v3 - 150028.083412663*v4 + 7168008.42971613*v5
414	}
415}
416
417// FastLinearRgb is much faster than and almost as accurate as LinearRgb.
418// BUT it is important to NOTE that they only produce good results for valid inputs r,g,b in [0,1].
419func FastLinearRgb(r, g, b float64) Color {
420	return Color{delinearize_fast(r), delinearize_fast(g), delinearize_fast(b)}
421}
422
423// XyzToLinearRgb converts from CIE XYZ-space to Linear RGB space.
424func XyzToLinearRgb(x, y, z float64) (r, g, b float64) {
425	r = 3.2404542*x - 1.5371385*y - 0.4985314*z
426	g = -0.9692660*x + 1.8760108*y + 0.0415560*z
427	b = 0.0556434*x - 0.2040259*y + 1.0572252*z
428	return
429}
430
431func LinearRgbToXyz(r, g, b float64) (x, y, z float64) {
432	x = 0.4124564*r + 0.3575761*g + 0.1804375*b
433	y = 0.2126729*r + 0.7151522*g + 0.0721750*b
434	z = 0.0193339*r + 0.1191920*g + 0.9503041*b
435	return
436}
437
438/// XYZ ///
439///////////
440// http://www.sjbrown.co.uk/2004/05/14/gamma-correct-rendering/
441
442func (col Color) Xyz() (x, y, z float64) {
443	return LinearRgbToXyz(col.LinearRgb())
444}
445
446func Xyz(x, y, z float64) Color {
447	return LinearRgb(XyzToLinearRgb(x, y, z))
448}
449
450/// xyY ///
451///////////
452// http://www.brucelindbloom.com/Eqn_XYZ_to_xyY.html
453
454// Well, the name is bad, since it's xyY but Golang needs me to start with a
455// capital letter to make the method public.
456func XyzToXyy(X, Y, Z float64) (x, y, Yout float64) {
457	return XyzToXyyWhiteRef(X, Y, Z, D65)
458}
459
460func XyzToXyyWhiteRef(X, Y, Z float64, wref [3]float64) (x, y, Yout float64) {
461	Yout = Y
462	N := X + Y + Z
463	if math.Abs(N) < 1e-14 {
464		// When we have black, Bruce Lindbloom recommends to use
465		// the reference white's chromacity for x and y.
466		x = wref[0] / (wref[0] + wref[1] + wref[2])
467		y = wref[1] / (wref[0] + wref[1] + wref[2])
468	} else {
469		x = X / N
470		y = Y / N
471	}
472	return
473}
474
475func XyyToXyz(x, y, Y float64) (X, Yout, Z float64) {
476	Yout = Y
477
478	if -1e-14 < y && y < 1e-14 {
479		X = 0.0
480		Z = 0.0
481	} else {
482		X = Y / y * x
483		Z = Y / y * (1.0 - x - y)
484	}
485
486	return
487}
488
489// Converts the given color to CIE xyY space using D65 as reference white.
490// (Note that the reference white is only used for black input.)
491// x, y and Y are in [0..1]
492func (col Color) Xyy() (x, y, Y float64) {
493	return XyzToXyy(col.Xyz())
494}
495
496// Converts the given color to CIE xyY space, taking into account
497// a given reference white. (i.e. the monitor's white)
498// (Note that the reference white is only used for black input.)
499// x, y and Y are in [0..1]
500func (col Color) XyyWhiteRef(wref [3]float64) (x, y, Y float64) {
501	X, Y2, Z := col.Xyz()
502	return XyzToXyyWhiteRef(X, Y2, Z, wref)
503}
504
505// Generates a color by using data given in CIE xyY space.
506// x, y and Y are in [0..1]
507func Xyy(x, y, Y float64) Color {
508	return Xyz(XyyToXyz(x, y, Y))
509}
510
511/// L*a*b* ///
512//////////////
513// http://en.wikipedia.org/wiki/Lab_color_space#CIELAB-CIEXYZ_conversions
514// For L*a*b*, we need to L*a*b*<->XYZ->RGB and the first one is device dependent.
515
516func lab_f(t float64) float64 {
517	if t > 6.0/29.0*6.0/29.0*6.0/29.0 {
518		return math.Cbrt(t)
519	}
520	return t/3.0*29.0/6.0*29.0/6.0 + 4.0/29.0
521}
522
523func XyzToLab(x, y, z float64) (l, a, b float64) {
524	// Use D65 white as reference point by default.
525	// http://www.fredmiranda.com/forum/topic/1035332
526	// http://en.wikipedia.org/wiki/Standard_illuminant
527	return XyzToLabWhiteRef(x, y, z, D65)
528}
529
530func XyzToLabWhiteRef(x, y, z float64, wref [3]float64) (l, a, b float64) {
531	fy := lab_f(y / wref[1])
532	l = 1.16*fy - 0.16
533	a = 5.0 * (lab_f(x/wref[0]) - fy)
534	b = 2.0 * (fy - lab_f(z/wref[2]))
535	return
536}
537
538func lab_finv(t float64) float64 {
539	if t > 6.0/29.0 {
540		return t * t * t
541	}
542	return 3.0 * 6.0 / 29.0 * 6.0 / 29.0 * (t - 4.0/29.0)
543}
544
545func LabToXyz(l, a, b float64) (x, y, z float64) {
546	// D65 white (see above).
547	return LabToXyzWhiteRef(l, a, b, D65)
548}
549
550func LabToXyzWhiteRef(l, a, b float64, wref [3]float64) (x, y, z float64) {
551	l2 := (l + 0.16) / 1.16
552	x = wref[0] * lab_finv(l2+a/5.0)
553	y = wref[1] * lab_finv(l2)
554	z = wref[2] * lab_finv(l2-b/2.0)
555	return
556}
557
558// Converts the given color to CIE L*a*b* space using D65 as reference white.
559func (col Color) Lab() (l, a, b float64) {
560	return XyzToLab(col.Xyz())
561}
562
563// Converts the given color to CIE L*a*b* space, taking into account
564// a given reference white. (i.e. the monitor's white)
565func (col Color) LabWhiteRef(wref [3]float64) (l, a, b float64) {
566	x, y, z := col.Xyz()
567	return XyzToLabWhiteRef(x, y, z, wref)
568}
569
570// Generates a color by using data given in CIE L*a*b* space using D65 as reference white.
571// WARNING: many combinations of `l`, `a`, and `b` values do not have corresponding
572//          valid RGB values, check the FAQ in the README if you're unsure.
573func Lab(l, a, b float64) Color {
574	return Xyz(LabToXyz(l, a, b))
575}
576
577// Generates a color by using data given in CIE L*a*b* space, taking
578// into account a given reference white. (i.e. the monitor's white)
579func LabWhiteRef(l, a, b float64, wref [3]float64) Color {
580	return Xyz(LabToXyzWhiteRef(l, a, b, wref))
581}
582
583// DistanceLab is a good measure of visual similarity between two colors!
584// A result of 0 would mean identical colors, while a result of 1 or higher
585// means the colors differ a lot.
586func (c1 Color) DistanceLab(c2 Color) float64 {
587	l1, a1, b1 := c1.Lab()
588	l2, a2, b2 := c2.Lab()
589	return math.Sqrt(sq(l1-l2) + sq(a1-a2) + sq(b1-b2))
590}
591
592// That's actually the same, but I don't want to break code.
593func (c1 Color) DistanceCIE76(c2 Color) float64 {
594	return c1.DistanceLab(c2)
595}
596
597// Uses the CIE94 formula to calculate color distance. More accurate than
598// DistanceLab, but also more work.
599func (cl Color) DistanceCIE94(cr Color) float64 {
600	l1, a1, b1 := cl.Lab()
601	l2, a2, b2 := cr.Lab()
602
603	// NOTE: Since all those formulas expect L,a,b values 100x larger than we
604	//       have them in this library, we either need to adjust all constants
605	//       in the formula, or convert the ranges of L,a,b before, and then
606	//       scale the distances down again. The latter is less error-prone.
607	l1, a1, b1 = l1*100.0, a1*100.0, b1*100.0
608	l2, a2, b2 = l2*100.0, a2*100.0, b2*100.0
609
610	kl := 1.0 // 2.0 for textiles
611	kc := 1.0
612	kh := 1.0
613	k1 := 0.045 // 0.048 for textiles
614	k2 := 0.015 // 0.014 for textiles.
615
616	deltaL := l1 - l2
617	c1 := math.Sqrt(sq(a1) + sq(b1))
618	c2 := math.Sqrt(sq(a2) + sq(b2))
619	deltaCab := c1 - c2
620
621	// Not taking Sqrt here for stability, and it's unnecessary.
622	deltaHab2 := sq(a1-a2) + sq(b1-b2) - sq(deltaCab)
623	sl := 1.0
624	sc := 1.0 + k1*c1
625	sh := 1.0 + k2*c1
626
627	vL2 := sq(deltaL / (kl * sl))
628	vC2 := sq(deltaCab / (kc * sc))
629	vH2 := deltaHab2 / sq(kh*sh)
630
631	return math.Sqrt(vL2+vC2+vH2) * 0.01 // See above.
632}
633
634// DistanceCIEDE2000 uses the Delta E 2000 formula to calculate color
635// distance. It is more expensive but more accurate than both DistanceLab
636// and DistanceCIE94.
637func (cl Color) DistanceCIEDE2000(cr Color) float64 {
638	return cl.DistanceCIEDE2000klch(cr, 1.0, 1.0, 1.0)
639}
640
641// DistanceCIEDE2000klch uses the Delta E 2000 formula with custom values
642// for the weighting factors kL, kC, and kH.
643func (cl Color) DistanceCIEDE2000klch(cr Color, kl, kc, kh float64) float64 {
644	l1, a1, b1 := cl.Lab()
645	l2, a2, b2 := cr.Lab()
646
647	// As with CIE94, we scale up the ranges of L,a,b beforehand and scale
648	// them down again afterwards.
649	l1, a1, b1 = l1*100.0, a1*100.0, b1*100.0
650	l2, a2, b2 = l2*100.0, a2*100.0, b2*100.0
651
652	cab1 := math.Sqrt(sq(a1) + sq(b1))
653	cab2 := math.Sqrt(sq(a2) + sq(b2))
654	cabmean := (cab1 + cab2) / 2
655
656	g := 0.5 * (1 - math.Sqrt(math.Pow(cabmean, 7)/(math.Pow(cabmean, 7)+math.Pow(25, 7))))
657	ap1 := (1 + g) * a1
658	ap2 := (1 + g) * a2
659	cp1 := math.Sqrt(sq(ap1) + sq(b1))
660	cp2 := math.Sqrt(sq(ap2) + sq(b2))
661
662	hp1 := 0.0
663	if b1 != ap1 || ap1 != 0 {
664		hp1 = math.Atan2(b1, ap1)
665		if hp1 < 0 {
666			hp1 += math.Pi * 2
667		}
668		hp1 *= 180 / math.Pi
669	}
670	hp2 := 0.0
671	if b2 != ap2 || ap2 != 0 {
672		hp2 = math.Atan2(b2, ap2)
673		if hp2 < 0 {
674			hp2 += math.Pi * 2
675		}
676		hp2 *= 180 / math.Pi
677	}
678
679	deltaLp := l2 - l1
680	deltaCp := cp2 - cp1
681	dhp := 0.0
682	cpProduct := cp1 * cp2
683	if cpProduct != 0 {
684		dhp = hp2 - hp1
685		if dhp > 180 {
686			dhp -= 360
687		} else if dhp < -180 {
688			dhp += 360
689		}
690	}
691	deltaHp := 2 * math.Sqrt(cpProduct) * math.Sin(dhp/2*math.Pi/180)
692
693	lpmean := (l1 + l2) / 2
694	cpmean := (cp1 + cp2) / 2
695	hpmean := hp1 + hp2
696	if cpProduct != 0 {
697		hpmean /= 2
698		if math.Abs(hp1-hp2) > 180 {
699			if hp1+hp2 < 360 {
700				hpmean += 180
701			} else {
702				hpmean -= 180
703			}
704		}
705	}
706
707	t := 1 - 0.17*math.Cos((hpmean-30)*math.Pi/180) + 0.24*math.Cos(2*hpmean*math.Pi/180) + 0.32*math.Cos((3*hpmean+6)*math.Pi/180) - 0.2*math.Cos((4*hpmean-63)*math.Pi/180)
708	deltaTheta := 30 * math.Exp(-sq((hpmean-275)/25))
709	rc := 2 * math.Sqrt(math.Pow(cpmean, 7)/(math.Pow(cpmean, 7)+math.Pow(25, 7)))
710	sl := 1 + (0.015*sq(lpmean-50))/math.Sqrt(20+sq(lpmean-50))
711	sc := 1 + 0.045*cpmean
712	sh := 1 + 0.015*cpmean*t
713	rt := -math.Sin(2*deltaTheta*math.Pi/180) * rc
714
715	return math.Sqrt(sq(deltaLp/(kl*sl))+sq(deltaCp/(kc*sc))+sq(deltaHp/(kh*sh))+rt*(deltaCp/(kc*sc))*(deltaHp/(kh*sh))) * 0.01
716}
717
718// BlendLab blends two colors in the L*a*b* color-space, which should result in a smoother blend.
719// t == 0 results in c1, t == 1 results in c2
720func (c1 Color) BlendLab(c2 Color, t float64) Color {
721	l1, a1, b1 := c1.Lab()
722	l2, a2, b2 := c2.Lab()
723	return Lab(l1+t*(l2-l1),
724		a1+t*(a2-a1),
725		b1+t*(b2-b1))
726}
727
728/// L*u*v* ///
729//////////////
730// http://en.wikipedia.org/wiki/CIELUV#XYZ_.E2.86.92_CIELUV_and_CIELUV_.E2.86.92_XYZ_conversions
731// For L*u*v*, we need to L*u*v*<->XYZ<->RGB and the first one is device dependent.
732
733func XyzToLuv(x, y, z float64) (l, a, b float64) {
734	// Use D65 white as reference point by default.
735	// http://www.fredmiranda.com/forum/topic/1035332
736	// http://en.wikipedia.org/wiki/Standard_illuminant
737	return XyzToLuvWhiteRef(x, y, z, D65)
738}
739
740func XyzToLuvWhiteRef(x, y, z float64, wref [3]float64) (l, u, v float64) {
741	if y/wref[1] <= 6.0/29.0*6.0/29.0*6.0/29.0 {
742		l = y / wref[1] * 29.0 / 3.0 * 29.0 / 3.0 * 29.0 / 3.0
743	} else {
744		l = 1.16*math.Cbrt(y/wref[1]) - 0.16
745	}
746	ubis, vbis := xyz_to_uv(x, y, z)
747	un, vn := xyz_to_uv(wref[0], wref[1], wref[2])
748	u = 13.0 * l * (ubis - un)
749	v = 13.0 * l * (vbis - vn)
750	return
751}
752
753// For this part, we do as R's graphics.hcl does, not as wikipedia does.
754// Or is it the same?
755func xyz_to_uv(x, y, z float64) (u, v float64) {
756	denom := x + 15.0*y + 3.0*z
757	if denom == 0.0 {
758		u, v = 0.0, 0.0
759	} else {
760		u = 4.0 * x / denom
761		v = 9.0 * y / denom
762	}
763	return
764}
765
766func LuvToXyz(l, u, v float64) (x, y, z float64) {
767	// D65 white (see above).
768	return LuvToXyzWhiteRef(l, u, v, D65)
769}
770
771func LuvToXyzWhiteRef(l, u, v float64, wref [3]float64) (x, y, z float64) {
772	//y = wref[1] * lab_finv((l + 0.16) / 1.16)
773	if l <= 0.08 {
774		y = wref[1] * l * 100.0 * 3.0 / 29.0 * 3.0 / 29.0 * 3.0 / 29.0
775	} else {
776		y = wref[1] * cub((l+0.16)/1.16)
777	}
778	un, vn := xyz_to_uv(wref[0], wref[1], wref[2])
779	if l != 0.0 {
780		ubis := u/(13.0*l) + un
781		vbis := v/(13.0*l) + vn
782		x = y * 9.0 * ubis / (4.0 * vbis)
783		z = y * (12.0 - 3.0*ubis - 20.0*vbis) / (4.0 * vbis)
784	} else {
785		x, y = 0.0, 0.0
786	}
787	return
788}
789
790// Converts the given color to CIE L*u*v* space using D65 as reference white.
791// L* is in [0..1] and both u* and v* are in about [-1..1]
792func (col Color) Luv() (l, u, v float64) {
793	return XyzToLuv(col.Xyz())
794}
795
796// Converts the given color to CIE L*u*v* space, taking into account
797// a given reference white. (i.e. the monitor's white)
798// L* is in [0..1] and both u* and v* are in about [-1..1]
799func (col Color) LuvWhiteRef(wref [3]float64) (l, u, v float64) {
800	x, y, z := col.Xyz()
801	return XyzToLuvWhiteRef(x, y, z, wref)
802}
803
804// Generates a color by using data given in CIE L*u*v* space using D65 as reference white.
805// L* is in [0..1] and both u* and v* are in about [-1..1]
806// WARNING: many combinations of `l`, `a`, and `b` values do not have corresponding
807//          valid RGB values, check the FAQ in the README if you're unsure.
808func Luv(l, u, v float64) Color {
809	return Xyz(LuvToXyz(l, u, v))
810}
811
812// Generates a color by using data given in CIE L*u*v* space, taking
813// into account a given reference white. (i.e. the monitor's white)
814// L* is in [0..1] and both u* and v* are in about [-1..1]
815func LuvWhiteRef(l, u, v float64, wref [3]float64) Color {
816	return Xyz(LuvToXyzWhiteRef(l, u, v, wref))
817}
818
819// DistanceLuv is a good measure of visual similarity between two colors!
820// A result of 0 would mean identical colors, while a result of 1 or higher
821// means the colors differ a lot.
822func (c1 Color) DistanceLuv(c2 Color) float64 {
823	l1, u1, v1 := c1.Luv()
824	l2, u2, v2 := c2.Luv()
825	return math.Sqrt(sq(l1-l2) + sq(u1-u2) + sq(v1-v2))
826}
827
828// BlendLuv blends two colors in the CIE-L*u*v* color-space, which should result in a smoother blend.
829// t == 0 results in c1, t == 1 results in c2
830func (c1 Color) BlendLuv(c2 Color, t float64) Color {
831	l1, u1, v1 := c1.Luv()
832	l2, u2, v2 := c2.Luv()
833	return Luv(l1+t*(l2-l1),
834		u1+t*(u2-u1),
835		v1+t*(v2-v1))
836}
837
838/// HCL ///
839///////////
840// HCL is nothing else than L*a*b* in cylindrical coordinates!
841// (this was wrong on English wikipedia, I fixed it, let's hope the fix stays.)
842// But it is widely popular since it is a "correct HSV"
843// http://www.hunterlab.com/appnotes/an09_96a.pdf
844
845// Converts the given color to HCL space using D65 as reference white.
846// H values are in [0..360], C and L values are in [0..1] although C can overshoot 1.0
847func (col Color) Hcl() (h, c, l float64) {
848	return col.HclWhiteRef(D65)
849}
850
851func LabToHcl(L, a, b float64) (h, c, l float64) {
852	// Oops, floating point workaround necessary if a ~= b and both are very small (i.e. almost zero).
853	if math.Abs(b-a) > 1e-4 && math.Abs(a) > 1e-4 {
854		h = math.Mod(57.29577951308232087721*math.Atan2(b, a)+360.0, 360.0) // Rad2Deg
855	} else {
856		h = 0.0
857	}
858	c = math.Sqrt(sq(a) + sq(b))
859	l = L
860	return
861}
862
863// Converts the given color to HCL space, taking into account
864// a given reference white. (i.e. the monitor's white)
865// H values are in [0..360], C and L values are in [0..1]
866func (col Color) HclWhiteRef(wref [3]float64) (h, c, l float64) {
867	L, a, b := col.LabWhiteRef(wref)
868	return LabToHcl(L, a, b)
869}
870
871// Generates a color by using data given in HCL space using D65 as reference white.
872// H values are in [0..360], C and L values are in [0..1]
873// WARNING: many combinations of `l`, `a`, and `b` values do not have corresponding
874//          valid RGB values, check the FAQ in the README if you're unsure.
875func Hcl(h, c, l float64) Color {
876	return HclWhiteRef(h, c, l, D65)
877}
878
879func HclToLab(h, c, l float64) (L, a, b float64) {
880	H := 0.01745329251994329576 * h // Deg2Rad
881	a = c * math.Cos(H)
882	b = c * math.Sin(H)
883	L = l
884	return
885}
886
887// Generates a color by using data given in HCL space, taking
888// into account a given reference white. (i.e. the monitor's white)
889// H values are in [0..360], C and L values are in [0..1]
890func HclWhiteRef(h, c, l float64, wref [3]float64) Color {
891	L, a, b := HclToLab(h, c, l)
892	return LabWhiteRef(L, a, b, wref)
893}
894
895// BlendHcl blends two colors in the CIE-L*C*h° color-space, which should result in a smoother blend.
896// t == 0 results in c1, t == 1 results in c2
897func (col1 Color) BlendHcl(col2 Color, t float64) Color {
898	h1, c1, l1 := col1.Hcl()
899	h2, c2, l2 := col2.Hcl()
900
901	// We know that h are both in [0..360]
902	return Hcl(interp_angle(h1, h2, t), c1+t*(c2-c1), l1+t*(l2-l1))
903}
904