(*
Copyright (c) Jean-Baptiste Rouquier

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 2
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.
**********************************************************************

This is an answer to
http://coherentgraphics.blogspot.com/2008/04/storing-colours-in-31-bits-part-one.html
"
For rendering vector graphics scenes, colours are usually stored with what we call "premultiplied" or "associated" alpha. For instance, opaque dark red is stored as:
0.5, 0, 0, 0.5 (R * A, G * A, B * A, A)
instead of
0.5, 0, 0, 1 (R, G, B, A).
This originally had to do with making compositing algorithms fast (fewer multiplications), but it has other advantages - for instance there is only a single value for the clear colour (0, 0, 0, 0) instead of many (x, y, z, 0).

We usually use 8 bits for each component, packing them into a 32 bit word.

Now, 32 bits can store 2^32 colours. In fact, in the premultiplied scheme, many bit patterns are unused (when R, G or B is > A). It turns out there are only slightly more than 2^30 unique premultiplied colours. In other words, with a suitable mapping, we should be able to store them in OCaml's 31 bit integers. This is important so we can store them in native arrays unboxed, for example.

Such a mapping (togther with a discussion of all this) is in Jim Blinn's book "Dirty Pixels", Chapter 20. Unfortunately, it's too slow for practical use. Can you think of a fast one?
"




(UTF8 symbols ahead, ask me for another version if it is unreadable. "≤" is "less than or equal")

<optional>
The volume of the continuous 4D pyramid {(x,y,z,t) | x,y,z≤t, 0≤x,y,z,t≤1} is 1/4.
Unfortunately, we are in a discrete space, and the boundaries are part of the volume.
Faces, edges and vertices cause the total volume of the discrete pyramid {(r,g,b,a) | r,g,b≤a, 0≤r,g,b,a≤255} to be slightly more than 256^4.
</optional>

Assume for now that we just need to encode {(r,g,b,a) | r,g,b<a, 0≤r,g,b,a≤255}, that is, a is always the unique maximum of {r,g,b,a}.
Swapping a with another coordinate is an injective function that maps our pyramid on another one.
Since there are 4 possible positions for a, this allows us to encode 2 bits.

The idea is to remove the least significant bit (LSB) from r g b a and concatenate the 28 remaining bits.
We also concatenate two of te LSB. The two remaining LSB are encoded with the previous trick (position of a).
This uses 30 bits.
<optional>
As seen before, if we ignore the boundaries of the pyramid, the volume is not more than 1/4 the whole hypercube, that's why 30 bits are enough.
</optional>

Now, remove the assumption. If a is not the only maximum after cutting the LSB, we use another scheme.
The previous scheme used only 30 bits, the last on is used to flag that we use another scheme.
Since we have either r=a or g=a or b=a,
we just need to remember one equality that holds, and encode three of the four channels.
In this scheme, we concatenate
* the flag mask_equality to say that we're using this scheme
* the four LSB
* three flags to say which equality holds (two flags would be enough)
* two padding bits
* three channels

Resulting code is longer than what I hoped.
Jean-Baptiste Rouquier
*)

(* ********************************************************************** *)
let mask_equality  = 0b1000000000000000000000000000000

let mask_r_lsb    = 0b0100000000000000000000000000000
let mask_g_lsb    = 0b0010000000000000000000000000000
let mask_channel3 = 0b0001111111000000000000000000000
let mask_channel2 = 0b0000000000111111100000000000000
let mask_channel1 = 0b0000000000000000011111110000000
let mask_channel0 = 0b0000000000000000000000001111111
	(* used only in the equality scheme: *)
let mask_b_lsb    = 0b0001000000000000000000000000000
let mask_a_lsb    = 0b0000100000000000000000000000000
let mask_r_eq_a   = 0b0000010000000000000000000000000
let mask_g_eq_a   = 0b0000001000000000000000000000000
let mask_b_eq_a   = 0b0000000100000000000000000000000


(** Concatenate four 7-bits integers *)
let concat r g b a =
	(r lsl 21) lor (g lsl 14) lor (b lsl 7) lor a

(** return the index of the maximum int among four. *)
let index_max4 a b c d =
	if a>b then
		if c>d then
			if a>c then 0 else 2
		else
			if a>d then 0 else 3
	else
		if c>d then
			if b>c then 1 else 2
		else
			if b>d then 1 else 3

let colour_of_rgba r g b a =
	let split i = i lsr 1, i land 1 <> 0 in
	let r, r_lsb = split r in
	let g, g_lsb = split g in
	let b, b_lsb = split b in
	let a, a_lsb = split a in
	if r <> a && g <> a && b <> a then
		(if r_lsb then mask_r_lsb else 0) lor
		(if g_lsb then mask_g_lsb else 0) lor
		(if b_lsb then
			 if a_lsb then
				 concat r g b a
			 else
				 concat r g a b
		 else
			 if a_lsb then
				 concat r a b g
			 else
				 concat a g b r
		)
	else
		mask_equality lor
		(if r_lsb then mask_r_lsb else 0) lor
		(if g_lsb then mask_g_lsb else 0) lor
		(if b_lsb then mask_b_lsb else 0) lor
		(if a_lsb then mask_a_lsb else 0) lor
		(if r = a then mask_r_eq_a else 0) lor
		(if g = a then mask_g_eq_a else 0) lor
		(if b = a then mask_b_eq_a else 0) lor
		(if      r=a then   concat 0 g b a
		 else if g=a then   concat 0 r b a
		 else (assert(b=a); concat 0 r g a)
		)

let rgba_of_colour c =
	let unsplit i lsb =
		i lsl 1 lor (if lsb then 1 else 0) in
	let r = ref 0 in
	let g = ref 0 in
	let b = ref 0 in
	let a = ref 0 in
	let r_lsb = ref (mask_r_lsb land c <> 0) in
	let g_lsb = ref (mask_g_lsb land c <> 0) in
	let b_lsb = ref false in
	let a_lsb = ref false in
	if c land mask_equality = 0 then
		let channel3 = (c land mask_channel3) lsr 21 in
		let channel2 = (c land mask_channel2) lsr 14 in
		let channel1 = (c land mask_channel1) lsr 7 in
		let channel0 = (c land mask_channel0) in
		(match index_max4 channel3 channel2 channel1 channel0 with
		 |3 -> b_lsb:=true ; a_lsb:=true ; r:=channel3; g:=channel2; b:=channel1; a:=channel0
		 |2 -> b_lsb:=true ; a_lsb:=false; r:=channel3; g:=channel2; a:=channel1; b:=channel0
		 |1 -> b_lsb:=false; a_lsb:=true ; r:=channel3; a:=channel2; b:=channel1; g:=channel0
		 |n -> assert (n=0);
		       b_lsb:=false; a_lsb:=false; a:=channel3; g:=channel2; b:=channel1; r:=channel0)
	else (
		b_lsb := mask_b_lsb land c <> 0;
		a_lsb := mask_a_lsb land c <> 0;
		let channel2 = (c land mask_channel2) lsr 14 in
		let channel1 = (c land mask_channel1) lsr 7 in
		let channel0 = (c land mask_channel0) in
		a := channel0;
		if            c land mask_r_eq_a <> 0  then (r:= !a; g:=channel2; b:=channel1)
		else       if c land mask_g_eq_a <> 0  then (g:= !a; r:=channel2; b:=channel1)
		else (assert (c land mask_b_eq_a <> 0);     (b:= !a; r:=channel2; g:=channel1));
	);
	unsplit !r !r_lsb,
	unsplit !g !g_lsb,
	unsplit !b !b_lsb,
	unsplit !a !a_lsb
;;

(* Exhaustive test: *)
let _ =
	for a = 0 to 255 do
		Printf.printf "\r%i%!" a; (* testing a takes time \Theta(a^3) *)
		for r = 0 to a do
		for g = 0 to a do
		for b = 0 to a do
			if not (rgba_of_colour (colour_of_rgba r g b a) = (r,g,b,a)) then (
				Printf.printf "Failure for r=%i, g=%i, b=%i, a=%i\n%!" r g b a; exit 1)
		done
		done
		done
	done;
	Printf.printf "\nOK!\n%!"