Deuteranomaly simulation with nip2-cli

colorblind.nip2

main = lut_deut_sim_xyz;

lut_deut_sim_srgb = colour_transform_to Image_type.sRGB lut_deut_sim_xyz;
lut_deut_sim_xyz = recomb deut_sim_xyz lut_neu_xyz;

lut_neu_xyz = colour_transform_to Image_type.XYZ lut_neu_srgb;
# https://github.com/obsproject/obs-studio/blob/master/plugins/obs-filters/data/LUTs/original.png
lut_neu_srgb = Image_file "lut-neutral.png";

deut_sim_xyz = lms_to_xyz * deut_sim_lms * xyz_to_lms;
deut_sim_lms = Matrix [[0.54, 0.46, 0],
                       [0.46, 0.54, 0],
                       [0,    0,    1.0]];

# https://en.wikipedia.org/wiki/LMS_color_space#CIECAM02
xyz_to_lms = Matrix [[ 0.7328, 0.4296, -0.1624],
                     [-0.7036, 1.6975,  0.0061],
                     [ 0.0030, 0.0136,  0.9834]];
lms_to_xyz = xyz_to_lms ** -1;

ImageToCubeLut.hs

{-# LANGUAGE FlexibleContexts #-}

module Main where

import Data.Foldable as F
import Data.List (intercalate)
import Data.Massiv.Array as A
import Data.Massiv.Array.IO
import Graphics.ColorSpace
import System.Environment (getArgs)
import System.IO

main :: IO ()
main = do
  args <- getArgs
  F.forM_ args $ \file -> do
    img <- readImage file :: IO (Image S RGB Word8)
    let imgFloat = (A.map . fmap) pxFloat img
    writeLut (file ++ ".cube") (imageToLut3D imgFloat)

pxFloat :: (Real a, Bounded a) => a -> Double
pxFloat n = realToFrac n / realToFrac (maxBound `asTypeOf` n)

writeLut :: Source r Ix3 (Pixel RGB Double)
         => FilePath
         -> Array r Ix3 (Pixel RGB Double)
         -> IO ()
writeLut file lut = withFile file WriteMode $ \h -> do
  hPutStrLn h "LUT_3D_SIZE 64"
  A.forM_ (lut3DToCube lut) $ \px -> do
    hPutStrLn h . intercalate " " . Prelude.map show . F.toList $ px

imageToLut3D :: Source r Ix2 e => Array r Ix2 e -> Array D Ix3 e
imageToLut3D = backpermute (Ix3 64 64 64) rgbCoord
  where
    rgbCoord (Ix3 r g b) = Ix2 y x
      where
        x = 64 * tileX + r
        y = 64 * tileY + g
        (tileY, tileX) = b `divMod` 8

lut3DToCube :: Source r Ix3 e => Array r Ix3 e -> Array D Ix1 e
lut3DToCube = backpermute (Ix1 (64 * 64 * 64)) flatIx
  where
    flatIx (Ix1 i) = Ix3 r g b
      where
        (bg, r) = i `divMod` 64
        (b, g) = bg `divMod` 64