From e8ed4000f0903520fc43ad7904452a94c14fce6a Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Thu, 2 Jun 2016 21:24:19 +0200 Subject: [PATCH 01/30] ColorConvert: the basic stuff --- mathics/builtin/graphics.py | 145 +++++++++++++++++++++++++++++++++++- 1 file changed, 144 insertions(+), 1 deletion(-) diff --git a/mathics/builtin/graphics.py b/mathics/builtin/graphics.py index f0877dc431..b296a4ef88 100644 --- a/mathics/builtin/graphics.py +++ b/mathics/builtin/graphics.py @@ -10,6 +10,7 @@ from __future__ import division from math import floor, ceil, log10 +import math import json import base64 from six.moves import map @@ -183,7 +184,7 @@ class Graphics(Builtin): In 'TeXForm', 'Graphics' produces Asymptote figures: >> Graphics[Circle[]] // TeXForm - = + = . \begin{asy} . size(5.8556cm, 5.8333cm); . draw(ellipse((175,175),175,175), rgb(0, 0, 0)+linewidth(0.66667)); @@ -350,6 +351,115 @@ class RGBColor(_Color): def to_rgba(self): return self.components + def to_xyza(self): + components = [max(0, min(1, c)) for c in self.components] + + # inverse sRGB companding. see http://www.brucelindbloom.com/Eqn_RGB_to_XYZ.html + r, g, b = (math.pow(((c + 0.055) / 1.055), 2.4) + if c > 0.04045 else c / 12.92 for c in components[:3]) + + # use rRGB D50 conversion like MMA. see http://www.brucelindbloom.com/Eqn_RGB_XYZ_Matrix.html + return [0.4360747 * r + 0.3850649 * g + 0.1430804 * b, + 0.2225045 * r + 0.7168786 * g + 0.0606169 * b, + 0.0139322 * r + 0.0971045 * g + 0.7141733 * b] + components[3:] + + def to_laba(self): + return XYZColor(components=self.to_xyza()).to_laba() + + def to_lcha(self): + return XYZColor(components=self.to_xyza()).to_lcha() + + +class LABColor(_Color): + components_sizes = [3, 4] + default_components = [0, 0, 0, 1] + + def to_rgba(self): + return XYZColor(components=self.to_xyza()).to_rgba() + + def to_xyza(self): + components = [max(0, min(1, c)) for c in self.components] + + # see https://en.wikipedia.org/wiki/Lab_color_space´ + def inv_f(t): + if t > (6. / 29.): + return math.pow(t, 3) + else: + return 3. * (36. / 841.) * (t - (4. / 29.)) + + l, a, b = components[:3] + l0 = (l + 0.16) / 1.16 + x, y, z = [inv_f(l0 + a / 5.), inv_f(l0), inv_f(l0 - b / 2.)] + + # D50 white; taken from http://www.easyrgb.com/index.php?X=MATH&H=15#text15 + xyz_ref_white = (0.96422, 1.0, 0.82521) + + return [xyz * w for xyz, w in zip((x, y, z), xyz_ref_white)] + components[3:] + + def to_laba(self): + return self.components + + +class LCHColor(_Color): + components_sizes = [3, 4] + default_components = [0, 0, 0, 1] + + def to_rgba(self): + return XYZColor(components=self.to_xyza()).to_rgba() + + def to_xyza(self): + raise NotImplementedError + + def to_lcha(self): + return self.components + + +class XYZColor(_Color): + components_sizes = [3, 4] + default_components = [0, 0, 0, 1] + + def to_rgba(self): + # FIXME this is still somewhat inexact compared to MMA + + components = [max(0, min(1, c)) for c in self.components] + + # the inverse matrix of the one in RGBColor.to_xyza() + x, y, z = components[:3] + r, g, b = [x * 3.13386 + y * -1.61687 + z * -0.490615, + x * -0.978769 + y * 1.91614 + z * 0.0334541, + x * 0.0719452 + y * -0.228991 + z * 1.40524] + + return [1.055 * math.pow(c, 1. / 2.4) - 0.055 if c > 0.0031308 + else c * 12.92 for c in (r, g, b)] + components[3:] + + def to_laba(self): + components = [max(0, min(1, c)) for c in self.components] + + # D50 white; taken from http://www.easyrgb.com/index.php?X=MATH&H=15#text15 + xyz_ref_white = (0.96422, 1.0, 0.82521) + + # computation of x, y, z; see https://en.wikipedia.org/wiki/Lab_color_space + components = [xyz / w for xyz, w in zip(components, xyz_ref_white)] + + t0 = 0.0088564516790356308172 # math.pow(6. / 29., 3) + a = 7.7870370370370 # (1. / 3.) * math.pow(29. / 6., 2) + x, y, z = (math.pow(t, 1. / 3.) if t > t0 + else a * t + (4. / 29.) for t in components) + + # MMA scales by 1/100 + return [(1.16 * y) - 0.16, 5. * (x - y), 2. * (y - z)] + components[3:] + + def to_lcha(self): + # see https://en.wikipedia.org/wiki/Lab_color_space + # see http://www.brucelindbloom.com/Eqn_Lab_to_LCH.html + laba = self.to_laba() + l, a, b = laba[:3] + h = math.atan2(b, a) + if h < 0: + h += 2 * math.pi + h /= 2 * math.pi # MMA specific + return [l, math.sqrt(a * a + b * b), h] + laba[3:] + class CMYKColor(_Color): """ @@ -466,6 +576,36 @@ def to_rgba(self): return (g, g, g, self.components[1]) +class ColorConvert(Builtin): + _convert = { + 'LAB': lambda color: LABColor(components=color.to_laba()), + 'LCH': lambda color: LCHColor(components=color.to_lcha()), + 'RGB': lambda color: RGBColor(components=color.to_rgba()), + 'XYZ': lambda color: XYZColor(components=color.to_xyza()) + } + + def apply(self, color, colorspace, evaluation): + 'ColorConvert[color_, colorspace_String]' + py_color = _Color.create(color) + if not isinstance(py_color, _Color): + return + + convert = ColorConvert._convert.get(colorspace.get_string_value()) + if not convert: + return + + converted_color = convert(py_color) + return Expression(converted_color.get_name(), *converted_color.components) + + +def _euclidean_distance(a, b): + return math.sqrt(sum((x1 - x2) * (x1 - x2) for x1, x2 in zip(a, b))) + + +def _component_distance(a, b, i): + return abs(a[i] - b[i]) + + class _Size(_GraphicsElement): def init(self, graphics, item=None, value=None): super(_Size, self).init(graphics, item) @@ -2364,6 +2504,9 @@ class Large(Builtin): styles = system_symbols_dict({ 'RGBColor': RGBColor, + 'XYZColor': XYZColor, + 'LABColor': LABColor, + 'LCHColor': LCHColor, 'CMYKColor': CMYKColor, 'Hue': Hue, 'GrayLevel': GrayLevel, From 3d806b8908c447a551bb6ede9f1ed10519c84125 Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Thu, 2 Jun 2016 22:51:06 +0200 Subject: [PATCH 02/30] ColorConvert: GrayLevel and more --- mathics/builtin/graphics.py | 106 +++++++++++++++++++++++++++--------- 1 file changed, 80 insertions(+), 26 deletions(-) diff --git a/mathics/builtin/graphics.py b/mathics/builtin/graphics.py index b296a4ef88..7f048000fe 100644 --- a/mathics/builtin/graphics.py +++ b/mathics/builtin/graphics.py @@ -326,6 +326,36 @@ def to_js(self): def to_expr(self): return Expression(self.get_name(), *self.components) + def to_cmyka(self): + rgba = self.to_rgba() + r, g, b = rgba[:3] + k = 1 - max(r, g, b) + w = 1 - k + return [(1 - r - k) / w, (1 - g - k) / w, (1 - b - k) / w] + rgba[3:] + + def to_ga(self): + # see https://en.wikipedia.org/wiki/Grayscale + rgba = self.to_rgba() + r, g, b = rgba[:3] + y = 0.299 * r + 0.587 * g + 0.114 * b # Y of Y'UV + return [y] + rgba[3:] + + def to_hsb(self): + raise NotImplementedError + + def to_lcha(self): + # see http://www.brucelindbloom.com/Eqn_Lab_to_LCH.html + laba = self.to_laba() + l, a, b = laba[:3] + h = math.atan2(b, a) + if h < 0: + h += 2 * math.pi + h /= 2 * math.pi # MMA specific + return [l, math.sqrt(a * a + b * b), h] + laba[3:] + + +def _clip(components): + return [max(0, min(1, c)) for c in components] class RGBColor(_Color): """ @@ -352,23 +382,20 @@ def to_rgba(self): return self.components def to_xyza(self): - components = [max(0, min(1, c)) for c in self.components] + components = _clip(self.components) # inverse sRGB companding. see http://www.brucelindbloom.com/Eqn_RGB_to_XYZ.html r, g, b = (math.pow(((c + 0.055) / 1.055), 2.4) if c > 0.04045 else c / 12.92 for c in components[:3]) # use rRGB D50 conversion like MMA. see http://www.brucelindbloom.com/Eqn_RGB_XYZ_Matrix.html - return [0.4360747 * r + 0.3850649 * g + 0.1430804 * b, + return _clip([0.4360747 * r + 0.3850649 * g + 0.1430804 * b, 0.2225045 * r + 0.7168786 * g + 0.0606169 * b, - 0.0139322 * r + 0.0971045 * g + 0.7141733 * b] + components[3:] + 0.0139322 * r + 0.0971045 * g + 0.7141733 * b]) + components[3:] def to_laba(self): return XYZColor(components=self.to_xyza()).to_laba() - def to_lcha(self): - return XYZColor(components=self.to_xyza()).to_lcha() - class LABColor(_Color): components_sizes = [3, 4] @@ -378,7 +405,7 @@ def to_rgba(self): return XYZColor(components=self.to_xyza()).to_rgba() def to_xyza(self): - components = [max(0, min(1, c)) for c in self.components] + components = self.components # see https://en.wikipedia.org/wiki/Lab_color_space´ def inv_f(t): @@ -394,7 +421,7 @@ def inv_f(t): # D50 white; taken from http://www.easyrgb.com/index.php?X=MATH&H=15#text15 xyz_ref_white = (0.96422, 1.0, 0.82521) - return [xyz * w for xyz, w in zip((x, y, z), xyz_ref_white)] + components[3:] + return _clip([xyz * w for xyz, w in zip((x, y, z), xyz_ref_white)]) + components[3:] def to_laba(self): return self.components @@ -408,7 +435,15 @@ def to_rgba(self): return XYZColor(components=self.to_xyza()).to_rgba() def to_xyza(self): - raise NotImplementedError + return LABColor(components=self.to_laba()).to_xyza() + + def to_laba(self): + lcha = self.to_lcha() + l, c, h = lcha[:3] + h *= 2 * math.pi # MMA specific + a = c * math.cos(h) + b = c * math.sin(h) + return [l, a, b] + lcha[3:] def to_lcha(self): return self.components @@ -419,9 +454,7 @@ class XYZColor(_Color): default_components = [0, 0, 0, 1] def to_rgba(self): - # FIXME this is still somewhat inexact compared to MMA - - components = [max(0, min(1, c)) for c in self.components] + components = _clip(self.components) # the inverse matrix of the one in RGBColor.to_xyza() x, y, z = components[:3] @@ -429,11 +462,14 @@ def to_rgba(self): x * -0.978769 + y * 1.91614 + z * 0.0334541, x * 0.0719452 + y * -0.228991 + z * 1.40524] - return [1.055 * math.pow(c, 1. / 2.4) - 0.055 if c > 0.0031308 - else c * 12.92 for c in (r, g, b)] + components[3:] + return _clip([1.055 * math.pow(c, 1. / 2.4) - 0.055 if c > 0.0031308 + else c * 12.92 for c in (r, g, b)]) + components[3:] + + def to_xyza(self): + return self.components def to_laba(self): - components = [max(0, min(1, c)) for c in self.components] + components = _clip(self.components) # D50 white; taken from http://www.easyrgb.com/index.php?X=MATH&H=15#text15 xyz_ref_white = (0.96422, 1.0, 0.82521) @@ -449,17 +485,6 @@ def to_laba(self): # MMA scales by 1/100 return [(1.16 * y) - 0.16, 5. * (x - y), 2. * (y - z)] + components[3:] - def to_lcha(self): - # see https://en.wikipedia.org/wiki/Lab_color_space - # see http://www.brucelindbloom.com/Eqn_Lab_to_LCH.html - laba = self.to_laba() - l, a, b = laba[:3] - h = math.atan2(b, a) - if h < 0: - h += 2 * math.pi - h /= 2 * math.pi # MMA specific - return [l, math.sqrt(a * a + b * b), h] + laba[3:] - class CMYKColor(_Color): """ @@ -476,6 +501,9 @@ class CMYKColor(_Color): components_sizes = [3, 4, 5] default_components = [0, 0, 0, 0, 1] + def to_cmyka(self): + return self.components + def to_rgba(self): k = self.components[3] k_ = 1 - k @@ -484,6 +512,15 @@ def to_rgba(self): rgb = (1 - cmy[0], 1 - cmy[1], 1 - cmy[2]) return rgb + (c[4],) + def to_xyza(self): + return RGBColor(components=self.to_rgba()).to_xyza() + + def to_laba(self): + return RGBColor(components=self.to_rgba()).to_laba() + + def to_lcha(self): + return RGBColor(components=self.to_rgba()).to_lcha() + class Hue(_Color): """ @@ -557,6 +594,8 @@ def trans(t): result = tuple([trans(list(map(t))) for t in rgb]) + (self.components[3],) return result + def to_hsb(self): + return self.components class GrayLevel(_Color): """ @@ -575,9 +614,24 @@ def to_rgba(self): g = self.components[0] return (g, g, g, self.components[1]) + def to_ga(self): + return self.components + + def to_xyza(self): + return RGBColor(components=self.to_rgba()).to_xyza() + + def to_laba(self): + return RGBColor(components=self.to_rgba()).to_laba() + + def to_lcha(self): + return RGBColor(components=self.to_rgba()).to_lcha() + class ColorConvert(Builtin): _convert = { + 'CMYK': lambda color: CMYKColor(components=color.to_cmyk()), + 'Grayscale': lambda color: GrayLevel(components=color.to_ga()), + 'HSB': lambda color: Hue(components=color.to_hsb()), 'LAB': lambda color: LABColor(components=color.to_laba()), 'LCH': lambda color: LCHColor(components=color.to_lcha()), 'RGB': lambda color: RGBColor(components=color.to_rgba()), From 3aabe9bf122d03617e219d5efd400e0d9e0dbc48 Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Sat, 4 Jun 2016 10:04:40 +0200 Subject: [PATCH 03/30] work in progress, cleanup: LUVColors, ref white --- mathics/builtin/graphics.py | 271 ++++++++++++++++++++++++++++++------ 1 file changed, 228 insertions(+), 43 deletions(-) diff --git a/mathics/builtin/graphics.py b/mathics/builtin/graphics.py index 7f048000fe..de952aac11 100644 --- a/mathics/builtin/graphics.py +++ b/mathics/builtin/graphics.py @@ -23,7 +23,7 @@ from mathics.builtin.options import options_to_rules from mathics.core.expression import ( Expression, Integer, Real, String, Symbol, strip_context, - system_symbols, system_symbols_dict) + system_symbols, system_symbols_dict, from_python) class CoordinatesError(BoxConstructError): @@ -164,6 +164,34 @@ def _data_and_options(leaves, defined_options): return data, options +class _PerfectReflectingDiffuser: + def __init__(self, x1, x2): + scale = 1.0 / x2 + self.xyz = (x1 * scale, 1.0, (1.0 - x1 - x2) * scale) + + q_r = self.xyz[0] + 15. * self.xyz[1] + 3. * self.xyz[2] + self.u_r = 4. * self.xyz[0] / q_r + self.v_r = 9. * self.xyz[1] / q_r + + +# MMA's reference white is a # D50, 2 degrees diffuser; for the +# values, see https://en.wikipedia.org/wiki/Standard_illuminant + +_ref_white = _PerfectReflectingDiffuser(0.34567, 0.35850) + + +def _clip(components): # clip to [0, 1] + return [max(0, min(1, c)) for c in components] + + +def _euclidean_distance(a, b): + return math.sqrt(sum((x1 - x2) * (x1 - x2) for x1, x2 in zip(a, b))) + + +def _component_distance(a, b, i): + return abs(a[i] - b[i]) + + class Graphics(Builtin): r"""
@@ -288,9 +316,13 @@ def init(self, item=None, components=None): components = [value.round_to_float() for value in leaves] if None in components or not all(0 <= c <= 1 for c in components): raise ColorError - if len(components) < len(self.default_components): - components.extend(self.default_components[ - len(components):]) + # the following lines always extend to the maximum available + # default_components, so RGBColor[0, 0, 0] will _always_ + # become RGBColor[0, 0, 0, 1]. does not seem the right thing + # to do in this general context. poke1024 + # if len(components) < len(self.default_components): + # components.extend(self.default_components[ + # len(components):]) self.components = components else: raise ColorError @@ -340,8 +372,11 @@ def to_ga(self): y = 0.299 * r + 0.587 * g + 0.114 * b # Y of Y'UV return [y] + rgba[3:] - def to_hsb(self): - raise NotImplementedError + def to_hsba(self): + return RGBColor(components=self.to_rgba()).to_hsba() + + def to_laba(self): + return XYZColor(components=self.to_xyza()).to_laba() def to_lcha(self): # see http://www.brucelindbloom.com/Eqn_Lab_to_LCH.html @@ -353,9 +388,15 @@ def to_lcha(self): h /= 2 * math.pi # MMA specific return [l, math.sqrt(a * a + b * b), h] + laba[3:] + def to_luva(self): + return XYZColor(components=self.to_xyza()).to_luva() + + def to_rgba(self): + return XYZColor(components=self.to_xyza()).to_rgba() + + def to_xyza(self): + raise NotImplementedError -def _clip(components): - return [max(0, min(1, c)) for c in components] class RGBColor(_Color): """ @@ -393,17 +434,46 @@ def to_xyza(self): 0.2225045 * r + 0.7168786 * g + 0.0606169 * b, 0.0139322 * r + 0.0971045 * g + 0.7141733 * b]) + components[3:] - def to_laba(self): - return XYZColor(components=self.to_xyza()).to_laba() + def to_hsba(self): + # see https://en.wikipedia.org/wiki/HSB_color_space. HSB is also known as HSV. + + components = _clip(self.components) + + r, g, b = components[:3] + m1 = max(r, g, b) + m0 = min(r, g, b) + c = m1 - m0 + + if c < 1e-5: # FIXME + h = 0 + elif m1 == r: + h = ((g - b) / c) % 6 + elif m1 == g: + h = (b - r) / c + 2 + else: # m1 == b + h = (r - g) / c + 4 + h = (h * 60.) / 360. + if h < 0.: + h += 1. + + b = m1 + s = c / b + + return [h, s, b] + components[3:] class LABColor(_Color): + """ +
+
'LABColor[$l$, $a$, $b$]' +
represents a color with the specified lightness, red/green and yellow/blue + components in the CIE 1976 L*a*b* (CIELAB) color space. +
+ """ + components_sizes = [3, 4] default_components = [0, 0, 0, 1] - def to_rgba(self): - return XYZColor(components=self.to_xyza()).to_rgba() - def to_xyza(self): components = self.components @@ -418,22 +488,24 @@ def inv_f(t): l0 = (l + 0.16) / 1.16 x, y, z = [inv_f(l0 + a / 5.), inv_f(l0), inv_f(l0 - b / 2.)] - # D50 white; taken from http://www.easyrgb.com/index.php?X=MATH&H=15#text15 - xyz_ref_white = (0.96422, 1.0, 0.82521) - - return _clip([xyz * w for xyz, w in zip((x, y, z), xyz_ref_white)]) + components[3:] + return _clip([xyz * w for xyz, w in zip((x, y, z), _ref_white.xyz)]) + components[3:] def to_laba(self): return self.components class LCHColor(_Color): + """ +
+
'LCHColor[$l$, $c$, $h$]' +
represents a color with the specified lightness, chroma and hue + components in the CIELCh CIELab cube color space. +
+ """ + components_sizes = [3, 4] default_components = [0, 0, 0, 1] - def to_rgba(self): - return XYZColor(components=self.to_xyza()).to_rgba() - def to_xyza(self): return LABColor(components=self.to_laba()).to_xyza() @@ -449,15 +521,52 @@ def to_lcha(self): return self.components +class LUVColor(_Color): + """ +
+
'LCHColor[$l$, $u$, $v$]' +
represents a color with the specified components in the CIE 1976 L*u*v* (CIELUV) color space. +
+ """ + + components_sizes = [3, 4] + default_components = [0, 0, 0, 1] + + def to_xyza(self): + lum, u, v = self.components[:3] + + u_0 = u / (13. * lum) + _ref_white.u_r + v_0 = v / (13. * lum) + _ref_white.v_r + + lum *= 100.0 # MMA specific + + if lum <= 8.: + y = _ref_white.xyz[1] * lum * math.pow(3. / 29., 3) + else: + y = _ref_white.xyz[1] * math.pow((lum + 16.) / 116., 3) + + x = y * (9. * u_0) / (4. * v_0) + z = y * (12. - 3. * u_0 - 20. * v_0) / (4. * v_0) + + return _clip([x, y, z]) + self.components[3:] + + class XYZColor(_Color): + """ +
+
'XYZColor[$x$, $y$, $z$]' +
represents a color with the specified components in the CIE 1931 XYZ color space. +
+ """ + components_sizes = [3, 4] default_components = [0, 0, 0, 1] def to_rgba(self): components = _clip(self.components) - # the inverse matrix of the one in RGBColor.to_xyza() x, y, z = components[:3] + # multiply with the inverse matrix of the one in RGBColor.to_xyza() r, g, b = [x * 3.13386 + y * -1.61687 + z * -0.490615, x * -0.978769 + y * 1.91614 + z * 0.0334541, x * 0.0719452 + y * -0.228991 + z * 1.40524] @@ -471,11 +580,8 @@ def to_xyza(self): def to_laba(self): components = _clip(self.components) - # D50 white; taken from http://www.easyrgb.com/index.php?X=MATH&H=15#text15 - xyz_ref_white = (0.96422, 1.0, 0.82521) - # computation of x, y, z; see https://en.wikipedia.org/wiki/Lab_color_space - components = [xyz / w for xyz, w in zip(components, xyz_ref_white)] + components = [xyz / w for xyz, w in zip(components, _ref_white.xyz)] t0 = 0.0088564516790356308172 # math.pow(6. / 29., 3) a = 7.7870370370370 # (1. / 3.) * math.pow(29. / 6., 2) @@ -485,6 +591,27 @@ def to_laba(self): # MMA scales by 1/100 return [(1.16 * y) - 0.16, 5. * (x - y), 2. * (y - z)] + components[3:] + def to_luva(self): + # see http://www.brucelindbloom.com/Eqn_XYZ_to_Luv.html + # and https://en.wikipedia.org/wiki/CIELUV + + components = _clip(self.components) + + x_orig, y_orig, z_orig = components[:3] + y = y_orig / _ref_white.xyz[1] + + lum = 116. * math.pow(y, 1. / 3.) - 16. if y > 0.008856 else 903.3 * y + + q_0 = x_orig + 15. * y_orig + 3. * z_orig + u_0 = 4. * x_orig / q_0 + v_0 = 9. * y_orig / q_0 + + lum /= 100.0 # MMA specific + u = 13. * lum * (u_0 - _ref_white.u_r) + v = 13. * lum * (v_0 - _ref_white.v_r) + + return [lum, u, v] + components[3:] + class CMYKColor(_Color): """ @@ -515,12 +642,6 @@ def to_rgba(self): def to_xyza(self): return RGBColor(components=self.to_rgba()).to_xyza() - def to_laba(self): - return RGBColor(components=self.to_rgba()).to_laba() - - def to_lcha(self): - return RGBColor(components=self.to_rgba()).to_lcha() - class Hue(_Color): """ @@ -594,9 +715,13 @@ def trans(t): result = tuple([trans(list(map(t))) for t in rgb]) + (self.components[3],) return result - def to_hsb(self): + def to_xyza(self): + return RGBColor(components=self.to_rgba()).to_xyza() + + def to_hsba(self): return self.components + class GrayLevel(_Color): """
@@ -620,20 +745,33 @@ def to_ga(self): def to_xyza(self): return RGBColor(components=self.to_rgba()).to_xyza() - def to_laba(self): - return RGBColor(components=self.to_rgba()).to_laba() - def to_lcha(self): - return RGBColor(components=self.to_rgba()).to_lcha() +class ColorConvert(Builtin): + """ +
+
'ColorConvert[$c$, $colspace$]' +
returns the representation of color $c$ in the color space $colspace$. +
+ Valid values for $colspace$ are: + + CMYK + Grayscale + HSB + LAB + LCH + LUV + RGB + XYZ + """ -class ColorConvert(Builtin): _convert = { - 'CMYK': lambda color: CMYKColor(components=color.to_cmyk()), + 'CMYK': lambda color: CMYKColor(components=color.to_cmyka()), 'Grayscale': lambda color: GrayLevel(components=color.to_ga()), - 'HSB': lambda color: Hue(components=color.to_hsb()), + 'HSB': lambda color: Hue(components=color.to_hsba()), 'LAB': lambda color: LABColor(components=color.to_laba()), 'LCH': lambda color: LCHColor(components=color.to_lcha()), + 'LUV': lambda color: LUVColor(components=color.to_luva()), 'RGB': lambda color: RGBColor(components=color.to_rgba()), 'XYZ': lambda color: XYZColor(components=color.to_xyza()) } @@ -652,12 +790,58 @@ def apply(self, color, colorspace, evaluation): return Expression(converted_color.get_name(), *converted_color.components) -def _euclidean_distance(a, b): - return math.sqrt(sum((x1 - x2) * (x1 - x2) for x1, x2 in zip(a, b))) +class ColorDistance(Builtin): + """ +
+
'ColorDistance[$c1$, $c2$]' +
returns a measure of color distance between the colors $c1$ and $c2$. +
'ColorDistance[$list$, $c2$]' +
returns a list of color distances between the colors in $list$ and $c2$. +
+ The optional parameter "Distance" specified the method used to measure the color + distance. Available options are: -def _component_distance(a, b, i): - return abs(a[i] - b[i]) + CIE76 + CIE74 + DeltaL + DeltaC + DeltaH + """ + + options = { + 'Distance': 'CIE76', + } + + rules = { + 'ColorDistance[l_List, c2_]': 'ColorDistance[#, c2]& /@ l', + } + + _distances = { + "CIE76": lambda c1, c2: _euclidean_distance(c1.to_laba()[:3], c2.to_laba()[:3]), + "CIE94": lambda c1, c2: _euclidean_distance(c1.to_lcha()[:3], c2.to_lcha()[:3]), + "DeltaL": lambda c1, c2: _component_distance(c1.to_lcha(), c2.to_lcha(), 0), + "DeltaC": lambda c1, c2: _component_distance(c1.to_lcha(), c2.to_lcha(), 1), + "DeltaH": lambda c1, c2: _component_distance(c1.to_lcha(), c2.to_lcha(), 2), + } + + def apply(self, c1, c2, evaluation, options): + 'ColorDistance[c1_, c2_, OptionsPattern[ColorDistance]]' + distance = options.get('System`Distance') + if isinstance(distance, String): + compute = ColorDistance._distances.get(distance.get_string_value()) + if not compute: + evaluation.message('ColorDistance') + return + else: + compute = lambda a, b: Expression(distance, a.to_laba(), b.to_laba()) + + if c1.get_head_name() == 'System`List' and isinstance(c2, _Color): + if any(not isinstance(c, _Color) for c in c1.leaves): + return # fail + return Expression('List', *[from_python(compute(c, c2)) for c in c1.leaves]) + elif isinstance(c1, _Color) and isinstance(c2, _Color): + return from_python(compute(c1, c2)) class _Size(_GraphicsElement): @@ -2561,6 +2745,7 @@ class Large(Builtin): 'XYZColor': XYZColor, 'LABColor': LABColor, 'LCHColor': LCHColor, + 'LUVColor': LUVColor, 'CMYKColor': CMYKColor, 'Hue': Hue, 'GrayLevel': GrayLevel, From 39a9a788d67fc15012be706c5b1592de85311d05 Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Sat, 4 Jun 2016 13:32:28 +0200 Subject: [PATCH 04/30] color inversion tests working now, docs --- mathics/builtin/graphics.py | 77 +++++++++++++++++++++---------------- test/test_color.py | 77 +++++++++++++++++++++++++++++++++++++ 2 files changed, 121 insertions(+), 33 deletions(-) create mode 100644 test/test_color.py diff --git a/mathics/builtin/graphics.py b/mathics/builtin/graphics.py index de952aac11..ba528f064a 100644 --- a/mathics/builtin/graphics.py +++ b/mathics/builtin/graphics.py @@ -180,8 +180,8 @@ def __init__(self, x1, x2): _ref_white = _PerfectReflectingDiffuser(0.34567, 0.35850) -def _clip(components): # clip to [0, 1] - return [max(0, min(1, c)) for c in components] +def _clip(*components): # clip to [0, 1] + return tuple(max(0, min(1, c)) for c in components) def _euclidean_distance(a, b): @@ -313,16 +313,23 @@ def init(self, item=None, components=None): if item is not None: leaves = item.leaves if len(leaves) in self.components_sizes: + # we must not clip here; we copy the components, without clipping, + # e.g. RGBColor[-1, 0, 0] stays RGBColor[-1, 0, 0]. this is especially + # important for color spaces like LAB that have negative components. + components = [value.round_to_float() for value in leaves] - if None in components or not all(0 <= c <= 1 for c in components): + if None in components: raise ColorError + # the following lines always extend to the maximum available # default_components, so RGBColor[0, 0, 0] will _always_ # become RGBColor[0, 0, 0, 1]. does not seem the right thing # to do in this general context. poke1024 + # if len(components) < len(self.default_components): # components.extend(self.default_components[ # len(components):]) + self.components = components else: raise ColorError @@ -362,15 +369,15 @@ def to_cmyka(self): rgba = self.to_rgba() r, g, b = rgba[:3] k = 1 - max(r, g, b) - w = 1 - k - return [(1 - r - k) / w, (1 - g - k) / w, (1 - b - k) / w] + rgba[3:] + k_ = 1 - k + return ((1 - r - k) / k_, (1 - g - k) / k_, (1 - b - k) / k_, k) + tuple(rgba[3:]) def to_ga(self): # see https://en.wikipedia.org/wiki/Grayscale rgba = self.to_rgba() r, g, b = rgba[:3] y = 0.299 * r + 0.587 * g + 0.114 * b # Y of Y'UV - return [y] + rgba[3:] + return (y,) + tuple(rgba[3:]) def to_hsba(self): return RGBColor(components=self.to_rgba()).to_hsba() @@ -386,7 +393,7 @@ def to_lcha(self): if h < 0: h += 2 * math.pi h /= 2 * math.pi # MMA specific - return [l, math.sqrt(a * a + b * b), h] + laba[3:] + return (l, math.sqrt(a * a + b * b), h) + tuple(laba[3:]) def to_luva(self): return XYZColor(components=self.to_xyza()).to_luva() @@ -423,21 +430,21 @@ def to_rgba(self): return self.components def to_xyza(self): - components = _clip(self.components) + components = _clip(*self.components) # inverse sRGB companding. see http://www.brucelindbloom.com/Eqn_RGB_to_XYZ.html r, g, b = (math.pow(((c + 0.055) / 1.055), 2.4) if c > 0.04045 else c / 12.92 for c in components[:3]) # use rRGB D50 conversion like MMA. see http://www.brucelindbloom.com/Eqn_RGB_XYZ_Matrix.html - return _clip([0.4360747 * r + 0.3850649 * g + 0.1430804 * b, + return _clip(0.4360747 * r + 0.3850649 * g + 0.1430804 * b, 0.2225045 * r + 0.7168786 * g + 0.0606169 * b, - 0.0139322 * r + 0.0971045 * g + 0.7141733 * b]) + components[3:] + 0.0139322 * r + 0.0971045 * g + 0.7141733 * b) + tuple(components[3:]) def to_hsba(self): # see https://en.wikipedia.org/wiki/HSB_color_space. HSB is also known as HSV. - components = _clip(self.components) + components = _clip(*self.components) r, g, b = components[:3] m1 = max(r, g, b) @@ -459,7 +466,7 @@ def to_hsba(self): b = m1 s = c / b - return [h, s, b] + components[3:] + return (h, s, b) + tuple(components[3:]) class LABColor(_Color): @@ -477,7 +484,7 @@ class LABColor(_Color): def to_xyza(self): components = self.components - # see https://en.wikipedia.org/wiki/Lab_color_space´ + # see https://en.wikipedia.org/wiki/Lab_color_space def inv_f(t): if t > (6. / 29.): return math.pow(t, 3) @@ -486,9 +493,9 @@ def inv_f(t): l, a, b = components[:3] l0 = (l + 0.16) / 1.16 - x, y, z = [inv_f(l0 + a / 5.), inv_f(l0), inv_f(l0 - b / 2.)] + xyz = [inv_f(l0 + a / 5.), inv_f(l0), inv_f(l0 - b / 2.)] - return _clip([xyz * w for xyz, w in zip((x, y, z), _ref_white.xyz)]) + components[3:] + return _clip(*[c * w for c, w in zip(xyz, _ref_white.xyz)]) + tuple(components[3:]) def to_laba(self): return self.components @@ -515,7 +522,7 @@ def to_laba(self): h *= 2 * math.pi # MMA specific a = c * math.cos(h) b = c * math.sin(h) - return [l, a, b] + lcha[3:] + return (l, a, b) + tuple(lcha[3:]) def to_lcha(self): return self.components @@ -548,7 +555,7 @@ def to_xyza(self): x = y * (9. * u_0) / (4. * v_0) z = y * (12. - 3. * u_0 - 20. * v_0) / (4. * v_0) - return _clip([x, y, z]) + self.components[3:] + return _clip(x, y, z) + tuple(self.components[3:]) class XYZColor(_Color): @@ -562,8 +569,11 @@ class XYZColor(_Color): components_sizes = [3, 4] default_components = [0, 0, 0, 1] + _lab_t0 = math.pow(6. / 29., 3) + _lab_a = math.pow(29. / 6., 2) + def to_rgba(self): - components = _clip(self.components) + components = _clip(*self.components) x, y, z = components[:3] # multiply with the inverse matrix of the one in RGBColor.to_xyza() @@ -571,31 +581,29 @@ def to_rgba(self): x * -0.978769 + y * 1.91614 + z * 0.0334541, x * 0.0719452 + y * -0.228991 + z * 1.40524] - return _clip([1.055 * math.pow(c, 1. / 2.4) - 0.055 if c > 0.0031308 - else c * 12.92 for c in (r, g, b)]) + components[3:] + return _clip(*[1.055 * math.pow(c, 1. / 2.4) - 0.055 if c > 0.0031308 + else c * 12.92 for c in (r, g, b)]) + tuple(components[3:]) def to_xyza(self): return self.components def to_laba(self): - components = _clip(self.components) + components = _clip(*self.components) # computation of x, y, z; see https://en.wikipedia.org/wiki/Lab_color_space - components = [xyz / w for xyz, w in zip(components, _ref_white.xyz)] + components = [c / w for c, w in zip(components, _ref_white.xyz)] - t0 = 0.0088564516790356308172 # math.pow(6. / 29., 3) - a = 7.7870370370370 # (1. / 3.) * math.pow(29. / 6., 2) - x, y, z = (math.pow(t, 1. / 3.) if t > t0 - else a * t + (4. / 29.) for t in components) + x, y, z = (math.pow(t, 1. / 3.) if t > XYZColor._lab_t0 + else XYZColor._lab_a * t + (4. / 29.) for t in components) # MMA scales by 1/100 - return [(1.16 * y) - 0.16, 5. * (x - y), 2. * (y - z)] + components[3:] + return ((1.16 * y) - 0.16, 5. * (x - y), 2. * (y - z)) + tuple(components[3:]) def to_luva(self): # see http://www.brucelindbloom.com/Eqn_XYZ_to_Luv.html # and https://en.wikipedia.org/wiki/CIELUV - components = _clip(self.components) + components = _clip(*self.components) x_orig, y_orig, z_orig = components[:3] y = y_orig / _ref_white.xyz[1] @@ -610,7 +618,7 @@ def to_luva(self): u = 13. * lum * (u_0 - _ref_white.u_r) v = 13. * lum * (v_0 - _ref_white.v_r) - return [lum, u, v] + components[3:] + return (lum, u, v) + tuple(components[3:]) class CMYKColor(_Color): @@ -632,12 +640,12 @@ def to_cmyka(self): return self.components def to_rgba(self): - k = self.components[3] + k = self.components[3] if len(self.components) >= 4 else 0 k_ = 1 - k c = self.components cmy = [c[0] * k_ + k, c[1] * k_ + k, c[2] * k_ + k] rgb = (1 - cmy[0], 1 - cmy[1], 1 - cmy[2]) - return rgb + (c[4],) + return rgb + tuple(c[4:]) def to_xyza(self): return RGBColor(components=self.to_rgba()).to_xyza() @@ -683,7 +691,7 @@ def to_rgba(self): 4: (t, p, v), 5: (v, p, q), }[i] - return rgb + (self.components[3],) + return rgb + tuple(self.components[3:]) def hsl_to_rgba(self): h, s, l = self.components[:3] @@ -778,7 +786,10 @@ class ColorConvert(Builtin): def apply(self, color, colorspace, evaluation): 'ColorConvert[color_, colorspace_String]' - py_color = _Color.create(color) + try: + py_color = _Color.create(color) + except ColorError: + return if not isinstance(py_color, _Color): return diff --git a/test/test_color.py b/test/test_color.py new file mode 100644 index 0000000000..3430893b04 --- /dev/null +++ b/test/test_color.py @@ -0,0 +1,77 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +from __future__ import absolute_import +from __future__ import unicode_literals + +import os +import sys +import pexpect +import unittest +from six.moves import range +from mathics.core.expression import Expression, Integer, Rational, Symbol +from mathics.core.definitions import Definitions +from mathics.core.evaluation import Evaluation + + +class ColorTest(unittest.TestCase): + def setUp(self): + definitions = Definitions(add_builtin=True) + self.evaluation = Evaluation(definitions, format='xml') + + def testInverseConversions(self): + spaces = ("CMYK", "HSB", "LAB", "LCH", "LUV", "RGB", "XYZ") + places = 4 + for i, from_space in enumerate(spaces): + for to_space in spaces[i + 1:]: + try: + if from_space == 'HSB': + construct_name = 'Hue' + else: + construct_name = from_space + 'Color' + + original = (0.5, 0.1, 0.2) # carefully chosen components + # that give useful transformations along all color spaces + + # now calculate from_space -> to_space -> from_space + inverted = [c.to_python() for c in Expression('ColorConvert', + Expression('ColorConvert', + Expression(construct_name, *original), + to_space), + from_space).evaluate(self.evaluation).leaves] + if from_space == 'CMYK': # if cmyk, cmyk -> cmy + k = inverted[3] + inverted = [c * (1 - k) + k for c in inverted[:3]] + self.assertEqual(len(original), len(inverted)) + for x, y in zip(original, inverted): + self.assertAlmostEqual(x, y, places) + except: + print('test failed for %s -> %s -> %s' % + (from_space, to_space, from_space)) + raise + + def testConversionsFromRGB(self): + pass + + def testConversionsFromXYZ(self): + self._checkConversion("XYZ", (0.5, 0.5, 0.5), "LAB", (0.7606, 0.0484, -0.1049)) + + def _checkConversion(self, from_space, from_components, to_space, to_components): + places = 3 + + if from_space == 'HSB': + construct_name = 'Hue' + else: + construct_name = from_space + 'Color' + + components = [c.to_python() for c in Expression('ColorConvert', Expression(construct_name, *from_components), + to_space).evaluate(self.evaluation).leaves] + self.assertEqual(len(components), len(to_components)) + for x, y in zip(components, to_components): + self.assertAlmostEqual(x, y, places) + + def testDistance(self): + pass + +if __name__ == "__main__": + unittest.main() From 3c9262fefd0cc19fd0e01ed136a9f06b6f128b9f Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Sat, 4 Jun 2016 15:14:00 +0200 Subject: [PATCH 05/30] ColorDistance works now, docs --- mathics/builtin/graphics.py | 92 ++++++++++++++++++++++++------------- test/test_color.py | 4 ++ 2 files changed, 65 insertions(+), 31 deletions(-) diff --git a/mathics/builtin/graphics.py b/mathics/builtin/graphics.py index ba528f064a..b5ae2a9a9b 100644 --- a/mathics/builtin/graphics.py +++ b/mathics/builtin/graphics.py @@ -763,14 +763,14 @@ class ColorConvert(Builtin): Valid values for $colspace$ are: - CMYK - Grayscale - HSB - LAB - LCH - LUV - RGB - XYZ + CMYK: convert to CMYKColor + Grayscale: convert to GrayLevel + HSB: convert to Hue + LAB: concert to LABColor + LCH: convert to LCHColor + LUV: convert to LUVColor + RGB: convert to RGBColor + XYZ: convert to XYZColor """ _convert = { @@ -784,17 +784,22 @@ class ColorConvert(Builtin): 'XYZ': lambda color: XYZColor(components=color.to_xyza()) } + messages = { + 'ccvinput': '`` should be a color.', + 'imgcstype': '`` is not a valid color space.', + } + def apply(self, color, colorspace, evaluation): 'ColorConvert[color_, colorspace_String]' try: py_color = _Color.create(color) except ColorError: - return - if not isinstance(py_color, _Color): + evaluation.message('ColorConvert', 'ccvinput', color) return convert = ColorConvert._convert.get(colorspace.get_string_value()) if not convert: + evaluation.message('ColorConvert', 'imgcstype', colorspace) return converted_color = convert(py_color) @@ -810,22 +815,27 @@ class ColorDistance(Builtin):
returns a list of color distances between the colors in $list$ and $c2$.
- The optional parameter "Distance" specified the method used to measure the color + The option DistanceFunction specifies the method used to measure the color distance. Available options are: - CIE76 - CIE74 - DeltaL - DeltaC - DeltaH + CIE76: euclidean distance in the LABColor space + CIE94: euclidean distance in the LCHColor space + DeltaL: difference in the L component of LCHColor + DeltaC: difference in the C component of LCHColor + DeltaH: difference in the H component of LCHColor + + >> N[ColorDistance[Magenta, Green], 5] + = 2.2507 """ options = { - 'Distance': 'CIE76', + 'DistanceFunction': '"CIE76"', } - rules = { - 'ColorDistance[l_List, c2_]': 'ColorDistance[#, c2]& /@ l', + messages = { + 'invdist': '`` is not a valid color distance function.', + 'invarg': '`1` and `2` should be two colors or a color and a lists of colors or ' + + 'two lists of colors of the same length.' } _distances = { @@ -838,21 +848,41 @@ class ColorDistance(Builtin): def apply(self, c1, c2, evaluation, options): 'ColorDistance[c1_, c2_, OptionsPattern[ColorDistance]]' - distance = options.get('System`Distance') - if isinstance(distance, String): - compute = ColorDistance._distances.get(distance.get_string_value()) + distance_function = options.get('System`DistanceFunction') + if isinstance(distance_function, String): + compute = ColorDistance._distances.get(distance_function.get_string_value()) if not compute: - evaluation.message('ColorDistance') + evaluation.message('ColorDistance', 'invdist', distance_function) return else: - compute = lambda a, b: Expression(distance, a.to_laba(), b.to_laba()) - - if c1.get_head_name() == 'System`List' and isinstance(c2, _Color): - if any(not isinstance(c, _Color) for c in c1.leaves): - return # fail - return Expression('List', *[from_python(compute(c, c2)) for c in c1.leaves]) - elif isinstance(c1, _Color) and isinstance(c2, _Color): - return from_python(compute(c1, c2)) + def compute(a, b): + Expression(distance_function, a.to_laba(), b.to_laba()) + + def distance(a, b): + try: + py_a = _Color.create(a) + py_b = _Color.create(b) + except ColorError: + evaluation.message('ColorDistance', 'invarg', a, b) + raise + return from_python(compute(py_a, py_b)) + + try: + if c1.get_head_name() == 'System`List': + if c2.get_head_name() == 'System`List': + if len(c1.leaves) != len(c2.leaves): + evaluation.message('ColorDistance', 'invarg', c1, c2) + return + else: + return Expression('List', *[distance(a, b) for a, b in zip(c1.leaves, c2.leaves)]) + else: + return Expression('List', *[distance(c, c2) for c in c1.leaves]) + elif c2.get_head_name() == 'System`List': + return Expression('List', *[distance(c1, c) for c in c2.leaves]) + else: + return distance(c1, c2) + except ColorError: + return class _Size(_GraphicsElement): diff --git a/test/test_color.py b/test/test_color.py index 3430893b04..78e6d0b830 100644 --- a/test/test_color.py +++ b/test/test_color.py @@ -20,6 +20,10 @@ def setUp(self): self.evaluation = Evaluation(definitions, format='xml') def testInverseConversions(self): + # check that a conversion A -> B -> A restores the original + # components. this tests color space transformations and their + # inverse transformations. + spaces = ("CMYK", "HSB", "LAB", "LCH", "LUV", "RGB", "XYZ") places = 4 for i, from_space in enumerate(spaces): From 4d00bb9cd2fa7b1216c6dce7055a17ef114e2c71 Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Sat, 4 Jun 2016 15:16:39 +0200 Subject: [PATCH 06/30] smaller tolerance for to_hsba --- mathics/builtin/graphics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mathics/builtin/graphics.py b/mathics/builtin/graphics.py index b5ae2a9a9b..5ec4820b95 100644 --- a/mathics/builtin/graphics.py +++ b/mathics/builtin/graphics.py @@ -451,7 +451,7 @@ def to_hsba(self): m0 = min(r, g, b) c = m1 - m0 - if c < 1e-5: # FIXME + if c < 1e-15: h = 0 elif m1 == r: h = ((g - b) / c) % 6 From 11ec0990c1bbfb2198af25b26bccc74a081a60bb Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Sat, 4 Jun 2016 15:28:54 +0200 Subject: [PATCH 07/30] fixes test cases --- mathics/builtin/graphics.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/mathics/builtin/graphics.py b/mathics/builtin/graphics.py index 5ec4820b95..a13fcab8cd 100644 --- a/mathics/builtin/graphics.py +++ b/mathics/builtin/graphics.py @@ -745,7 +745,7 @@ class GrayLevel(_Color): def to_rgba(self): g = self.components[0] - return (g, g, g, self.components[1]) + return (g, g, g) + tuple(self.components[1:]) def to_ga(self): return self.components @@ -2441,11 +2441,11 @@ class Blend(Builtin):
>> Blend[{Red, Blue}] - = RGBColor[0.5, 0., 0.5, 1.] + = RGBColor[0.5, 0., 0.5] >> Blend[{Red, Blue}, 0.3] - = RGBColor[0.7, 0., 0.3, 1.] + = RGBColor[0.7, 0., 0.3] >> Blend[{Red, Blue, Green}, 0.75] - = RGBColor[0., 0.5, 0.5, 1.] + = RGBColor[0., 0.5, 0.5] >> Graphics[Table[{Blend[{Red, Green, Blue}, x], Rectangle[{10 x, 0}]}, {x, 0, 1, 1/10}]] = -Graphics- From 04bfbdc3371204de84ca6468976f9e083f416689 Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Sat, 4 Jun 2016 16:20:13 +0200 Subject: [PATCH 08/30] adapted to rgba without alpha --- mathics/builtin/graphics.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/mathics/builtin/graphics.py b/mathics/builtin/graphics.py index a13fcab8cd..071c3b1840 100644 --- a/mathics/builtin/graphics.py +++ b/mathics/builtin/graphics.py @@ -350,14 +350,16 @@ def create_as_style(klass, graphics, item): def to_css(self): rgba = self.to_rgba() + alpha = rgba[3] if len(rgba) > 3 else 1. return (r'rgb(%f%%, %f%%, %f%%)' % ( - rgba[0] * 100, rgba[1] * 100, rgba[2] * 100), rgba[3]) + rgba[0] * 100, rgba[1] * 100, rgba[2] * 100), alpha) def to_asy(self): rgba = self.to_rgba() + alpha = rgba[3] if len(rgba) > 3 else 1. return (r'rgb(%s, %s, %s)' % ( asy_number(rgba[0]), asy_number(rgba[1]), asy_number(rgba[2])), - rgba[3]) + alpha) def to_js(self): return self.to_rgba() From 9c6557e50091fdc2bb50000f5e6f81d7cad98479 Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Sat, 4 Jun 2016 16:28:12 +0200 Subject: [PATCH 09/30] fixes to test cases --- mathics/builtin/graphics.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mathics/builtin/graphics.py b/mathics/builtin/graphics.py index 071c3b1840..3aa9ffd94c 100644 --- a/mathics/builtin/graphics.py +++ b/mathics/builtin/graphics.py @@ -212,7 +212,7 @@ class Graphics(Builtin): In 'TeXForm', 'Graphics' produces Asymptote figures: >> Graphics[Circle[]] // TeXForm - = + = . \begin{asy} . size(5.8556cm, 5.8333cm); . draw(ellipse((175,175),175,175), rgb(0, 0, 0)+linewidth(0.66667)); @@ -2550,7 +2550,7 @@ class Lighter(Builtin): >> Lighter[Orange, 1/4] - = RGBColor[1., 0.625, 0.25, 1.] + = RGBColor[1., 0.625, 0.25] >> Graphics[{Lighter[Orange, 1/4], Disk[]}] = -Graphics- >> Graphics[Table[{Lighter[Orange, x], Disk[{12x, 0}]}, {x, 0, 1, 1/6}]] From 8e9ced7ef85b006412273e63c6e3e341cf392a58 Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Sat, 4 Jun 2016 16:32:41 +0200 Subject: [PATCH 10/30] one more xyz test case --- test/test_color.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/test/test_color.py b/test/test_color.py index 78e6d0b830..1123d929d7 100644 --- a/test/test_color.py +++ b/test/test_color.py @@ -54,11 +54,9 @@ def testInverseConversions(self): (from_space, to_space, from_space)) raise - def testConversionsFromRGB(self): - pass - def testConversionsFromXYZ(self): self._checkConversion("XYZ", (0.5, 0.5, 0.5), "LAB", (0.7606, 0.0484, -0.1049)) + self._checkConversion("XYZ", (0.5, 0.5, 0.5), "RGB", (0.7440, 0.7257, 0.8118)) def _checkConversion(self, from_space, from_components, to_space, to_components): places = 3 @@ -74,8 +72,5 @@ def _checkConversion(self, from_space, from_components, to_space, to_components) for x, y in zip(components, to_components): self.assertAlmostEqual(x, y, places) - def testDistance(self): - pass - if __name__ == "__main__": unittest.main() From c6bd775d7ba1cea408aa23b60c2ffba6727cb58f Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Sat, 4 Jun 2016 19:59:43 +0200 Subject: [PATCH 11/30] rgb-xyz transforms are exact now --- mathics/builtin/graphics.py | 15 +++++++++------ test/test_color.py | 15 +++++++++++---- 2 files changed, 20 insertions(+), 10 deletions(-) diff --git a/mathics/builtin/graphics.py b/mathics/builtin/graphics.py index 3aa9ffd94c..24dfb9b1af 100644 --- a/mathics/builtin/graphics.py +++ b/mathics/builtin/graphics.py @@ -439,9 +439,10 @@ def to_xyza(self): if c > 0.04045 else c / 12.92 for c in components[:3]) # use rRGB D50 conversion like MMA. see http://www.brucelindbloom.com/Eqn_RGB_XYZ_Matrix.html - return _clip(0.4360747 * r + 0.3850649 * g + 0.1430804 * b, - 0.2225045 * r + 0.7168786 * g + 0.0606169 * b, - 0.0139322 * r + 0.0971045 * g + 0.7141733 * b) + tuple(components[3:]) + # MMA seems to round matrix values to six significant digits. we do the same. + return _clip(0.436075 * r + 0.385065 * g + 0.14308 * b, + 0.222504 * r + 0.716879 * g + 0.0606169 * b, + 0.0139322 * r + 0.0971045 * g + 0.714173 * b) + tuple(components[3:]) def to_hsba(self): # see https://en.wikipedia.org/wiki/HSB_color_space. HSB is also known as HSV. @@ -578,10 +579,12 @@ def to_rgba(self): components = _clip(*self.components) x, y, z = components[:3] - # multiply with the inverse matrix of the one in RGBColor.to_xyza() + + # for matrix, see http://www.brucelindbloom.com/Eqn_RGB_XYZ_Matrix.html + # MMA seems to round matrix values to six significant digits. we do the same. r, g, b = [x * 3.13386 + y * -1.61687 + z * -0.490615, - x * -0.978769 + y * 1.91614 + z * 0.0334541, - x * 0.0719452 + y * -0.228991 + z * 1.40524] + x * -0.978768 + y * 1.91614 + z * 0.033454, + x * 0.0719453 + y * -0.228991 + z * 1.40524] return _clip(*[1.055 * math.pow(c, 1. / 2.4) - 0.055 if c > 0.0031308 else c * 12.92 for c in (r, g, b)]) + tuple(components[3:]) diff --git a/test/test_color.py b/test/test_color.py index 1123d929d7..ea91ede91b 100644 --- a/test/test_color.py +++ b/test/test_color.py @@ -54,12 +54,19 @@ def testInverseConversions(self): (from_space, to_space, from_space)) raise - def testConversionsFromXYZ(self): - self._checkConversion("XYZ", (0.5, 0.5, 0.5), "LAB", (0.7606, 0.0484, -0.1049)) - self._checkConversion("XYZ", (0.5, 0.5, 0.5), "RGB", (0.7440, 0.7257, 0.8118)) + def testConversions(self): + self._checkConversion("RGB", (0.5, 0.5, 0.5), + "XYZ", (.20638274847577825, 0.2140411190781185, 0.17662882532500096)) + self._checkConversion("RGB", (0.4, 0.2, 0.3), + "XYZ", (0.08116707006828128, 0.05773536343816594, 0.057371054671583044)) + + self._checkConversion("XYZ", (0.5, 0.5, 0.5), + "RGB", (0.743976775016277, 0.7256665017497576, 0.8118438573490818)) + self._checkConversion("XYZ", (0.4, 0.2, 0.3), + "RGB", (0.8977592548573999, 0.022700440000000155, 0.6685886356522144)) def _checkConversion(self, from_space, from_components, to_space, to_components): - places = 3 + places = 12 if from_space == 'HSB': construct_name = 'Hue' From 7ed2f75e6799b3421a52da2c08842e41177a50ab Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Sun, 5 Jun 2016 11:38:42 +0200 Subject: [PATCH 12/30] bug fixes and better accurary in the LAB conversions --- mathics/builtin/graphics.py | 34 +++++++++++++---------- test/test_color.py | 55 +++++++++++++++++++------------------ 2 files changed, 47 insertions(+), 42 deletions(-) diff --git a/mathics/builtin/graphics.py b/mathics/builtin/graphics.py index 24dfb9b1af..b68ab32584 100644 --- a/mathics/builtin/graphics.py +++ b/mathics/builtin/graphics.py @@ -166,8 +166,11 @@ def _data_and_options(leaves, defined_options): class _PerfectReflectingDiffuser: def __init__(self, x1, x2): - scale = 1.0 / x2 - self.xyz = (x1 * scale, 1.0, (1.0 - x1 - x2) * scale) + # MMA seems to use the following constants, and not the + # values derived via calculation (commented out below) + self.xyz = (0.96422, 1., 0.82521) + #scale = 1.0 / x2 + #self.xyz = (x1 * scale, 1.0, (1.0 - x1 - x2) * scale) q_r = self.xyz[0] + 15. * self.xyz[1] + 3. * self.xyz[2] self.u_r = 4. * self.xyz[0] / q_r @@ -487,18 +490,21 @@ class LABColor(_Color): def to_xyza(self): components = self.components - # see https://en.wikipedia.org/wiki/Lab_color_space + # see http://www.brucelindbloom.com/Eqn_Lab_to_XYZ.html def inv_f(t): - if t > (6. / 29.): + if t > 0.008856: return math.pow(t, 3) else: - return 3. * (36. / 841.) * (t - (4. / 29.)) + return (116 * t - 16) / 903.3 l, a, b = components[:3] - l0 = (l + 0.16) / 1.16 - xyz = [inv_f(l0 + a / 5.), inv_f(l0), inv_f(l0 - b / 2.)] + f_y = (l * 100. + 16.) / 116. - return _clip(*[c * w for c, w in zip(xyz, _ref_white.xyz)]) + tuple(components[3:]) + x = inv_f(a / 5. + f_y) + y = math.pow(f_y, 3) if (l * 100.0) > 903.3 * 0.008856 else (l * 100.0) / 903.3 + z = inv_f(f_y - b / 2.) + + return tuple(c * w for c, w in zip((x, y, z), _ref_white.xyz)) + tuple(components[3:]) def to_laba(self): return self.components @@ -572,9 +578,6 @@ class XYZColor(_Color): components_sizes = [3, 4] default_components = [0, 0, 0, 1] - _lab_t0 = math.pow(6. / 29., 3) - _lab_a = math.pow(29. / 6., 2) - def to_rgba(self): components = _clip(*self.components) @@ -593,13 +596,14 @@ def to_xyza(self): return self.components def to_laba(self): - components = _clip(*self.components) + components = self.components[:3] + + # see http://www.brucelindbloom.com/Eqn_XYZ_to_Lab.html - # computation of x, y, z; see https://en.wikipedia.org/wiki/Lab_color_space components = [c / w for c, w in zip(components, _ref_white.xyz)] - x, y, z = (math.pow(t, 1. / 3.) if t > XYZColor._lab_t0 - else XYZColor._lab_a * t + (4. / 29.) for t in components) + x, y, z = (math.pow(t, 0.33333333) if t > 0.008856 + else (903.3 * t + 16.) / 116. for t in components) # MMA scales by 1/100 return ((1.16 * y) - 0.16, 5. * (x - y), 2. * (y - z)) + tuple(components[3:]) diff --git a/test/test_color.py b/test/test_color.py index ea91ede91b..73fc4c9214 100644 --- a/test/test_color.py +++ b/test/test_color.py @@ -25,34 +25,32 @@ def testInverseConversions(self): # inverse transformations. spaces = ("CMYK", "HSB", "LAB", "LCH", "LUV", "RGB", "XYZ") - places = 4 - for i, from_space in enumerate(spaces): - for to_space in spaces[i + 1:]: - try: - if from_space == 'HSB': - construct_name = 'Hue' - else: - construct_name = from_space + 'Color' + places = 3 + for original in ((0.5, 0.1, 0.2), (0.9, 0.1, 0.1)): + for i, from_space in enumerate(spaces): + for to_space in spaces[i + 1:]: + try: + if from_space == 'HSB': + construct_name = 'Hue' + else: + construct_name = from_space + 'Color' - original = (0.5, 0.1, 0.2) # carefully chosen components - # that give useful transformations along all color spaces - - # now calculate from_space -> to_space -> from_space - inverted = [c.to_python() for c in Expression('ColorConvert', - Expression('ColorConvert', - Expression(construct_name, *original), - to_space), - from_space).evaluate(self.evaluation).leaves] - if from_space == 'CMYK': # if cmyk, cmyk -> cmy - k = inverted[3] - inverted = [c * (1 - k) + k for c in inverted[:3]] - self.assertEqual(len(original), len(inverted)) - for x, y in zip(original, inverted): - self.assertAlmostEqual(x, y, places) - except: - print('test failed for %s -> %s -> %s' % - (from_space, to_space, from_space)) - raise + # now calculate from_space -> to_space -> from_space + inverted = [c.to_python() for c in Expression('ColorConvert', + Expression('ColorConvert', + Expression(construct_name, *original), + to_space), + from_space).evaluate(self.evaluation).leaves] + if from_space == 'CMYK': # if cmyk, cmyk -> cmy + k = inverted[3] + inverted = [c * (1 - k) + k for c in inverted[:3]] + self.assertEqual(len(original), len(inverted)) + for x, y in zip(original, inverted): + self.assertAlmostEqual(x, y, places) + except: + print('test failed for %s(%f, %f, %f) -> %s -> %s' % + (from_space, original[0], original[1], original[2], to_space, from_space)) + raise def testConversions(self): self._checkConversion("RGB", (0.5, 0.5, 0.5), @@ -65,6 +63,9 @@ def testConversions(self): self._checkConversion("XYZ", (0.4, 0.2, 0.3), "RGB", (0.8977592548573999, 0.022700440000000155, 0.6685886356522144)) + self._checkConversion("XYZ", (0.1, 0.1, 0.1), + "LAB", (0.3784243088316416, 0.02835852516741566, -0.06139347105996884)) + def _checkConversion(self, from_space, from_components, to_space, to_components): places = 12 From 5a12ef6f6acdee7bd5c13caffb8c5120bf699922 Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Tue, 14 Jun 2016 21:11:34 +0200 Subject: [PATCH 13/30] work in progress: color conversion with and without numpy --- mathics/builtin/colors.py | 330 +++++++++++++++++++++++++++++++ mathics/builtin/graphics.py | 289 ++++----------------------- mathics/builtin/numpy_utils.py | 170 ++++++++++++++-- mathics/builtin/randomnumbers.py | 4 +- 4 files changed, 528 insertions(+), 265 deletions(-) create mode 100644 mathics/builtin/colors.py diff --git a/mathics/builtin/colors.py b/mathics/builtin/colors.py new file mode 100644 index 0000000000..48f3ca854e --- /dev/null +++ b/mathics/builtin/colors.py @@ -0,0 +1,330 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +from math import pi +from mathics.builtin.numpy_utils import stack, unstack, concat, array, clip, conditional, switch, choose +from mathics.builtin.numpy_utils import sqrt, floor, mod, cos, sin, arctan2, minimum, maximum, dot_t + +# use rRGB D50 conversion like MMA. see http://www.brucelindbloom.com/Eqn_RGB_XYZ_Matrix.html +# MMA seems to round matrix values to six significant digits. we do the same. + +xyz_from_rgb = [ + [0.436075, 0.385065, 0.14308], + [0.222504, 0.716879, 0.0606169], + [0.0139322, 0.0971045, 0.714173] +] + +# for matrix, see http://www.brucelindbloom.com/Eqn_RGB_XYZ_Matrix.html +# MMA seems to round matrix values to six significant digits. we do the same. +rgb_from_xyz = [ + [3.13386, -1.61687, -0.490615], + [-0.978768, 1.91614, 0.033454], + [0.0719453, -0.228991, 1.40524] +] + + +class _PerfectReflectingDiffuser: + def __init__(self, x1, x2): + # MMA seems to use the following constants, and not the + # values derived via calculation (commented out below) + self.xyz = [0.96422, 1., 0.82521] + # scale = 1.0 / x2 + # self.xyz = (x1 * scale, 1.0, (1.0 - x1 - x2) * scale) + + q_r = self.xyz[0] + 15. * self.xyz[1] + 3. * self.xyz[2] + self.u_r = 4. * self.xyz[0] / q_r + self.v_r = 9. * self.xyz[1] / q_r + + +# MMA's reference white is a # D50, 2 degrees diffuser; for the +# values, see https://en.wikipedia.org/wiki/Standard_illuminant + +_ref_white = _PerfectReflectingDiffuser(0.34567, 0.35850) + + +def rgb_to_grayscale(pixels): + # see https://en.wikipedia.org/wiki/Grayscale + components = unstack(pixels) + r, g, b = components[:3] + y = 0.299 * r + 0.587 * g + 0.114 * b # Y of Y'UV + return stack(y, *components[3:]) + + +def grayscale_to_rgb(pixels): + components = unstack(pixels) + g = components[0] + return stack(g, g, g, *components[1:]) + + +def rgb_to_xyz(pixels): + components = unstack(pixels) + rgb = components[:3] + + x, y, z = conditional(rgb, lambda t: t > 0.04045, + lambda t: ((t + 0.055) / 1.055) ** 2.4, + lambda t: t / 12.92) + xyz = dot_t(stack(x, y, z), xyz_from_rgb) + + return concat(clip(xyz, 0, 1), components[3:]) + + +def rgb_to_hsb(pixels): + # see https://en.wikipedia.org/wiki/HSB_color_space. HSB is also known as HSV. + components = unstack(pixels) + r, g, b = clip(components[:3], 0, 1) + + m1 = maximum(r, g, b) + m0 = minimum(r, g, b) + c = m1 - m0 + + h = switch( + c < 1e-15, + lambda s: 0, + m1 == r, + lambda s: mod(((s(g) - s(b)) / s(c)), 6.), + m1 == g, + lambda s: (s(b) - s(r)) / s(c) + 2., + m1 == b, + lambda s: (s(r) - s(g)) / s(c) + 4. + ) + + h = conditional(h * (60. / 360.), lambda t: t < 0., + lambda t: t + 1., + lambda t: t) + + return stack(h, c / m1, m1, *components[3:]) + + +def hsb_to_rgb(pixels): + components = unstack(pixels) + h, s, v = components[:3] + + i = floor(6 * h) + f = 6 * h - i + i = mod(i, 6) + p = v * (1 - s) + q = v * (1 - f * s) + t = v * (1 - (1 - f) * s) + + r, g, b = choose(i, + (v, t, p), + (q, v, p), + (p, v, t), + (p, q, v), + (t, p, v), + (v, p, q)) + + return stack(r, g, b, *components[3:]) + + +def cmyk_to_rgb(pixels): + c = unstack(pixels) + + if len(c) >= 4: + k = c[3] + else: + k = 0 + k_ = 1 - k + + cmy = [(x * k_ + k) for x in c[:3]] + rgb = [1 - x for x in cmy] + + return stack(*rgb, *c[4:]) + + +def rgb_to_cmyk(pixels): + c = unstack(pixels) + + r, g, b = c[:3] + k = 1 - maximum(r, g, b) + k_ = 1 - k + + return stack((1 - r - k) / k_, (1 - g - k) / k_, (1 - b - k) / k_, k, *c[3:]) + + +def xyz_to_rgb(pixels): + components = unstack(pixels) + x, y, z = clip(components[:3], 0, 1) + xyz = dot_t(stack(x, y, z), rgb_from_xyz) + + rgb = conditional(xyz, lambda t: t > 0.0031308, + lambda t: 1.055 * (t ** (1. / 2.4)) - 0.055, + lambda t: t * 12.92) + + return concat(clip(rgb, 0, 1), components[3:]) + + +def xyz_to_lab(pixels): + # see http://www.brucelindbloom.com/Eqn_XYZ_to_Lab.html + components = unstack(pixels) + xyz = components[:3] + + xyz = array([u / v for u, v in zip(xyz, _ref_white.xyz)]) + xyz = conditional(xyz, lambda c: c > 0.008856, + lambda c: c ** 0.33333333, + lambda c: (903.3 * c + 16.) / 116.) + x, y, z = xyz + + # MMA scales by 1/100 + return stack((1.16 * y) - 0.16, 5. * (x - y), 2. * (y - z), *components[3:]) + + +def xyz_to_luv(pixels): + # see http://www.brucelindbloom.com/Eqn_XYZ_to_Luv.html + # and https://en.wikipedia.org/wiki/CIELUV + + components = unstack(pixels) + xyz = clip(components[:3], 0, 1) + + x_orig, y_orig, z_orig = xyz + y = y_orig / _ref_white.xyz[1] + + lum = conditional(y, lambda y: y > 0.008856, + lambda y: 116. * (y ** (1. / 3.)) - 16, + lambda y: 903.3 * y) + + q_0 = x_orig + 15. * y_orig + 3. * z_orig + u_0 = 4. * x_orig / q_0 + v_0 = 9. * y_orig / q_0 + + lum /= 100.0 # MMA specific + + u = 13. * lum * (u_0 - _ref_white.u_r) + v = 13. * lum * (v_0 - _ref_white.v_r) + + return stack(lum, u, v, *components[3:]) + + +def luv_to_xyz(pixels): + components = unstack(pixels) + lum, u, v = components[:3] + + u_0 = u / (13. * lum) + _ref_white.u_r + v_0 = v / (13. * lum) + _ref_white.v_r + + lum *= 100.0 # MMA specific + y = conditional(lum, lambda y: lum <= 8., + lambda y: _ref_white.xyz[1] * lum * ((3. / 29.) ** 3), + lambda y: _ref_white.xyz[1] * (((lum + 16.) / 116.) ** 3)) + x = y * (9. * u_0) / (4. * v_0) + z = y * (12. - 3. * u_0 - 20. * v_0) / (4. * v_0) + + return stack(x, y, z, *components[3:]) + + +def lch_to_lab(pixels): + components = unstack(pixels) + l, c, h = components[:3] + h *= 2 * pi # MMA specific + return stack(l, c * cos(h), c * sin(h), *components[3:]) + + +def lab_to_lch(pixels): + components = unstack(pixels) + l, a, b = components[:3] + h = arctan2(b, a) + h = conditional(h, lambda t: t < 0, + lambda t: t + 2 * pi, lambda t: t) + h /= 2 * pi # MMA specific + return stack(l, sqrt(a * a + b * b), h, *components[3:]) + + +def lab_to_xyz(pixels): + components = unstack(pixels) + + # see http://www.brucelindbloom.com/Eqn_Lab_to_XYZ.html + l, a, b = components[:3] + f_y = (l * 100. + 16.) / 116. + + x = conditional(a / 5. + f_y, lambda t: t > 0.008856, + lambda t: t ** 3, + lambda t: (116 * t - 16) / 903.3) + y = conditional(l * 100.0, lambda t: t > 903.3 * 0.008856, + lambda t: ((t + 16.) / 116.) ** 3, + lambda t: t / 903.3) + z = conditional(f_y - b / 2., lambda t: t > 0.008856, + lambda t: t ** 3, + lambda t: (116 * t - 16) / 903.3) + + x, y, z = [u * v for u, v in zip((x, y, z), _ref_white.xyz)] + + return stack(x, y, z, *components[3:]) + + +spaces = ('rgb', 'lab', 'luv') + +functions = ( + rgb_to_grayscale, + grayscale_to_rgb, + rgb_to_xyz, + rgb_to_hsb, + hsb_to_rgb, + cmyk_to_rgb, + rgb_to_cmyk, + xyz_to_rgb, + xyz_to_lab, + xyz_to_luv, + luv_to_xyz, + lch_to_lab, + lab_to_lch, + lab_to_xyz, +) + +_rewrites = { # see http://www.brucelindbloom.com/Math.html + 'XYZ>LCH': ('XYZ', 'LAB', 'LCH'), + 'LAB>LUV': ('LAB', 'XYZ', 'LUV'), + 'LAB>RGB': ('LAB', 'XYZ', 'RGB'), + 'LCH>XYZ': ('LCH', 'LAB', 'XYZ'), + 'LCH>LUV': ('LCH', 'LAB', 'XYZ', 'LUV'), + 'LCH>RGB': ('LCH', 'LAB', 'XYZ', 'RGB'), + 'LUV>LAB': ('LUV', 'XYZ', 'LAB'), + 'LUV>LCH': ('LUV', 'XYZ', 'LAB', 'LCH'), + 'LUV>RGB': ('LUV', 'XYZ', 'RGB'), + 'RGB>LAB': ('RGB', 'XYZ', 'LAB'), + 'RGB>LCH': ('RGB', 'XYZ', 'LAB', 'LCH'), + 'RGB>LUV': ('RGB', 'XYZ', 'LUV'), +} + +_conversions = { + 'Grayscale>RGB': grayscale_to_rgb, + 'RGB>Grayscale': rgb_to_grayscale, + 'CMYK>RGB': cmyk_to_rgb, + 'RGB>CMYK': rgb_to_cmyk, + 'RGB>HSB': rgb_to_hsb, + 'HSB>RGB': hsb_to_rgb, + 'XYZ>LAB': xyz_to_lab, + 'XYZ>LUV': xyz_to_luv, + 'XYZ>RGB': xyz_to_rgb, + 'LAB>XYZ': lab_to_xyz, + 'LAB>LCH': lab_to_lch, + 'LCH>LAB': lch_to_lab, + 'LUV>XYZ': luv_to_xyz, + 'RGB>XYZ': rgb_to_xyz, +} + +_rgb_transitions = set(['Grayscale', 'CMYK', 'HSB']) + + +def convert(components, src, dst): + if src == dst: + return components + + if src in _rgb_transitions or dst in _rgb_transitions: + flow = (src, 'RGB', dst) + else: + flow = (src, dst) + + new_flow = [flow[0]] + for s, d in (flow[i:i + 2] for i in range(len(flow) - 1)): + r = _rewrites.get('%s>%s' % (s, d)) + new_flow.extend(r[1:] if r else [d]) + + for s, d in (new_flow[i:i + 2] for i in range(len(new_flow) - 1)): + if s == d: + continue + func = _conversions.get('%s>%s' % (s, d)) + if not func: + return None + components = func(components) + + return components diff --git a/mathics/builtin/graphics.py b/mathics/builtin/graphics.py index b68ab32584..6a78be20c2 100644 --- a/mathics/builtin/graphics.py +++ b/mathics/builtin/graphics.py @@ -24,6 +24,7 @@ from mathics.core.expression import ( Expression, Integer, Real, String, Symbol, strip_context, system_symbols, system_symbols_dict, from_python) +from mathics.builtin.colors import convert as convert_color class CoordinatesError(BoxConstructError): @@ -370,44 +371,14 @@ def to_js(self): def to_expr(self): return Expression(self.get_name(), *self.components) - def to_cmyka(self): - rgba = self.to_rgba() - r, g, b = rgba[:3] - k = 1 - max(r, g, b) - k_ = 1 - k - return ((1 - r - k) / k_, (1 - g - k) / k_, (1 - b - k) / k_, k) + tuple(rgba[3:]) - - def to_ga(self): - # see https://en.wikipedia.org/wiki/Grayscale - rgba = self.to_rgba() - r, g, b = rgba[:3] - y = 0.299 * r + 0.587 * g + 0.114 * b # Y of Y'UV - return (y,) + tuple(rgba[3:]) - - def to_hsba(self): - return RGBColor(components=self.to_rgba()).to_hsba() - - def to_laba(self): - return XYZColor(components=self.to_xyza()).to_laba() - - def to_lcha(self): - # see http://www.brucelindbloom.com/Eqn_Lab_to_LCH.html - laba = self.to_laba() - l, a, b = laba[:3] - h = math.atan2(b, a) - if h < 0: - h += 2 * math.pi - h /= 2 * math.pi # MMA specific - return (l, math.sqrt(a * a + b * b), h) + tuple(laba[3:]) - - def to_luva(self): - return XYZColor(components=self.to_xyza()).to_luva() - def to_rgba(self): - return XYZColor(components=self.to_xyza()).to_rgba() + return self.to_color_space("RGB") - def to_xyza(self): - raise NotImplementedError + def to_color_space(self, color_space): + components = convert_color(self.components, self.color_space, color_space) + if components is None: + raise ValueError('cannot convert from color space %s to %s.' % (self.color_space, color_space)) + return components class RGBColor(_Color): @@ -428,52 +399,13 @@ class RGBColor(_Color): = StyleBox[GraphicsBox[...], ...] """ + color_space = 'RGB' components_sizes = [3, 4] default_components = [0, 0, 0, 1] def to_rgba(self): return self.components - def to_xyza(self): - components = _clip(*self.components) - - # inverse sRGB companding. see http://www.brucelindbloom.com/Eqn_RGB_to_XYZ.html - r, g, b = (math.pow(((c + 0.055) / 1.055), 2.4) - if c > 0.04045 else c / 12.92 for c in components[:3]) - - # use rRGB D50 conversion like MMA. see http://www.brucelindbloom.com/Eqn_RGB_XYZ_Matrix.html - # MMA seems to round matrix values to six significant digits. we do the same. - return _clip(0.436075 * r + 0.385065 * g + 0.14308 * b, - 0.222504 * r + 0.716879 * g + 0.0606169 * b, - 0.0139322 * r + 0.0971045 * g + 0.714173 * b) + tuple(components[3:]) - - def to_hsba(self): - # see https://en.wikipedia.org/wiki/HSB_color_space. HSB is also known as HSV. - - components = _clip(*self.components) - - r, g, b = components[:3] - m1 = max(r, g, b) - m0 = min(r, g, b) - c = m1 - m0 - - if c < 1e-15: - h = 0 - elif m1 == r: - h = ((g - b) / c) % 6 - elif m1 == g: - h = (b - r) / c + 2 - else: # m1 == b - h = (r - g) / c + 4 - h = (h * 60.) / 360. - if h < 0.: - h += 1. - - b = m1 - s = c / b - - return (h, s, b) + tuple(components[3:]) - class LABColor(_Color): """ @@ -484,31 +416,10 @@ class LABColor(_Color): """ + color_space = 'LAB' components_sizes = [3, 4] default_components = [0, 0, 0, 1] - def to_xyza(self): - components = self.components - - # see http://www.brucelindbloom.com/Eqn_Lab_to_XYZ.html - def inv_f(t): - if t > 0.008856: - return math.pow(t, 3) - else: - return (116 * t - 16) / 903.3 - - l, a, b = components[:3] - f_y = (l * 100. + 16.) / 116. - - x = inv_f(a / 5. + f_y) - y = math.pow(f_y, 3) if (l * 100.0) > 903.3 * 0.008856 else (l * 100.0) / 903.3 - z = inv_f(f_y - b / 2.) - - return tuple(c * w for c, w in zip((x, y, z), _ref_white.xyz)) + tuple(components[3:]) - - def to_laba(self): - return self.components - class LCHColor(_Color): """ @@ -519,23 +430,10 @@ class LCHColor(_Color): """ + color_space = 'LCH' components_sizes = [3, 4] default_components = [0, 0, 0, 1] - def to_xyza(self): - return LABColor(components=self.to_laba()).to_xyza() - - def to_laba(self): - lcha = self.to_lcha() - l, c, h = lcha[:3] - h *= 2 * math.pi # MMA specific - a = c * math.cos(h) - b = c * math.sin(h) - return (l, a, b) + tuple(lcha[3:]) - - def to_lcha(self): - return self.components - class LUVColor(_Color): """ @@ -545,27 +443,10 @@ class LUVColor(_Color): """ + color_space = 'LUV' components_sizes = [3, 4] default_components = [0, 0, 0, 1] - def to_xyza(self): - lum, u, v = self.components[:3] - - u_0 = u / (13. * lum) + _ref_white.u_r - v_0 = v / (13. * lum) + _ref_white.v_r - - lum *= 100.0 # MMA specific - - if lum <= 8.: - y = _ref_white.xyz[1] * lum * math.pow(3. / 29., 3) - else: - y = _ref_white.xyz[1] * math.pow((lum + 16.) / 116., 3) - - x = y * (9. * u_0) / (4. * v_0) - z = y * (12. - 3. * u_0 - 20. * v_0) / (4. * v_0) - - return _clip(x, y, z) + tuple(self.components[3:]) - class XYZColor(_Color): """ @@ -575,60 +456,10 @@ class XYZColor(_Color): """ + color_space = 'XYZ' components_sizes = [3, 4] default_components = [0, 0, 0, 1] - def to_rgba(self): - components = _clip(*self.components) - - x, y, z = components[:3] - - # for matrix, see http://www.brucelindbloom.com/Eqn_RGB_XYZ_Matrix.html - # MMA seems to round matrix values to six significant digits. we do the same. - r, g, b = [x * 3.13386 + y * -1.61687 + z * -0.490615, - x * -0.978768 + y * 1.91614 + z * 0.033454, - x * 0.0719453 + y * -0.228991 + z * 1.40524] - - return _clip(*[1.055 * math.pow(c, 1. / 2.4) - 0.055 if c > 0.0031308 - else c * 12.92 for c in (r, g, b)]) + tuple(components[3:]) - - def to_xyza(self): - return self.components - - def to_laba(self): - components = self.components[:3] - - # see http://www.brucelindbloom.com/Eqn_XYZ_to_Lab.html - - components = [c / w for c, w in zip(components, _ref_white.xyz)] - - x, y, z = (math.pow(t, 0.33333333) if t > 0.008856 - else (903.3 * t + 16.) / 116. for t in components) - - # MMA scales by 1/100 - return ((1.16 * y) - 0.16, 5. * (x - y), 2. * (y - z)) + tuple(components[3:]) - - def to_luva(self): - # see http://www.brucelindbloom.com/Eqn_XYZ_to_Luv.html - # and https://en.wikipedia.org/wiki/CIELUV - - components = _clip(*self.components) - - x_orig, y_orig, z_orig = components[:3] - y = y_orig / _ref_white.xyz[1] - - lum = 116. * math.pow(y, 1. / 3.) - 16. if y > 0.008856 else 903.3 * y - - q_0 = x_orig + 15. * y_orig + 3. * z_orig - u_0 = 4. * x_orig / q_0 - v_0 = 9. * y_orig / q_0 - - lum /= 100.0 # MMA specific - u = 13. * lum * (u_0 - _ref_white.u_r) - v = 13. * lum * (v_0 - _ref_white.v_r) - - return (lum, u, v) + tuple(components[3:]) - class CMYKColor(_Color): """ @@ -642,23 +473,10 @@ class CMYKColor(_Color): = -Graphics- """ + color_space = 'CMYK' components_sizes = [3, 4, 5] default_components = [0, 0, 0, 0, 1] - def to_cmyka(self): - return self.components - - def to_rgba(self): - k = self.components[3] if len(self.components) >= 4 else 0 - k_ = 1 - k - c = self.components - cmy = [c[0] * k_ + k, c[1] * k_ + k, c[2] * k_ + k] - rgb = (1 - cmy[0], 1 - cmy[1], 1 - cmy[2]) - return rgb + tuple(c[4:]) - - def to_xyza(self): - return RGBColor(components=self.to_rgba()).to_xyza() - class Hue(_Color): """ @@ -680,28 +498,10 @@ class Hue(_Color): = -Graphics- """ + color_space = 'HSB' components_sizes = [1, 2, 3, 4] default_components = [0, 1, 1, 1] - def to_rgba(self): - h, s, v = self.components[:3] - i = floor(6 * h) - f = 6 * h - i - i = i % 6 - p = v * (1 - s) - q = v * (1 - f * s) - t = v * (1 - (1 - f) * s) - - rgb = { - 0: (v, t, p), - 1: (q, v, p), - 2: (p, v, t), - 3: (p, q, v), - 4: (t, p, v), - 5: (v, p, q), - }[i] - return rgb + tuple(self.components[3:]) - def hsl_to_rgba(self): h, s, l = self.components[:3] if l < 0.5: @@ -732,12 +532,6 @@ def trans(t): result = tuple([trans(list(map(t))) for t in rgb]) + (self.components[3],) return result - def to_xyza(self): - return RGBColor(components=self.to_rgba()).to_xyza() - - def to_hsba(self): - return self.components - class GrayLevel(_Color): """ @@ -749,19 +543,11 @@ class GrayLevel(_Color):
represents a shade of gray specified by $g$ with opacity $a$. """ + + color_space = 'Grayscale' components_sizes = [1, 2] default_components = [0, 1] - def to_rgba(self): - g = self.components[0] - return (g, g, g) + tuple(self.components[1:]) - - def to_ga(self): - return self.components - - def to_xyza(self): - return RGBColor(components=self.to_rgba()).to_xyza() - class ColorConvert(Builtin): """ @@ -782,17 +568,6 @@ class ColorConvert(Builtin): XYZ: convert to XYZColor """ - _convert = { - 'CMYK': lambda color: CMYKColor(components=color.to_cmyka()), - 'Grayscale': lambda color: GrayLevel(components=color.to_ga()), - 'HSB': lambda color: Hue(components=color.to_hsba()), - 'LAB': lambda color: LABColor(components=color.to_laba()), - 'LCH': lambda color: LCHColor(components=color.to_lcha()), - 'LUV': lambda color: LUVColor(components=color.to_luva()), - 'RGB': lambda color: RGBColor(components=color.to_rgba()), - 'XYZ': lambda color: XYZColor(components=color.to_xyza()) - } - messages = { 'ccvinput': '`` should be a color.', 'imgcstype': '`` is not a valid color space.', @@ -806,13 +581,23 @@ def apply(self, color, colorspace, evaluation): evaluation.message('ColorConvert', 'ccvinput', color) return - convert = ColorConvert._convert.get(colorspace.get_string_value()) - if not convert: + py_colorspace = colorspace.get_string_value() + converted_components = convert_color(py_color.components, + py_color.color_space, + py_colorspace) + + if converted_components is None: evaluation.message('ColorConvert', 'imgcstype', colorspace) return - converted_color = convert(py_color) - return Expression(converted_color.get_name(), *converted_color.components) + if py_colorspace == 'Grayscale': + converted_color_name = 'GrayLevel' + elif py_colorspace == 'HSB': + converted_color_name = 'Hue' + else: + converted_color_name = py_colorspace + 'Color' + + return Expression(converted_color_name, *converted_components) class ColorDistance(Builtin): @@ -848,11 +633,11 @@ class ColorDistance(Builtin): } _distances = { - "CIE76": lambda c1, c2: _euclidean_distance(c1.to_laba()[:3], c2.to_laba()[:3]), - "CIE94": lambda c1, c2: _euclidean_distance(c1.to_lcha()[:3], c2.to_lcha()[:3]), - "DeltaL": lambda c1, c2: _component_distance(c1.to_lcha(), c2.to_lcha(), 0), - "DeltaC": lambda c1, c2: _component_distance(c1.to_lcha(), c2.to_lcha(), 1), - "DeltaH": lambda c1, c2: _component_distance(c1.to_lcha(), c2.to_lcha(), 2), + "CIE76": lambda c1, c2: _euclidean_distance(c1.to_color_space('LAB')[:3], c2.to_color_space('LAB')[:3]), + "CIE94": lambda c1, c2: _euclidean_distance(c1.to_color_space('LCH')[:3], c2.to_color_space('LCH')[:3]), + "DeltaL": lambda c1, c2: _component_distance(c1.to_color_space('LCH'), c2.to_color_space('LCH'), 0), + "DeltaC": lambda c1, c2: _component_distance(c1.to_color_space('LCH'), c2.to_color_space('LCH'), 1), + "DeltaH": lambda c1, c2: _component_distance(c1.to_color_space('LCH'), c2.to_color_space('LCH'), 2), } def apply(self, c1, c2, evaluation, options): @@ -865,7 +650,9 @@ def apply(self, c1, c2, evaluation, options): return else: def compute(a, b): - Expression(distance_function, a.to_laba(), b.to_laba()) + Expression(distance_function, + a.to_color_space('LAB'), + b.to_color_space('LAB')) def distance(a, b): try: diff --git a/mathics/builtin/numpy_utils.py b/mathics/builtin/numpy_utils.py index ee3bfa431c..100d7b8ae1 100644 --- a/mathics/builtin/numpy_utils.py +++ b/mathics/builtin/numpy_utils.py @@ -7,6 +7,8 @@ from mathics.core.expression import Expression from itertools import chain +from functools import reduce +from math import sin as sinf, cos as cosf, sqrt as sqrtf, atan2 as atan2f, floor as floorf try: import numpy @@ -31,15 +33,6 @@ def py_instantiate_elements(a, new_element, d=1): leaves = [py_instantiate_elements(e, new_element, d) for e in a] return Expression('List', *leaves) - -def py_stack_along_inner_axis(a): - # explanation of functionality: see numpy version of _stack_array above - - if not isinstance(a[0], list): - return list(chain(a)) - else: - return [py_stack_along_inner_axis([x[i] for x in a]) for i in range(len(a[0]))] - if _numpy: def instantiate_elements(a, new_element, d=1): # given a numpy array 'a' and a python element constructor 'new_element', generate a python array of the @@ -52,7 +45,14 @@ def instantiate_elements(a, new_element, d=1): leaves = [instantiate_elements(e, new_element, d) for e in a] return Expression('List', *leaves) - def stack_along_inner_axis(a): + def array(a): + return numpy.array(a) + + def unstack(a): + a = array(a) + return array(numpy.split(a, a.shape[-1], axis=-1)) + + def stack(*a): # numpy.stack with axis=-1 stacks arrays along the most inner axis: # e.g. numpy.stack([ [1, 2], [3, 4] ], axis=-1) @@ -61,10 +61,156 @@ def stack_along_inner_axis(a): # e.g. numpy.stack([ [[1, 2], [3, 4]], [[4, 5], [6, 7]] ], axis=-1) # gives: array([[[1, 4], [2, 5]], [[3, 6], [4, 7]]]) - return numpy.stack(a, axis=-1) + a = array(a) + b = numpy.stack(a, axis=-1) + + if a.shape[-1] == 1 and b.shape[0] > 0: # e.g. [[a], [b], [c]] + b = b[0] # makes stack(unstack(x)) == x + return b + + def concat(*a): + a = [x for x in a if x.shape[0]] # skip empty + return numpy.concatenate(a, axis=-1) + + def conditional(a, cond, t, f): + b = array(a)[:] + mask = cond(a) + b[mask] = t(b[mask]) + b[~mask] = f(b[~mask]) + return b + + def switch(*a): + assert a and len(a) % 2 == 0 + b = numpy.ndarray(a[0].shape) + # we apply the rules in reversed order, so that the first rules + # always make the last changes, which corresponds to the unreversed + # processing in the non-vectorized (non-numpy) implementation. + for mask, v in reversed([a[i:i + 2] for i in range(0, len(a), 2)]): + b[mask] = v(lambda x: x[mask]) + return b + + def choose(i, *options): + assert 0 < len(options) < 256 + i_int = i.astype(numpy.uint8) + dim = len(options[0]) + return [numpy.choose(i_int, [o[d] for o in options]) for d in range(dim)] + + def clip(a, t0, t1): + return numpy.clip(array(a), t0, t1) + + def dot_t(u, v): + return numpy.dot(array(u), array(v).T) + + def mod(a, b): + return numpy.mod(a, b) + + def sin(a): + return numpy.sin(array(a)) + + def cos(a): + return numpy.cos(array(a)) + + def arctan2(y, x): + return numpy.arctan2(array(y), array(x)) + + def sqrt(a): + return numpy.sqrt(array(a)) + + def floor(a): + return numpy.floor(array(a)) + + def maximum(*a): + return reduce(numpy.maximum, [array(x) for x in a]) + + def minimum(*a): + return reduce(numpy.minimum, [array(x) for x in a]) else: # If numpy is not available, we define the following fallbacks that are useful for implementing a similar # logic in pure python without numpy. They obviously work on regular python array though, not numpy arrays. instantiate_elements = py_instantiate_elements - stack_along_inner_axis = py_stack_along_inner_axis + + def array(a): + return a + + def _is_bottom(a): + return any(not isinstance(x, list) for x in a) + + def unstack(a): + if not a: + return [] + + def length(b): + return max(length(x) for x in b) if not _is_bottom(b) else len(b) + + def split(b, i): + if not _is_bottom(b): + return [split(x, i) for x in b] + else: + return b[i] + + return [split(a, i) for i in range(length(a))] + + def stack(*a): + if _is_bottom(a): + return list(chain(a)) + else: + return [stack(*[x[i] for x in a]) for i in range(len(a[0]))] + + def concat(a, b): + return stack(*(unstack(a) + unstack(b))) + + def _apply(a, f): + if isinstance(a, (list, tuple)): + return [_apply(t, f) for t in a] + else: + return f(a) + + def _apply_n(*p, f): + if isinstance(p[0], list): + return [_apply_n(*q, f=f) for q in zip(*p)] + else: + return f(*p) + + def conditional(a, cond, t, f): + return _apply(a, lambda x: t(x) if cond(x) else f(x)) + + def switch(*a): + assert a and len(a) % 2 == 0 + for cond, v in (a[i:i + 2] for i in range(0, len(a), 2)): + if cond: + return v(lambda x: x) + raise ValueError('no matching case in switch') + + def choose(i, *options): + return options[i] + + def clip(a, t0, t1): + return _apply(a, lambda x: max(t0, min(t1, x))) + + def dot_t(u, v): + return [sum(x * y for x, y in zip(u, r)) for r in v] + + def mod(a, b): + return _apply_n(a, b, f=lambda x, y: x % y) + + def sin(a): + return _apply(a, lambda t: sinf(t)) + + def cos(a): + return _apply(a, lambda t: cosf(t)) + + def arctan2(y, x): + return _apply_n(y, x, f=atan2f) + + def sqrt(a): + return _apply(a, lambda t: sqrtf(t)) + + def floor(a): + return _apply(a, lambda t: floorf(t)) + + def maximum(*a): + return _apply_n(*a, f=lambda *x: max(*x)) + + def minimum(*a): + return _apply_n(*a, f=lambda *x: min(*x)) diff --git a/mathics/builtin/randomnumbers.py b/mathics/builtin/randomnumbers.py index 1446d473b6..07bb411320 100644 --- a/mathics/builtin/randomnumbers.py +++ b/mathics/builtin/randomnumbers.py @@ -20,7 +20,7 @@ from functools import reduce from mathics.builtin.base import Builtin -from mathics.builtin.numpy_utils import instantiate_elements, stack_along_inner_axis +from mathics.builtin.numpy_utils import instantiate_elements, stack from mathics.core.expression import (Integer, String, Symbol, Real, Expression, Complex) @@ -502,7 +502,7 @@ def apply_list(self, zmin, zmax, ns, evaluation): real = rand.randreal(min_value.real, max_value.real, py_ns) imag = rand.randreal(min_value.imag, max_value.imag, py_ns) return instantiate_elements( - stack_along_inner_axis([real, imag]), + stack(real, imag), lambda c: Complex(Real(c[0]), Real(c[1])), d=2) From 5b2821b4014b67d0359ac797c9584584427f72df Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Tue, 14 Jun 2016 22:56:16 +0200 Subject: [PATCH 14/30] should fix pre-Python 3.5 --- mathics/builtin/colors.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/mathics/builtin/colors.py b/mathics/builtin/colors.py index 48f3ca854e..f4c03cd314 100644 --- a/mathics/builtin/colors.py +++ b/mathics/builtin/colors.py @@ -65,7 +65,7 @@ def rgb_to_xyz(pixels): lambda t: t / 12.92) xyz = dot_t(stack(x, y, z), xyz_from_rgb) - return concat(clip(xyz, 0, 1), components[3:]) + return clip(concat(xyz, components[3:]), 0, 1) def rgb_to_hsb(pixels): @@ -127,9 +127,9 @@ def cmyk_to_rgb(pixels): k_ = 1 - k cmy = [(x * k_ + k) for x in c[:3]] - rgb = [1 - x for x in cmy] + r, g, b = [1 - x for x in cmy] - return stack(*rgb, *c[4:]) + return concat(stack(r, g, b), c[4:]) def rgb_to_cmyk(pixels): From 00bb6f3d802f501d25477793a955ff34b6915a11 Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Tue, 14 Jun 2016 23:07:56 +0200 Subject: [PATCH 15/30] simpler calculation of conversion flows --- mathics/builtin/colors.py | 30 ++++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/mathics/builtin/colors.py b/mathics/builtin/colors.py index f4c03cd314..6f5013e456 100644 --- a/mathics/builtin/colors.py +++ b/mathics/builtin/colors.py @@ -4,6 +4,7 @@ from math import pi from mathics.builtin.numpy_utils import stack, unstack, concat, array, clip, conditional, switch, choose from mathics.builtin.numpy_utils import sqrt, floor, mod, cos, sin, arctan2, minimum, maximum, dot_t +from itertools import chain # use rRGB D50 conversion like MMA. see http://www.brucelindbloom.com/Eqn_RGB_XYZ_Matrix.html # MMA seems to round matrix values to six significant digits. we do the same. @@ -270,7 +271,7 @@ def lab_to_xyz(pixels): lab_to_xyz, ) -_rewrites = { # see http://www.brucelindbloom.com/Math.html +_flows = { # see http://www.brucelindbloom.com/Math.html 'XYZ>LCH': ('XYZ', 'LAB', 'LCH'), 'LAB>LUV': ('LAB', 'XYZ', 'LUV'), 'LAB>RGB': ('LAB', 'XYZ', 'RGB'), @@ -285,6 +286,18 @@ def lab_to_xyz(pixels): 'RGB>LUV': ('RGB', 'XYZ', 'LUV'), } + +_rgb_flows = set(['Grayscale', 'CMYK', 'HSB']) + + +def _flow(src, dst): + if (src in _rgb_flows and dst != 'RGB') or (dst in _rgb_flows and src != 'RGB'): + return list(chain(_flow(src, 'RGB'), _flow('RGB', dst))) + else: + r = _flows.get('%s>%s' % (src, dst)) + return list(r) if r else [src, dst] + + _conversions = { 'Grayscale>RGB': grayscale_to_rgb, 'RGB>Grayscale': rgb_to_grayscale, @@ -302,24 +315,13 @@ def lab_to_xyz(pixels): 'RGB>XYZ': rgb_to_xyz, } -_rgb_transitions = set(['Grayscale', 'CMYK', 'HSB']) - def convert(components, src, dst): if src == dst: return components - if src in _rgb_transitions or dst in _rgb_transitions: - flow = (src, 'RGB', dst) - else: - flow = (src, dst) - - new_flow = [flow[0]] - for s, d in (flow[i:i + 2] for i in range(len(flow) - 1)): - r = _rewrites.get('%s>%s' % (s, d)) - new_flow.extend(r[1:] if r else [d]) - - for s, d in (new_flow[i:i + 2] for i in range(len(new_flow) - 1)): + flows = _flow(src, dst) + for s, d in (flows[i:i + 2] for i in range(len(flows) - 1)): if s == d: continue func = _conversions.get('%s>%s' % (s, d)) From bd5de4cc4df3babef08b53ca38f1c529cb766b82 Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Tue, 14 Jun 2016 23:19:44 +0200 Subject: [PATCH 16/30] more Python 2 fixes --- mathics/builtin/numpy_utils.py | 33 +++++++++++++++++++-------------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/mathics/builtin/numpy_utils.py b/mathics/builtin/numpy_utils.py index 100d7b8ae1..26a1a3c669 100644 --- a/mathics/builtin/numpy_utils.py +++ b/mathics/builtin/numpy_utils.py @@ -9,6 +9,7 @@ from itertools import chain from functools import reduce from math import sin as sinf, cos as cosf, sqrt as sqrtf, atan2 as atan2f, floor as floorf +import operator try: import numpy @@ -160,20 +161,22 @@ def stack(*a): def concat(a, b): return stack(*(unstack(a) + unstack(b))) - def _apply(a, f): + def _apply(f, a): if isinstance(a, (list, tuple)): - return [_apply(t, f) for t in a] + return [_apply(f, t) for t in a] else: return f(a) - def _apply_n(*p, f): + def _apply_n(f, *p): if isinstance(p[0], list): - return [_apply_n(*q, f=f) for q in zip(*p)] + return [_apply_n(f, *q) for q in zip(*p)] else: return f(*p) def conditional(a, cond, t, f): - return _apply(a, lambda x: t(x) if cond(x) else f(x)) + def _eval(x): + return t(x) if cond(x) else f(x) + return _apply(_eval, a) def switch(*a): assert a and len(a) % 2 == 0 @@ -186,31 +189,33 @@ def choose(i, *options): return options[i] def clip(a, t0, t1): - return _apply(a, lambda x: max(t0, min(t1, x))) + def _eval(x): + return max(t0, min(t1, x)) + return _apply(_eval, a) def dot_t(u, v): return [sum(x * y for x, y in zip(u, r)) for r in v] def mod(a, b): - return _apply_n(a, b, f=lambda x, y: x % y) + return _apply_n(operator.mod, a, b) def sin(a): - return _apply(a, lambda t: sinf(t)) + return _apply(sinf, a) def cos(a): - return _apply(a, lambda t: cosf(t)) + return _apply(cosf, a) def arctan2(y, x): - return _apply_n(y, x, f=atan2f) + return _apply_n(atan2f, y, x) def sqrt(a): - return _apply(a, lambda t: sqrtf(t)) + return _apply(sqrtf, a) def floor(a): - return _apply(a, lambda t: floorf(t)) + return _apply(floorf, a) def maximum(*a): - return _apply_n(*a, f=lambda *x: max(*x)) + return _apply_n(max, *a) def minimum(*a): - return _apply_n(*a, f=lambda *x: min(*x)) + return _apply_n(min, *a) From f550a9c4ac7e9c4f7dfaccfefdf6b7cbd7288e54 Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Tue, 14 Jun 2016 23:46:34 +0200 Subject: [PATCH 17/30] fixes PyPy --- mathics/builtin/numpy_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mathics/builtin/numpy_utils.py b/mathics/builtin/numpy_utils.py index 26a1a3c669..e21090711c 100644 --- a/mathics/builtin/numpy_utils.py +++ b/mathics/builtin/numpy_utils.py @@ -186,7 +186,7 @@ def switch(*a): raise ValueError('no matching case in switch') def choose(i, *options): - return options[i] + return options[int(i)] # int cast needed for PyPy def clip(a, t0, t1): def _eval(x): From 910f912fb21715518f8bca21c6bd5eeae382d6dc Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Wed, 15 Jun 2016 07:37:17 +0200 Subject: [PATCH 18/30] bug fixes and test cases for image color space conversions --- mathics/builtin/colors.py | 46 ++++++++++------------------------ mathics/builtin/numpy_utils.py | 31 ++++++++++++++++++----- test/test_color.py | 24 ++++++++++++++++++ 3 files changed, 62 insertions(+), 39 deletions(-) diff --git a/mathics/builtin/colors.py b/mathics/builtin/colors.py index 6f5013e456..fad8b1edd4 100644 --- a/mathics/builtin/colors.py +++ b/mathics/builtin/colors.py @@ -2,7 +2,7 @@ # -*- coding: utf-8 -*- from math import pi -from mathics.builtin.numpy_utils import stack, unstack, concat, array, clip, conditional, switch, choose +from mathics.builtin.numpy_utils import vectorize, stack, unstack, concat, array, clip, conditional, switch, choose from mathics.builtin.numpy_utils import sqrt, floor, mod, cos, sin, arctan2, minimum, maximum, dot_t from itertools import chain @@ -180,9 +180,9 @@ def xyz_to_luv(pixels): x_orig, y_orig, z_orig = xyz y = y_orig / _ref_white.xyz[1] - lum = conditional(y, lambda y: y > 0.008856, - lambda y: 116. * (y ** (1. / 3.)) - 16, - lambda y: 903.3 * y) + lum = conditional(y, lambda t: t > 0.008856, + lambda t: 116. * (t ** (1. / 3.)) - 16, + lambda t: 903.3 * t) q_0 = x_orig + 15. * y_orig + 3. * z_orig u_0 = 4. * x_orig / q_0 @@ -204,9 +204,9 @@ def luv_to_xyz(pixels): v_0 = v / (13. * lum) + _ref_white.v_r lum *= 100.0 # MMA specific - y = conditional(lum, lambda y: lum <= 8., - lambda y: _ref_white.xyz[1] * lum * ((3. / 29.) ** 3), - lambda y: _ref_white.xyz[1] * (((lum + 16.) / 116.) ** 3)) + y = conditional(lum, lambda t: t <= 8., + lambda t: _ref_white.xyz[1] * t * ((3. / 29.) ** 3), + lambda t: _ref_white.xyz[1] * (((t + 16.) / 116.) ** 3)) x = y * (9. * u_0) / (4. * v_0) z = y * (12. - 3. * u_0 - 20. * v_0) / (4. * v_0) @@ -223,10 +223,9 @@ def lch_to_lab(pixels): def lab_to_lch(pixels): components = unstack(pixels) l, a, b = components[:3] - h = arctan2(b, a) - h = conditional(h, lambda t: t < 0, - lambda t: t + 2 * pi, lambda t: t) - h /= 2 * pi # MMA specific + h = conditional(arctan2(b, a), lambda t: t < 0., + lambda t: t + 2. * pi, lambda t: t) + h /= 2. * pi # MMA specific return stack(l, sqrt(a * a + b * b), h, *components[3:]) @@ -252,25 +251,6 @@ def lab_to_xyz(pixels): return stack(x, y, z, *components[3:]) -spaces = ('rgb', 'lab', 'luv') - -functions = ( - rgb_to_grayscale, - grayscale_to_rgb, - rgb_to_xyz, - rgb_to_hsb, - hsb_to_rgb, - cmyk_to_rgb, - rgb_to_cmyk, - xyz_to_rgb, - xyz_to_lab, - xyz_to_luv, - luv_to_xyz, - lch_to_lab, - lab_to_lch, - lab_to_xyz, -) - _flows = { # see http://www.brucelindbloom.com/Math.html 'XYZ>LCH': ('XYZ', 'LAB', 'LCH'), 'LAB>LUV': ('LAB', 'XYZ', 'LUV'), @@ -298,7 +278,7 @@ def _flow(src, dst): return list(r) if r else [src, dst] -_conversions = { +conversions = { 'Grayscale>RGB': grayscale_to_rgb, 'RGB>Grayscale': rgb_to_grayscale, 'CMYK>RGB': cmyk_to_rgb, @@ -324,9 +304,9 @@ def convert(components, src, dst): for s, d in (flows[i:i + 2] for i in range(len(flows) - 1)): if s == d: continue - func = _conversions.get('%s>%s' % (s, d)) + func = conversions.get('%s>%s' % (s, d)) if not func: return None - components = func(components) + components = vectorize(func, components) return components diff --git a/mathics/builtin/numpy_utils.py b/mathics/builtin/numpy_utils.py index e21090711c..1355ec99a6 100644 --- a/mathics/builtin/numpy_utils.py +++ b/mathics/builtin/numpy_utils.py @@ -49,9 +49,22 @@ def instantiate_elements(a, new_element, d=1): def array(a): return numpy.array(a) + def vectorize(f, a): + # if a is a flat array, we embed the components in an array, i.e. + # [[a, b, c]], so that unstack will yield [a], [b], [c], and we operate on + # arrays of length 1 insteaf of scalars. + a = array(a) + if len(a.shape) == 1: + return f([a])[0] + else: + return f(a) + def unstack(a): a = array(a) - return array(numpy.split(a, a.shape[-1], axis=-1)) + b = array(numpy.split(a, a.shape[-1], axis=-1)) + if b.shape[-1] == 1: + b = b.reshape(b.shape[:-1]) + return b def stack(*a): # numpy.stack with axis=-1 stacks arrays along the most inner axis: @@ -64,9 +77,6 @@ def stack(*a): a = array(a) b = numpy.stack(a, axis=-1) - - if a.shape[-1] == 1 and b.shape[0] > 0: # e.g. [[a], [b], [c]] - b = b[0] # makes stack(unstack(x)) == x return b def concat(*a): @@ -131,11 +141,20 @@ def minimum(*a): instantiate_elements = py_instantiate_elements + def _is_bottom(a): + return any(not isinstance(x, list) for x in a) + def array(a): return a - def _is_bottom(a): - return any(not isinstance(x, list) for x in a) + def vectorize(f, a): + # we work on a scalar level, i.e. we always pass flat arrays like + # [a, b, c, ...] to f, and unstack will resolve this + # into a, b, c there. + if _is_bottom(a): + return f(a) + else: + return [vectorize(f, x) for x in a] def unstack(a): if not a: diff --git a/test/test_color.py b/test/test_color.py index 73fc4c9214..4e0a0ce215 100644 --- a/test/test_color.py +++ b/test/test_color.py @@ -9,9 +9,12 @@ import pexpect import unittest from six.moves import range +from random import random from mathics.core.expression import Expression, Integer, Rational, Symbol from mathics.core.definitions import Definitions from mathics.core.evaluation import Evaluation +import mathics.builtin.colors as colors +from mathics.builtin.numpy_utils import array, vectorize class ColorTest(unittest.TestCase): @@ -66,6 +69,12 @@ def testConversions(self): self._checkConversion("XYZ", (0.1, 0.1, 0.1), "LAB", (0.3784243088316416, 0.02835852516741566, -0.06139347105996884)) + def testImageConversions(self): + # test that f([x, y, ...]) = [f(x), f(y), ...] for rectangular image arrays. + + for convert in colors.conversions.values(): + self._checkImageConversion(4, convert) + def _checkConversion(self, from_space, from_components, to_space, to_components): places = 12 @@ -80,5 +89,20 @@ def _checkConversion(self, from_space, from_components, to_space, to_components) for x, y in zip(components, to_components): self.assertAlmostEqual(x, y, places) + def _checkImageConversion(self, size, convert): + pixels = [[random(), random(), random()] for _ in range(size * size)] + refs = [vectorize(convert, p) for p in pixels] + + image = [[pixels[x * size + y] for y in range(size)] for x in range(size)] + image = vectorize(convert, array(image)) + + for x in range(size): + for y in range(size): + p1 = image[x][y] + p2 = refs[x * size + y] + self.assertEqual(len(p1), len(p2)) + for a, b in zip(p1, p2): + self.assertAlmostEqual(a, b, 12) + if __name__ == "__main__": unittest.main() From ca4b4e7628f6fadf80b07222ff4f2bea0270022c Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Wed, 15 Jun 2016 19:32:57 +0200 Subject: [PATCH 19/30] new numpy utils cleanup --- mathics/builtin/colors.py | 6 +-- mathics/builtin/numpy_utils.py | 82 ++++++++++++++++++++++++---------- test/test_color.py | 6 +-- 3 files changed, 64 insertions(+), 30 deletions(-) diff --git a/mathics/builtin/colors.py b/mathics/builtin/colors.py index fad8b1edd4..729eeadf4e 100644 --- a/mathics/builtin/colors.py +++ b/mathics/builtin/colors.py @@ -2,7 +2,7 @@ # -*- coding: utf-8 -*- from math import pi -from mathics.builtin.numpy_utils import vectorize, stack, unstack, concat, array, clip, conditional, switch, choose +from mathics.builtin.numpy_utils import vectorized, stack, unstack, concat, array, clip, conditional, compose, choose from mathics.builtin.numpy_utils import sqrt, floor, mod, cos, sin, arctan2, minimum, maximum, dot_t from itertools import chain @@ -78,7 +78,7 @@ def rgb_to_hsb(pixels): m0 = minimum(r, g, b) c = m1 - m0 - h = switch( + h = compose( c < 1e-15, lambda s: 0, m1 == r, @@ -307,6 +307,6 @@ def convert(components, src, dst): func = conversions.get('%s>%s' % (s, d)) if not func: return None - components = vectorize(func, components) + components = vectorized(func, components, 1) return components diff --git a/mathics/builtin/numpy_utils.py b/mathics/builtin/numpy_utils.py index 1355ec99a6..99e2277b7a 100644 --- a/mathics/builtin/numpy_utils.py +++ b/mathics/builtin/numpy_utils.py @@ -49,13 +49,14 @@ def instantiate_elements(a, new_element, d=1): def array(a): return numpy.array(a) - def vectorize(f, a): - # if a is a flat array, we embed the components in an array, i.e. - # [[a, b, c]], so that unstack will yield [a], [b], [c], and we operate on - # arrays of length 1 insteaf of scalars. + def vectorized(f, a, depth): + # we ignore depth, as numpy arrays will handle the vectorized cases. a = array(a) if len(a.shape) == 1: - return f([a])[0] + # if a is a flat array, we embed the components in an array, i.e. + # [[a, b, c]], so that unstack will yield [a], [b], [c], and we + # operate on arrays of length 1 instead of scalars. + return f(array([a]))[0] else: return f(a) @@ -80,6 +81,7 @@ def stack(*a): return b def concat(*a): + a = [array(x) for x in a] a = [x for x in a if x.shape[0]] # skip empty return numpy.concatenate(a, axis=-1) @@ -90,7 +92,7 @@ def conditional(a, cond, t, f): b[~mask] = f(b[~mask]) return b - def switch(*a): + def compose(*a): assert a and len(a) % 2 == 0 b = numpy.ndarray(a[0].shape) # we apply the rules in reversed order, so that the first rules @@ -100,12 +102,6 @@ def switch(*a): b[mask] = v(lambda x: x[mask]) return b - def choose(i, *options): - assert 0 < len(options) < 256 - i_int = i.astype(numpy.uint8) - dim = len(options[0]) - return [numpy.choose(i_int, [o[d] for o in options]) for d in range(dim)] - def clip(a, t0, t1): return numpy.clip(array(a), t0, t1) @@ -135,6 +131,20 @@ def maximum(*a): def minimum(*a): return reduce(numpy.minimum, [array(x) for x in a]) + + def choose(i, *options): + assert options + dim = len(options[0]) + columns = [[o[d] for o in options] for d in range(dim)] + if isinstance(i, (int, float)): + return [column[int(i)] for column in columns] # int cast needed for PyPy + else: + assert len(options) < 256 + i_int = array(i).astype(numpy.uint8) + return [numpy.choose(i_int, column) for column in columns] + + def allclose(a, b): + return numpy.allclose(array(a), array(b)) else: # If numpy is not available, we define the following fallbacks that are useful for implementing a similar # logic in pure python without numpy. They obviously work on regular python array though, not numpy arrays. @@ -147,14 +157,16 @@ def _is_bottom(a): def array(a): return a - def vectorize(f, a): - # we work on a scalar level, i.e. we always pass flat arrays like - # [a, b, c, ...] to f, and unstack will resolve this - # into a, b, c there. + def vectorized(f, a, depth): + # depth == 0 means: f will get only scalars. depth == 1 means: f will + # get lists of scalars. if _is_bottom(a): - return f(a) + if depth == 0: + return [f(x) for x in a] + else: + return f(a) else: - return [vectorize(f, x) for x in a] + return [vectorized(f, x, depth) for x in a] def unstack(a): if not a: @@ -197,15 +209,12 @@ def _eval(x): return t(x) if cond(x) else f(x) return _apply(_eval, a) - def switch(*a): + def compose(*a): assert a and len(a) % 2 == 0 for cond, v in (a[i:i + 2] for i in range(0, len(a), 2)): if cond: return v(lambda x: x) - raise ValueError('no matching case in switch') - - def choose(i, *options): - return options[int(i)] # int cast needed for PyPy + raise ValueError('no matching case in compose') def clip(a, t0, t1): def _eval(x): @@ -213,7 +222,10 @@ def _eval(x): return _apply(_eval, a) def dot_t(u, v): - return [sum(x * y for x, y in zip(u, r)) for r in v] + if not isinstance(v[0], list): + return sum(x * y for x, y in zip(u, v)) + else: + return [sum(x * y for x, y in zip(u, r)) for r in v] def mod(a, b): return _apply_n(operator.mod, a, b) @@ -238,3 +250,25 @@ def maximum(*a): def minimum(*a): return _apply_n(min, *a) + + def _choose_descend(i, options): + if isinstance(i, (int, float)): + return options[int(i)] # int cast needed for PyPy + else: + return [_choose_descend(next_i, [o[k] for o in options]) for k, next_i in enumerate(i)] + + def choose(i, *options): + assert options + dim = len(options[0]) + columns = [[o[d] for o in options] for d in range(dim)] + return [_choose_descend(i, column) for column in columns] + + def allclose(a, b): + if isinstance(a, list) and isinstance(b, list): + if len(a) != len(b): + return False + return all(allclose(x, y) for x, y in zip(a, b)) + elif isinstance(a, list) or isinstance(b, list): + return False + else: + return abs(a - b) < 1e-12 diff --git a/test/test_color.py b/test/test_color.py index 4e0a0ce215..15e9b9afc4 100644 --- a/test/test_color.py +++ b/test/test_color.py @@ -14,7 +14,7 @@ from mathics.core.definitions import Definitions from mathics.core.evaluation import Evaluation import mathics.builtin.colors as colors -from mathics.builtin.numpy_utils import array, vectorize +from mathics.builtin.numpy_utils import array, vectorized class ColorTest(unittest.TestCase): @@ -91,10 +91,10 @@ def _checkConversion(self, from_space, from_components, to_space, to_components) def _checkImageConversion(self, size, convert): pixels = [[random(), random(), random()] for _ in range(size * size)] - refs = [vectorize(convert, p) for p in pixels] + refs = [vectorized(convert, p, 1) for p in pixels] image = [[pixels[x * size + y] for y in range(size)] for x in range(size)] - image = vectorize(convert, array(image)) + image = vectorized(convert, array(image), 1) for x in range(size): for y in range(size): From 6f5495056dc28ff70349797500346b8428c0ef65 Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Wed, 15 Jun 2016 19:34:19 +0200 Subject: [PATCH 20/30] test cases for (new) numpy utils --- test/test_numpy_utils.py | 169 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 169 insertions(+) create mode 100644 test/test_numpy_utils.py diff --git a/test/test_numpy_utils.py b/test/test_numpy_utils.py new file mode 100644 index 0000000000..2c51fa8d2e --- /dev/null +++ b/test/test_numpy_utils.py @@ -0,0 +1,169 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +from __future__ import absolute_import +from __future__ import unicode_literals + +from mathics.builtin.numpy_utils import allclose, stack, unstack, concat, conditional, compose, clip, array, choose +from mathics.builtin.numpy_utils import vectorized, minimum, maximum, sin, dot_t, mod, floor, sqrt +import sys +import unittest + + +class NumpyUtils(unittest.TestCase): + def testUnstack(self): + # flat lists remain flat lists. + self.assertEqualArrays( + unstack([1, 2, 3]), + [1, 2, 3]) + + # lists of lists get unstacked. + self.assertEqualArrays( + unstack([[1, 2], [3, 4]]), + [[1, 3], [2, 4]]) + + # matrices stay matrices, e.g. each r, g, b + # components is split into grayscale images. + self.assertEqualArrays( + unstack([ + [[1, 2, 3], [4, 5, 6]], + [[7, 8, 9], [10, 11, 12]] + ]), + [ + [[1, 4], [7, 10]], + [[2, 5], [8, 11]], + [[3, 6], [9, 12]] + ]) + + def testStackUnstackIdentity(self): + a = [[[1, 2, 3], [4, 5, 6]], [[7, 8, 9], [10, 11, 12]]] + b = [1, 2, 3] + c = [[1, 2], [3, 4]] + + for m in (a, b, c): + self.assertEqualArrays(stack(*unstack(m)), m) + + def testConcatSimple(self): + # concat concatenates arrays. + self.assertEqualArrays(concat([1], [2]), [1, 2]) + + def testConcatComplex(self): + # concat concatenates the most inner axis. + a = [[[1, 2], [4, 5]], [[7, 8], [10, 11]]] + b = [[[3], [6]], [[9], [12]]] + c = [[[1, 2, 3], [4, 5, 6]], [[7, 8, 9], [10, 11, 12]]] + self.assertEqualArrays(concat(a, b), c) + + def testConditional(self): + # conditional operates on each element independently. + a = array([[[0.1, 0.6], [1.8, 0.4]], [[-0.1, -0.8], [1.1, 0.5]]]) + a = conditional(a, lambda t: t > 0.5, # if + lambda t: t + 1., # true + lambda t: -t) # false + self.assertEqualArrays(a, [[[-0.1, 1.6], [2.8, -0.4]], [[0.1, 0.8], [2.1, -0.5]]]) + + def testCompose(self): + a = array([[[1, 2], [4, 5]], [[7, 8], [10, 11]]]) + + def f(a): + return compose( + a > 10, + lambda s: s(a * 10 + 1), + a > 3, + lambda s: s(a * 10), + a < 3, + lambda s: -1) + + a = vectorized(f, a, 0) + self.assertEqualArrays(a, [[[-1, -1], [40, 50]], [[70, 80], [100, 111]]]) + + def testChooseSimple(self): + # select a single value from a list of values. + options = [ + [0, 1, 2], + [3, 4, 5], + [6, 7, 8] + ] + + for i in range(len(options)): + self.assertEqual(choose(i, *options), options[i]) + + def testChooseComplex(self): + def m(i): + return [[10 * i + 0, 10 * i + 1], + [10 * i + 2, 10 * i + 3]] + + selector = [[0, 1], + [1, 2]] + + a = choose(selector, + (m(1), m(2), m(3)), + (m(4), m(5), m(6)), + (m(7), m(8), m(9))) + + # choose() operates on column after column of the options matrix, i.e. + # in the example above, the selector is first applied to (m(1), m(4), m(7)), + # then to (m(2), m(5), m(8)), and so on. let's calculate the result for the + # first column: + + # selector (integers indicate from which option to take a value): + # [0 1] + # [1 2] + + # option #0 / column 1 (i.e. m(1)): + # [10 11] + # [12 13] + + # option #1 / column 1 (i.e. m(4)): + # [40 41] + # [42 43] + + # option #2 / column 1 (i.e. m(7)): + # [70 71] + # [72 73] + + # choose() now picks the right value from each options depending on selector: + + # [10 41] + # [42 73] + + self.assertEqual(len(a), 3) + self.assertEqualArrays(a[0], [[10, 41], [42, 73]]) + self.assertEqualArrays(a[1], [[20, 51], [52, 83]]) + self.assertEqualArrays(a[2], [[30, 61], [62, 93]]) + + def testClip(self): + a = array([[[-0.1, 0.6], [-1.8, -0.4]], [[0.1, 0.8], [1.1, 0.5]]]) + a = clip(a, 0, 1) + self.assertEqualArrays(a, [[[0., 0.6], [0., 0.]], [[0.1, 0.8], [1., 0.5]]]) + + def testDot(self): + self.assertEqual(dot_t([1, 2, 3], [4, 5, 6]), 32) + self.assertEqualArrays(dot_t([1, 2, 3], [[4, 5, 6], [7, 8, 9], [10, 11, 12]]), [32, 50, 68]) + + def testMod(self): + self.assertEqualArrays(mod(array([[10, 20], [30, 40]]), [[7, 7], [7, 7]]), [[3, 6], [2, 5]]) + + def testMaximum(self): + self.assertEqualArrays(maximum( + [[1, 2], [3, 4]], + [[-1, 4], [-8, 5]], + [[8, -4], [1, 10]]), [[8, 4], [3, 10]]) + + def testMinimum(self): + self.assertEqualArrays(minimum( + [[1, 2], [3, 4]], + [[-1, 4], [-8, 5]], + [[8, -4], [1, 10]]), [[-1, -4], [-8, 4]]) + + def testFloor(self): + self.assertEqualArrays(floor([[1.2, 5.8], [-1.2, 3.5]]), [[1., 5.], [-2., 3.]]) + + def testSqrt(self): + self.assertEqualArrays(sqrt([[9, 100], [25, 16]]), [[3, 10], [5, 4]]) + + def assertEqualArrays(self, a, b): + self.assertEqual(allclose(a, b), True) + +if __name__ == '__main__': + unittest.main() From 5f380ab00846aba451aad484a1a6ba0e69ffb4c6 Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Thu, 16 Jun 2016 07:39:05 +0200 Subject: [PATCH 21/30] moved numpy_utils into a new "numpy" package --- mathics/builtin/colors.py | 7 +- mathics/builtin/numpy/__init__.py | 18 ++ mathics/builtin/numpy/with_numpy.py | 144 ++++++++++ mathics/builtin/numpy/without_numpy.py | 182 +++++++++++++ mathics/builtin/numpy_utils.py | 274 -------------------- mathics/builtin/randomnumbers.py | 2 +- test/test_color.py | 13 +- test/{test_numpy_utils.py => test_numpy.py} | 8 +- 8 files changed, 359 insertions(+), 289 deletions(-) create mode 100644 mathics/builtin/numpy/__init__.py create mode 100644 mathics/builtin/numpy/with_numpy.py create mode 100644 mathics/builtin/numpy/without_numpy.py delete mode 100644 mathics/builtin/numpy_utils.py rename test/{test_numpy_utils.py => test_numpy.py} (95%) diff --git a/mathics/builtin/colors.py b/mathics/builtin/colors.py index 729eeadf4e..8823fdf0fa 100644 --- a/mathics/builtin/colors.py +++ b/mathics/builtin/colors.py @@ -1,10 +1,11 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -from math import pi -from mathics.builtin.numpy_utils import vectorized, stack, unstack, concat, array, clip, conditional, compose, choose -from mathics.builtin.numpy_utils import sqrt, floor, mod, cos, sin, arctan2, minimum, maximum, dot_t from itertools import chain +from math import pi + +from mathics.builtin.numpy import sqrt, floor, mod, cos, sin, arctan2, minimum, maximum, dot_t +from mathics.builtin.numpy import vectorized, stack, unstack, concat, array, clip, conditional, compose, choose # use rRGB D50 conversion like MMA. see http://www.brucelindbloom.com/Eqn_RGB_XYZ_Matrix.html # MMA seems to round matrix values to six significant digits. we do the same. diff --git a/mathics/builtin/numpy/__init__.py b/mathics/builtin/numpy/__init__.py new file mode 100644 index 0000000000..03cb4f0fe7 --- /dev/null +++ b/mathics/builtin/numpy/__init__.py @@ -0,0 +1,18 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +from __future__ import unicode_literals +from __future__ import absolute_import + +try: + import numpy + from mathics.builtin.numpy import with_numpy as numpy_layer +except ImportError: + from mathics.builtin.numpy import without_numpy as numpy_layer + +import types + +for s in dir(numpy_layer): + f = numpy_layer.__dict__.get(s) + if isinstance(f, types.FunctionType): + globals()[s] = getattr(numpy_layer, s) diff --git a/mathics/builtin/numpy/with_numpy.py b/mathics/builtin/numpy/with_numpy.py new file mode 100644 index 0000000000..a02c31aae0 --- /dev/null +++ b/mathics/builtin/numpy/with_numpy.py @@ -0,0 +1,144 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +""" +A couple of helper functions for doing numpy-like stuff with numpy. +""" + +from mathics.core.expression import Expression +from functools import reduce +import numpy + + +def is_numpy_available(): + return True + + +def instantiate_elements(a, new_element, d=1): + # given a numpy array 'a' and a python element constructor 'new_element', generate a python array of the + # same shape as 'a' with python elements constructed through 'new_element'. 'new_element' will get called + # if an array of dimension 'd' is reached. + + if len(a.shape) == d: + leaves = [new_element(x) for x in a] + else: + leaves = [instantiate_elements(e, new_element, d) for e in a] + return Expression('List', *leaves) + + +def array(a): + return numpy.array(a) + + +def vectorized(f, a, depth): + # we ignore depth, as numpy arrays will handle the vectorized cases. + a = array(a) + if len(a.shape) == 1: + # if a is a flat array, we embed the components in an array, i.e. + # [[a, b, c]], so that unstack will yield [a], [b], [c], and we + # operate on arrays of length 1 instead of scalars. + return f(array([a]))[0] + else: + return f(a) + + +def unstack(a): + a = array(a) + b = array(numpy.split(a, a.shape[-1], axis=-1)) + if b.shape[-1] == 1: + b = b.reshape(b.shape[:-1]) + return b + +def stack(*a): + # numpy.stack with axis=-1 stacks arrays along the most inner axis: + + # e.g. numpy.stack([ [1, 2], [3, 4] ], axis=-1) + # gives: array([ [1, 3], [2, 4] ]) + + # e.g. numpy.stack([ [[1, 2], [3, 4]], [[4, 5], [6, 7]] ], axis=-1) + # gives: array([[[1, 4], [2, 5]], [[3, 6], [4, 7]]]) + + a = array(a) + b = numpy.stack(a, axis=-1) + return b + + +def concat(*a): + a = [array(x) for x in a] + a = [x for x in a if x.shape[0]] # skip empty + return numpy.concatenate(a, axis=-1) + + +def conditional(a, cond, t, f): + b = array(a)[:] + mask = cond(a) + b[mask] = t(b[mask]) + b[~mask] = f(b[~mask]) + return b + + +def compose(*a): + assert a and len(a) % 2 == 0 + b = numpy.ndarray(a[0].shape) + # we apply the rules in reversed order, so that the first rules + # always make the last changes, which corresponds to the unreversed + # processing in the non-vectorized (non-numpy) implementation. + for mask, v in reversed([a[i:i + 2] for i in range(0, len(a), 2)]): + b[mask] = v(lambda x: x[mask]) + return b + + +def clip(a, t0, t1): + return numpy.clip(array(a), t0, t1) + + +def dot_t(u, v): + return numpy.dot(array(u), array(v).T) + + +def mod(a, b): + return numpy.mod(a, b) + + +def sin(a): + return numpy.sin(array(a)) + + +def cos(a): + return numpy.cos(array(a)) + + +def arctan2(y, x): + return numpy.arctan2(array(y), array(x)) + + +def sqrt(a): + return numpy.sqrt(array(a)) + + +def floor(a): + return numpy.floor(array(a)) + + +def maximum(*a): + return reduce(numpy.maximum, [array(x) for x in a]) + + +def minimum(*a): + return reduce(numpy.minimum, [array(x) for x in a]) + + +def choose(i, *options): + assert options + dim = len(options[0]) + columns = [[o[d] for o in options] for d in range(dim)] + if isinstance(i, (int, float)): + return [column[int(i)] for column in columns] # int cast needed for PyPy + else: + assert len(options) < 256 + i_int = array(i).astype(numpy.uint8) + return [numpy.choose(i_int, column) for column in columns] + + +def allclose(a, b): + return numpy.allclose(array(a), array(b)) diff --git a/mathics/builtin/numpy/without_numpy.py b/mathics/builtin/numpy/without_numpy.py new file mode 100644 index 0000000000..3f07c8a923 --- /dev/null +++ b/mathics/builtin/numpy/without_numpy.py @@ -0,0 +1,182 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +""" +A couple of helper functions for doing numpy-like stuff without numpy. +""" + +from mathics.core.expression import Expression +from itertools import chain +from math import sin as sinf, cos as cosf, sqrt as sqrtf, atan2 as atan2f, floor as floorf +import operator + +# If numpy is not available, we define the following fallbacks that are useful for implementing a similar +# logic in pure python without numpy. They obviously work on regular python array though, not numpy arrays. + + +def is_numpy_available(): + return False + + +def instantiate_elements(a, new_element, d=1): + # given a python array 'a' and a python element constructor 'new_element', generate a python array of the + # same shape as 'a' with python elements constructed through 'new_element'. 'new_element' will get called + # if an array of dimension 'd' is reached. + + e = a[0] + depth = 1 + while depth <= d and isinstance(e, list): + e = e[0] + depth += 1 + if d == depth: + leaves = [new_element(x) for x in a] + else: + leaves = [instantiate_elements(e, new_element, d) for e in a] + return Expression('List', *leaves) + + +def _is_bottom(a): + return any(not isinstance(x, list) for x in a) + + +def array(a): + return a + + +def vectorized(f, a, depth): + # depth == 0 means: f will get only scalars. depth == 1 means: f will + # get lists of scalars. + if _is_bottom(a): + if depth == 0: + return [f(x) for x in a] + else: + return f(a) + else: + return [vectorized(f, x, depth) for x in a] + + +def unstack(a): + if not a: + return [] + + def length(b): + return max(length(x) for x in b) if not _is_bottom(b) else len(b) + + def split(b, i): + if not _is_bottom(b): + return [split(x, i) for x in b] + else: + return b[i] + + return [split(a, i) for i in range(length(a))] + + +def stack(*a): + if _is_bottom(a): + return list(chain(a)) + else: + return [stack(*[x[i] for x in a]) for i in range(len(a[0]))] + + +def concat(a, b): + return stack(*(unstack(a) + unstack(b))) + + +def _apply(f, a): + if isinstance(a, (list, tuple)): + return [_apply(f, t) for t in a] + else: + return f(a) + + +def _apply_n(f, *p): + if isinstance(p[0], list): + return [_apply_n(f, *q) for q in zip(*p)] + else: + return f(*p) + + +def conditional(a, cond, t, f): + def _eval(x): + return t(x) if cond(x) else f(x) + + return _apply(_eval, a) + + +def compose(*a): + assert a and len(a) % 2 == 0 + for cond, v in (a[i:i + 2] for i in range(0, len(a), 2)): + if cond: + return v(lambda x: x) + raise ValueError('no matching case in compose') + + +def clip(a, t0, t1): + def _eval(x): + return max(t0, min(t1, x)) + + return _apply(_eval, a) + + +def dot_t(u, v): + if not isinstance(v[0], list): + return sum(x * y for x, y in zip(u, v)) + else: + return [sum(x * y for x, y in zip(u, r)) for r in v] + + +def mod(a, b): + return _apply_n(operator.mod, a, b) + + +def sin(a): + return _apply(sinf, a) + + +def cos(a): + return _apply(cosf, a) + + +def arctan2(y, x): + return _apply_n(atan2f, y, x) + + +def sqrt(a): + return _apply(sqrtf, a) + + +def floor(a): + return _apply(floorf, a) + + +def maximum(*a): + return _apply_n(max, *a) + + +def minimum(*a): + return _apply_n(min, *a) + + +def _choose_descend(i, options): + if isinstance(i, (int, float)): + return options[int(i)] # int cast needed for PyPy + else: + return [_choose_descend(next_i, [o[k] for o in options]) for k, next_i in enumerate(i)] + + +def choose(i, *options): + assert options + dim = len(options[0]) + columns = [[o[d] for o in options] for d in range(dim)] + return [_choose_descend(i, column) for column in columns] + + +def allclose(a, b): + if isinstance(a, list) and isinstance(b, list): + if len(a) != len(b): + return False + return all(allclose(x, y) for x, y in zip(a, b)) + elif isinstance(a, list) or isinstance(b, list): + return False + else: + return abs(a - b) < 1e-12 diff --git a/mathics/builtin/numpy_utils.py b/mathics/builtin/numpy_utils.py deleted file mode 100644 index 99e2277b7a..0000000000 --- a/mathics/builtin/numpy_utils.py +++ /dev/null @@ -1,274 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- - -""" -A couple of helper functions for working with numpy (and a couple of fallbacks, if numpy is unavailable) -""" - -from mathics.core.expression import Expression -from itertools import chain -from functools import reduce -from math import sin as sinf, cos as cosf, sqrt as sqrtf, atan2 as atan2f, floor as floorf -import operator - -try: - import numpy - _numpy = True -except ImportError: # no numpy? - _numpy = False - - -def py_instantiate_elements(a, new_element, d=1): - # given a python array 'a' and a python element constructor 'new_element', generate a python array of the - # same shape as 'a' with python elements constructed through 'new_element'. 'new_element' will get called - # if an array of dimension 'd' is reached. - - e = a[0] - depth = 1 - while depth <= d and isinstance(e, list): - e = e[0] - depth += 1 - if d == depth: - leaves = [new_element(x) for x in a] - else: - leaves = [py_instantiate_elements(e, new_element, d) for e in a] - return Expression('List', *leaves) - -if _numpy: - def instantiate_elements(a, new_element, d=1): - # given a numpy array 'a' and a python element constructor 'new_element', generate a python array of the - # same shape as 'a' with python elements constructed through 'new_element'. 'new_element' will get called - # if an array of dimension 'd' is reached. - - if len(a.shape) == d: - leaves = [new_element(x) for x in a] - else: - leaves = [instantiate_elements(e, new_element, d) for e in a] - return Expression('List', *leaves) - - def array(a): - return numpy.array(a) - - def vectorized(f, a, depth): - # we ignore depth, as numpy arrays will handle the vectorized cases. - a = array(a) - if len(a.shape) == 1: - # if a is a flat array, we embed the components in an array, i.e. - # [[a, b, c]], so that unstack will yield [a], [b], [c], and we - # operate on arrays of length 1 instead of scalars. - return f(array([a]))[0] - else: - return f(a) - - def unstack(a): - a = array(a) - b = array(numpy.split(a, a.shape[-1], axis=-1)) - if b.shape[-1] == 1: - b = b.reshape(b.shape[:-1]) - return b - - def stack(*a): - # numpy.stack with axis=-1 stacks arrays along the most inner axis: - - # e.g. numpy.stack([ [1, 2], [3, 4] ], axis=-1) - # gives: array([ [1, 3], [2, 4] ]) - - # e.g. numpy.stack([ [[1, 2], [3, 4]], [[4, 5], [6, 7]] ], axis=-1) - # gives: array([[[1, 4], [2, 5]], [[3, 6], [4, 7]]]) - - a = array(a) - b = numpy.stack(a, axis=-1) - return b - - def concat(*a): - a = [array(x) for x in a] - a = [x for x in a if x.shape[0]] # skip empty - return numpy.concatenate(a, axis=-1) - - def conditional(a, cond, t, f): - b = array(a)[:] - mask = cond(a) - b[mask] = t(b[mask]) - b[~mask] = f(b[~mask]) - return b - - def compose(*a): - assert a and len(a) % 2 == 0 - b = numpy.ndarray(a[0].shape) - # we apply the rules in reversed order, so that the first rules - # always make the last changes, which corresponds to the unreversed - # processing in the non-vectorized (non-numpy) implementation. - for mask, v in reversed([a[i:i + 2] for i in range(0, len(a), 2)]): - b[mask] = v(lambda x: x[mask]) - return b - - def clip(a, t0, t1): - return numpy.clip(array(a), t0, t1) - - def dot_t(u, v): - return numpy.dot(array(u), array(v).T) - - def mod(a, b): - return numpy.mod(a, b) - - def sin(a): - return numpy.sin(array(a)) - - def cos(a): - return numpy.cos(array(a)) - - def arctan2(y, x): - return numpy.arctan2(array(y), array(x)) - - def sqrt(a): - return numpy.sqrt(array(a)) - - def floor(a): - return numpy.floor(array(a)) - - def maximum(*a): - return reduce(numpy.maximum, [array(x) for x in a]) - - def minimum(*a): - return reduce(numpy.minimum, [array(x) for x in a]) - - def choose(i, *options): - assert options - dim = len(options[0]) - columns = [[o[d] for o in options] for d in range(dim)] - if isinstance(i, (int, float)): - return [column[int(i)] for column in columns] # int cast needed for PyPy - else: - assert len(options) < 256 - i_int = array(i).astype(numpy.uint8) - return [numpy.choose(i_int, column) for column in columns] - - def allclose(a, b): - return numpy.allclose(array(a), array(b)) -else: - # If numpy is not available, we define the following fallbacks that are useful for implementing a similar - # logic in pure python without numpy. They obviously work on regular python array though, not numpy arrays. - - instantiate_elements = py_instantiate_elements - - def _is_bottom(a): - return any(not isinstance(x, list) for x in a) - - def array(a): - return a - - def vectorized(f, a, depth): - # depth == 0 means: f will get only scalars. depth == 1 means: f will - # get lists of scalars. - if _is_bottom(a): - if depth == 0: - return [f(x) for x in a] - else: - return f(a) - else: - return [vectorized(f, x, depth) for x in a] - - def unstack(a): - if not a: - return [] - - def length(b): - return max(length(x) for x in b) if not _is_bottom(b) else len(b) - - def split(b, i): - if not _is_bottom(b): - return [split(x, i) for x in b] - else: - return b[i] - - return [split(a, i) for i in range(length(a))] - - def stack(*a): - if _is_bottom(a): - return list(chain(a)) - else: - return [stack(*[x[i] for x in a]) for i in range(len(a[0]))] - - def concat(a, b): - return stack(*(unstack(a) + unstack(b))) - - def _apply(f, a): - if isinstance(a, (list, tuple)): - return [_apply(f, t) for t in a] - else: - return f(a) - - def _apply_n(f, *p): - if isinstance(p[0], list): - return [_apply_n(f, *q) for q in zip(*p)] - else: - return f(*p) - - def conditional(a, cond, t, f): - def _eval(x): - return t(x) if cond(x) else f(x) - return _apply(_eval, a) - - def compose(*a): - assert a and len(a) % 2 == 0 - for cond, v in (a[i:i + 2] for i in range(0, len(a), 2)): - if cond: - return v(lambda x: x) - raise ValueError('no matching case in compose') - - def clip(a, t0, t1): - def _eval(x): - return max(t0, min(t1, x)) - return _apply(_eval, a) - - def dot_t(u, v): - if not isinstance(v[0], list): - return sum(x * y for x, y in zip(u, v)) - else: - return [sum(x * y for x, y in zip(u, r)) for r in v] - - def mod(a, b): - return _apply_n(operator.mod, a, b) - - def sin(a): - return _apply(sinf, a) - - def cos(a): - return _apply(cosf, a) - - def arctan2(y, x): - return _apply_n(atan2f, y, x) - - def sqrt(a): - return _apply(sqrtf, a) - - def floor(a): - return _apply(floorf, a) - - def maximum(*a): - return _apply_n(max, *a) - - def minimum(*a): - return _apply_n(min, *a) - - def _choose_descend(i, options): - if isinstance(i, (int, float)): - return options[int(i)] # int cast needed for PyPy - else: - return [_choose_descend(next_i, [o[k] for o in options]) for k, next_i in enumerate(i)] - - def choose(i, *options): - assert options - dim = len(options[0]) - columns = [[o[d] for o in options] for d in range(dim)] - return [_choose_descend(i, column) for column in columns] - - def allclose(a, b): - if isinstance(a, list) and isinstance(b, list): - if len(a) != len(b): - return False - return all(allclose(x, y) for x, y in zip(a, b)) - elif isinstance(a, list) or isinstance(b, list): - return False - else: - return abs(a - b) < 1e-12 diff --git a/mathics/builtin/randomnumbers.py b/mathics/builtin/randomnumbers.py index 07bb411320..ed1ee65185 100644 --- a/mathics/builtin/randomnumbers.py +++ b/mathics/builtin/randomnumbers.py @@ -20,7 +20,7 @@ from functools import reduce from mathics.builtin.base import Builtin -from mathics.builtin.numpy_utils import instantiate_elements, stack +from mathics.builtin.numpy import instantiate_elements, stack from mathics.core.expression import (Integer, String, Symbol, Real, Expression, Complex) diff --git a/test/test_color.py b/test/test_color.py index 15e9b9afc4..5dbaa15871 100644 --- a/test/test_color.py +++ b/test/test_color.py @@ -4,17 +4,16 @@ from __future__ import absolute_import from __future__ import unicode_literals -import os -import sys -import pexpect import unittest -from six.moves import range from random import random -from mathics.core.expression import Expression, Integer, Rational, Symbol + +from six.moves import range + +import mathics.builtin.colors as colors +from mathics.builtin.numpy.with_numpy import array, vectorized from mathics.core.definitions import Definitions from mathics.core.evaluation import Evaluation -import mathics.builtin.colors as colors -from mathics.builtin.numpy_utils import array, vectorized +from mathics.core.expression import Expression class ColorTest(unittest.TestCase): diff --git a/test/test_numpy_utils.py b/test/test_numpy.py similarity index 95% rename from test/test_numpy_utils.py rename to test/test_numpy.py index 2c51fa8d2e..73264f0760 100644 --- a/test/test_numpy_utils.py +++ b/test/test_numpy.py @@ -4,13 +4,13 @@ from __future__ import absolute_import from __future__ import unicode_literals -from mathics.builtin.numpy_utils import allclose, stack, unstack, concat, conditional, compose, clip, array, choose -from mathics.builtin.numpy_utils import vectorized, minimum, maximum, sin, dot_t, mod, floor, sqrt -import sys import unittest +from mathics.builtin.numpy.with_numpy import stack, unstack, concat, conditional, compose, clip, array, choose +from mathics.builtin.numpy.with_numpy import vectorized, minimum, maximum, dot_t, mod, floor, sqrt, allclose -class NumpyUtils(unittest.TestCase): + +class Numpy(unittest.TestCase): def testUnstack(self): # flat lists remain flat lists. self.assertEqualArrays( From 38b93ee808cd59b25d0d28c937ea219c83bb94e4 Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Thu, 16 Jun 2016 07:58:22 +0200 Subject: [PATCH 22/30] renamed package "numpy" back to "numpy_utils" --- mathics/builtin/colors.py | 4 ++-- mathics/builtin/{numpy => numpy_utils}/__init__.py | 4 ++-- mathics/builtin/{numpy => numpy_utils}/with_numpy.py | 0 mathics/builtin/{numpy => numpy_utils}/without_numpy.py | 0 mathics/builtin/randomnumbers.py | 2 +- test/test_color.py | 2 +- test/{test_numpy.py => test_numpy_utils.py} | 4 ++-- 7 files changed, 8 insertions(+), 8 deletions(-) rename mathics/builtin/{numpy => numpy_utils}/__init__.py (70%) rename mathics/builtin/{numpy => numpy_utils}/with_numpy.py (100%) rename mathics/builtin/{numpy => numpy_utils}/without_numpy.py (100%) rename test/{test_numpy.py => test_numpy_utils.py} (95%) diff --git a/mathics/builtin/colors.py b/mathics/builtin/colors.py index 8823fdf0fa..9104e1987e 100644 --- a/mathics/builtin/colors.py +++ b/mathics/builtin/colors.py @@ -4,8 +4,8 @@ from itertools import chain from math import pi -from mathics.builtin.numpy import sqrt, floor, mod, cos, sin, arctan2, minimum, maximum, dot_t -from mathics.builtin.numpy import vectorized, stack, unstack, concat, array, clip, conditional, compose, choose +from mathics.builtin.numpy_utils import sqrt, floor, mod, cos, sin, arctan2, minimum, maximum, dot_t +from mathics.builtin.numpy_utils import vectorized, stack, unstack, concat, array, clip, conditional, compose, choose # use rRGB D50 conversion like MMA. see http://www.brucelindbloom.com/Eqn_RGB_XYZ_Matrix.html # MMA seems to round matrix values to six significant digits. we do the same. diff --git a/mathics/builtin/numpy/__init__.py b/mathics/builtin/numpy_utils/__init__.py similarity index 70% rename from mathics/builtin/numpy/__init__.py rename to mathics/builtin/numpy_utils/__init__.py index 03cb4f0fe7..9377871133 100644 --- a/mathics/builtin/numpy/__init__.py +++ b/mathics/builtin/numpy_utils/__init__.py @@ -6,9 +6,9 @@ try: import numpy - from mathics.builtin.numpy import with_numpy as numpy_layer + from mathics.builtin.numpy_utils import with_numpy as numpy_layer except ImportError: - from mathics.builtin.numpy import without_numpy as numpy_layer + from mathics.builtin.numpy_utils import without_numpy as numpy_layer import types diff --git a/mathics/builtin/numpy/with_numpy.py b/mathics/builtin/numpy_utils/with_numpy.py similarity index 100% rename from mathics/builtin/numpy/with_numpy.py rename to mathics/builtin/numpy_utils/with_numpy.py diff --git a/mathics/builtin/numpy/without_numpy.py b/mathics/builtin/numpy_utils/without_numpy.py similarity index 100% rename from mathics/builtin/numpy/without_numpy.py rename to mathics/builtin/numpy_utils/without_numpy.py diff --git a/mathics/builtin/randomnumbers.py b/mathics/builtin/randomnumbers.py index ed1ee65185..07bb411320 100644 --- a/mathics/builtin/randomnumbers.py +++ b/mathics/builtin/randomnumbers.py @@ -20,7 +20,7 @@ from functools import reduce from mathics.builtin.base import Builtin -from mathics.builtin.numpy import instantiate_elements, stack +from mathics.builtin.numpy_utils import instantiate_elements, stack from mathics.core.expression import (Integer, String, Symbol, Real, Expression, Complex) diff --git a/test/test_color.py b/test/test_color.py index 5dbaa15871..21bc70df7b 100644 --- a/test/test_color.py +++ b/test/test_color.py @@ -10,7 +10,7 @@ from six.moves import range import mathics.builtin.colors as colors -from mathics.builtin.numpy.with_numpy import array, vectorized +from mathics.builtin.numpy_utils import array, vectorized from mathics.core.definitions import Definitions from mathics.core.evaluation import Evaluation from mathics.core.expression import Expression diff --git a/test/test_numpy.py b/test/test_numpy_utils.py similarity index 95% rename from test/test_numpy.py rename to test/test_numpy_utils.py index 73264f0760..f5053efcab 100644 --- a/test/test_numpy.py +++ b/test/test_numpy_utils.py @@ -6,8 +6,8 @@ import unittest -from mathics.builtin.numpy.with_numpy import stack, unstack, concat, conditional, compose, clip, array, choose -from mathics.builtin.numpy.with_numpy import vectorized, minimum, maximum, dot_t, mod, floor, sqrt, allclose +from mathics.builtin.numpy_utils import stack, unstack, concat, conditional, compose, clip, array, choose +from mathics.builtin.numpy_utils import vectorized, minimum, maximum, dot_t, mod, floor, sqrt, allclose class Numpy(unittest.TestCase): From dbf13ec05e7c7d6180573e1faaec438662386fa3 Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Thu, 16 Jun 2016 18:06:28 +0200 Subject: [PATCH 23/30] added "numpy_utils" to setup.py --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 91eb915adc..cf17010869 100644 --- a/setup.py +++ b/setup.py @@ -159,7 +159,7 @@ def run(self): 'mathics.algorithm', 'mathics.core', 'mathics.core.parser', - 'mathics.builtin', 'mathics.builtin.pymimesniffer', + 'mathics.builtin', 'mathics.builtin.pymimesniffer', 'mathics.builtin.numpy_utils', 'mathics.doc', 'mathics.web', 'mathics.web.templatetags' ], From eb916b6a52ca45110f52b71f1e963fdeec2a4c37 Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Thu, 16 Jun 2016 21:37:52 +0200 Subject: [PATCH 24/30] work in progress: bug fixes, test cases --- mathics/builtin/colors.py | 86 ++++---- mathics/builtin/numpy_utils/with_numpy.py | 13 +- mathics/builtin/numpy_utils/without_numpy.py | 10 + test/test_color.py | 203 ++++++++++++++++++- 4 files changed, 267 insertions(+), 45 deletions(-) diff --git a/mathics/builtin/colors.py b/mathics/builtin/colors.py index 9104e1987e..d77c1275a1 100644 --- a/mathics/builtin/colors.py +++ b/mathics/builtin/colors.py @@ -4,8 +4,9 @@ from itertools import chain from math import pi -from mathics.builtin.numpy_utils import sqrt, floor, mod, cos, sin, arctan2, minimum, maximum, dot_t +from mathics.builtin.numpy_utils import sqrt, floor, mod, cos, sin, arctan2, minimum, maximum, dot_t, errstate from mathics.builtin.numpy_utils import vectorized, stack, unstack, concat, array, clip, conditional, compose, choose +from mathics.builtin.numpy_utils import constant # use rRGB D50 conversion like MMA. see http://www.brucelindbloom.com/Eqn_RGB_XYZ_Matrix.html # MMA seems to round matrix values to six significant digits. we do the same. @@ -79,22 +80,26 @@ def rgb_to_hsb(pixels): m0 = minimum(r, g, b) c = m1 - m0 - h = compose( - c < 1e-15, - lambda s: 0, - m1 == r, - lambda s: mod(((s(g) - s(b)) / s(c)), 6.), - m1 == g, - lambda s: (s(b) - s(r)) / s(c) + 2., - m1 == b, - lambda s: (s(r) - s(g)) / s(c) + 4. - ) + with errstate(divide='ignore', invalid='ignore'): + eps = 1e-15 + + h = compose( + c < eps, + lambda mask: 0, + m1 == r, + lambda mask: mod(((mask(g) - mask(b)) / mask(c)), 6.), + m1 == g, + lambda mask: (mask(b) - mask(r)) / mask(c) + 2., + m1 == b, + lambda mask: (mask(r) - mask(g)) / mask(c) + 4.) + + s = compose(m1 < eps, lambda mask: 0, m1 >= eps, lambda mask: mask(c / m1)) h = conditional(h * (60. / 360.), lambda t: t < 0., lambda t: t + 1., lambda t: t) - return stack(h, c / m1, m1, *components[3:]) + return stack(h, s, m1, *components[3:]) def hsb_to_rgb(pixels): @@ -138,10 +143,15 @@ def rgb_to_cmyk(pixels): c = unstack(pixels) r, g, b = c[:3] - k = 1 - maximum(r, g, b) - k_ = 1 - k + k_ = maximum(r, g, b) + k = 1 - k_ - return stack((1 - r - k) / k_, (1 - g - k) / k_, (1 - b - k) / k_, k, *c[3:]) + eps = 1e-15 + with errstate(divide='ignore', invalid='ignore'): + c, m, y = [compose(k_ < eps, lambda s: s(constant(0, a)), + k_ >= eps, lambda s: (1 - s(a) - s(k)) / s(k_)) for a in (r, g, b)] + + return stack(c, m, y, k, *c[3:]) def xyz_to_rgb(pixels): @@ -190,7 +200,6 @@ def xyz_to_luv(pixels): v_0 = 9. * y_orig / q_0 lum /= 100.0 # MMA specific - u = 13. * lum * (u_0 - _ref_white.u_r) v = 13. * lum * (v_0 - _ref_white.v_r) @@ -199,19 +208,23 @@ def xyz_to_luv(pixels): def luv_to_xyz(pixels): components = unstack(pixels) - lum, u, v = components[:3] + cie_l, cie_u, cie_v = components[:3] + + # see http://www.easyrgb.com/index.php?X=MATH&H=17#text17 + u = cie_u / (13. * cie_l) + _ref_white.u_r + v = cie_v / (13. * cie_l) + _ref_white.v_r - u_0 = u / (13. * lum) + _ref_white.u_r - v_0 = v / (13. * lum) + _ref_white.v_r + cie_l_100 = cie_l * 100.0 # MMA specific + y = conditional((cie_l_100 + 16.) / 116., lambda t: t > 0.2068930, # 0.008856 ** (1/3) + lambda t: t ** 3, + lambda t: (t - 16. / 116.) / 7.787) - lum *= 100.0 # MMA specific - y = conditional(lum, lambda t: t <= 8., - lambda t: _ref_white.xyz[1] * t * ((3. / 29.) ** 3), - lambda t: _ref_white.xyz[1] * (((t + 16.) / 116.) ** 3)) - x = y * (9. * u_0) / (4. * v_0) - z = y * (12. - 3. * u_0 - 20. * v_0) / (4. * v_0) + x = -(9 * y * u) / ((u - 4) * v - u * v) + z = (9 * y - (15 * v * y) - (v * x)) / (3 * v) - return stack(x, y, z, *components[3:]) + x, y, z = [compose(cie_l < 0.078, lambda s: s(constant(0, a)), + cie_l >= 0.078, lambda s: s(a)) for a in (x, y, z)] + return clip(stack(x, y, z, *components[3:]), 0, 1) def lch_to_lab(pixels): @@ -232,24 +245,21 @@ def lab_to_lch(pixels): def lab_to_xyz(pixels): components = unstack(pixels) - - # see http://www.brucelindbloom.com/Eqn_Lab_to_XYZ.html l, a, b = components[:3] + + # see http://www.easyrgb.com/index.php?X=MATH&H=08#text8 f_y = (l * 100. + 16.) / 116. - x = conditional(a / 5. + f_y, lambda t: t > 0.008856, - lambda t: t ** 3, - lambda t: (116 * t - 16) / 903.3) - y = conditional(l * 100.0, lambda t: t > 903.3 * 0.008856, - lambda t: ((t + 16.) / 116.) ** 3, - lambda t: t / 903.3) - z = conditional(f_y - b / 2., lambda t: t > 0.008856, - lambda t: t ** 3, - lambda t: (116 * t - 16) / 903.3) + xyz = stack(a / 5. + f_y, f_y, f_y - b / 2.) + + xyz = conditional(xyz, lambda t: t > 0.2068930, # 0.008856 ** (1/3) + lambda t: t ** 3, + lambda t: (t - 16. / 116.) / 7.787) + x, y, z = unstack(xyz) x, y, z = [u * v for u, v in zip((x, y, z), _ref_white.xyz)] - return stack(x, y, z, *components[3:]) + return clip(stack(x, y, z, *components[3:]), 0, 1) _flows = { # see http://www.brucelindbloom.com/Math.html diff --git a/mathics/builtin/numpy_utils/with_numpy.py b/mathics/builtin/numpy_utils/with_numpy.py index a02c31aae0..ad68a3d1f6 100644 --- a/mathics/builtin/numpy_utils/with_numpy.py +++ b/mathics/builtin/numpy_utils/with_numpy.py @@ -49,6 +49,7 @@ def unstack(a): b = b.reshape(b.shape[:-1]) return b + def stack(*a): # numpy.stack with axis=-1 stacks arrays along the most inner axis: @@ -58,7 +59,7 @@ def stack(*a): # e.g. numpy.stack([ [[1, 2], [3, 4]], [[4, 5], [6, 7]] ], axis=-1) # gives: array([[[1, 4], [2, 5]], [[3, 6], [4, 7]]]) - a = array(a) + a = [array(x) for x in a] b = numpy.stack(a, axis=-1) return b @@ -142,3 +143,13 @@ def choose(i, *options): def allclose(a, b): return numpy.allclose(array(a), array(b)) + + +def errstate(**kwargs): + return numpy.errstate(**kwargs) + + +def constant(x, a): + b = numpy.ndarray(a.shape) + b.fill(x) + return b diff --git a/mathics/builtin/numpy_utils/without_numpy.py b/mathics/builtin/numpy_utils/without_numpy.py index 3f07c8a923..83da0a0b27 100644 --- a/mathics/builtin/numpy_utils/without_numpy.py +++ b/mathics/builtin/numpy_utils/without_numpy.py @@ -7,6 +7,7 @@ from mathics.core.expression import Expression from itertools import chain +from contextlib import contextmanager from math import sin as sinf, cos as cosf, sqrt as sqrtf, atan2 as atan2f, floor as floorf import operator @@ -180,3 +181,12 @@ def allclose(a, b): return False else: return abs(a - b) < 1e-12 + + +@contextmanager +def errstate(**kwargs): + yield + + +def constant(x, a): + return x diff --git a/test/test_color.py b/test/test_color.py index 21bc70df7b..93c57d7398 100644 --- a/test/test_color.py +++ b/test/test_color.py @@ -16,6 +16,176 @@ from mathics.core.expression import Expression +_color_tests = [[ + [ + [0.1], + [0.1, 0.1, 0.1], + [0., 0., 0., 0.9], + [0., 0., 0.1], + [0.009664, 0.010023, 0.008271], + [0.090104, 0., 0.], + [0.090104, 0., 0.129196], + [0.090104, 0., 0.] + ], + [ + [0.1815], + [0.1, 0.2, 0.3], + [0.666667, 0.333333, 0., 0.7], + [0.583333, 0.666667, 0.3], + [0.027597, 0.030402, 0.05566], + [0.202041, -0.031078, -0.189912], + [0.202041, 0.192438, 0.724184], + [0.202041, -0.103717, -0.177331] + ], + [ + [0.8185], + [0.9, 0.8, 0.7], + [0.1, 0.2, 0.3, 0.], + [0.083333, 0.222222, 0.9], + [0.639982, 0.635229, 0.389546], + [0.837168, 0.063343, 0.161993], + [0.837168, 0.173937, 0.190676], + [0.837168, 0.181115, 0.176395] + ], + [ + [0.279072], + [0.3, 0.276, 0.24], + [0., 0.08, 0.2, 0.7], + [0.1, 0.2, 0.3], + [0.062498, 0.063527, 0.040573], + [0.302854, 0.013415, 0.065333], + [0.302854, 0.066696, 0.217769], + [0.302854, 0.042072, 0.057975] + ], + [ + [0.414652], + [0., 0.57972, 0.652246], + [1., 0.111193, 0., 0.347754], + [0.518532, 1., 0.652246], + [0.1, 0.2, 0.3], + [0.518372, -0.574865, -0.257804], + [0.518372, 0.630026, 0.567095], + [0.518372, -0.735612, -0.256568] + ], + [ + [0.094788], + [0.213536, 0.05271, 0.], + [0., 0.753156, 1., 0.786464], + [0.041141, 1., 0.213536], + [0.017769, 0.01126, 0.], + [0.1, 0.2, 0.3], + [0.1, 0.360555, 0.156416], + [0.1, 0.223077, 0.071257] + ], + [ + [0.098639], + [0.102999, 0.115575, 0.], + [0.10882, 0., 1., 0.884425], + [0.184803, 1., 0.115575], + [0.009158, 0.01126, 0.], + [0.1, -0.061803, 0.190211], + [0.1, 0.2, 0.3], + [0.1, -0.004455, 0.105387] + ], + [ + [0.102653], + [0.159623, 0.09357, 0.], + [0., 0.41381, 1., 0.840377], + [0.097698, 1., 0.159623], + [0.012794, 0.01126, 0.], + [0.1, 0.063027, 0.172414], + [0.1, 0.183573, 0.194222], + [0.1, 0.2, 0.3] + ] +], +[ + [ + [0.001], + [0.001, 0.001, 0.001], + [0., 0., 0., 0.999], + [0., 0., 0.001], + [0.000075, 0.000077, 0.000064], + [0.000699, 0., 0.], + [0.000699, 0., 0.129196], + [0.000699, 0., 0.] + ], + [ + [0.001815], + [0.001, 0.002, 0.003], + [0.666667, 0.333333, 0., 0.997], + [0.583333, 0.666667, 0.003], + [0.000127, 0.000142, 0.000182], + [0.001285, -0.000428, -0.001218], + [0.001285, 0.001291, 0.696235], + [0.001285, -0.00048, -0.000532] + ], + [ + [0.998185], + [0.999, 0.998, 0.997], + [0.001, 0.002, 0.003, 0.], + [0.083333, 0.002002, 0.999], + [0.960504, 0.995824, 0.819873], + [0.998383, 0.000539, 0.001533], + [0.998383, 0.001625, 0.196175], + [0.998383, 0.001674, 0.001856] + ], + [ + [0.002996], + [0.003, 0.002994, 0.002994], + [0., 0.001988, 0.002, 0.997], + [0.001, 0.002, 0.003], + [0.000224, 0.000232, 0.000191], + [0.002094, 4.e-6, 2.e-6], + [0.002094, 4.e-6, 0.056043], + [0.002094, 3.e-6, 0.] + ], + [ + [0.027967], + [0., 0.038164, 0.048809], + [1., 0.218102, 0., 0.951191], + [0.53635, 1., 0.048809], + [0.001, 0.002, 0.003], + [0.018066, -0.03749, -0.02547], + [0.018066, 0.045324, 0.594977], + [0.018066, -0.025637, -0.008942] + ], + [ + [0.001649], + [0.004015, 0.000765, 0.], + [0., 0.80957, 1., 0.995985], + [0.031738, 1., 0.004015], + [0.000156, 0.000111, 0.], + [0.001, 0.002, 0.003], + [0.001, 0.003606, 0.156416], + [0.001, 0.001754, 0.000784] + ], + [ + [0.001667], + [0.003448, 0.000803, 0.00144], + [0., 0.767033, 0.582337, 0.996552], + [0.959868, 0.767033, 0.003448], + [0.000156, 0.000111, 0.000089], + [0.001, 0.002, 0.000038], + [0.001, 0.002, 0.003], + [0.001, 0.001178, -0.000132] + ], + [ + [0.], + [0., 0., 0.], + [0., 0., 0., 1.], + [0., 0., 0.], + [0., 0., 0.], + [0., 0., 0.], + [0., 0., 0.], + [0.001, 0.002, 0.003] + ] +]] + + +def _color_string(head, components): + return '%s[%s]' % (head, ', '.join(['%.6f' % c for c in components])) + + class ColorTest(unittest.TestCase): def setUp(self): definitions = Definitions(add_builtin=True) @@ -39,10 +209,10 @@ def testInverseConversions(self): # now calculate from_space -> to_space -> from_space inverted = [c.to_python() for c in Expression('ColorConvert', - Expression('ColorConvert', - Expression(construct_name, *original), - to_space), - from_space).evaluate(self.evaluation).leaves] + Expression('ColorConvert', + Expression(construct_name, *original), + to_space), + from_space).evaluate(self.evaluation).leaves] if from_space == 'CMYK': # if cmyk, cmyk -> cmy k = inverted[3] inverted = [c * (1 - k) + k for c in inverted[:3]] @@ -50,8 +220,8 @@ def testInverseConversions(self): for x, y in zip(original, inverted): self.assertAlmostEqual(x, y, places) except: - print('test failed for %s(%f, %f, %f) -> %s -> %s' % - (from_space, original[0], original[1], original[2], to_space, from_space)) + print('test failed for %s -> %s -> %s' % + (_color_string(from_space, original), to_space, from_space)) raise def testConversions(self): @@ -68,6 +238,27 @@ def testConversions(self): self._checkConversion("XYZ", (0.1, 0.1, 0.1), "LAB", (0.3784243088316416, 0.02835852516741566, -0.06139347105996884)) + values = [[0.1, 0.2, 0.3, 0.0], [0.001, 0.002, 0.003, 0.0]] + spaces = ["Grayscale", "RGB", "CMYK", "HSB", "XYZ", "LAB", "LCH", "LUV"] + + for components, t1 in zip(values, _color_tests): + for src, t2 in zip(spaces, t1): + for dst, expected in zip(spaces, t2): + if src == "Grayscale": + c = components[:1] + elif src == 'CMYK': + c = components + else: + c = components[:3] + result = colors.convert(c, src, dst) + res_name = _color_string(src, c) + dst_name = _color_string(dst, result) + exp_name = _color_string(dst, expected) + self.assertIsNotNone(result, '%s -> %s != %s' % (res_name, dst_name, exp_name)) + self.assertEqual(len(result), len(expected), '%s -> %s != %s' % (res_name, dst_name, exp_name)) + for a, b in zip(result, expected): + self.assertAlmostEqual(a, b, 2, '%s -> %s != %s' % (res_name, dst_name, exp_name)) + def testImageConversions(self): # test that f([x, y, ...]) = [f(x), f(y), ...] for rectangular image arrays. From cd045717f21e75e62d9c6fabf92208486b3dbb1d Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Thu, 16 Jun 2016 21:45:21 +0200 Subject: [PATCH 25/30] bug fixes --- mathics/builtin/colors.py | 6 +++--- test/test_color.py | 7 ++++--- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/mathics/builtin/colors.py b/mathics/builtin/colors.py index d77c1275a1..a4cec9c1b9 100644 --- a/mathics/builtin/colors.py +++ b/mathics/builtin/colors.py @@ -140,9 +140,9 @@ def cmyk_to_rgb(pixels): def rgb_to_cmyk(pixels): - c = unstack(pixels) + components = unstack(pixels) - r, g, b = c[:3] + r, g, b = components[:3] k_ = maximum(r, g, b) k = 1 - k_ @@ -151,7 +151,7 @@ def rgb_to_cmyk(pixels): c, m, y = [compose(k_ < eps, lambda s: s(constant(0, a)), k_ >= eps, lambda s: (1 - s(a) - s(k)) / s(k_)) for a in (r, g, b)] - return stack(c, m, y, k, *c[3:]) + return stack(c, m, y, k, *components[3:]) def xyz_to_rgb(pixels): diff --git a/test/test_color.py b/test/test_color.py index 93c57d7398..3daec4ad18 100644 --- a/test/test_color.py +++ b/test/test_color.py @@ -257,13 +257,14 @@ def testConversions(self): self.assertIsNotNone(result, '%s -> %s != %s' % (res_name, dst_name, exp_name)) self.assertEqual(len(result), len(expected), '%s -> %s != %s' % (res_name, dst_name, exp_name)) for a, b in zip(result, expected): - self.assertAlmostEqual(a, b, 2, '%s -> %s != %s' % (res_name, dst_name, exp_name)) + self.assertAlmostEqual(a, b, 4, '%s -> %s != %s' % (res_name, dst_name, exp_name)) def testImageConversions(self): # test that f([x, y, ...]) = [f(x), f(y), ...] for rectangular image arrays. - for convert in colors.conversions.values(): - self._checkImageConversion(4, convert) + for name, convert in colors.conversions.items(): + if name.find('CMYK') < 0: + self._checkImageConversion(4, convert) def _checkConversion(self, from_space, from_components, to_space, to_components): places = 12 From b24e62805f1f8c42351dfa64faaeda8707955e06 Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Thu, 23 Jun 2016 18:43:29 +0200 Subject: [PATCH 26/30] major cleanup --- mathics/builtin/colors.py | 310 +++++++++------- mathics/builtin/numpy_utils/__init__.py | 35 +- mathics/builtin/numpy_utils/with_numpy.py | 365 ++++++++++++++++--- mathics/builtin/numpy_utils/without_numpy.py | 158 ++++---- test/test_color.py | 8 +- test/test_numpy_utils.py | 54 +-- 6 files changed, 646 insertions(+), 284 deletions(-) mode change 100644 => 100755 mathics/builtin/numpy_utils/__init__.py mode change 100644 => 100755 mathics/builtin/numpy_utils/with_numpy.py mode change 100644 => 100755 mathics/builtin/numpy_utils/without_numpy.py mode change 100644 => 100755 test/test_color.py diff --git a/mathics/builtin/colors.py b/mathics/builtin/colors.py index a4cec9c1b9..f6a91a0401 100644 --- a/mathics/builtin/colors.py +++ b/mathics/builtin/colors.py @@ -4,14 +4,43 @@ from itertools import chain from math import pi -from mathics.builtin.numpy_utils import sqrt, floor, mod, cos, sin, arctan2, minimum, maximum, dot_t, errstate -from mathics.builtin.numpy_utils import vectorized, stack, unstack, concat, array, clip, conditional, compose, choose -from mathics.builtin.numpy_utils import constant +from mathics.builtin.numpy_utils import sqrt, floor, mod, cos, sin, arctan2, minimum, maximum, dot_t +from mathics.builtin.numpy_utils import stack, unstack, array, clip, conditional, choose, stacked + +# in the long run, we might want to implement these functions using Compile[]. until Compile[] is available for all +# Mathics configurations, this implementation works on Python and PyPy and with or without numpy. + + +def _clip1(t): + return clip(t, 0, 1) + + +class _PerfectReflectingDiffuser: + def __init__(self, *x): # accepts (x1, y2) or (x, y, z) + assert 2 <= len(x) <= 3 + + if len(x) == 3: + self.xyz = x + else: + x1, x2 = x + scale = 1.0 / x2 + self.xyz = (x1 * scale, 1.0, (1.0 - x1 - x2) * scale) + + q_r = self.xyz[0] + 15. * self.xyz[1] + 3. * self.xyz[2] + self.u_r = 4. * self.xyz[0] / q_r + self.v_r = 9. * self.xyz[1] / q_r + + +# MMA's reference white is a D50, 2 degrees diffuser; the standard values for this configuration +# can be found at https://en.wikipedia.org/wiki/Standard_illuminant, MMA seems to use slightly +# different values though. + +_ref_white = _PerfectReflectingDiffuser(0.96422, 1., 0.82521) # use rRGB D50 conversion like MMA. see http://www.brucelindbloom.com/Eqn_RGB_XYZ_Matrix.html # MMA seems to round matrix values to six significant digits. we do the same. -xyz_from_rgb = [ +_xyz_from_rgb = [ [0.436075, 0.385065, 0.14308], [0.222504, 0.716879, 0.0606169], [0.0139322, 0.0971045, 0.714173] @@ -19,93 +48,83 @@ # for matrix, see http://www.brucelindbloom.com/Eqn_RGB_XYZ_Matrix.html # MMA seems to round matrix values to six significant digits. we do the same. -rgb_from_xyz = [ + +_rgb_from_xyz = [ [3.13386, -1.61687, -0.490615], [-0.978768, 1.91614, 0.033454], [0.0719453, -0.228991, 1.40524] ] -class _PerfectReflectingDiffuser: - def __init__(self, x1, x2): - # MMA seems to use the following constants, and not the - # values derived via calculation (commented out below) - self.xyz = [0.96422, 1., 0.82521] - # scale = 1.0 / x2 - # self.xyz = (x1 * scale, 1.0, (1.0 - x1 - x2) * scale) +def rgb_to_grayscale(r, g, b, *rest): + # see https://en.wikipedia.org/wiki/Grayscale + y = 0.299 * r + 0.587 * g + 0.114 * b # Y of Y'UV + return (y,) + rest - q_r = self.xyz[0] + 15. * self.xyz[1] + 3. * self.xyz[2] - self.u_r = 4. * self.xyz[0] / q_r - self.v_r = 9. * self.xyz[1] / q_r +def grayscale_to_rgb(g, *rest): + return (g, g, g) + rest -# MMA's reference white is a # D50, 2 degrees diffuser; for the -# values, see https://en.wikipedia.org/wiki/Standard_illuminant -_ref_white = _PerfectReflectingDiffuser(0.34567, 0.35850) +@conditional +def _inverse_compand_srgb(x): + if x > 0.04045: + return ((x + 0.055) / 1.055) ** 2.4 + else: + return x / 12.92 -def rgb_to_grayscale(pixels): - # see https://en.wikipedia.org/wiki/Grayscale - components = unstack(pixels) - r, g, b = components[:3] - y = 0.299 * r + 0.587 * g + 0.114 * b # Y of Y'UV - return stack(y, *components[3:]) +def rgb_to_xyz(r, g, b, *rest): + r, g, b = map(_inverse_compand_srgb, (r, g, b)) + x, y, z = unstack(_clip1(dot_t(stack(r, g, b), _xyz_from_rgb))) + return map(_clip1, (x, y, z) + rest) -def grayscale_to_rgb(pixels): - components = unstack(pixels) - g = components[0] - return stack(g, g, g, *components[1:]) +@conditional +def _compute_hsb_s(m1, c, eps): + if m1 < eps: + return 0 + else: + return c / m1 -def rgb_to_xyz(pixels): - components = unstack(pixels) - rgb = components[:3] +@conditional +def _compute_hsb_h(m1, c, eps, r, g, b): + if c < eps: + return 0 + elif m1 == r: + return mod(((g - b) / c), 6.) + elif m1 == g: + return (b - r) / c + 2. + elif m1 == b: + return (r - g) / c + 4. - x, y, z = conditional(rgb, lambda t: t > 0.04045, - lambda t: ((t + 0.055) / 1.055) ** 2.4, - lambda t: t / 12.92) - xyz = dot_t(stack(x, y, z), xyz_from_rgb) - return clip(concat(xyz, components[3:]), 0, 1) +@conditional +def _wrap_hsb_h(h): + if h < 0.: + return h + 1. + else: + return h -def rgb_to_hsb(pixels): +def rgb_to_hsb(r, g, b, *rest): # see https://en.wikipedia.org/wiki/HSB_color_space. HSB is also known as HSV. - components = unstack(pixels) - r, g, b = clip(components[:3], 0, 1) + r, g, b = _clip1((r, g, b)) m1 = maximum(r, g, b) m0 = minimum(r, g, b) c = m1 - m0 + eps = 1e-15 - with errstate(divide='ignore', invalid='ignore'): - eps = 1e-15 - - h = compose( - c < eps, - lambda mask: 0, - m1 == r, - lambda mask: mod(((mask(g) - mask(b)) / mask(c)), 6.), - m1 == g, - lambda mask: (mask(b) - mask(r)) / mask(c) + 2., - m1 == b, - lambda mask: (mask(r) - mask(g)) / mask(c) + 4.) - - s = compose(m1 < eps, lambda mask: 0, m1 >= eps, lambda mask: mask(c / m1)) - - h = conditional(h * (60. / 360.), lambda t: t < 0., - lambda t: t + 1., - lambda t: t) - - return stack(h, s, m1, *components[3:]) + h = _compute_hsb_h(m1, c, eps, r, g, b) + h = _wrap_hsb_h(h * (60. / 360.)) + s = _compute_hsb_s(m1, c, eps) + return (h, s, m1) + rest -def hsb_to_rgb(pixels): - components = unstack(pixels) - h, s, v = components[:3] +def hsb_to_rgb(h, s, v, *rest): i = floor(6 * h) f = 6 * h - i i = mod(i, 6) @@ -121,79 +140,88 @@ def hsb_to_rgb(pixels): (t, p, v), (v, p, q)) - return stack(r, g, b, *components[3:]) + return (r, g, b) + rest -def cmyk_to_rgb(pixels): - c = unstack(pixels) - - if len(c) >= 4: - k = c[3] +def cmyk_to_rgb(c, m, y, *rest): + if rest: + k = rest[0] + rest = rest[1:] else: k = 0 k_ = 1 - k - cmy = [(x * k_ + k) for x in c[:3]] + cmy = [(x * k_ + k) for x in (c, m, y)] r, g, b = [1 - x for x in cmy] - return concat(stack(r, g, b), c[4:]) + return (r, g, b) + rest + +@conditional +def _scale_rgb_to_cmyk(t, k, k_, eps): + if k_ < eps: + return 0 + else: + return (1 - t - k) / k_ -def rgb_to_cmyk(pixels): - components = unstack(pixels) - r, g, b = components[:3] +def rgb_to_cmyk(r, g, b, *rest): k_ = maximum(r, g, b) k = 1 - k_ + c, m, y = map(lambda t: _scale_rgb_to_cmyk(t, k, k_, 1e-15), (r, g, b)) + return (c, m, y, k) + rest - eps = 1e-15 - with errstate(divide='ignore', invalid='ignore'): - c, m, y = [compose(k_ < eps, lambda s: s(constant(0, a)), - k_ >= eps, lambda s: (1 - s(a) - s(k)) / s(k_)) for a in (r, g, b)] - return stack(c, m, y, k, *components[3:]) +@conditional +def _compand_srgb(t): + if t > 0.0031308: + return 1.055 * (t ** (1. / 2.4)) - 0.055 + else: + return t * 12.92 -def xyz_to_rgb(pixels): - components = unstack(pixels) - x, y, z = clip(components[:3], 0, 1) - xyz = dot_t(stack(x, y, z), rgb_from_xyz) +def xyz_to_rgb(x, y, z, *rest): + x, y, z = map(_clip1, (x, y, z)) + r, g, b = unstack(dot_t(stack(x, y, z), _rgb_from_xyz)) + r, g, b = map(_clip1, (r, g, b)) + r, g, b = map(_compand_srgb, (r, g, b)) + return map(_clip1, (r, g, b) + rest) - rgb = conditional(xyz, lambda t: t > 0.0031308, - lambda t: 1.055 * (t ** (1. / 2.4)) - 0.055, - lambda t: t * 12.92) - return concat(clip(rgb, 0, 1), components[3:]) +@conditional +def _scale_xyz_to_lab(t): + if t > 0.008856: + return t ** 0.33333333 # MMA specific + else: + return (903.3 * t + 16.) / 116. -def xyz_to_lab(pixels): +def xyz_to_lab(x, y, z, *rest): # see http://www.brucelindbloom.com/Eqn_XYZ_to_Lab.html - components = unstack(pixels) - xyz = components[:3] - - xyz = array([u / v for u, v in zip(xyz, _ref_white.xyz)]) - xyz = conditional(xyz, lambda c: c > 0.008856, - lambda c: c ** 0.33333333, - lambda c: (903.3 * c + 16.) / 116.) - x, y, z = xyz + xyz = array([u / v for u, v in zip([x, y, z], _ref_white.xyz)]) + x, y, z = map(_scale_xyz_to_lab, xyz) # MMA scales by 1/100 - return stack((1.16 * y) - 0.16, 5. * (x - y), 2. * (y - z), *components[3:]) + return ((1.16 * y) - 0.16, 5. * (x - y), 2. * (y - z)) + rest + + +@conditional +def _xyz_to_luv_scale_y(y): + if y > 0.008856: + return 116. * (y ** (1. / 3.)) - 16 + else: + return 903.3 * y -def xyz_to_luv(pixels): +def xyz_to_luv(x, y, z, *rest): # see http://www.brucelindbloom.com/Eqn_XYZ_to_Luv.html # and https://en.wikipedia.org/wiki/CIELUV - - components = unstack(pixels) - xyz = clip(components[:3], 0, 1) + xyz = _clip1((x, y, z)) x_orig, y_orig, z_orig = xyz y = y_orig / _ref_white.xyz[1] - lum = conditional(y, lambda t: t > 0.008856, - lambda t: 116. * (t ** (1. / 3.)) - 16, - lambda t: 903.3 * t) + lum = _xyz_to_luv_scale_y(y) q_0 = x_orig + 15. * y_orig + 3. * z_orig u_0 = 4. * x_orig / q_0 @@ -203,63 +231,73 @@ def xyz_to_luv(pixels): u = 13. * lum * (u_0 - _ref_white.u_r) v = 13. * lum * (v_0 - _ref_white.v_r) - return stack(lum, u, v, *components[3:]) + return (lum, u, v) + rest + +@conditional +def _scale_lab_to_xyz(t): + if t > 0.2068930: # 0.008856 ** (1/3) + return t ** 3 + else: + return (t - 16. / 116.) / 7.787 -def luv_to_xyz(pixels): - components = unstack(pixels) - cie_l, cie_u, cie_v = components[:3] +@conditional +def _luv_to_xyz_clip_zero(cie_l_is_zero, t): + if cie_l_is_zero: + return 0 + else: + return t + + +def luv_to_xyz(cie_l, cie_u, cie_v, *rest): # see http://www.easyrgb.com/index.php?X=MATH&H=17#text17 u = cie_u / (13. * cie_l) + _ref_white.u_r v = cie_v / (13. * cie_l) + _ref_white.v_r cie_l_100 = cie_l * 100.0 # MMA specific - y = conditional((cie_l_100 + 16.) / 116., lambda t: t > 0.2068930, # 0.008856 ** (1/3) - lambda t: t ** 3, - lambda t: (t - 16. / 116.) / 7.787) + y = (cie_l_100 + 16.) / 116. + y = _scale_lab_to_xyz(y) x = -(9 * y * u) / ((u - 4) * v - u * v) z = (9 * y - (15 * v * y) - (v * x)) / (3 * v) - x, y, z = [compose(cie_l < 0.078, lambda s: s(constant(0, a)), - cie_l >= 0.078, lambda s: s(a)) for a in (x, y, z)] - return clip(stack(x, y, z, *components[3:]), 0, 1) + cie_l_is_zero = cie_l < 0.078 + x, y, z = map(lambda t: _luv_to_xyz_clip_zero(cie_l_is_zero, t), (x, y, z)) + x, y, z = clip([x, y, z], 0, 1) + + return (x, y, z) + rest -def lch_to_lab(pixels): - components = unstack(pixels) - l, c, h = components[:3] +def lch_to_lab(l, c, h, *rest): h *= 2 * pi # MMA specific - return stack(l, c * cos(h), c * sin(h), *components[3:]) + return (l, c * cos(h), c * sin(h)) + rest -def lab_to_lch(pixels): - components = unstack(pixels) - l, a, b = components[:3] - h = conditional(arctan2(b, a), lambda t: t < 0., - lambda t: t + 2. * pi, lambda t: t) - h /= 2. * pi # MMA specific - return stack(l, sqrt(a * a + b * b), h, *components[3:]) +@conditional +def _wrap_lch_h(h, pi2): + if h < 0.: + return h + pi2 + else: + return h + +def lab_to_lch(l, a, b, *rest): + h = _wrap_lch_h(arctan2(b, a), 2. * pi) + h /= 2. * pi # MMA specific + return (l, sqrt(a * a + b * b), h) + rest -def lab_to_xyz(pixels): - components = unstack(pixels) - l, a, b = components[:3] +def lab_to_xyz(l, a, b, *rest): # see http://www.easyrgb.com/index.php?X=MATH&H=08#text8 f_y = (l * 100. + 16.) / 116. + x, y, z = a / 5. + f_y, f_y, f_y - b / 2. + x, y, z = map(_scale_lab_to_xyz, (x, y, z)) - xyz = stack(a / 5. + f_y, f_y, f_y - b / 2.) - - xyz = conditional(xyz, lambda t: t > 0.2068930, # 0.008856 ** (1/3) - lambda t: t ** 3, - lambda t: (t - 16. / 116.) / 7.787) - - x, y, z = unstack(xyz) x, y, z = [u * v for u, v in zip((x, y, z), _ref_white.xyz)] + x, y, z = clip([x, y, z], 0, 1) - return clip(stack(x, y, z, *components[3:]), 0, 1) + return (x, y, z) + rest _flows = { # see http://www.brucelindbloom.com/Math.html @@ -318,6 +356,6 @@ def convert(components, src, dst): func = conversions.get('%s>%s' % (s, d)) if not func: return None - components = vectorized(func, components, 1) + components = stacked(func, components) return components diff --git a/mathics/builtin/numpy_utils/__init__.py b/mathics/builtin/numpy_utils/__init__.py old mode 100644 new mode 100755 index 9377871133..ad5d4d876c --- a/mathics/builtin/numpy_utils/__init__.py +++ b/mathics/builtin/numpy_utils/__init__.py @@ -10,9 +10,34 @@ except ImportError: from mathics.builtin.numpy_utils import without_numpy as numpy_layer -import types +# we explicitly list all imported symbols so IDEs as PyCharm can properly +# do their code intelligence. + +array = numpy_layer.array + +stack = numpy_layer.stack +unstack = numpy_layer.unstack +concat = numpy_layer.concat +stacked = numpy_layer.stacked +vectorize = numpy_layer.vectorize + +conditional = numpy_layer.conditional +choose = numpy_layer.choose + +clip = numpy_layer.clip +sqrt = numpy_layer.sqrt +floor = numpy_layer.floor +mod = numpy_layer.mod +cos = numpy_layer.cos +sin = numpy_layer.sin +arctan2 = numpy_layer.arctan2 + +minimum = numpy_layer.minimum +maximum = numpy_layer.maximum +dot_t = numpy_layer.dot_t + +is_numpy_available = numpy_layer.is_numpy_available +allclose = numpy_layer.allclose +errstate = numpy_layer.errstate +instantiate_elements = numpy_layer.instantiate_elements -for s in dir(numpy_layer): - f = numpy_layer.__dict__.get(s) - if isinstance(f, types.FunctionType): - globals()[s] = getattr(numpy_layer, s) diff --git a/mathics/builtin/numpy_utils/with_numpy.py b/mathics/builtin/numpy_utils/with_numpy.py old mode 100644 new mode 100755 index ad68a3d1f6..8b9cca7903 --- a/mathics/builtin/numpy_utils/with_numpy.py +++ b/mathics/builtin/numpy_utils/with_numpy.py @@ -8,40 +8,36 @@ from mathics.core.expression import Expression from functools import reduce import numpy +import ast +import inspect +import sys -def is_numpy_available(): - return True +# +# INTERNAL FUNCTIONS +# +def _promote(x, shape): + if isinstance(x, (int, float)): + data = numpy.ndarray(shape) + data.fill(x) + return data + else: + return x -def instantiate_elements(a, new_element, d=1): - # given a numpy array 'a' and a python element constructor 'new_element', generate a python array of the - # same shape as 'a' with python elements constructed through 'new_element'. 'new_element' will get called - # if an array of dimension 'd' is reached. - if len(a.shape) == d: - leaves = [new_element(x) for x in a] - else: - leaves = [instantiate_elements(e, new_element, d) for e in a] - return Expression('List', *leaves) +def _is_scalar(x): + return not isinstance(x, (list, tuple, numpy.ndarray)) + +# +# ARRAY CREATION AND REORGANIZATION: STACK, UNSTACK, CONCAT, ... +# def array(a): return numpy.array(a) -def vectorized(f, a, depth): - # we ignore depth, as numpy arrays will handle the vectorized cases. - a = array(a) - if len(a.shape) == 1: - # if a is a flat array, we embed the components in an array, i.e. - # [[a, b, c]], so that unstack will yield [a], [b], [c], and we - # operate on arrays of length 1 instead of scalars. - return f(array([a]))[0] - else: - return f(a) - - def unstack(a): a = array(a) b = array(numpy.split(a, a.shape[-1], axis=-1)) @@ -64,29 +60,32 @@ def stack(*a): return b +def stacked(f, a): + a = array(a) + unwrap = False + if len(a.shape) == 1: + a = array([a]) + unwrap = True + components = unstack(a) + result = f(*components) + result = stack(*result) + if unwrap: + result = result[0] + return result + + def concat(*a): a = [array(x) for x in a] a = [x for x in a if x.shape[0]] # skip empty return numpy.concatenate(a, axis=-1) -def conditional(a, cond, t, f): - b = array(a)[:] - mask = cond(a) - b[mask] = t(b[mask]) - b[~mask] = f(b[~mask]) - return b - +def vectorize(a, depth, f): + return f(a) -def compose(*a): - assert a and len(a) % 2 == 0 - b = numpy.ndarray(a[0].shape) - # we apply the rules in reversed order, so that the first rules - # always make the last changes, which corresponds to the unreversed - # processing in the non-vectorized (non-numpy) implementation. - for mask, v in reversed([a[i:i + 2] for i in range(0, len(a), 2)]): - b[mask] = v(lambda x: x[mask]) - return b +# +# MATHEMATICAL OPERATIONS +# def clip(a, t0, t1): @@ -129,6 +128,87 @@ def minimum(*a): return reduce(numpy.minimum, [array(x) for x in a]) +# +# PUBLIC HELPER FUNCTIONS +# + +def is_numpy_available(): + return True + + +def allclose(a, b): + return numpy.allclose(array(a), array(b)) + + +def errstate(**kwargs): + return numpy.errstate(**kwargs) + + +def instantiate_elements(a, new_element, d=1): + # given a numpy array 'a' and a python element constructor 'new_element', generate a python array of the + # same shape as 'a' with python elements constructed through 'new_element'. 'new_element' will get called + # if an array of dimension 'd' is reached. + + if len(a.shape) == d: + leaves = [new_element(x) for x in a] + else: + leaves = [instantiate_elements(e, new_element, d) for e in a] + return Expression('List', *leaves) + + +# +# CONDITIONALS AND PROGRAM FLOW +# + +# @conditional is an annotation that basically invoked a mini compiler in order to compile numpy +# conditional expressions from code that would run as regular Python code in the no-numpy case. +# the idea is to have readable Python code that states what happens in both numpy and no-numpy +# cases. the alternative would be using wrapper functions that are kind of ugly and hard to read +# for the numpy case. +# +# as an example, take the following code: +# +# @conditional +# def f(a): +# if a > 10: +# return a + 1.5 +# elif a > 9: +# return 7 +# else: +# return 2 +# +# if numpy is not available, f() will just take one scalar and work as expected. if numpy is +# available, @conditional will recompile the function into something that allows "a" (or any +# other parameters of f) to be numpy arrays. +# +# this is necessary, as for numpy arrays, there is no single execution branch in conditionals +# as above, since the conditional has to be evaluated for each element of the numpy array. +# internally, in the numpy case, the function above will be transformed to something like: +# +# a[a > 10] = a + 1.5 +# ... +# +# in general, and to make this transformation as simple as possible, @conditional expects the +# function that is annotated with it to be of the following restrained form: +# +# if [binary comparisons or variable]: +# return [expression1] +# elif [binary comparisons or variable]: # optional +# return [expression2] +# else: # optional +# return [expression3] +# +# all relevant rules for @conditional functions are: +# - all "if" branches must exit immediately with "return". +# - "if"s must rely on simple binary comparisons, e.g. "b < 4" or "4 > b", or variables +# - the occurence of "elif" is optional, as is the occurence of "else" +# - if "else" is not provided, the provided "if" cases must cover all possible cases, +# otherwise there will be undefined results. +# - code in @conditional must not reference global variables. +# +# if a function does not adhere to the rules described above, a MalformedConditional is thrown when +# constructing the @vectorized function. + def choose(i, *options): assert options dim = len(options[0]) @@ -141,15 +221,206 @@ def choose(i, *options): return [numpy.choose(i_int, column) for column in columns] -def allclose(a, b): - return numpy.allclose(array(a), array(b)) +_else_case_id = 'else_case' -def errstate(**kwargs): - return numpy.errstate(**kwargs) +def _numpy_conditional(shape_id, *paths): + # called during runtime when we actually evaluate a conditional. + if _is_scalar(shape_id): + for test, args, comp in paths: + if test == _else_case_id or test(): + return comp(*args) + assert False # one case must be true -def constant(x, a): - b = numpy.ndarray(a.shape) - b.fill(x) - return b + shape = shape_id.shape + result = numpy.ndarray(shape) + + has_else = paths[-1][0] == _else_case_id + if_paths = paths[:-1] if has_else else paths + masks = [path[0]() for path in if_paths] + + def efficient_masks(): + rest = ~masks[0] + yield masks[0] + for mask in masks[1:]: + yield rest & mask + rest &= ~mask + + if has_else: + if len(paths) == 2: # simple case: just 1 if, and 1 else + masks.append(~masks[0]) + else: # n ifs, and 1 else + masks.append(numpy.ones(masks[0].shape, dtype=bool)) + masks = list(efficient_masks()) + else: + masks = list(efficient_masks()) + + # process in reverse order in order to reflect order that is written + # down in Python code. for example, writing + # + # if a > 5: + # b = True + # elif a > 4: + # b = False + # + # in Python, will always set b to True, if a > 5. if we evaluated the + # statements above in non-reversed order for numpy, for some element + # a > 5, we would first set b = True, and then b = False, which would + # be wrong. + + for mask, path in zip(reversed(masks), reversed(paths)): + test, args, comp = path + result[mask] = comp(*[x if _is_scalar(x) else x[mask] for x in args]) + return result + + +class MalformedConditional(Exception): + def __init__(self, func, node, error): + Exception.__init__(self, 'in function %s in line %d: %s' % (func.__name__, node.lineno, error)) + + +class _NameCollector(ast.NodeVisitor): + def __init__(self): + self.names = set() + + def visit_Name(self, node): + assert isinstance(node.ctx, ast.Load) + self.names.add(node.id) + + +def _create_ast_lambda(names, body): + if sys.version_info >= (3, 0): # change in AST structure for Python 3 + args = [ast.arg(arg=name, annotation=None) for name in names] + else: + args = [ast.Name(id=name, ctx=ast.Load()) for name in names] + + return ast.Lambda(args=ast.arguments( + args=args, vararg=None, kwonlyargs=[], kw_defaults=[], kwarg=None, defaults=[]), body=body) + + +def _expression_lambda(node): + # convert some expression, e.g. b + a * 2 + b * 7, into a lambda with with all values used in the + # expression as parameters, e.g. lambda a, b: b + a * 2 + b * 7 + + collector = _NameCollector() + collector.visit(node) + names = list(collector.names) + elements = [ast.Name(id=name, ctx=ast.Load()) for name in names] + + # return (1) a tuple of all the variable names in the expression, (2) a lambda that takes values + # for these variables as parameters and evaluates the expression. + return ast.Tuple(elts=elements, ctx=ast.Load()), _create_ast_lambda(names, node) + + +class _ConditionalTransformer(ast.NodeTransformer): + def __init__(self, f): + self._func = f + + def transform(self): + tree = ast.parse(inspect.getsource(self._func)) + self._expect(tree, 'function not a Module', type(tree), ast.Module) + self._expect(tree, 'Module body too large', len(tree.body), 1) + self._expect(tree, 'FunctionDef not found', type(tree.body[0]), ast.FunctionDef) + + func_def = tree.body[0] + func_name = func_def.name + self._expect(func_def, 'FunctionDef body too large', len(func_def.body), 1) + self._expect(func_def, 'function must start with "if"', type(func_def.body[0]), ast.If) + + tree = self.visit(tree) + tree = ast.fix_missing_locations(tree) + code = compile(tree, ''% func_name, 'exec') + + data = {} + eval(code, globals(), data) + return data[func_name] + + def _expect(self, node, error, value, expected): + if value != expected: + raise MalformedConditional(self._func, node, "%s (%s != %s)" % (error, expected, value)) + + def visit_FunctionDef(self, node): + assert len(node.decorator_list) == 1 # we expect that we are the only decorator + assert isinstance(node.decorator_list[0], ast.Name) + assert node.decorator_list[0].id == 'conditional' + return ast.FunctionDef(name=node.name, args=node.args, + body=[self.visit(x) for x in node.body], decorator_list=[]) + + def visit_If(self, node): + blocks = [] + tests = [] + shapes = [] + + while True: + body = node.body + self._expect(node, '"if" code body must contain exactly 1 element', len(body), 1) + blocks.append(body[0]) + + test = node.test + tests.append(test) + + if type(test) == ast.Name: + shapes.append(test) + elif type(test) == ast.Compare: + if isinstance(test.left, ast.Name): + shapes.append(test.left) + elif isinstance(test.right, ast.Name): + shapes.append(test.right) + else: + MalformedConditional(self._func, test, 'expected variable in comparison') + else: + self._expect(test, 'expected single comparison or name', type(test), ast.Compare) + + or_elses = node.orelse + if not or_elses: + break + + self._expect(node, '"else" code body must contain 1 element', len(or_elses), 1) + or_else = or_elses[0] + + if isinstance(or_else, ast.If): + node = or_else + else: + blocks.append(or_else) + tests.append(_else_case_id) + break + + for block in blocks: + self._expect(block, '"if" blocks must exit with "return"', type(block), ast.Return) + + # now build a call to _numpy_conditional() using cond_args as arguments, that will handle + # the runtime evaluation of the conditional. + + cond_args = [shapes[0]] + for test, value in zip(tests, (block.value for block in blocks)): + elements = [] + if test == _else_case_id: + elements.append(ast.Str(s=test)) + else: + elements.append(_create_ast_lambda([], test)) + + elements.extend(_expression_lambda(value)) + cond_args.append(ast.Tuple(elts=elements, ctx=ast.Load())) + + return ast.Return(value=ast.Call( + func=ast.Name(id='_numpy_conditional', ctx=ast.Load()), keywords=[], args=cond_args)) + + +def conditional(*args, **kwargs): + if len(args) == 1 and callable(args[0]): + f = args[0] # @conditional without arguments? + else: + return lambda f: conditional(f) # with arguments + + if not inspect.isfunction(f): + raise Exception('@conditional can only be applied to functions') + + transformer = _ConditionalTransformer(f) + f_transformed = transformer.transform() + + def wrapper(*a): + a = [numpy.array(x) if isinstance(x, (list, tuple)) else x for x in a] + return f_transformed(*a) + + return wrapper diff --git a/mathics/builtin/numpy_utils/without_numpy.py b/mathics/builtin/numpy_utils/without_numpy.py old mode 100644 new mode 100755 index 83da0a0b27..4f6974b4ae --- a/mathics/builtin/numpy_utils/without_numpy.py +++ b/mathics/builtin/numpy_utils/without_numpy.py @@ -10,52 +10,42 @@ from contextlib import contextmanager from math import sin as sinf, cos as cosf, sqrt as sqrtf, atan2 as atan2f, floor as floorf import operator +import inspect # If numpy is not available, we define the following fallbacks that are useful for implementing a similar # logic in pure python without numpy. They obviously work on regular python array though, not numpy arrays. -def is_numpy_available(): - return False +# +# INTERNAL FUNCTIONS +# +def _is_bottom(a): + return any(not isinstance(x, list) for x in a) -def instantiate_elements(a, new_element, d=1): - # given a python array 'a' and a python element constructor 'new_element', generate a python array of the - # same shape as 'a' with python elements constructed through 'new_element'. 'new_element' will get called - # if an array of dimension 'd' is reached. - e = a[0] - depth = 1 - while depth <= d and isinstance(e, list): - e = e[0] - depth += 1 - if d == depth: - leaves = [new_element(x) for x in a] +def _apply(f, a): + if isinstance(a, (list, tuple)): + return [_apply(f, t) for t in a] else: - leaves = [instantiate_elements(e, new_element, d) for e in a] - return Expression('List', *leaves) + return f(a) -def _is_bottom(a): - return any(not isinstance(x, list) for x in a) +def _apply_n(f, *p): + if isinstance(p[0], list): + return [_apply_n(f, *q) for q in zip(*p)] + else: + return f(*p) + +# +# ARRAY CREATION AND REORGANIZATION: STACK, UNSTACK, CONCAT, ... +# def array(a): return a -def vectorized(f, a, depth): - # depth == 0 means: f will get only scalars. depth == 1 means: f will - # get lists of scalars. - if _is_bottom(a): - if depth == 0: - return [f(x) for x in a] - else: - return f(a) - else: - return [vectorized(f, x, depth) for x in a] - - def unstack(a): if not a: return [] @@ -79,38 +69,31 @@ def stack(*a): return [stack(*[x[i] for x in a]) for i in range(len(a[0]))] +def stacked(f, a): + components = unstack(a) + result = f(*components) + return stack(*result) + + def concat(a, b): return stack(*(unstack(a) + unstack(b))) -def _apply(f, a): - if isinstance(a, (list, tuple)): - return [_apply(f, t) for t in a] +def vectorize(a, depth, f): + # depth == 0 means: f will get only scalars. depth == 1 means: f will + # get lists of scalars. + if _is_bottom(a): + if depth == 0: + return [f(x) for x in a] + else: + return f(a) else: - return f(a) - + return [vectorize(x, f, depth) for x in a] -def _apply_n(f, *p): - if isinstance(p[0], list): - return [_apply_n(f, *q) for q in zip(*p)] - else: - return f(*p) - - -def conditional(a, cond, t, f): - def _eval(x): - return t(x) if cond(x) else f(x) - - return _apply(_eval, a) - - -def compose(*a): - assert a and len(a) % 2 == 0 - for cond, v in (a[i:i + 2] for i in range(0, len(a), 2)): - if cond: - return v(lambda x: x) - raise ValueError('no matching case in compose') +# +# MATHEMATICAL OPERATIONS +# def clip(a, t0, t1): def _eval(x): @@ -158,18 +141,12 @@ def minimum(*a): return _apply_n(min, *a) -def _choose_descend(i, options): - if isinstance(i, (int, float)): - return options[int(i)] # int cast needed for PyPy - else: - return [_choose_descend(next_i, [o[k] for o in options]) for k, next_i in enumerate(i)] +# +# PUBLIC HELPER FUNCTIONS +# - -def choose(i, *options): - assert options - dim = len(options[0]) - columns = [[o[d] for o in options] for d in range(dim)] - return [_choose_descend(i, column) for column in columns] +def is_numpy_available(): + return False def allclose(a, b): @@ -188,5 +165,52 @@ def errstate(**kwargs): yield -def constant(x, a): - return x +def instantiate_elements(a, new_element, d=1): + # given a python array 'a' and a python element constructor 'new_element', generate a python array of the + # same shape as 'a' with python elements constructed through 'new_element'. 'new_element' will get called + # if an array of dimension 'd' is reached. + + e = a[0] + depth = 1 + while depth <= d and isinstance(e, list): + e = e[0] + depth += 1 + if d == depth: + leaves = [new_element(x) for x in a] + else: + leaves = [instantiate_elements(e, new_element, d) for e in a] + return Expression('List', *leaves) + + +# +# CONDITIONALS AND PROGRAM FLOW +# + + +def _choose_descend(i, options): + if isinstance(i, (int, float)): + return options[int(i)] # int cast needed for PyPy + else: + return [_choose_descend(next_i, [o[k] for o in options]) for k, next_i in enumerate(i)] + + +def choose(i, *options): + assert options + dim = len(options[0]) + columns = [[o[d] for o in options] for d in range(dim)] + return [_choose_descend(i, column) for column in columns] + + +def conditional(*args): # essentially a noop + if len(args) == 1 and callable(args[0]): + f = args[0] # @conditional without arguments? + else: + return lambda f: conditional(f) # with arguments + + if not inspect.isfunction(f): + raise Exception('@conditional can only be applied to functions') + + def wrapper(*a): + return f(*a) + + return wrapper diff --git a/test/test_color.py b/test/test_color.py old mode 100644 new mode 100755 index 3daec4ad18..ad7704b3f9 --- a/test/test_color.py +++ b/test/test_color.py @@ -10,7 +10,7 @@ from six.moves import range import mathics.builtin.colors as colors -from mathics.builtin.numpy_utils import array, vectorized +from mathics.builtin.numpy_utils import array, stacked, vectorize from mathics.core.definitions import Definitions from mathics.core.evaluation import Evaluation from mathics.core.expression import Expression @@ -264,7 +264,7 @@ def testImageConversions(self): for name, convert in colors.conversions.items(): if name.find('CMYK') < 0: - self._checkImageConversion(4, convert) + self._checkImageConversion(4, lambda p: vectorize(p, 1, lambda q: stacked(convert, q))) def _checkConversion(self, from_space, from_components, to_space, to_components): places = 12 @@ -282,10 +282,10 @@ def _checkConversion(self, from_space, from_components, to_space, to_components) def _checkImageConversion(self, size, convert): pixels = [[random(), random(), random()] for _ in range(size * size)] - refs = [vectorized(convert, p, 1) for p in pixels] + refs = [convert(p) for p in pixels] image = [[pixels[x * size + y] for y in range(size)] for x in range(size)] - image = vectorized(convert, array(image), 1) + image = convert(array(image)) for x in range(size): for y in range(size): diff --git a/test/test_numpy_utils.py b/test/test_numpy_utils.py index f5053efcab..e8a5320fbf 100644 --- a/test/test_numpy_utils.py +++ b/test/test_numpy_utils.py @@ -6,8 +6,25 @@ import unittest -from mathics.builtin.numpy_utils import stack, unstack, concat, conditional, compose, clip, array, choose -from mathics.builtin.numpy_utils import vectorized, minimum, maximum, dot_t, mod, floor, sqrt, allclose +from mathics.builtin.numpy_utils import stack, unstack, concat, vectorize, conditional, clip, array, choose +from mathics.builtin.numpy_utils import minimum, maximum, dot_t, mod, floor, sqrt, allclose + + +@conditional +def _test_simple_conditional(t): + if t > 0.5: + return t + 1. + else: + return -t + +@conditional +def _test_complex_conditional(t): + if t > 10: + return t * 10 + 1 + elif t > 3: + return t * 10 + elif t <= 3: + return -1 class Numpy(unittest.TestCase): @@ -54,29 +71,6 @@ def testConcatComplex(self): c = [[[1, 2, 3], [4, 5, 6]], [[7, 8, 9], [10, 11, 12]]] self.assertEqualArrays(concat(a, b), c) - def testConditional(self): - # conditional operates on each element independently. - a = array([[[0.1, 0.6], [1.8, 0.4]], [[-0.1, -0.8], [1.1, 0.5]]]) - a = conditional(a, lambda t: t > 0.5, # if - lambda t: t + 1., # true - lambda t: -t) # false - self.assertEqualArrays(a, [[[-0.1, 1.6], [2.8, -0.4]], [[0.1, 0.8], [2.1, -0.5]]]) - - def testCompose(self): - a = array([[[1, 2], [4, 5]], [[7, 8], [10, 11]]]) - - def f(a): - return compose( - a > 10, - lambda s: s(a * 10 + 1), - a > 3, - lambda s: s(a * 10), - a < 3, - lambda s: -1) - - a = vectorized(f, a, 0) - self.assertEqualArrays(a, [[[-1, -1], [40, 50]], [[70, 80], [100, 111]]]) - def testChooseSimple(self): # select a single value from a list of values. options = [ @@ -162,6 +156,16 @@ def testFloor(self): def testSqrt(self): self.assertEqualArrays(sqrt([[9, 100], [25, 16]]), [[3, 10], [5, 4]]) + def testSimpleConditional(self): + a = array([[[0.1, 0.6], [1.8, 0.4]], [[-0.1, -0.8], [1.1, 0.5]]]) + a = vectorize(a, 0, _test_simple_conditional) + self.assertEqualArrays(a, [[[-0.1, 1.6], [2.8, -0.4]], [[0.1, 0.8], [2.1, -0.5]]]) + + def testConditionalComplex(self): + a = array([[[1, 2], [4, 5]], [[7, 8], [10, 11]]]) + a = vectorize(a, 0, _test_complex_conditional) + self.assertEqualArrays(a, [[[-1, -1], [40, 50]], [[70, 80], [100, 111]]]) + def assertEqualArrays(self, a, b): self.assertEqual(allclose(a, b), True) From bb6c45adc197ee05a2b000e5fdf31ea973077e6c Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Thu, 23 Jun 2016 18:51:06 +0200 Subject: [PATCH 27/30] fixes identation --- mathics/builtin/colors.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mathics/builtin/colors.py b/mathics/builtin/colors.py index f6a91a0401..8bf33bfe9e 100644 --- a/mathics/builtin/colors.py +++ b/mathics/builtin/colors.py @@ -160,9 +160,9 @@ def cmyk_to_rgb(c, m, y, *rest): @conditional def _scale_rgb_to_cmyk(t, k, k_, eps): if k_ < eps: - return 0 + return 0 else: - return (1 - t - k) / k_ + return (1 - t - k) / k_ def rgb_to_cmyk(r, g, b, *rest): From eb9ddab4ba9d032588fc40e6c4616338d7eb672b Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Wed, 24 Aug 2016 10:25:35 +0200 Subject: [PATCH 28/30] rebase, central color conversion, more robust test cases --- mathics/builtin/colors.py | 14 ++++- mathics/builtin/graphics.py | 82 ++++--------------------- mathics/builtin/image.py | 118 +++++++++++++++++++++++------------- test/test_color.py | 35 +++++++---- 4 files changed, 123 insertions(+), 126 deletions(-) diff --git a/mathics/builtin/colors.py b/mathics/builtin/colors.py index 8bf33bfe9e..84319748d4 100644 --- a/mathics/builtin/colors.py +++ b/mathics/builtin/colors.py @@ -56,6 +56,15 @@ def __init__(self, *x): # accepts (x1, y2) or (x, y, z) ] +def omit_alpha(*c): + if len(c) == 2: + return (c,) + elif len(c) == 4: + return c[:3] + else: + return c + + def rgb_to_grayscale(r, g, b, *rest): # see https://en.wikipedia.org/wiki/Grayscale y = 0.299 * r + 0.587 * g + 0.114 * b # Y of Y'UV @@ -345,7 +354,10 @@ def _flow(src, dst): } -def convert(components, src, dst): +def convert(components, src, dst, preserve_alpha=True): + if not preserve_alpha: + components = stacked(omit_alpha, components) + if src == dst: return components diff --git a/mathics/builtin/graphics.py b/mathics/builtin/graphics.py index 6a78be20c2..2b89a6e092 100644 --- a/mathics/builtin/graphics.py +++ b/mathics/builtin/graphics.py @@ -165,29 +165,6 @@ def _data_and_options(leaves, defined_options): return data, options -class _PerfectReflectingDiffuser: - def __init__(self, x1, x2): - # MMA seems to use the following constants, and not the - # values derived via calculation (commented out below) - self.xyz = (0.96422, 1., 0.82521) - #scale = 1.0 / x2 - #self.xyz = (x1 * scale, 1.0, (1.0 - x1 - x2) * scale) - - q_r = self.xyz[0] + 15. * self.xyz[1] + 3. * self.xyz[2] - self.u_r = 4. * self.xyz[0] / q_r - self.v_r = 9. * self.xyz[1] / q_r - - -# MMA's reference white is a # D50, 2 degrees diffuser; for the -# values, see https://en.wikipedia.org/wiki/Standard_illuminant - -_ref_white = _PerfectReflectingDiffuser(0.34567, 0.35850) - - -def _clip(*components): # clip to [0, 1] - return tuple(max(0, min(1, c)) for c in components) - - def _euclidean_distance(a, b): return math.sqrt(sum((x1 - x2) * (x1 - x2) for x1, x2 in zip(a, b))) @@ -549,55 +526,22 @@ class GrayLevel(_Color): default_components = [0, 1] -class ColorConvert(Builtin): - """ -
-
'ColorConvert[$c$, $colspace$]' -
returns the representation of color $c$ in the color space $colspace$. -
- - Valid values for $colspace$ are: +def expression_to_color(color): + try: + return _Color.create(color) + except ColorError: + return None - CMYK: convert to CMYKColor - Grayscale: convert to GrayLevel - HSB: convert to Hue - LAB: concert to LABColor - LCH: convert to LCHColor - LUV: convert to LUVColor - RGB: convert to RGBColor - XYZ: convert to XYZColor - """ - messages = { - 'ccvinput': '`` should be a color.', - 'imgcstype': '`` is not a valid color space.', - } - - def apply(self, color, colorspace, evaluation): - 'ColorConvert[color_, colorspace_String]' - try: - py_color = _Color.create(color) - except ColorError: - evaluation.message('ColorConvert', 'ccvinput', color) - return - - py_colorspace = colorspace.get_string_value() - converted_components = convert_color(py_color.components, - py_color.color_space, - py_colorspace) - - if converted_components is None: - evaluation.message('ColorConvert', 'imgcstype', colorspace) - return - - if py_colorspace == 'Grayscale': - converted_color_name = 'GrayLevel' - elif py_colorspace == 'HSB': - converted_color_name = 'Hue' - else: - converted_color_name = py_colorspace + 'Color' +def color_to_expression(components, colorspace): + if colorspace == 'Grayscale': + converted_color_name = 'GrayLevel' + elif colorspace == 'HSB': + converted_color_name = 'Hue' + else: + converted_color_name = colorspace + 'Color' - return Expression(converted_color_name, *converted_components) + return Expression(converted_color_name, *components) class ColorDistance(Builtin): diff --git a/mathics/builtin/image.py b/mathics/builtin/image.py index 0ecce9c992..c6c1c55b17 100644 --- a/mathics/builtin/image.py +++ b/mathics/builtin/image.py @@ -10,6 +10,7 @@ Builtin, AtomBuiltin, Test, BoxConstruct, String) from mathics.core.expression import ( Atom, Expression, Integer, Rational, Real, Symbol, from_python) +from mathics.builtin.colors import convert as convert_color import six import base64 @@ -54,33 +55,6 @@ from io import BytesIO -if _enabled: - _color_space_conversions = { - 'RGB2Grayscale': skimage.color.rgb2gray, - 'Grayscale2RGB': skimage.color.gray2rgb, - - 'HSV2RGB': skimage.color.hsv2rgb, - 'RGB2HSV': skimage.color.rgb2hsv, - - 'LAB2LCH': skimage.color.lab2lch, - 'LCH2LAB': skimage.color.lch2lab, - - 'LAB2RGB': skimage.color.lab2rgb, - 'LAB2XYZ': skimage.color.lab2xyz, - - 'LUV2RGB': skimage.color.luv2rgb, - 'LUV2XYZ': skimage.color.luv2xyz, - - 'RGB2LAB': skimage.color.rgb2lab, - 'RGB2LUV': skimage.color.rgb2luv, - 'RGB2XYZ': skimage.color.rgb2xyz, - - 'XYZ2LAB': skimage.color.xyz2lab, - 'XYZ2LUV': skimage.color.xyz2luv, - 'XYZ2RGB': skimage.color.xyz2rgb, - } - - class _ImageBuiltin(Builtin): requires = _image_requires @@ -843,7 +817,7 @@ class _MorphologyFilter(_ImageBuiltin): def apply(self, image, k, evaluation): '%(name)s[image_Image, k_?MatrixQ]' if image.color_space != 'Grayscale': - image = image.color_convert('Grayscale') + image = image.grayscale() evaluation.message(self.name, 'grayscale') f = getattr(skimage.morphology, self.get_name(True).lower()) img = f(image.pixels, numpy.array(k.to_python())) @@ -923,15 +897,61 @@ def apply(self, image, evaluation): class ColorConvert(_ImageBuiltin): - def apply(self, image, colorspace, evaluation): - 'ColorConvert[image_Image, colorspace_String]' - return image.color_convert(colorspace.get_string_value()) + """ +
+
'ColorConvert[$c$, $colspace$]' +
returns the representation of $c$ in the color space $colspace$. $c$ + may be a color or an image. +
+ + Valid values for $colspace$ are: + + CMYK: convert to CMYKColor + Grayscale: convert to GrayLevel + HSB: convert to Hue + LAB: concert to LABColor + LCH: convert to LCHColor + LUV: convert to LUVColor + RGB: convert to RGBColor + XYZ: convert to XYZColor + """ + + messages = { + 'ccvinput': '`` should be a color.', + 'imgcstype': '`` is not a valid color space.', + } + + def apply(self, input, colorspace, evaluation): + 'ColorConvert[input_, colorspace_String]' + + if isinstance(input, Image): + return input.color_convert(colorspace.get_string_value()) + else: + from mathics.builtin.graphics import expression_to_color, color_to_expression + + py_color = expression_to_color(input) + if py_color is None: + evaluation.message('ColorConvert', 'ccvinput', input) + return + + py_colorspace = colorspace.get_string_value() + converted_components = convert_color( + py_color.components, py_color.color_space, py_colorspace) + + if converted_components is None: + evaluation.message('ColorConvert', 'imgcstype', colorspace) + return + + return color_to_expression(converted_components, py_colorspace) class ColorQuantize(_ImageBuiltin): def apply(self, image, n, evaluation): 'ColorQuantize[image_Image, n_Integer]' - pixels = skimage.img_as_ubyte(image.color_convert('RGB').pixels) + converted = image.color_convert('RGB') + if converted is None: + return + pixels = skimage.img_as_ubyte(converted.pixels) im = PIL.Image.fromarray(pixels).quantize(n.to_python()) im = im.convert('RGB') return Image(numpy.array(im), 'RGB') @@ -1086,7 +1106,11 @@ def apply(self, image, n, evaluation): class PixelValue(_ImageBuiltin): def apply(self, image, x, y, evaluation): 'PixelValue[image_Image, {x_?RealNumberQ, y_?RealNumberQ}]' - return Real(image.pixels[int(y.to_python() - 1), int(x.to_python() - 1)]) + pixel = image.pixels[int(y.to_python() - 1), int(x.to_python() - 1)] + if isinstance(pixel, (numpy.ndarray, numpy.generic, list)): + return Expression('List', *[Real(float(x)) for x in list(pixel)]) + else: + return Real(float(pixel)) class PixelValuePositions(_ImageBuiltin): @@ -1214,25 +1238,31 @@ def __init__(self, pixels, color_space, **kwargs): def pil(self): return PIL.Image.fromarray(self.pixels) - def color_convert(self, to_color_space): + def color_convert(self, to_color_space, preserve_alpha=True): if to_color_space == self.color_space: return self else: - conversion = '%s2%s' % (self.color_space, to_color_space) - if conversion in _color_space_conversions: - return Image(_color_space_conversions[conversion](self.pixels), to_color_space) - else: - raise ValueError('cannot convert from color space %s to %s' % (self.color_space, to_color_space)) + pixels = skimage.img_as_float(self.pixels) + converted = convert_color(pixels, self.color_space, to_color_space, preserve_alpha) + if converted is None: + return None + return Image(converted, to_color_space) def grayscale(self): return self.color_convert('Grayscale') def atom_to_boxes(self, f, evaluation): try: - if self.color_space == 'Grayscale': + if self.color_space == 'Grayscale' and self.pixels.shape[2] == 1: pixels = self.pixels + # convert_color gives us [[[a], [b], ...]], but + # skimage.io.imsave only wants [[a, b, ...]], so: + pixels = pixels.reshape(pixels.shape[:2]) else: - pixels = self.color_convert('RGB').pixels + pixels = self.color_convert('RGB', False).pixels + + if pixels is None: + raise ValueError('could not convert pixels to RGB.') if pixels.dtype == numpy.bool: pixels = skimage.img_as_ubyte(pixels) @@ -1267,8 +1297,10 @@ def atom_to_boxes(self, f, evaluation): encoded = 'data:image/png;base64,' + encoded return Expression('ImageBox', String(encoded), Integer(scaled_width), Integer(scaled_height)) - except: - return Symbol("$Failed") + except Exception as e: + evaluation.print_out('Image processing failed: %s' % str(e)) + return String('') + def __str__(self): return '-Image-' diff --git a/test/test_color.py b/test/test_color.py index ad7704b3f9..5db1cd1ee0 100755 --- a/test/test_color.py +++ b/test/test_color.py @@ -196,28 +196,36 @@ def testInverseConversions(self): # components. this tests color space transformations and their # inverse transformations. + def space_to_head(name): + if name == 'HSB': + return 'System`Hue' + else: + return 'System`%sColor' % name + spaces = ("CMYK", "HSB", "LAB", "LCH", "LUV", "RGB", "XYZ") places = 3 for original in ((0.5, 0.1, 0.2), (0.9, 0.1, 0.1)): for i, from_space in enumerate(spaces): for to_space in spaces[i + 1:]: try: - if from_space == 'HSB': - construct_name = 'Hue' - else: - construct_name = from_space + 'Color' + construct_name = space_to_head(from_space) + source_color = Expression(construct_name, *original) # now calculate from_space -> to_space -> from_space - inverted = [c.to_python() for c in Expression('ColorConvert', - Expression('ColorConvert', - Expression(construct_name, *original), - to_space), - from_space).evaluate(self.evaluation).leaves] + target_color = Expression('ColorConvert', source_color, to_space).evaluate(self.evaluation) + self.assertEqual(target_color.get_head_name(), space_to_head(to_space)) + + checked_color = Expression('ColorConvert', target_color, from_space).evaluate(self.evaluation) + self.assertEqual(checked_color.get_head_name(), source_color.get_head_name()) + + checked_components = [c.to_python() for c in checked_color.leaves] if from_space == 'CMYK': # if cmyk, cmyk -> cmy - k = inverted[3] - inverted = [c * (1 - k) + k for c in inverted[:3]] - self.assertEqual(len(original), len(inverted)) - for x, y in zip(original, inverted): + k = checked_components[3] + checked_components = [c * (1 - k) + k for c in checked_components[:3]] + + self.assertEqual(len(original), len(checked_components)) + + for x, y in zip(original, checked_components): self.assertAlmostEqual(x, y, places) except: print('test failed for %s -> %s -> %s' % @@ -244,6 +252,7 @@ def testConversions(self): for components, t1 in zip(values, _color_tests): for src, t2 in zip(spaces, t1): for dst, expected in zip(spaces, t2): + components = list(components) if src == "Grayscale": c = components[:1] elif src == 'CMYK': From 33ec0c5f32f9c21dad106d09e8f193090967a004 Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Thu, 25 Aug 2016 20:41:55 +0200 Subject: [PATCH 29/30] reduced color convert test case accuracy, make ColorConvert[] usable in non-scipy installations --- mathics/builtin/image.py | 2 +- test/test_color.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/mathics/builtin/image.py b/mathics/builtin/image.py index c6c1c55b17..074307edf3 100644 --- a/mathics/builtin/image.py +++ b/mathics/builtin/image.py @@ -896,7 +896,7 @@ def apply(self, image, evaluation): return String(image.color_space) -class ColorConvert(_ImageBuiltin): +class ColorConvert(Builtin): """
'ColorConvert[$c$, $colspace$]' diff --git a/test/test_color.py b/test/test_color.py index 5db1cd1ee0..bccd0a2cce 100755 --- a/test/test_color.py +++ b/test/test_color.py @@ -276,7 +276,7 @@ def testImageConversions(self): self._checkImageConversion(4, lambda p: vectorize(p, 1, lambda q: stacked(convert, q))) def _checkConversion(self, from_space, from_components, to_space, to_components): - places = 12 + places = 6 if from_space == 'HSB': construct_name = 'Hue' From d3ef32acab314fc2c145d890b0a23e57ac9dd2f8 Mon Sep 17 00:00:00 2001 From: Bernhard Liebl Date: Thu, 1 Sep 2016 09:19:38 +0200 Subject: [PATCH 30/30] cleaner color conversion paths --- mathics/builtin/colors.py | 105 +++++++++++++++++++++++++++----------- 1 file changed, 74 insertions(+), 31 deletions(-) diff --git a/mathics/builtin/colors.py b/mathics/builtin/colors.py index 84319748d4..76f33f454b 100644 --- a/mathics/builtin/colors.py +++ b/mathics/builtin/colors.py @@ -308,33 +308,75 @@ def lab_to_xyz(l, a, b, *rest): return (x, y, z) + rest - -_flows = { # see http://www.brucelindbloom.com/Math.html - 'XYZ>LCH': ('XYZ', 'LAB', 'LCH'), - 'LAB>LUV': ('LAB', 'XYZ', 'LUV'), - 'LAB>RGB': ('LAB', 'XYZ', 'RGB'), - 'LCH>XYZ': ('LCH', 'LAB', 'XYZ'), - 'LCH>LUV': ('LCH', 'LAB', 'XYZ', 'LUV'), - 'LCH>RGB': ('LCH', 'LAB', 'XYZ', 'RGB'), - 'LUV>LAB': ('LUV', 'XYZ', 'LAB'), - 'LUV>LCH': ('LUV', 'XYZ', 'LAB', 'LCH'), - 'LUV>RGB': ('LUV', 'XYZ', 'RGB'), - 'RGB>LAB': ('RGB', 'XYZ', 'LAB'), - 'RGB>LCH': ('RGB', 'XYZ', 'LAB', 'LCH'), - 'RGB>LUV': ('RGB', 'XYZ', 'LUV'), -} - - -_rgb_flows = set(['Grayscale', 'CMYK', 'HSB']) - - -def _flow(src, dst): - if (src in _rgb_flows and dst != 'RGB') or (dst in _rgb_flows and src != 'RGB'): - return list(chain(_flow(src, 'RGB'), _flow('RGB', dst))) - else: - r = _flows.get('%s>%s' % (src, dst)) - return list(r) if r else [src, dst] - +# for an overview of color conversions see http://www.brucelindbloom.com/Math.html + +# the following table was computed by starting with the allowed hard +# coded conversions from "conversions" and then finding the shortest +# paths: + +# g = Graph[{"Grayscale" -> "RGB", "RGB" -> "Grayscale", "CMYK" -> "RGB", "RGB" -> "CMYK", "RGB" -> "HSB", +# "HSB" -> "RGB", "XYZ" -> "LAB", "XYZ" -> "LUV", "XYZ" -> "RGB", "LAB" -> "XYZ", "LAB" -> "LCH", +# "LCH" -> "LAB", "LUV" -> "XYZ", "RGB" -> "XYZ"}] +# s = FindShortestPath[g, All, All]; {#, s @@ #} & /@ Permutations[{ +# "Grayscale", "RGB", "CMYK", "HSB", "XYZ", "LAB", "LUV", "LCH"}, {2}] // CForm + +_paths = dict(( + (("Grayscale","RGB"),("Grayscale","RGB")), + (("Grayscale","CMYK"),("Grayscale","RGB","CMYK")), + (("Grayscale","HSB"),("Grayscale","RGB","HSB")), + (("Grayscale","XYZ"),("Grayscale","RGB","XYZ")), + (("Grayscale","LAB"),("Grayscale","RGB","XYZ","LAB")), + (("Grayscale","LUV"),("Grayscale","RGB","XYZ","LUV")), + (("Grayscale","LCH"),("Grayscale","RGB","XYZ","LAB","LCH")), + (("RGB","Grayscale"),("RGB","Grayscale")), + (("RGB","CMYK"),("RGB","CMYK")), + (("RGB","HSB"),("RGB","HSB")), + (("RGB","XYZ"),("RGB","XYZ")), + (("RGB","LAB"),("RGB","XYZ","LAB")), + (("RGB","LUV"),("RGB","XYZ","LUV")), + (("RGB","LCH"),("RGB","XYZ","LAB","LCH")), + (("CMYK","Grayscale"),("CMYK","RGB","Grayscale")), + (("CMYK","RGB"),("CMYK","RGB")), + (("CMYK","HSB"),("CMYK","RGB","HSB")), + (("CMYK","XYZ"),("CMYK","RGB","XYZ")), + (("CMYK","LAB"),("CMYK","RGB","XYZ","LAB")), + (("CMYK","LUV"),("CMYK","RGB","XYZ","LUV")), + (("CMYK","LCH"),("CMYK","RGB","XYZ","LAB","LCH")), + (("HSB","Grayscale"),("HSB","RGB","Grayscale")), + (("HSB","RGB"),("HSB","RGB")), + (("HSB","CMYK"),("HSB","RGB","CMYK")), + (("HSB","XYZ"),("HSB","RGB","XYZ")), + (("HSB","LAB"),("HSB","RGB","XYZ","LAB")), + (("HSB","LUV"),("HSB","RGB","XYZ","LUV")), + (("HSB","LCH"),("HSB","RGB","XYZ","LAB","LCH")), + (("XYZ","Grayscale"),("XYZ","RGB","Grayscale")), + (("XYZ","RGB"),("XYZ","RGB")), + (("XYZ","CMYK"),("XYZ","RGB","CMYK")), + (("XYZ","HSB"),("XYZ","RGB","HSB")), + (("XYZ","LAB"),("XYZ","LAB")), + (("XYZ","LUV"),("XYZ","LUV")), + (("XYZ","LCH"),("XYZ","LAB","LCH")), + (("LAB","Grayscale"),("LAB","XYZ","RGB","Grayscale")), + (("LAB","RGB"),("LAB","XYZ","RGB")), + (("LAB","CMYK"),("LAB","XYZ","RGB","CMYK")), + (("LAB","HSB"),("LAB","XYZ","RGB","HSB")), + (("LAB","XYZ"),("LAB","XYZ")), + (("LAB","LUV"),("LAB","XYZ","LUV")), + (("LAB","LCH"),("LAB","LCH")), + (("LUV","Grayscale"),("LUV","XYZ","RGB","Grayscale")), + (("LUV","RGB"),("LUV","XYZ","RGB")), + (("LUV","CMYK"),("LUV","XYZ","RGB","CMYK")), + (("LUV","HSB"),("LUV","XYZ","RGB","HSB")), + (("LUV","XYZ"),("LUV","XYZ")), + (("LUV","LAB"),("LUV","XYZ","LAB")), + (("LUV","LCH"),("LUV","XYZ","LAB","LCH")), + (("LCH","Grayscale"),("LCH","LAB","XYZ","RGB","Grayscale")), + (("LCH","RGB"),("LCH","LAB","XYZ","RGB")), + (("LCH","CMYK"),("LCH","LAB","XYZ","RGB","CMYK")), + (("LCH","HSB"),("LCH","LAB","XYZ","RGB","HSB")), + (("LCH","XYZ"),("LCH","LAB","XYZ")), + (("LCH","LAB"),("LCH","LAB")), + (("LCH","LUV"),("LCH","LAB","XYZ","LUV")))) conversions = { 'Grayscale>RGB': grayscale_to_rgb, @@ -361,10 +403,11 @@ def convert(components, src, dst, preserve_alpha=True): if src == dst: return components - flows = _flow(src, dst) - for s, d in (flows[i:i + 2] for i in range(len(flows) - 1)): - if s == d: - continue + path = _paths.get((src, dst), None) + if path is None: + return None + + for s, d in zip(path[:-1], path[1:]): func = conversions.get('%s>%s' % (s, d)) if not func: return None