{-
Explore CIELab colors.

Copyright 2017 Ken Takusagawa

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or (at
your option) any later version.

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
-}

{-# LANGUAGE LambdaCase, ScopedTypeVariables #-}
module Main where {

import Data.Word(Word8);
import qualified Data.Colour.CIE.Illuminant;
import Data.Colour.CIE(Chromaticity,lightness,cieLABView,cieLAB);
import qualified Data.Colour.SRGB.Linear as Linear;
import Data.Colour.SRGB(sRGB,channelRed,channelGreen,channelBlue,Colour,sRGB24,RGB,toSRGB24,toSRGB);
import Data.Tuple.Curry(uncurryN);
import Control.Monad(join,msum,MonadPlus,mzero,guard);
import System.Environment(getArgs);
import Data.List(genericTake,genericIndex,genericLength,genericDrop,group,sortBy,intersperse,genericReplicate);
import Control.Exception(assert);
import Data.Maybe(fromJust);
import Codec.Picture(writePng,encodePng);
import Codec.Picture.Types(PixelRGBF(..),MutableImage,createMutableImage,Pixel,writePixel,PixelRGB8(..),Image,freezeImage,pixelMap);
import Control.Monad.ST(ST,runST);
import Control.Monad.Primitive(PrimState,PrimMonad);
import System.Random(randomRIO,randomIO,setStdGen,mkStdGen);
import Debug.Trace(trace);
import qualified Data.ByteString.Lazy as ByteString;
import qualified Data.Colour.Names as Color;
import Control.Arrow((***));

main :: IO();
main = getArgs >>= \case {
["band"] -> outppm_band rgbcolors;
-- uses the hardcoded lightness
["palette",n] -> outppm_band
$ map (get_color . nphi)
$ take (read n)
$ enumFrom 0;
["plot"] -> mapM_ putStrLn $ map flow rgbcolors; -- gets the shape of the butterfly for plotting in octave
["test-juicy"] -> bar;
["butterfly"] -> many_pixel; -- most colorful butterfly

-- the edge of the gamut, by angle on a given lightness plane
["lab",width,l] -> outppm_band
$ map (ab_find $ Lightness $ read l)
$ range01 $ read width;

-- perceptually uniform
["pband",l] -> outppm_band
$ all_edge $ Lightness $ read l;
-- at a fixed (hardcoded) lightness
["pbandf"] -> outppm_band $ all_colors;
["manylab",width,height] -> many_lab (read width) (read height);
["edgelength",l] -> print $ length $ all_edge $ Lightness $ read l;
["goline-scaled",n] -> mapM_ print $ scaled_pband (line_to_light $ read n) numColumns;

-- calculates the points on the perimeter, at perceptually equidistant fenceposts, in AB space
-- join (***) is an idiom for mapTuple;
["perim",n] -> mapM_ print $ map (join (***) (realToFrac::Double->Float)) $ all_ab $ line_to_light $ read n;

["slice",l] -> slice $ Lightness $ read l;
["lof",r,g,b] -> print $ lightness_of (read r)(read g)(read b); -- 0..1
["rslice",l] -> radial_slice $ Lightness $ read l;
["islice",n] -> radial_slice $ genericIndex lightness_index ((read n)::Integer);
["prainbow"] -> outppm_band all_rainbow;
["urainbow"] -> outppm_band uncorrected_rainbow;
["urainbown",n] -> outppm_band $ uncorrected_rainbow_of_width $ read n;
["read-goline"] -> all_read_goline;
["frame", w, h] -> out_rainbow_frame (read w) (read h);
["circle",shift] -> out_rainbow_circle $ read shift; -- 0..6
["read-perim"] -> all_read_perim;
["rainbow-shift",width,shift] -> rainbow_shift_plot (read width) (read shift);
["prainbow-disc", shift] -> out_perceptually_uniform_rainbow_circle $ read shift;
_ -> undefined;
};

-- this is to avoid warnings that trace is unused
dummy :: ();
dummy = trace "" ();

three_range :: [(Word8, Word8, Word8)];
three_range = let {
c :: Word8;
c = maxBound;
a :: [Word8];
a = enumFromTo 0 c;
b :: [Word8];
b = reverse a;
} in tail -- remove the circular duplicate
$ map head $ group -- remove duplicates
                                -- ends at 00c
$  do {x :: Word8 <- a;return (0,x,c);} -- 0cc
++ do {x :: Word8 <- b;return (0,c,x);} -- 0c0
++ do {x :: Word8 <- a;return (x,c,0);} -- cc0
++ do {x :: Word8 <- b;return (c,x,0);} -- c00
++ do {x :: Word8 <- a;return (c,0,x);} -- c0c
++ do {x :: Word8 <- b;return (x,0,c);} -- 00c
;

continuous3 :: Double -> (Double,Double,Double);
continuous3 x = let {
pieces :: [Double -> (Double, Double, Double)];
pieces =
[\t -> (0  ,t  ,1)
,\t -> (0  ,1  ,1-t)
,\t -> (t  ,1  ,0)
,\t -> (1  ,1-t,0)
,\t -> (1  ,0  ,t)
,\t -> (1-t,0  ,1)
];
npieces :: [(Double, Double -> (Double, Double, Double))];
npieces = zip (map fromInteger [(0::Integer)..]) pieces;
} in assert (0 <= x) $ if (6 < x) then continuous3 $ x-6 else
fromJust $ msum (map (eval_piece x) npieces);

eval_piece :: (MonadPlus m) => Double -> (Double,Double -> b) -> m b;
eval_piece x (z,f) = if x <= z+1 then
return $ f $ x - z else mzero;

type C = Colour Double;

tupleToColor :: (Word8,Word8,Word8) -> C;
tupleToColor = uncurryN sRGB24;

rgbcolors :: [C];
rgbcolors = map tupleToColor three_range;

white_point :: Chromaticity Double;
white_point=Data.Colour.CIE.Illuminant.d65;

-- the darkest color is solid blue
-- the lightest color is solid yellow
darkest :: Lightness;
darkest = Lightness $ minimum $ map (lightness white_point) rgbcolors;

-- see note below about out-of-gamut RGB colors
scale_color1 :: C -> C;
scale_color1 c = let {
(l,a,b) = cieLABView white_point c;
} in assert (darkest <= Lightness l) $ cieLAB_l darkest a b;

in01 :: Double -> Bool;
in01 x = 0<=x && x<=1;

scale_color_M :: forall m . (MonadPlus m) => C -> m C;
scale_color_M c1 = let {
c :: C;
c = scale_color1 c1;
} in do {
mapM_ (guard . in01) $ rgblist $ Linear.toRGB c;
return c;
};

-- a failure
step1fail :: [C];
step1fail = map scale_color1 rgbcolors;

circular_differences :: [a] -> [(a,a)];
circular_differences l = zip l (tail l ++ l);

color_difference :: C -> C -> Double;
color_difference x y = distance2 (list3 $ cieLABView white_point x) (list3 $ cieLABView white_point y);

list3 :: (a,a,a) -> [a];
list3 (x,y,z) = [x,y,z];

mapTuple3 :: (a -> b) -> (a,a,a) -> (b,b,b);
mapTuple3 f (x,y,z) = (f x, f y, f z);

distance2 :: (Num a) => [a] -> [a] -> a;
distance2 p q = sum $ map (\x -> x*x) $ zipWith (-) p q;

sortWith :: (Ord b) => (a -> b) -> [a] -> [a];
sortWith f = sortBy (\x y -> compare (f x) (f y));

-- fastest speed 0,72,196
-- slowest speed 7,0,255
speeds :: [C] -> [(Double,RGB Word8)];
speeds cs = do {
(x,y) <- circular_differences cs;
return (color_difference x y,toSRGB24 x);
};

fastest :: Double;
fastest = maximum $ map (uncurry color_difference) $circular_differences rgbcolors;
-- step1;

-- should be a small number
check_continuous :: Double;
check_continuous = let {
f :: Double -> (Word8,Word8,Word8) -> Double;
f n www = distance2 (map (* 255) $ list3 $ continuous3 $ get_t n)
(map fromIntegral $ list3 www);
get_t :: Double -> Double;
get_t n = 6*n/(fromIntegral $ length three_range);
} in maximum $ zipWith f [1..] three_range; -- start from 1 because tail in three_range

epsilon :: Double;
epsilon = 1e-11;

root_bisect :: (Double -> Double) -> Double -> Double -> Double;
root_bisect f left right =
let {mid = (left+right)/2;
z = f mid;
} in -- trace (show mid) $
if (abs (left - right))<epsilon then right else -- prefer too large a step
if z<0 then root_bisect f mid right
else if 0<z then root_bisect f left mid
else mid;

continuous3rgb :: Double -> C;
continuous3rgb t = let {
(x,y,z) = continuous3 t;
} in sRGB x y z;

desired_difference :: Double;
-- fastest is about 0.585; (for 1st try failure)
-- in the saturated rainbow, fastest is about 0.656;
-- desired_difference = fastest;
desired_difference = 0.1;

-- a failure
continuous_from_failure :: Double -> Double -> Double ;
continuous_from_failure x y = let {
f :: Double -> C;
f = scale_color1 . continuous3rgb }
in step_the_desired_distance (f x) (f y);


tstep :: Double;
tstep = 1/255; -- magic, = length three_range / 6

root_find_from :: (Double -> Double) -> Double -> Double;
root_find_from f left = let { right = left + tstep*0.1 } in
if (f right>0) then root_bisect f left right
else root_find_from f right;

-- next_t :: Double -> Double;
-- next_t = next_x continuous_from;

next_x :: (Double -> Double -> Double) -> Double -> Double;
next_x distance t = root_find_from (distance t) t;

-- all_t :: [Double];
-- all_t = takeWhile (<6) $ iterate next_t 0;

-- distance from the last color to the first color, offset fastest
-- edge_effect_error :: Double;
-- edge_effect_error = (continuous_from 0 $ last all_t)/desired_difference;

all_colors :: [C];
-- This fails terribly due to weirdness of Lab colors.
-- Taking a color, bringing down its lightness, results in
-- out-of-gamut RGB colors all the time, at the lightness of blue (darkest).
-- I do not fully understand why;
-- all_colors = map (scale_color . continuous3rgb) all_t;
all_colors = all_edge
 $ line_to_light 546; -- discovered to have the longest perimeter
-- (Lightness 50) ;
-- $ Lightness darkest ;

unlist3 :: [a] -> (a,a,a);
unlist3 [x,y,z] = (x,y,z);
unlist3 _ = undefined;

quantize :: C -> C;
quantize = (uncurryN sRGB24). unlist3. rgblist. toSRGB24;

rgblist :: RGB a -> [a];
rgblist c = [channelRed c,channelGreen c,channelBlue c];

outcolor :: C -> String;
outcolor = concat . intersperse " " . map show . rgblist . toSRGB24;

phi :: Double;
phi = 0.5*(1+sqrt(5));

frac :: Double -> Double;
frac x = x - fromInteger( floor x);

nphi :: Integer -> Double;
nphi n = frac $ (fromIntegral n)*phi;

get_color_index :: Double -> Integer;
get_color_index x = if x==1 then pred $ genericLength all_colors
else (floor $ x*(fromIntegral $ length all_colors));

get_color :: Double -> C;
get_color x = genericIndex all_colors $ get_color_index x;

flow :: C -> String;
flow c = let {
(_,a,b) = cieLABView white_point c;
} in (show a)++" "++(show b);

foo :: forall s . ST s (Image PixelRGBF);
foo = do {
pic :: MutableImage s PixelRGBF <- createMutableImage 100 100 $ PixelRGBF 0 0 0;
writePixel pic 50 50 $ PixelRGBF 0 0.8 0;
freezeImage pic;
};

bar2 :: IO();
bar2 = writePng "/mit/bitbucket/www/ss1.png" $ pixelMap quantize_pixel $ runST foo;

-- origin is top left
bar :: IO();
bar = do {
pic :: MutableImage (PrimState IO) PixelRGBF <- createMutableImage 100 100 $ PixelRGBF 0 0 0;
writePixel pic 20 60 $ PixelRGBF 0 0.8 0;
freezeImage pic >>= writePng "/mit/bitbucket/www/ss1.png" . pixelMap quantize_pixel;
};

addpoint :: IO ((Double,Double),PixelRGBF);
addpoint = do {
c :: C <- randomRIO (0,6) >>= return . continuous3rgb;
let {(_,a,b) = cieLABView white_point c};
t :: Double <- randomIO >>= return . sqrt; --sqrt pushes toward edge
let {v :: [Double];v = map (realToFrac .(*t)) [a,b]};
return ((v!!0,v!!1),fpixel c);
};

image_width :: Integer;
image_width = 4800;

image_height :: Integer;
image_height = image_width;

-- usually 1000000
num_rays :: Integer;
num_rays = 20000;

offset :: Double;
offset = 128;

scale_location :: (Double,Double) -> (Double,Double);
scale_location (x,y) = ((x+offset)/(2*offset)*fromIntegral image_width
,(1-(y+offset)/(2*offset))*fromIntegral image_height
);

one_add_pixel :: MutableImage (PrimState IO) PixelRGBF -> IO();
one_add_pixel pic = do {
(p,c) <- addpoint;
let {(x,y) = scale_location p;};
writePixel pic (round x) (round y) c;
};

many_pixel :: IO ();
many_pixel = do {
setStdGen $ mkStdGen 1;
_ :: Double <- randomIO;  -- avoid the first sample;
pic :: MutableImage (PrimState IO) PixelRGBF <- standard_image $ PixelRGBF 0 0 0;
-- replicateM_ 10000000 $ one_add_pixel pic;
genericReplicateM_ num_rays $ one_rainbow_ray pic;
outpic pic;
-- freezeImage pic >>= writePng "/mit/kenta/0/ss2.png" . pixelMap quantize_pixel;
};

fpixel :: C -> PixelRGBF;
fpixel = uncurryN PixelRGBF . unlist3 . map realToFrac . rgblist . toSRGB;
-- supposedly realToFrac should be optimized via RULES to use GHC.Float.double2float

quantize_pixel :: PixelRGBF -> PixelRGB8;
quantize_pixel (PixelRGBF x y z) = let {
q :: Float -> Word8;
q t
| t>1 = 0
| t==1 = 255
| t<0 = 255
| True = floor $ t*256;
}
in uncurryN PixelRGB8 $ mapTuple3 q (x,y,z);

draw_line :: (Pixel px,PrimMonad m) => MutableImage (PrimState m) px -> (Double,Double) -> (Double,Double) -> px -> m();
draw_line pic a b c = mapM_
(\(p,q) -> writePixel pic (fromIntegral p) (fromIntegral q) c) $ line_tuples a b;

line_tuples :: (Double,Double) -> (Double,Double) -> [(Integer,Integer)];
line_tuples (ax,ay) (bx,by) = let {
x = bx-ax;
y = by-ay;
maxt = max (abs x)(abs y);
} in do {
t :: Integer <- enumFromTo 0 $ round maxt;
let { d :: Double = (fromIntegral t)/maxt};
return (round $ ax+d*x,round $ ay+d*y);
};
{-
in if x < y
then map Tuple.swap $ line_tuples (swap a) (swap b)
else do {
t::Integer <- enumFromTo 0 $ round x;
-- this aint bresenham
return (t, round $ y*(fromIntegral t)/x);
};
-}

cieLAB_l :: Lightness -> Double -> Double -> C;
cieLAB_l (Lightness l) = cieLAB white_point l;

ab_find2 :: Lightness -> Double -> (Double,Double);
ab_find2 l theta = let {
a = cos (2*pi*theta);
b = sin (2*pi*theta);
mkcolor :: Double -> C;
mkcolor t = cieLAB_l l (t*a) (t*b);
mktuple :: Double -> (Double,Double);
mktuple t = (t*a,t*b);
f :: Double -> Double;
-- doing it backward because we want the last point just inside, according to "right" in root_bisect
f t = if and $ map in01 $ rgblist $ toSRGB $ mkcolor t then 1 else -1;
far = 200
} in assert (f far<0)
$ assert (0 < f 0)
$ mktuple $ root_bisect f far 0;

ab_find :: Lightness -> Double -> C;
ab_find l theta = (uncurry $ cieLAB_l l) $ ab_find2 l theta;

range01 :: Integer -> [Double];
range01 width =
map (\(i::Double) -> i/(fromIntegral width))
$ genericTake width
$ enumFrom 0;

edge_from :: Lightness -> Double -> Double -> Double;
edge_from l x y = step_the_desired_distance (ab_find l x)(ab_find l y);

newtype Lightness = Lightness Double deriving (Show, Ord, Eq);

all_edge :: Lightness -> [C];
all_edge l = map (uncurry $ cieLAB_l l) $ all_ab l;
-- map (ab_find l) $ takeWhile (<1) $ iterate (next_x $ edge_from l ) 0;

all_ab :: Lightness -> [(Double,Double)];
all_ab l = map (ab_find2 l) $ takeWhile (<1) $ iterate (next_x $ edge_from l ) 0;

outppm_band :: [C] -> IO();
outppm_band l = mapM_ putStrLn $ ["P3"
,show $ length l
,"1","255"]
++ map outcolor l;

red_lightness :: Lightness;
red_lightness = Lightness $ lightness white_point $ sRGB 1 0 0;

many_lab :: Integer -> Integer -> IO();
many_lab width height = do {
 putStrLn "P3";
 print width;
 print height;
 putStrLn "255";
 sequence_ $ do {
 l :: Lightness <- range01 height >>= return . Lightness . (*100);
 return$ mapM_ putStrLn $ map outcolor $ map (ab_find l) $ range01 width;
}
};

length_band :: Lightness -> Integer;
length_band (Lightness 0) = 0;
length_band l = genericLength $ all_edge l;

scaled_pband :: Lightness -> Integer -> [C];
scaled_pband (Lightness 0) n = genericReplicate n Color.black;
scaled_pband l n = let {
band :: [C];
band = all_edge l;
index :: Double -> Integer;
index x = floor $ x*(fromIntegral $ length band);
} in trace (show l ++ " " ++ (show $ length band)) $ do {
t <- range01 n;
return $ genericIndex band $ index t;
};

numLines :: Integer;
numLines=1024; -- 0 and 1024 not used

line_to_light :: Integer -> Lightness;
line_to_light n = Lightness $ (fromInteger n)/(fromInteger numLines)*100;

numColumns :: Integer;
numColumns = 1024;

lightness_of :: Double -> Double -> Double -> Lightness;
lightness_of r g b = Lightness $ lightness white_point $ sRGB r g b;

lightness_plot :: Lightness -> Double -> Double -> C;
lightness_plot l a b = let {
c :: C;
c = cieLAB_l l a b;
} in if and $ map in01 $ rgblist $ toSRGB c then c else Color.black;

slice :: Lightness -> IO ();
slice l = do {
let {size :: Integer; size= image_width};
let {grow :: Double -> Double;
grow x = (realToFrac offset)*(x*2-1);
};
putStrLn "P3";
print size;print size;
putStrLn "255";
sequence_ $ do {
b <- range01 size >>= return . grow . (1-);
a <- range01 size >>= return . grow;
return $ putStrLn $ outcolor $ lightness_plot l a b;
};
};

one_radial_ray :: MutableImage (PrimState IO) PixelRGBF -> Lightness -> IO();
one_radial_ray pic l = do {
-- nphi might be better than randomIO
c :: C <- randomIO >>= return . ab_find l;
let {(_,a,b) = cieLABView white_point c;
q = scale_location (a,b);
};
if (sqrt(a*a+b*b)>200) then
  error  $ "one_radial_ray " ++ show (l,a,b);
else
  return ();
draw_line pic (scale_location (0,0)) q $ fpixel c;
};

outpic :: MutableImage (PrimState IO) PixelRGBF -> IO();
outpic pic = freezeImage pic >>= ByteString.putStr . encodePng . pixelMap quantize_pixel;

radial_slice :: Lightness -> IO();
radial_slice l = do {
setStdGen $ mkStdGen 1;
_ :: Double <- randomIO;  -- avoid the first sample;
let {bg :: Float = if l > Lightness 50 then 0 else 1};
pic :: MutableImage (PrimState IO) PixelRGBF <- standard_image $ PixelRGBF bg bg bg;
-- 10000 takes about 17 seconds compiled
-- 100000 is probably fine
-- 1000000 is overkill
genericReplicateM_ num_rays $ one_radial_ray pic l;
outpic pic;
};

-- not in Control.Monad!
genericReplicateM_ :: (Monad m, Integral i) => i -> m a -> m ();
genericReplicateM_ n x = sequence_ (genericReplicate n x);

one_rainbow_ray :: MutableImage (PrimState IO) PixelRGBF -> IO();
one_rainbow_ray pic = do {
c :: C <- randomRIO (0,6) >>= return . continuous3rgb;
let {(_,a,b) = cieLABView white_point c;
q = scale_location (a,b);
};
draw_line pic (scale_location (0,0)) q $ fpixel c;
};

rainbow_distance :: Double -> Double -> Double;
rainbow_distance p q = step_the_desired_distance (continuous3rainbow p)(continuous3rainbow q);

step_the_desired_distance :: C -> C -> Double;
step_the_desired_distance p q = color_difference p q - desired_difference;

-- perceptually uniform rainbow
all_rainbow :: [C];
all_rainbow = map continuous3rainbow $ takeWhile (<6) $ iterate (next_x rainbow_distance) 0;

rainbow_in_proper_order :: (a,a,a) -> (a,a,a);
rainbow_in_proper_order (x,y,z) = (z,y,x);

continuous3rainbow :: Double -> C;
continuous3rainbow t = uncurryN sRGB $ rainbow_in_proper_order $ continuous3 t;

uncorrected_rainbow :: [C];
uncorrected_rainbow = map tupleToColor $ map rainbow_in_proper_order three_range;

read_goline :: Integer -> IO [C];
read_goline n = readFile ("goline."++ show n)  >>= return . map read .lines;

read_and_emit :: Integer -> IO ();
read_and_emit n = read_goline n >>= mapM_ (putStrLn .outcolor);

-- emits the giant collection of perceptually scaled colors by lightness
all_read_goline :: IO ();
all_read_goline = do {
 putStrLn "P3\n1024 1023\n255";
 mapM_ read_and_emit [1.. (pred numLines)];
};

to_float_color :: C -> Colour Float;
to_float_color = uncurryN Linear.rgb . unlist3 . map realToFrac . rgblist . Linear.toRGB;

colors_by_lightness :: [C];
colors_by_lightness = [sRGB 0 0 0, sRGB 0 0 1, sRGB 1 0 0, sRGB 1 0 1,  sRGB 0 1 0, sRGB 0 1 1, sRGB 1 1 0, sRGB 1 1 1];

and_midpoints :: (Num a, Fractional a) => [a] -> [a];
and_midpoints [] = [];
and_midpoints [x] = [x];
and_midpoints (x: yrest@(y:_)) = x:((x+y)/2):(and_midpoints yrest);

lightness_index :: [Lightness];
lightness_index = map Lightness $ and_midpoints $ map (lightness white_point) colors_by_lightness;


distance_from_edge :: (Ord a, Num a) => a -> a -> a;
distance_from_edge t y = assert (0<=t) $ assert (t<y)
$ min t (y-1-t);

frame_distance :: (Ord a, Num a) => a -> a -> (a,a) -> a;
frame_distance x y (x1, y1) = min (distance_from_edge x x1)
(distance_from_edge y y1);

distance_color_scaling :: Double;
distance_color_scaling = 30;

distance_to_color :: Integer -> C;
distance_to_color i = let {z = fromIntegral i / distance_color_scaling } in
assert (z>=0) $ if z <6 then continuous3rainbow $ colorshift $ z else Color.black;

colorshift :: Double -> Double;
colorshift = (6*) . frac . (+ (1/6)) . (/6);

out_rainbow_frame :: Integer -> Integer -> IO ();
out_rainbow_frame width height = do {
putStrLn "P3";
print width;
print height;
putStrLn "255";
sequence_ $ do {
j :: Integer <- genericTake height $ enumFrom 0;
i :: Integer <- genericTake width $ enumFrom 0;
return $ putStrLn $ outcolor $ distance_to_color $ frame_distance i j (width,height);
};
};

standard_image :: (Pixel px, PrimMonad m) => px -> m (MutableImage (PrimState m) px);
standard_image bg = createMutableImage (fromIntegral image_width) (fromIntegral image_height) bg;

out_rainbow_circle :: Double -> IO();
out_rainbow_circle shift = do {
pic :: MutableImage (PrimState IO) PixelRGBF <- standard_image $ PixelRGBF 0 0 0;
sequence_ $ do {
t <- range01 num_rays;
let {
cx = fromIntegral image_width /2;
cy = fromIntegral image_height /2;
r = 0.8* min cx cy;
th = 2*pi*t;
};
return $ draw_line pic (cx,cy) (cx+r*sin th, cy+r*cos th) $ fpixel $ continuous3rainbow $ shift + t*6;
};
outpic pic;
};

tolab :: Lightness -> (Double, Double) -> C;
tolab l (a,b) = cieLAB_l l a b;

read_perim :: Integer -> IO [C];
read_perim n = readFile ("perim."++show n) >>= return . map (tolab (line_to_light n) . read) . lines;

smooth_index :: Integer -> Double -> Integer;
smooth_index n x = floor $ x * fromIntegral n;

c_rescale :: Integer -> [a] -> [a];
c_rescale desired l = do {
t :: Double <- range01 desired;
return $ genericIndex l $ smooth_index (genericLength l) t;
};

all_read_perim :: IO();
all_read_perim = do {
 putStrLn "P3";
 print numColumns; print numLines; putStrLn "255";
 genericReplicateM_ numColumns (putStrLn $ outcolor Color.black);
 mapM_ (\n -> read_perim n >>= mapM_ (putStrLn . outcolor) . c_rescale numColumns) [1 .. pred numLines];
};

uncorrected_rainbow_of_width :: Integer -> [C];
uncorrected_rainbow_of_width n = map (continuous3rainbow .  (*6)) $ range01 n;

circular_shift :: Integer -> [a] -> [a];
circular_shift n l = genericDrop n l ++ genericTake n l;

rainbow_shift_diff :: Integer -> [Double];
rainbow_shift_diff n = map sqrt $ zipWith color_difference all_rainbow (circular_shift n all_rainbow);

range_plot :: Integer -> [Double] -> IO();
range_plot width l = do {
putStrLn "P1";
print width;
print $ length l;
mapM_ (\x ->
let {
n1 :: Integer;
n1 = round $ fromIntegral width * x;
n2 = width - n1;
} in do {
genericReplicateM_ n1 $ putStrLn "1";
genericReplicateM_ n2 $ putStrLn "0";
}) l;
};

rainbow_shift_plot :: Integer -> Integer -> IO();
rainbow_shift_plot width n = let {
xs = rainbow_shift_diff n;
} in range_plot width $ map (/ (maximum xs)) xs;

out_generic_circle :: (Double -> C) -> IO();
out_generic_circle f = do {
pic :: MutableImage (PrimState IO) PixelRGBF <- standard_image $ PixelRGBF 0 0 0;
sequence_ $ do {
t <- range01 num_rays;
let {
cx = fromIntegral image_width /2;
cy = fromIntegral image_height /2;
r = 0.8* min cx cy;
th = 2*pi*t;
};
return $ draw_line pic (cx,cy) (cx+r*sin th, cy+r*cos th) $ fpixel $ f t;
};
outpic pic;
};

perceptually_uniform_rainbow :: Double -> C;
perceptually_uniform_rainbow t = index_by_float all_rainbow t;

index_by_float :: [a] -> Double -> a;
index_by_float band x = genericIndex band $ (floor $ (frac x)*(fromIntegral $ length band)::Integer);

out_perceptually_uniform_rainbow_circle :: Double -> IO ();
out_perceptually_uniform_rainbow_circle shift = out_generic_circle (\t -> perceptually_uniform_rainbow $ t-shift);

} -- end of program