{-

Color space conversion between RGB and Y"CbCr implemented with exact rational arithmetic

Copyright 2020 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/>.


Usage:

Run without arguments to print out the conversion to YCbCr, then back to RGB, of all 16777216 RGB colors.  Of interest might be lines whose final RGB value is the same as the starting RGB value.

The "halfoff fraction s" arguments prints colors which require an intermediate value "near" 1/2 to be rounded, where near is defined as the window s*fraction <= x < (s+1)*fraction.

This lists all the differences with cjpeg when distance from 1/2 is in the range 1/1024 .. 2/2024.  There are only 2 colors in that window (0 207 35 and 255 48 220) among the 8140 tests.  Their distance from 1/2 is 1/1000, with which we used to choose farenoughthreshold.

time ./a.out halfoff "1%1024" 1 > p2 && wc -l p2 && time for file in $(perl -nlwae 'print$F[1]' p2) ; do ppmmake $file 16 16 | cjpeg -q 100 -dct float | djpeg -dct float | pnmnoraw | head -4 | tail -1 | perl -nlwae 'print"$F[0] $F[1] $F[2]"' ; done > p3 && perl -nlwe 'die unless /^\[(\d+),(\d+),(\d+)/;print"$1 $2 $3"' p2 > p4 && diff p4 p3

-}

{-# LANGUAGE ScopedTypeVariables, LambdaCase, PackageImports #-}
module Main(main,trace_placeholder) where {
import System.Environment(getArgs);
import Control.Exception(assert);
import Debug.Trace(trace);
import Data.Function((&));
import Control.Category((>>>));
import Prelude hiding((.),(>>),round,exponent);
import qualified Prelude;
import System.IO(hSetBuffering,stdout,BufferMode(LineBuffering));
--import System.IO(stderr,hPutStrLn);
import qualified Data.List as List;
import qualified Control.Monad as Monad;
--import Control.Monad(guard);
--import qualified Data.Maybe as Maybe;
--import qualified Data.Map as Map; import Data.Map(Map);
--import qualified Data.Set as Set; import Data.Set(Set);
--import qualified Data.Bifunctor as Bifunctor;
--import qualified Data.Tuple as Tuple;
import Text.Printf(printf);

type Ii = Integer;

-- to avoid the redundancy warning
trace_placeholder :: ();
trace_placeholder = (trace,assert) & (id >>> error) "trace_placeholder";

type Endo a = a->a;

main :: IO();
main = getArgs >>= \case{
[] -> mapM_ showevolution run255;
("do1":rest) -> showevolution (map read rest);
("yuv1":rest) -> evolutionycbcr rgbround (map read rest);
-- sort -t, -k2n,2 -k3n yuv2rgb.out | ./a.out ysizes
["ysizes"] -> groupbychroma >>= mapM_ printsize;
["y4sum"] -> groupbychroma >>= (map (snd >>> List.genericLength >>> (id :: Endo Integer) >>> (\x -> x^(4 :: Integer))) >>> sum >>> print);
["matrix"] -> outputmatrices;
["halfoff",blocksize,startblock] -> do {
hSetBuffering stdout LineBuffering;
filter (potentiallyroundedwrong (read blocksize) (read startblock)) run255 & mapM_ showevolution;
};
_ -> undefined;
};


toycbcr_dot :: (Num a) => [a] -> [a] -> a;
toycbcr_dot (c:coeff) rgb = c+dot coeff rgb;
toycbcr_dot [] _ = error "toycbcr_dot";

dot :: (Num a) => [a] -> [a] -> a;
dot x y = sum $ zipwith_check_same_length (*) x y;

zipwith_check_same_length :: (a->b->c) -> [a] -> [b] -> [c];
zipwith_check_same_length _ [] [] = [];
zipwith_check_same_length f (a:arest) (b:brest) = f a b:zipwith_check_same_length f arest brest;
zipwith_check_same_length _ _ _ = error "not same length";

toycbcr_rational :: [Rational] -> [Rational];
toycbcr_rational rgb = map (\row -> toycbcr_dot row rgb) toycbcr_matrix;

torgb_dot :: (Num a) => [a] -> [a] -> a;
torgb_dot coeff ycbcr = dot (1:coeff) (offsetc ycbcr);
-- hopefully the repeated calls to offsetc get optimized

offsetc :: (Num a) => [a] -> [a];
offsetc [y,cb,cr] = [y,cb-128,cr-128];
offsetc _ = error "offsetc";

torgb_rational :: [Rational] -> [Rational];
torgb_rational ycbcr = map (\row -> torgb_dot row ycbcr) torgb_matrix;

toycbcr_matrix_approx :: [[Rational]];
toycbcr_matrix_approx = List.transpose
-- https://en.wikipedia.org/w/index.php?title=YCbCr&oldid=946790958
[[0    ,128      ,128]
,[0.299,-0.168736,0.5]
,[0.587,-0.331264,-0.418688]
,[0.114,0.5      ,-0.081312]
];

torgb_matrix_approx :: [[Rational]];
torgb_matrix_approx = List.transpose
[[0    ,-0.344136,1.772]
,[1.402,-0.714136,0]
];

-- Rec. ITU-T T.871 (05/2011)
toycbcr_matrix_exact :: [[Rational]];
toycbcr_matrix_exact =
[   [0, 0.299,  0.587,  0.114]
,128:([-0.299, -0.587,  0.886] & map (/1.772))
,128:([ 0.701, -0.587, -0.114] & map (/1.402))
];

torgb_matrix_exact :: [[Rational]];
torgb_matrix_exact =
[[0,                   1.402]
,[-0.114*1.772, -0.299*1.402] & map (/0.587)
,[       1.772, 0]
];

onlypartial :: Bool;
onlypartial = False;

partial :: [Ii];
partial = [0..16];

run255 :: [[Ii]];
run255 = Monad.replicateM 3
$ if onlypartial then partialrange else [0..255];

toycbcr_matrix :: [[Rational]];
toycbcr_matrix = toycbcr_matrix_exact;
--toycbcr_matrix = toycbcr_matrix_approx;

torgb_matrix :: [[Rational]];
torgb_matrix = torgb_matrix_exact;
--torgb_matrix = torgb_matrix_approx;

clamp :: Ii -> Ii;
clamp x = if x<0 then 0 else if x>255 then 255 else x;

roundtoeven :: Rational -> Integer;
roundtoeven = Prelude.round >>> clamp;

-- traditional rounding
roundingup :: Rational -> Integer;
roundingup x = floor (x+0.5) & clamp;

roundingdown :: Rational -> Integer;
roundingdown x = let {
tentative = roundingup x;
} in if (fromInteger tentative - x == 1/2)  -- amazingly no space necessary
then tentative-1 & clamp
else tentative & clamp;

partialrange :: [Ii];
partialrange = partial ++ (map (\x -> 255-x) partial & reverse);

todouble :: Rational -> Double;
todouble = fromRational;

printrational :: [Rational] -> String;
printrational l@[_,_,_] = printmanyrational l;
printrational _ = error "printrational";

printonerational :: Rational -> String;
printonerational x = printf "%.4f = {%s}" (todouble x) (show x);

printmanyrational :: [Rational] -> String;
printmanyrational l = "[" ++ (map printonerational l & List.intersperse "," & concat) ++ "]";

-- cjpeg
yuvround :: Rational -> Ii;
-- yuvround = roundingdown; -- closest to cjpeg
yuvround = parserounding "up"; -- what the spec says

-- djpeg
rgbround :: Rational -> Ii;
rgbround = parserounding "up";

showrgbhex :: [Ii] -> String;
showrgbhex [r,g,b] = printf "rgb:%02x/%02x/%02x" r g b;
showrgbhex e = error $"showrgbhex " ++ show e;

parserounding :: String -> Rational -> Ii;
parserounding "even" = roundtoeven;
parserounding "down" = roundingdown;
parserounding "up" = roundingup;
parserounding _ = error "parserounding";

showevolution :: [Ii] -> IO ();
showevolution = showevolution_rounding yuvround rgbround;

showevolution_rounding :: (Rational -> Ii) -> (Rational -> Ii) -> [Ii] -> IO ();
showevolution_rounding myyuvround myrgbround rgb1 = do {
putStr $ show rgb1 ++ " " ++ showrgbhex rgb1;
let {yuv1 :: [Rational]; yuv1 = map fromInteger rgb1 & toycbcr_rational };
putStr $ " " ++ printrational yuv1;
putStr $ " ishalfyuv="++(any ishalf yuv1 & show);
putStr $ " toocloseyuv="++(any tooclose yuv1 & show);
let {yuv2 :: [Ii]; yuv2 = map myyuvround yuv1};
evolutionycbcr myrgbround yuv2
};

evolutionycbcr :: (Rational -> Ii) -> [Ii] -> IO ();
evolutionycbcr myrgbround yuv2 = do {
putStr $ " ycbcr8 " ++ show yuv2;
let {rgb2 :: [Rational]; rgb2 = map fromInteger yuv2 & torgb_rational };
putStr $ " " ++ printrational rgb2;
putStr $ " ishalfrgb="++(any ishalf rgb2& show);
putStr $ " tooclosergb="++(any tooclose rgb2 & show);
putStr $ " " ++ (show $ map myrgbround rgb2);
putStrLn "";
};

-- expects sorted by 2nd and 3rd key
groupbychroma :: IO [([Ii],[Ii])];
groupbychroma = getContents >>= (lines >>> map (read >>> totuple) >>> List.groupBy f >>> map collate >>> return)
where {
totuple :: [Ii] -> ([Ii],Ii);
totuple [y,cb,cr] = ([cb,cr],y);
totuple _ = error "groupbychroma";
f :: ([Ii],Ii) -> ([Ii],Ii) -> Bool;
f x y = fst x == fst y;

collate :: [([Ii],Ii)] -> ([Ii],[Ii]);
collate x = (fst $ head x, map snd x);
};

printsize :: ([Ii],[Ii]) -> IO ();
printsize (chroma, ys) = putStrLn $ show chroma ++ " " ++ (show $ length ys);

ishalf :: Rational -> Bool;
ishalf x = distancefromhalf x == 0;

printlist :: [String] -> String;
printlist l = "[" ++ (List.intersperse "," l & concat) ++ "]";

outputmatrices :: IO ();
outputmatrices = do {
putStr "to_ycbcr_matrix ";
map printmanyrational toycbcr_matrix_exact & printlist & putStrLn;

putStr "to_rgb_matrix ";
map printmanyrational torgb_matrix_exact & printlist & putStrLn;

};

inseparatedwindow :: (Num a, Ord a) => a -> Integer -> a -> Bool;
inseparatedwindow blocksize startblock x = (blocksize * (fromInteger startblock) <= x) && (x < blocksize * (1+startblock & fromInteger));

distancefromhalf :: Rational -> Rational;
distancefromhalf x = (x - (x & floor & fromInteger)) - 1/2 & abs;

-- how far off from 1/2 can you get and still get rounded the wrong way?
nearx :: Rational -> Integer -> Rational -> Bool;
nearx blocksize  startblock x = inseparatedwindow blocksize startblock (distancefromhalf x);

calcpieces :: [Ii] -> ([Rational],[Rational],[Ii]);
calcpieces rgb1 = let
{ yuv1 :: [Rational]
; yuv1 = map fromInteger rgb1 & toycbcr_rational
; yuv2 :: [Ii]
; yuv2 = map yuvround yuv1
; rgb2 :: [Rational]
; rgb2 = map fromInteger yuv2 & torgb_rational
; rgb3 :: [Ii];
; rgb3 = map rgbround rgb2;
} in (yuv1,rgb2,rgb3);

nottooclose :: Rational -> Integer -> Rational -> Bool;
nottooclose blocksize startblock x = (blocksize * fromInteger startblock) < distancefromhalf x;

potentiallyroundedwrong :: Rational -> Integer -> [Ii] -> Bool;
potentiallyroundedwrong blocksize startblock rgb = let {(x,y,final)=calcpieces rgb} in (rgb==final) && (any (nearx blocksize startblock) $ x ++ y) && (all (nottooclose blocksize startblock) $ x ++ y);

farenoughthreshold :: Rational;
farenoughthreshold = 1/1000;  -- note 1/1000 is tight for errors in 0 207 35 and 255 48 220

tooclose :: Rational -> Bool;
tooclose x = distancefromhalf x <= farenoughthreshold;  -- equality happens at 0 207 35 and 255 48 220

} --end