-- Color calculations in Haskell -- Jan Behrens -- Version 0.3 (2008-06-29) -- -- This haskell program offers several functions to do calculations with -- color representations. Constants and equitations from literature are -- marked as such by a comment with a number in square brackets, referring -- to one of the following publications: -- -- [1] Color Appearance Models (Second Edition) by Mark D. Fairchild, 2005 -- ISBN 0-470-01216-1 / EAN 9-780470-01216-1 -- [2] A Standard Default Color Space for the Internet - sRGB -- Version 1.10, November 5, 1996 -- http://www.w3.org/Graphics/Color/sRGB -- [3] CIE 1931 Standard Colorimetric Observer xyz2(lambda) data -- http://www.cie.co.at/main/freepubs.html -- [4] WCS Data Archives from the World Color Survey -- Data from Berlin & Kay (1969 [1991]) & WCS Mapping Tables -- http://www.icsi.berkeley.edu/wcs/data.html -- - http://www.icsi.berkeley.edu/wcs/data/berlin-kay/BK-dict.txt -- - http://www.icsi.berkeley.edu/wcs/data/berlin-kay/BK-foci.txt -- - http://www.icsi.berkeley.edu/wcs/data/cnum-maps/cnum-vhcm-lab-new.txt -- -- -- I am allowing you to use my work under the following license conditions: -- -- Permission is hereby granted, free of charge, to any person obtaining a -- copy of this software and associated documentation files (the -- "Software"), to deal in the Software without restriction, including -- without limitation the rights to use, copy, modify, merge, publish, -- distribute, sublicense, and/or sell copies of the Software, and to -- permit persons to whom the Software is furnished to do so, subject to -- the following conditions: -- -- The above copyright notice and this permission notice shall be included -- in all copies or substantial portions of the Software. -- -- THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS -- OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF -- MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. -- IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY -- CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, -- TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE -- SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. -- {- This program has NO user interface, this means you have to use a haskell interpreter to calculate values by directly calling functions of this program. Here are some examples how to do this with "Hugs" or the "Glasgow Haskell Compiler" (one of them have to be installed on your system): $ hugs color.hs [...] Haskell 98 mode: Restart with command line option -98 to enable extensions Type :? for help Main> xyzTristimulusValues (Srgb (1,128/255,0)) -- calculate the XYZ values of an sRGB color (0.48959171484073,0.366983429681461,0.0450305716135768) Main> srgbValues (XyzTristimulus (0.6,0.5,0.4)) -- calculate the sRGB color of an XYZ-tristimulus (0.989580265800429,0.644610007668167,0.629571668699223) Main> toHexTriple (srgbValues (XyzTristimulus (0.6,0.5,0.4))) -- calculate the sRGB hex-triplet of an XYZ-tristimulus "#fca4a1" Main> lcwValues illuminantD65 (Srgb (0,1,0)) -- calculate the luminance, chroma and wavelength of sRGB green (0.7152,0.854549251536819,548.971196153777) Main> rectToPolar (cielabValues illuminantD65 (Srgb (0,1,0))) -- calculate the polar CIELab values of sRGB green (87.7370334735442,120.597837688152,136.391037447175) Main> srgbValues (RelativeToWhitepoint illuminantD65 (Cielab (polarToRect (50,30,10)))) -- calculate the sRGB values for a given polar CIELab color (0.665449434674922,0.38470490437161,0.436390789986915) Main> toHexTriple (srgbValues (RelativeToWhitepoint illuminantD65 (Cielab (polarToRect (80, 30, 90))))) -- convert the RAL DESIGN color 090-80-30 to a web color (sRGB hex-triplet) "#ddc48f" Main> toCiecam02Values monitorConditions (Srgb (1,0,0)) -- calculate the CIECAM02 color appearance values (lightness, chroma and hue) for sRGB red for typical monitor conditions (52.5614895146042,106.85923618289,32.2968736756162) Main> toHexTriple (srgbValues (fromCiecam02Values monitorConditions (50,0,0))) -- calculate the sRGB values for grey appearing to have 50% lightness "# Program error: argument out of range Main> :q [Leaving Hugs] $ # Hugs fails here due to some reason, but the Glasgow Haskell Compiler works better in this case: $ ghci color.hs [...] Loading package base ... linking ... done. [1 of 1] Compiling Main ( color.hs, interpreted ) Ok, modules loaded: Main. *Main> toHexTriple (srgbValues (fromCiecam02Values monitorConditions (50,0,0))) -- calculate the sRGB values for grey appearing to have 50% lightness "#868181" End of examples. -} ---------------------- -- libarary imports -- ---------------------- import Data.Map (Map) import qualified Data.Map as Map import Data.Maybe import Text.Printf (printf) -------------------------- -- physical color types -- -------------------------- data PhysicalColor = XyzTristimulus (Double, Double, Double) | LinearSrgb (Double, Double, Double) | Srgb (Double, Double, Double) | RelativeToWhitepoint PhysicalColor RelativeColor deriving (Show, Read) data RelativeColor = Lcw (Double, Double, Double) | -- Luminance, Chroma, Wavelength Cielab (Double, Double, Double) deriving (Show, Read) ------------------------- -- english color terms -- ------------------------- englishBlack = Cielab ( 15.60, -0.02, 0.02 ) -- J0 englishWhite = Cielab ((96.00 + 91.08) / 2, (- 0.06 - 0.05) / 2, 0.06 ) -- A0 B0 englishRed = Cielab ((41.22 + 30.77) / 2, ( 56.60 + 52.68) / 2, ( 40.99 + 34.06) / 2) -- G3 H3 englishGreen = Cielab ((51.57 + 41.22) / 2, (- 63.28 - 53.57) / 2, ( 28.95 + 23.28) / 2) -- F17 G17 englishYellow = Cielab ( 81.35, 7.28, 109.12 ) -- C9 englishBlue = Cielab ((51.57 + 41.22) / 2, (- 11.88 - 10.54) / 2, (- 38.56 - 39.65) / 2) -- (F27 G27) F28 G28 (F29 G29) englishBrown = Cielab ( 20.54, 13.28, 18.12 ) -- (I5) I6 (I7) englishPurple = Cielab ((30.77 + 20.54) / 2, ( 39.78 + 22.74) / 2, (- 27.99 - 16.65) / 2) -- H35 I35 englishPink = Cielab ((71.60 + 61.70) / 2, ( 40.04 + 40.07) / 2, (- 2.70 - 3.44) / 2) -- D38 E38 englishOrange = Cielab ( 51.57, 55.20, 68.32 ) -- F4 englishGrey = Cielab ( 51.57, -0.03, 0.04 ) -- F0 -- according to Berlin & Kay, 1969 [4] ----------------------------------- -- colors of certain wavelengths -- ----------------------------------- wavelengthColors :: Map Double PhysicalColor wavelengthColors = Map.fromList [ (380, XyzTristimulus (0.001368, 0.000039, 0.006450)), (385, XyzTristimulus (0.002236, 0.000064, 0.010550)), (390, XyzTristimulus (0.004243, 0.000120, 0.020050)), (395, XyzTristimulus (0.007650, 0.000217, 0.036210)), (400, XyzTristimulus (0.014310, 0.000396, 0.067850)), (405, XyzTristimulus (0.023190, 0.000640, 0.110200)), (410, XyzTristimulus (0.043510, 0.001210, 0.207400)), (415, XyzTristimulus (0.077630, 0.002180, 0.371300)), (420, XyzTristimulus (0.134380, 0.004000, 0.645600)), (425, XyzTristimulus (0.214770, 0.007300, 1.039050)), (430, XyzTristimulus (0.283900, 0.011600, 1.385600)), (435, XyzTristimulus (0.328500, 0.016840, 1.622960)), (440, XyzTristimulus (0.348280, 0.023000, 1.747060)), (445, XyzTristimulus (0.348060, 0.029800, 1.782600)), (450, XyzTristimulus (0.336200, 0.038000, 1.772110)), (455, XyzTristimulus (0.318700, 0.048000, 1.744100)), (460, XyzTristimulus (0.290800, 0.060000, 1.669200)), (465, XyzTristimulus (0.251100, 0.073900, 1.528100)), (470, XyzTristimulus (0.195360, 0.090980, 1.287640)), (475, XyzTristimulus (0.142100, 0.112600, 1.041900)), (480, XyzTristimulus (0.095640, 0.139020, 0.812950)), (485, XyzTristimulus (0.057950, 0.169300, 0.616200)), (490, XyzTristimulus (0.032010, 0.208020, 0.465180)), (495, XyzTristimulus (0.014700, 0.258600, 0.353300)), (500, XyzTristimulus (0.004900, 0.323000, 0.272000)), (505, XyzTristimulus (0.002400, 0.407300, 0.212300)), (510, XyzTristimulus (0.009300, 0.503000, 0.158200)), (515, XyzTristimulus (0.029100, 0.608200, 0.111700)), (520, XyzTristimulus (0.063270, 0.710000, 0.078250)), (525, XyzTristimulus (0.109600, 0.793200, 0.057250)), (530, XyzTristimulus (0.165500, 0.862000, 0.042160)), (535, XyzTristimulus (0.225750, 0.914850, 0.029840)), (540, XyzTristimulus (0.290400, 0.954000, 0.020300)), (545, XyzTristimulus (0.359700, 0.980300, 0.013400)), (550, XyzTristimulus (0.433450, 0.994950, 0.008750)), (555, XyzTristimulus (0.512050, 1.000000, 0.005750)), (560, XyzTristimulus (0.594500, 0.995000, 0.003900)), (565, XyzTristimulus (0.678400, 0.978600, 0.002750)), (570, XyzTristimulus (0.762100, 0.952000, 0.002100)), (575, XyzTristimulus (0.842500, 0.915400, 0.001800)), (580, XyzTristimulus (0.916300, 0.870000, 0.001650)), (585, XyzTristimulus (0.978600, 0.816300, 0.001400)), (590, XyzTristimulus (1.026300, 0.757000, 0.001100)), (595, XyzTristimulus (1.056700, 0.694900, 0.001000)), (600, XyzTristimulus (1.062200, 0.631000, 0.000800)), (605, XyzTristimulus (1.045600, 0.566800, 0.000600)), (610, XyzTristimulus (1.002600, 0.503000, 0.000340)), (615, XyzTristimulus (0.938400, 0.441200, 0.000240)), (620, XyzTristimulus (0.854450, 0.381000, 0.000190)), (625, XyzTristimulus (0.751400, 0.321000, 0.000100)), (630, XyzTristimulus (0.642400, 0.265000, 0.000050)), (635, XyzTristimulus (0.541900, 0.217000, 0.000030)), (640, XyzTristimulus (0.447900, 0.175000, 0.000020)), (645, XyzTristimulus (0.360800, 0.138200, 0.000010)), (650, XyzTristimulus (0.283500, 0.107000, 0.000000)), (655, XyzTristimulus (0.218700, 0.081600, 0.000000)), (660, XyzTristimulus (0.164900, 0.061000, 0.000000)), (665, XyzTristimulus (0.121200, 0.044580, 0.000000)), (670, XyzTristimulus (0.087400, 0.032000, 0.000000)), (675, XyzTristimulus (0.063600, 0.023200, 0.000000)), (680, XyzTristimulus (0.046770, 0.017000, 0.000000)), (685, XyzTristimulus (0.032900, 0.011920, 0.000000)), (690, XyzTristimulus (0.022700, 0.008210, 0.000000)), (695, XyzTristimulus (0.015840, 0.005723, 0.000000)), (700, XyzTristimulus (0.011359, 0.004102, 0.000000)), (705, XyzTristimulus (0.008111, 0.002929, 0.000000)), (710, XyzTristimulus (0.005790, 0.002091, 0.000000)), (715, XyzTristimulus (0.004109, 0.001484, 0.000000)), (720, XyzTristimulus (0.002899, 0.001047, 0.000000)), (725, XyzTristimulus (0.002049, 0.000740, 0.000000)), (730, XyzTristimulus (0.001440, 0.000520, 0.000000)), (735, XyzTristimulus (0.001000, 0.000361, 0.000000)), (740, XyzTristimulus (0.000690, 0.000249, 0.000000)), (745, XyzTristimulus (0.000476, 0.000172, 0.000000)), (750, XyzTristimulus (0.000332, 0.000120, 0.000000)), (755, XyzTristimulus (0.000235, 0.000085, 0.000000)), (760, XyzTristimulus (0.000166, 0.000060, 0.000000)), (765, XyzTristimulus (0.000117, 0.000042, 0.000000)), (770, XyzTristimulus (0.000083, 0.000030, 0.000000)), (775, XyzTristimulus (0.000059, 0.000021, 0.000000)), (780, XyzTristimulus (0.000042, 0.000015, 0.000000)) ] -- the standard observer with 2 degrees vision, as published in [3] -- also available with less precision as table 3.3 in [1] on page 74 -- wavelength range being "well" visible (and where data is good enough) minWellVisibleWavelength = 430 :: Double maxWellVisibleWavelength = 685 :: Double -- hugs complains, if the type 'Double' is not specified, why? -------------------------- -- standard illuminants -- -------------------------- illuminantA = XyzTristimulus (1.0985, 1, 0.3558) illuminantC = XyzTristimulus (0.9807, 1, 1.1823) illuminantD65 = XyzTristimulus (0.9595, 1, 1.0888) illuminantD50 = XyzTristimulus (0.9642, 1, 0.8249) illuminantF2 = XyzTristimulus (0.9920, 1, 0.6740) illuminantF8 = XyzTristimulus (0.9642, 1, 0.8246) illuminantF11 = XyzTristimulus (1.0096, 1, 0.6437) -- as shown in table 3.1 in [1] on page 62 illuminantE = XyzTristimulus (1, 1, 1) ------------------------------------------------------- -- functions for transforming tristimulus components -- ------------------------------------------------------- -- uses the given function to transform each component of a triple in the same way transformEachComponent :: (Double -> Double) -> (Double,Double,Double) -> (Double,Double,Double) transformEachComponent func (v1, v2, v3) = (func v1, func v2, func v3) -- transforms linear sRGB values to gamma-corrected values srgbGammaTransform = transformEachComponent ( \v -> if v <= 0.00304 then 12.92 * v else 1.055 * v ** (1.0 / 2.4) - 0.055 ) -- as in equitations 1.2a and 1.2b in [2] -- transforms gamma-corrected sRGB values to linear values srgbGammaTransformInverse = transformEachComponent ( \v -> if v <= 0.03928 then v / 12.92 else ((v + 0.055) / 1.055) ** 2.4 ) -- as in equitations 1.7a and 1.7b in [2] -- non-linear transformation used for normalized XYZ-values in CIELab cielabNonlinearTransform = transformEachComponent ( \v -> if v <= 0.008856 then 7.787 * v + (16/116) else v ** (1/3) ) -- as in equitation 10.4 in [1] on page 186 -- inverse of the non-linear transformation used for normalized XYZ-values in CIELab cielabNonlinearTransformInverse = transformEachComponent ( \v -> if v <= 0.008856 * 7.787 + (16/116) then (v - (16/116)) / 7.787 else v ** 3 ) -- inverse of equitation 10.4 in [1] on page 186 -- uses the given function to combine each component of two given triples to form a new triple combineEachComponent :: (Double -> Double -> Double) -> (Double,Double,Double) -> (Double,Double,Double) -> (Double,Double,Double) combineEachComponent func (v1, v2, v3) (w1, w2, w3) = (func v1 w1, func v2 w2, func v3 w3) -- linear conversion matrix for triples data LinearTransformation = LinearTransformation Double Double Double Double Double Double Double Double Double deriving (Show, Read) -- inverted matrix invertLT :: LinearTransformation -> LinearTransformation invertLT (LinearTransformation a11 a12 a13 a21 a22 a23 a31 a32 a33) = LinearTransformation (factor * (a22 * a33 - a23 * a32)) (factor * (a13 * a32 - a12 * a33)) (factor * (a12 * a23 - a13 * a22)) (factor * (a23 * a31 - a21 * a33)) (factor * (a11 * a33 - a13 * a31)) (factor * (a13 * a21 - a11 * a23)) (factor * (a21 * a32 - a22 * a31)) (factor * (a12 * a31 - a11 * a32)) (factor * (a11 * a22 - a12 * a21)) where factor = 1 / (a11 * a22 * a33 + a12 * a23 * a31 + a13 * a21 * a32 - a31 * a22 * a13 - a32 * a23 * a11 - a33 * a21 * a12 ) -- multiply a matrix with another matrix chainLT :: LinearTransformation -> LinearTransformation -> LinearTransformation chainLT (LinearTransformation a11 a12 a13 a21 a22 a23 a31 a32 a33) (LinearTransformation b11 b12 b13 b21 b22 b23 b31 b32 b33) = LinearTransformation (a11 * b11 + a12 * b21 + a13 * b31) (a11 * b12 + a12 * b22 + a13 * b32) (a11 * b13 + a12 * b23 + a13 * b33) (a21 * b11 + a22 * b21 + a23 * b31) (a21 * b12 + a22 * b22 + a23 * b32) (a21 * b13 + a22 * b23 + a23 * b33) (a31 * b11 + a32 * b21 + a33 * b31) (a31 * b12 + a32 * b22 + a33 * b32) (a31 * b13 + a32 * b23 + a33 * b33) -- multiply a matrix with a vector (triple) applyLT :: LinearTransformation -> (Double, Double, Double) -> (Double, Double, Double) applyLT (LinearTransformation a11 a12 a13 a21 a22 a23 a31 a32 a33) (v1, v2, v3) = (a11 * v1 + a12 * v2 + a13 * v3, a21 * v1 + a22 * v2 + a23 * v3, a31 * v1 + a32 * v2 + a33 * v3) -- transformation matrix to convert XYZ-values to linear sRGB-values xyzToSrgb = LinearTransformation ( 3.2410) (-1.5374) (-0.4986) (-0.9692) ( 1.8760) ( 0.0416) ( 0.0556) (-0.2040) ( 1.0570) -- as in equitation 1.1 in [2] -- transformation matrix to convert sRGB-values to XYZ-values srgbToXyz = LinearTransformation 0.4124 0.3576 0.1805 0.2126 0.7152 0.0722 0.0193 0.1192 0.9505 -- as in equitation 1.8 in [2] ciecam02Cat02 = LinearTransformation ( 0.7328) ( 0.4296) (-0.1624) (-0.7036) ( 1.6975) ( 0.0061) ( 0.0030) ( 0.0136) ( 0.9834) -- as in equitation 16.2 in [1] on page 268 ciecam02Cat02Inverse = invertLT ciecam02Cat02 ciecam02Hpe = LinearTransformation ( 0.38971) ( 0.68898) (-0.07868) (-0.22981) ( 1.18340) ( 0.04641) ( 0.0 ) ( 0.0 ) ( 1.0 ) -- as in equitation 16.13 in [1] on page 271 ciecam02Cat02ToHpe = ciecam02Hpe `chainLT` ciecam02Cat02Inverse ciecam02HpeToCat02 = invertLT ciecam02Cat02ToHpe ---------------------------------------------------------------------------- -- functions to force a physical color to be represented in a certain way -- ---------------------------------------------------------------------------- -- force a color to be represented as an XYZ-tristimulus xyzTristimulus :: PhysicalColor -> PhysicalColor xyzTristimulus color = XyzTristimulus (xyzTristimulusValues color) -- force a color to be represented as linear sRGB linearSrgb :: PhysicalColor -> PhysicalColor linearSrgb color = LinearSrgb (linearSrgbValues color) -- force a color to be represented as (gamma corrected) sRGB srgb :: PhysicalColor -> PhysicalColor srgb color = Srgb (srgbValues color) -- force a color to be represented by a given whitepoint and luminance, chroma and wavelength lcw :: PhysicalColor -> PhysicalColor -> PhysicalColor lcw white color = RelativeToWhitepoint white (Lcw $ lcwValues white $ color) -- force a color to be represented by a given whitepoint and CIELab values cielab :: PhysicalColor -> PhysicalColor -> PhysicalColor cielab white color = RelativeToWhitepoint white (Cielab $ cielabValues white $ color) --------------------------------------------------------- -- functions to extract values out of a physical color -- --------------------------------------------------------- -- returns the XYZ-values for a given color xyzTristimulusValues :: PhysicalColor -> (Double, Double, Double) xyzTristimulusValues (XyzTristimulus xyz) = xyz xyzTristimulusValues (LinearSrgb rgb) = applyLT srgbToXyz rgb xyzTristimulusValues (Srgb rgb') = xyzTristimulusValues $ linearSrgb $ Srgb rgb' xyzTristimulusValues (RelativeToWhitepoint white (Lcw (l, c, w))) = xyzTristimulusValues dimmedColor where pureColor = colorByWavelength w scaledColor = normalizeEnergy pureColor scaledWhite = normalizeEnergy white desaturatedColor = addRays (scaleEnergy c scaledColor) (scaleEnergy (1-c) scaledWhite) dimmedColor = scaleEnergy (l * linearLuminance white) desaturatedColor xyzTristimulusValues (RelativeToWhitepoint white (Cielab (l, a, b))) = xyz where x' = y' + a / 500 y' = (l + 16) / 116 z' = y' - b / 200 xyz = combineEachComponent (*) (cielabNonlinearTransformInverse (x', y', z')) (xyzTristimulusValues white) -- inverse of equitations 10.1 to 10.4 in [1] on page 186 -- returns the luminance for a given color (Y-value) -- (for the same spectral composition the luminance is linear to the amount of energy) linearLuminance :: PhysicalColor -> Double linearLuminance color = y where (x,y,z) = xyzTristimulusValues color -- Given a whitepoint, the following function calculates a hue of a physical color. -- The hue is NOT distortion free, but usable to be compared with other hues returned by this function. -- It is aligned in a way that all hues are positive and the wrap-around occurs at the purple colors (which have no wavelengths). -- TODO: give black and white a reasonable hue simpleHue :: PhysicalColor -> PhysicalColor -> Double simpleHue white color = if h' < 0 then h' + 2 * pi else h' where simpleHue' :: PhysicalColor -> PhysicalColor -> Double simpleHue' white (XyzTristimulus xyz) = atan2 (y' - z') (x' - y') where (x',y',z') = combineEachComponent (/) xyz (xyzTristimulusValues white) simpleHue' white color = simpleHue' white $ xyzTristimulus $ color redHue' = simpleHue' white $ colorByWavelength maxWellVisibleWavelength violetHue' = simpleHue' white $ colorByWavelength minWellVisibleWavelength offset = if violetHue' < redHue' then (violetHue' + redHue') / 2 else pi + (violetHue' + redHue') / 2 h' = align $ simpleHue' white color - offset align x = if x < 0 then align (x + 2 * pi) else if x >= 2 * pi then align (x - 2 * pi) else x -- Given a whitepoint, the following function calculates the chroma of a physical color. -- The chroma is NOT distortion free, but usable to be compared with other chroma values of colors having the same hue. simpleChroma :: PhysicalColor -> PhysicalColor -> Double simpleChroma white (XyzTristimulus xyz) = sqrt ((x'-y')^2 + (y'-z')^2) / y' where (x',y',z') = combineEachComponent (/) xyz (xyzTristimulusValues white) simpleChroma white color = simpleChroma white $ xyzTristimulus $ color -- limits for hues which can be mapped to a wavelength minWellVisibleHue white = simpleHue white $ colorByWavelength maxWellVisibleWavelength maxWellVisibleHue white = simpleHue white $ colorByWavelength minWellVisibleWavelength -- returns the linear sRGB-values for a given color linearSrgbValues :: PhysicalColor -> (Double, Double, Double) linearSrgbValues (XyzTristimulus xyz) = applyLT xyzToSrgb xyz linearSrgbValues (LinearSrgb rgb) = rgb linearSrgbValues (Srgb rgb') = srgbGammaTransformInverse rgb' linearSrgbValues color = linearSrgbValues $ xyzTristimulus $ color -- returns the (gamma-corrected) sRGB-values for a given color srgbValues :: PhysicalColor -> (Double, Double, Double) srgbValues (XyzTristimulus xyz) = srgbValues $ linearSrgb $ XyzTristimulus xyz srgbValues (LinearSrgb rgb) = srgbGammaTransform rgb srgbValues (Srgb rgb') = rgb' srgbValues color = srgbValues $ xyzTristimulus $ color -- returns the luminance, the chroma (relative to a monochromatic color) and a wavelength of a color using a given white-point -- purple colors are represented using negative chroma and the wavelength of the complement color lcwValues :: PhysicalColor -> PhysicalColor -> (Double, Double, Double) lcwValues white color = (l, c, w) where l = linearLuminance color / linearLuminance white c = (if complement then (-1) else 1) * simpleChroma white color / simpleChroma white pureColor (complement, w) = unpack $ wavelengthOfColor white $ color where unpack (Just wavelength) = (False, wavelength) unpack Nothing = (True, complementWavelength) Just complementWavelength = wavelengthOfColor white $ complementColor white $ color pureColor = colorByWavelength w -- returns the chroma relative to a monochromatic color using a given white-point linearChroma :: PhysicalColor -> PhysicalColor -> Double linearChroma white color = c where (l, c, w) = lcwValues white color -- returns the CIELab-values using a given white-point cielabValues :: PhysicalColor -> PhysicalColor -> (Double, Double, Double) cielabValues white (XyzTristimulus xyz) = (l, a, b) where (x', y', z') = cielabNonlinearTransform $ combineEachComponent (/) xyz (xyzTristimulusValues white) l = 116 * y' - 16 a = 500 * (x' - y') b = 200 * (y' - z') -- as in equitations 10.1 to 10.4 in [1] on page 186 cielabValues white color = cielabValues white $ xyzTristimulus $ color -------------------------------------- -- polar coordinate transformations -- -------------------------------------- -- transforms a triple using rectangular coordinates to a triple with polar coordinates rectToPolar :: (Double, Double, Double) -> (Double, Double, Double) rectToPolar (l, x, y) = (l, c, h) where c = sqrt (x^2 + y^2) h' = atan2 y x * 180 / pi h = if h' < 0 then h' + 360 else h' -- transforms a triple using polar coordinates to a triple using rectangular coordinates polarToRect :: (Double, Double, Double) -> (Double, Double, Double) polarToRect (l, c, h) = (l, x, y) where h' = h * pi / 180 x = c * cos h' y = c * sin h' ---------------------------------------- -- conversion from and to byte values -- ---------------------------------------- to8BitTriple :: (Double, Double, Double) -> (Int, Int, Int) to8BitTriple (v1, v2, v3) = (func v1, func v2, func v3) where func value = round (value * 255) toHexTriple :: (Double, Double, Double) -> String toHexTriple values = "#" ++ s1 ++ s2 ++ s3 where (v1, v2, v3) = to8BitTriple values (s1, s2, s3) = (format v1, format v2, format v3) format v = if v < 0 then "__" else if v > 255 then "^^" else printf "%02x" v from8BitTriple :: (Int, Int, Int) -> (Double, Double, Double) from8BitTriple (v1, v2, v3) = (func v1, func v2, func v3) where func value = fromIntegral value / 255 -- TODO: fromHexTriple ---------------------------------- -- scaling and mixing of colors -- ---------------------------------- scaleEnergy :: Double -> PhysicalColor -> PhysicalColor scaleEnergy factor (XyzTristimulus xyz) = XyzTristimulus $ transformEachComponent (factor *) $ xyz scaleEnergy factor (LinearSrgb rgb) = LinearSrgb $ transformEachComponent (factor *) $ rgb scaleEnergy factor (Srgb rgb') = scaleEnergy factor $ linearSrgb $ Srgb rgb' scaleEnergy factor color = scaleEnergy factor $ xyzTristimulus $ color normalizeEnergy :: PhysicalColor -> PhysicalColor normalizeEnergy color = scaleEnergy (1 / linearLuminance color) color addRays :: PhysicalColor -> PhysicalColor -> PhysicalColor addRays (XyzTristimulus xyz1) (XyzTristimulus xyz2) = XyzTristimulus $ combineEachComponent (+) xyz1 xyz2 addRays (LinearSrgb rgb1) (LinearSrgb rgb2) = LinearSrgb $ combineEachComponent (+) rgb1 rgb2 addRays (LinearSrgb rgb1) (Srgb rgb2') = addRays (LinearSrgb rgb1) (linearSrgb $ Srgb rgb2') addRays (Srgb rgb1') (LinearSrgb rgb2) = addRays (linearSrgb $ Srgb rgb1') (LinearSrgb rgb2) addRays (Srgb rgb1') (Srgb rgb2') = addRays (linearSrgb $ Srgb rgb1') (linearSrgb $ Srgb rgb2') addRays color1 color2 = addRays (xyzTristimulus color1) (xyzTristimulus color2) complementColor :: PhysicalColor -> PhysicalColor -> PhysicalColor complementColor white color = addRays white (scaleEnergy (-1) color) ------------------------------------------------- -- relationship between colors and wavelengths -- ------------------------------------------------- -- Map with wavelengths for a set of hues -- (made from the sane parts of the wavelengthColors Map) hueWavelengths :: PhysicalColor -> Map Double Double hueWavelengths white = Map.fromList $ map mappingFunc $ filter filterFunc $ Map.toList $ wavelengthColors where mappingFunc (wavelength, color) = (simpleHue white color, wavelength) filterFunc (wavelength, _) = wavelength >= minWellVisibleWavelength && wavelength <= maxWellVisibleWavelength -- returns a color for a given wavelength -- (with different lumiances for different wavelengths) colorByWavelength :: Double -> PhysicalColor colorByWavelength wavelength | wavelength < minWavelength || wavelength > maxWavelength = XyzTristimulus (0, 0, 0) | Map.member wavelength wavelengthColors = wavelengthColors Map.! wavelength | otherwise = interpolation where minWavelength = fst $ Map.findMin $ wavelengthColors maxWavelength = fst $ Map.findMax $ wavelengthColors (map1, map2) = Map.split wavelength wavelengthColors entry1 = Map.findMax map1 entry2 = Map.findMin map2 wavelength1 = fst entry1 wavelength2 = fst entry2 color1 = snd entry1 color2 = snd entry2 interpolation = addRays (scaleEnergy ((wavelength2 - wavelength) / (wavelength2 - wavelength1)) color1) (scaleEnergy ((wavelength - wavelength1) / (wavelength2 - wavelength1)) color2) -- uses a white-point to return Just the wavelength for a given color, or Nothing, if the color is purple wavelengthOfColor :: PhysicalColor -> PhysicalColor -> Maybe Double wavelengthOfColor white color | hasWavelength = Just wavelength | otherwise = Nothing where hue = simpleHue white color hasWavelength = hue >= minWellVisibleHue white && hue <= maxWellVisibleHue white hueMap = hueWavelengths white wavelength | Map.member hue hueMap = hueMap Map.! hue | otherwise = interpolation (map1, map2) = Map.split hue hueMap entry1 = Map.findMax map1 entry2 = Map.findMin map2 hue1 = fst entry1 hue2 = fst entry2 wavelength1 = snd entry1 wavelength2 = snd entry2 interpolation = (wavelength1 * (hue2 - hue) / (hue2 - hue1)) + (wavelength2 * (hue - hue1) / (hue2 - hue1)) -- wavelength of color according to CIECAM02 and given viewing conditions (compared with 99% (physically) saturated almost monochromatic light) -- TODO: This function is REALLY slow, because of adjusting the perceived chroma in small fixed with steps, instead of doing a binary search or something similar. perceivedWavelengthOfColor :: Ciecam02Condition -> PhysicalColor -> Maybe Double perceivedWavelengthOfColor viewingConditions color = wavelengthCalcFromChroma chroma where white = ciecam02White viewingConditions (j, chroma, h) = toCiecam02Values viewingConditions color wavelengthCalcFromChroma chroma' = if c >= 0.99 then maybeWl else if c >= 0 then wavelengthCalcFromChroma (chroma' + 0.1) else Nothing where c = linearChroma white $ fromCiecam02Values viewingConditions (j, chroma', h) maybeWl = wavelengthOfColor white $ fromCiecam02Values viewingConditions (j, chroma', h) ----------------------------------------- -- limitation of chroma for sRGB gamut -- ----------------------------------------- -- eventually reduces chroma to make the color fit in sRGB -- uses illuminantD65 limitToSrgb :: PhysicalColor -> PhysicalColor limitToSrgb color = LinearSrgb $ if l < 0 then(0, 0, 0) else if l > 1 then (1, 1, 1) else (s * r' + (1-s) * l, s * g' + (1-s) * l, s * b' + (1-s) * l) where l = linearLuminance color (r', g', b') = linearSrgbValues color limitMin v = if v < 0 then l / (l - v) else 1 limitMax v = if v > 1 then (1 - l) / (v - l) else 1 limit v = min (limitMin v) (limitMax v) s = minimum [limit r', limit g', limit b'] ------------------------------- -- special lists of interest -- ------------------------------- allSaturatedSrgbColors = map (Srgb . from8BitTriple) $ [ (255,0,255-n) | n <- [0..255] ] ++ [ (255,n,0) | n <- [1..255] ] ++ [ (255-n,255,0) | n <- [1..255] ] ++ [ (0,255,n) | n <- [1..255] ] ++ [ (0,255-n,255) | n <- [1..255] ] ++ [ (n,0,255) | n <- [1..255] ] allSaturatedSrgbColorsWithWavelength :: [(PhysicalColor, Maybe Double, Maybe Double)] allSaturatedSrgbColorsWithWavelength = [ (color, maybeWl color, maybeWl' color) | color <- allSaturatedSrgbColors, isJust (maybeWl color) || isJust (maybeWl' color) ] where maybeWl color = wavelengthOfColor illuminantD65 color maybeWl' color = perceivedWavelengthOfColor monitorConditions color allSaturatedSrgbColorsWithWavelengthInHtml = concat $ map formatter $ allSaturatedSrgbColorsWithWavelength where -- formatter :: (PhysicalColor, Double) -> String formatter (color, maybeWl, maybeWl') = "" ++ toHexTriple (srgbValues color) ++ "" ++ (if isJust maybeWl' then printf "%3.2f nm" (fromJust maybeWl') else "-") ++ "(" ++ (if isJust maybeWl then printf "%3.2f nm" (fromJust maybeWl) else "-") ++ ")\n" ---------------------- -- color appearance -- ---------------------- -- data type for viewing conditions data Ciecam02Condition = Ciecam02Condition { ciecam02White :: PhysicalColor, ciecam02RelativeBackgroundLuminance :: Double, ciecam02AdaptingLuminance :: Double, ciecam02Ambient :: Double, -- 0.0 = dark, 0.5 = dim, 1.0 = average ciecam02DiscountIlluminant :: Bool } deriving (Show, Read) -- predefined viewing conditions normalConditions = Ciecam02Condition illuminantD65 0.2 (0.2 * 1000 / pi) 1.0 True monitorConditions = Ciecam02Condition illuminantD65 0.2 (0.2 * 80) 0.5 False -- property derived from a viewing condition ciecam02ExponentialNonlinearityC :: Ciecam02Condition -> Double ciecam02ExponentialNonlinearityC viewingConditions | a >= 0.0 && a <= 0.5 = 0.525 + 2 * (0.59 - 0.525) * a | a >= 0.5 && a <= 1.0 = 0.59 + 2 * (0.69 - 0.59) * (a - 0.5) | otherwise = undefined where a = ciecam02Ambient viewingConditions -- interpolation from values in table 16.1 in [1] on page 267 -- property derived from a viewing condition ciecam02ChromaticInductionFactorNc :: Ciecam02Condition -> Double ciecam02ChromaticInductionFactorNc viewingConditions | a >= 0.0 && a <= 1.0 = 0.8 + 0.2 * a | otherwise = undefined where a = ciecam02Ambient viewingConditions -- interpolation from values in table 16.1 in [1] on page 267 -- property derived from a viewing condition ciecam02MaximumDegreeOfAdaptionF :: Ciecam02Condition -> Double ciecam02MaximumDegreeOfAdaptionF = ciecam02ChromaticInductionFactorNc -- property derived from a viewing condition ciecam02DegreeOfAdaptionD :: Ciecam02Condition -> Double ciecam02DegreeOfAdaptionD viewingConditions | discount = 1.0 | otherwise = f * (1 - (1 / 3.6 * exp (-(la + 42) / 92))) -- equitation 16.3 in [1] on page 268 where discount = ciecam02DiscountIlluminant viewingConditions f = ciecam02MaximumDegreeOfAdaptionF viewingConditions la = ciecam02AdaptingLuminance viewingConditions -- property derived from a viewing condition ciecam02LuminanceLevelAdaptionFactorFl :: Ciecam02Condition -> Double ciecam02LuminanceLevelAdaptionFactorFl viewingConditions = 0.2 * k^4 * 5 * la + 0.1 * (1 - k^4)^2 * (5 * la)**(1/3) -- equitation 16.8 in [1] on page 270 where la = ciecam02AdaptingLuminance viewingConditions k = 1 / (5 * la + 1) -- equitation 16.7 in [1] on page 270 -- property derived from a viewing condition ciecam02InductionFactorsNbbNcb :: Ciecam02Condition -> Double ciecam02InductionFactorsNbbNcb viewingConditions = 0.725 * (1 / n)**0.2 -- equitation 16.10 in [1] on page 270 where n = ciecam02RelativeBackgroundLuminance viewingConditions -- property derived from a viewing condition ciecam02BaseExponentialNonlinearityZ :: Ciecam02Condition -> Double ciecam02BaseExponentialNonlinearityZ viewingConditions = 1.48 + sqrt n -- equitation 16.11 in [1] on page 270 where n = ciecam02RelativeBackgroundLuminance viewingConditions -- converts a physical color to relative CIECAM02 color appearance values (lightness, chroma and hue) for given viewing conditions toCiecam02Values :: Ciecam02Condition -> PhysicalColor -> (Double, Double, Double) toCiecam02Values viewingConditions color = (j, chroma, if h < 0 then h / pi * 180 + 360 else h / pi * 180) where white' = ciecam02White viewingConditions -- unnormalized white white = normalizeEnergy white' -- normalized white d = ciecam02DegreeOfAdaptionD viewingConditions fl = ciecam02LuminanceLevelAdaptionFactorFl viewingConditions nbb = ciecam02InductionFactorsNbbNcb viewingConditions ncb = ciecam02InductionFactorsNbbNcb viewingConditions ec = ciecam02ExponentialNonlinearityC viewingConditions z = ciecam02BaseExponentialNonlinearityZ viewingConditions nc = ciecam02ChromaticInductionFactorNc viewingConditions n = ciecam02RelativeBackgroundLuminance viewingConditions nColor = scaleEnergy (1 / linearLuminance white') color -- normalized color rgb = applyLT ciecam02Cat02 (xyzTristimulusValues nColor) -- equitation 16.1 in [1] on page 268 rgbw = applyLT ciecam02Cat02 (xyzTristimulusValues white) -- according to equitation 16.1 in [1] on page 268 rgbc = combineEachComponent (\v1 v2 -> v1 * (d / v2 + 1 - d)) rgb rgbw -- equitations 16.4 to 16.6 in [1] on page 268 without factor 100 rgb' = applyLT ciecam02Cat02ToHpe rgbc -- equitation 16.12 in [1] on page 271 rgb'w = applyLT ciecam02Hpe (xyzTristimulusValues white) -- according to equitations 16.1 and 16.12 in [1] on pages 268 and 271 rgb'a = transformEachComponent (\v -> let t = (fl * v)**0.42 in 400 * t / (27.13 + t) + 0.1) rgb' -- equitations 16.15 to 16.17 in [1] on page 271, without factor 100 rgb'aw = transformEachComponent (\v -> let t = (fl * v)**0.42 in 400 * t / (27.13 + t) + 0.1) rgb'w -- according to equitations 16.15 to 16.17 in [1] on page 271, without factor 100 (r'a, g'a, b'a ) = rgb'a (r'aw, g'aw, b'aw) = rgb'aw a = r'a - (12/11) * g'a + (1/11) * b'a -- equitation 16.18 in [1] on page 271 b = (r'a + g'a - 2 * b'a) / 9 -- equitation 16.19 in [1] on page 271 h = atan2 b a -- equitation 16.20 in [1] on page 272 et = (cos (h + 2) + 3.8) / 4 -- equitation 16.21 in [1] on page 272, using radians instead of degrees aj = (2 * r'a + g'a + b'a / 20 - 0.305) * nbb -- equitation 16.23 in [1] on page 272 aw = (2 * r'aw + g'aw + b'aw / 20 - 0.305) * nbb -- according to equitation 16.23 in [1] on page 272 j = 100 * (aj / aw) ** (ec * z) -- equitation 16.24 in [1] on page 272 t = (50000/13) * nc * ncb * et * sqrt (a^2 + b^2) / (r'a + g'a + (21/20) * b'a) -- equitation 16.26 in [1] on page 273 chroma = t**0.9 * sqrt (j / 100) * (1.64 - 0.29**n)**0.73 -- equitation 16.27 in [1] on page 273 -- converts relative CIECAM02 color appearance values (lightness, chroma and hue) at given viewing conditions to a physical color fromCiecam02Values:: Ciecam02Condition -> (Double, Double, Double) -> PhysicalColor fromCiecam02Values viewingConditions (j, chroma, h') = scaleEnergy (linearLuminance white') (XyzTristimulus xyz) where white' = ciecam02White viewingConditions -- unnormalized white white = normalizeEnergy white' -- normalized white d = ciecam02DegreeOfAdaptionD viewingConditions fl = ciecam02LuminanceLevelAdaptionFactorFl viewingConditions nbb = ciecam02InductionFactorsNbbNcb viewingConditions ncb = ciecam02InductionFactorsNbbNcb viewingConditions ec = ciecam02ExponentialNonlinearityC viewingConditions z = ciecam02BaseExponentialNonlinearityZ viewingConditions nc = ciecam02ChromaticInductionFactorNc viewingConditions n = ciecam02RelativeBackgroundLuminance viewingConditions rgbw = applyLT ciecam02Cat02 (xyzTristimulusValues white) rgb'w = applyLT ciecam02Hpe (xyzTristimulusValues white) rgb'aw = transformEachComponent (\v -> let t = (fl * v)**0.42 in 400 * t / (27.13 + t) + 0.1) rgb'w -- according to equitations 16.15 to 16.17 in [1] on page 271, without factor 100 (r'aw, g'aw, b'aw) = rgb'aw aw = (2 * r'aw + g'aw + b'aw / 20 - 0.305) * nbb -- according to equitation 16.23 in [1] on page 272 t = (chroma / sqrt (j / 100) / (1.64 - 0.29**n)**0.73) ** (1/0.9) -- inverse of equitation 16.27 in [1] on page 273 h = h' / 180 * pi et = (cos (h + 2) + 3.8) / 4 -- equitation 16.21 in [1] on page 272, using radians instead of degrees aj = aw * (j / 100) ** (1 / (ec*z)) -- inverse of equitation 16.24 in [1] on page 272 ct = (50000/13) * nc * ncb * et / t -- help variable ca = cos h / ct -- help variable cb = sin h / ct -- help variable -- Following equitations are true: -- aj / nbb + 0.305 = 2 * r'a + g'a + 1/20 * b'a -- a/ca = r'a + g'a + 21/20 * b'a -- b/cb = r'a + g'a + 21/20 * b'a -- a = r'a - 12/11 * g'a + 1/11 * b'a -- b = 1/9 * r'a + 1/9 * g'a - 2/9 * b'a -- They can be transformed to: -- aj / nbb + 0.305 = 2 * r'a + g'a + 1/20 * b'a -- 0 = (ca - 1) * r'a + (ca + 12/11) * g'a + (21/20 * ca - 1/11) * b'a -- 0 = (cb - 1/9) * r'a + (cb - 1/9) * g'a + (21/20 * cb + 2/9) * b'a -- The system of linear equitations is solved by matrix inversion: matrix = LinearTransformation ( 2 ) ( 1 ) ( 1/20 ) (ca - 1 ) (ca + 12/11) (21/20 * ca - 1/11) (cb - 1/9) (cb - 1/9 ) (21/20 * cb + 2/9 ) vector = (aj / nbb + 0.305, 0, 0) rgb'a = applyLT (invertLT matrix) vector rgb' = transformEachComponent (\v -> ( (v-0.1) * 27.13 / (400 - (v-0.1)) )**(1/0.42) / fl) rgb'a -- inversion of equitations 16.15 to 16.17 in [1] on page 271, without factor 100 rgbc = applyLT ciecam02HpeToCat02 rgb' -- inversion of equitation 16.12 in [1] on page 271 rgb = combineEachComponent (\v1 v2 -> v1 / (d / v2 + 1 - d)) rgbc rgbw -- inversion of equitations 16.4 to 16.6 in [1] on page 268, without factor 100 xyz = applyLT ciecam02Cat02Inverse rgb -- inversion of equitation 16.1 in [1] on page 268 -- converts a CIECAM02 hue from 0-360 degrees to a value between 0 and 400, -- where 0 = red, 100 = yellow, 200 = green and 300 = blue ciecam02HueQuadrature :: Double -> Double ciecam02HueQuadrature h | h >= 0.00 && h < 20.14 = 400 - 100 * ( 20.14 - h) / (380.14 - 237.53) | h >= 20.14 && h < 90.00 = 0 + 100 * (h - 20.14) / ( 90.00 - 20.14) | h >= 90.00 && h < 164.25 = 100 + 100 * (h - 90.00) / (164.25 - 90.00) | h >= 164.25 && h < 237.53 = 200 + 100 * (h - 164.25) / (237.53 - 164.25) | h >= 237.53 && h < 360.00 = 300 + 100 * (h - 237.53) / (380.14 - 237.53) | otherwise = error "Hue out of range." -- TODO: corner cases -- linear interpolation of table 16.2 in [1] on page 272 -- inverse of ciecam02HueQuadrature ciecam02HueQuadratureInverse :: Double -> Double ciecam02HueQuadratureInverse h' = if h >= 360 then h - 360 else h where h | h' >= 0 && h' < 100 = 20.14 + (h' - 0) / 100 * ( 90.00 - 20.14) | h' >= 100 && h' < 200 = 90.00 + (h' - 100) / 100 * (164.25 - 90.00) | h' >= 200 && h' < 300 = 164.25 + (h' - 200) / 100 * (237.53 - 164.25) | h' >= 300 && h' < 400 = 237.53 + (h' - 300) / 100 * (380.14 - 237.53) | otherwise = error "Hue out of range." -- TODO: corner cases -- converts the hue of CIECAM02 color appearance values according to the ciecam02HueQuadrature function quadHueOfCiecam02Values :: (Double, Double, Double) -> (Double, Double, Double) quadHueOfCiecam02Values (j, c, h) = (j, c, ciecam02HueQuadrature h) -- inverse of quadHueOfCiecam02Values unquadHueOfCiecam02Values :: (Double, Double, Double) -> (Double, Double, Double) unquadHueOfCiecam02Values (j, c, h') = (j, c, ciecam02HueQuadratureInverse h')