// Title: CIELUV Color Solid Cube Isosurface
// Author: Michael Horvath, Christoph Lipka
// Website: http://isometricland.net
// Created: 2016-12-06
// Updated: 2016-12-06
// This file is licensed under the terms of the CC-GNU LGPL.
// http://www.gnu.org/licenses/lgpl-2.1.html
// Illuminant =  D65
// Observer   =  2° (1931)
// +kfi0 +kff15 +kc

#version 3.7;
#include "math.inc"
#include "colors.inc"
#declare ScaleAmount = 200;
#declare RotateAmount = 15;
#declare AxesScale = 8/8;


global_settings
{
	assumed_gamma	1.0
	adc_bailout	0.005
	max_trace_level	50
	ambient_light	1/3
/*
	ambient_light	0
	radiosity
	{
		pretrace_start	0.08
		pretrace_end	0.01
		count		50
		error_bound	0.1
		recursion_limit	1
		normal		on
		brightness	0.8
		always_sample	yes
		gray_threshold	0.8
		media		on
	}
*/
}

camera
{
//	#local cam_distance	= 40;	// render at 8192x8192
//	#local cam_planesize	= 4;	// render at 8192x8192
	#local cam_distance	= 17;
	#local cam_planesize	= 1.7;
//	#local cam_distance	= 20;	// render at 4096x4096
//	#local cam_planesize	= 2;	// render at 4096x4096
	#local cam_aspectratio	= image_width/image_height;
	orthographic
	location	-z * cam_distance
	direction	+z * cam_distance
	right		+x * cam_planesize*cam_aspectratio
	up		+y * cam_planesize
//	rotate		+x * asind(tand(30))
//	rotate		+x * 90
	rotate		-y * 045
//	rotate		+y * 180
	rotate		+y * RotateAmount
	translate	1/2
	scale		ScaleAmount
}

light_source
{
	<-32,+32,-32,>
	color rgb	1
	rotate		-y * 330
	rotate		+y * RotateAmount
	parallel
	point_at	0
	shadowless
	scale		ScaleAmount
}

light_source
{
	<-32,+32,-32,>
	color rgb	1
	rotate		-y * 90
	rotate		+y * RotateAmount
	parallel
	point_at	0
	shadowless	// was disabled originally
	scale		ScaleAmount
}

//background {srgb 1}


#declare XYZWhiteReference1 =  95.047;
#declare XYZWhiteReference2 = 100.000;
#declare XYZWhiteReference3 = 108.883;
#declare XYZEpsilon = 0.008856;
#declare XYZKappa = 903.3;
#declare uPrime = 4 * XYZWhiteReference1/(XYZWhiteReference1 + 15 * XYZWhiteReference2 + 3 * XYZWhiteReference3);
#declare vPrime = 9 * XYZWhiteReference2/(XYZWhiteReference1 + 15 * XYZWhiteReference2 + 3 * XYZWhiteReference3);


// utility
#declare selectRGBa = function(C,D) {select(C,-1,D)}
#declare selectRGBb = function(C,D) {select(C-1,D,D,-1)}
#declare linearRGBa = function(C) {C/12.92}
#declare linearRGBb = function(C) {pow((C + 0.055)/(1 + 0.055), 2.4)}
#declare linearRGBc = function(C) {select(C-0.04045,linearRGBa(C),linearRGBa(C),linearRGBb(C))}


// input H = between 0 and 360
// input C = between 0 and 128
// input L = between 0 and 100
// output L = between 0 and 100
// output U = between -128 and +128
// output V = between -128 and +128
#declare funcHCL2LUVa1 = function(H,C,L) {L}
#declare funcHCL2LUVa2 = function(H,C,L) {cosd(H) * C}
#declare funcHCL2LUVa3 = function(H,C,L) {sind(H) * C}


// input L = between 0 and 100
// input U = between -128 and +128
// input V = between -128 and +128
// output X: between 0 and 100
// output Y: between 0 and 100
// output Z: between 0 and 100
#declare funcLUV2XYZa1 = function(L,U,V) {1/3 * (52 * L/(U + 13 * L * uPrime) - 1)}		// a
#declare funcLUV2XYZb2 = function(L,U,V) {select(XYZKappa * XYZEpsilon - L, pow((L + 16)/116,3), L/XYZKappa)}		// Y
#declare funcLUV2XYZa2 = function(L,U,V) {-5 * funcLUV2XYZb2(L,U,V)}		// b
#declare funcLUV2XYZa3 = function(L,U,V) {-1/3}		// c
#declare funcLUV2XYZa4 = function(L,U,V) {funcLUV2XYZb2(L,U,V) * (39 * L/(V + 13 * L * vPrime) - 5)}		// d
#declare funcLUV2XYZb1 = function(L,U,V) {(funcLUV2XYZa4(L,U,V) - funcLUV2XYZa2(L,U,V))/(funcLUV2XYZa1(L,U,V) - funcLUV2XYZa3(L,U,V))}		// X
#declare funcLUV2XYZb3 = function(L,U,V) {funcLUV2XYZb1(L,U,V) * funcLUV2XYZa1(L,U,V) + funcLUV2XYZa2(L,U,V)}		// Z
#declare funcLUV2XYZc1 = function(L,U,V) {funcLUV2XYZb1(L,U,V) * 100}		// X
#declare funcLUV2XYZc2 = function(L,U,V) {funcLUV2XYZb2(L,U,V) * 100}		// Y
#declare funcLUV2XYZc3 = function(L,U,V) {funcLUV2XYZb3(L,U,V) * 100}		// Z


// input L = between 0 and 100
// input C = between 0 and 128
// input H = between 0 and 360
// output L = between 0 and 100
// output A = between -128 and +128
// output B = between -128 and +128
#declare funcLCH2LABa1 = function(L,C,H) {L}
#declare funcLCH2LABa2 = function(L,C,H) {cosd(H) * C}
#declare funcLCH2LABa3 = function(L,C,H) {sind(H) * C}


// input L = between 0 and 100
// input A = between -128 and +128
// input B = between -128 and +128
// output X: between 0 and 100
// output Y: between 0 and 100
// output Z: between 0 and 100
// Note that some of these functions are out of order.
#declare funcLAB2XYZa2 = function(L,A,B) {(L + 16) / 116}
#declare funcLAB2XYZa1 = function(L,A,B) {funcLAB2XYZa2(L,A,B) + A / 500}
#declare funcLAB2XYZa3 = function(L,A,B) {funcLAB2XYZa2(L,A,B) - B / 200}
#declare funcLAB2XYZb1 = function(L,A,B) {XYZWhiteReference1 * select(pow(funcLAB2XYZa1(L,A,B), 3) - XYZEpsilon, (funcLAB2XYZa1(L,A,B) - 16.0 / 116.0) / 7.787, (funcLAB2XYZa1(L,A,B) - 16.0 / 116.0) / 7.787, pow(funcLAB2XYZa1(L,A,B), 3))}
#declare funcLAB2XYZb2 = function(L,A,B) {XYZWhiteReference2 * select(L - XYZKappa * XYZEpsilon, L / XYZKappa, L / XYZKappa, pow(((L + 16) / 116), 3))}
#declare funcLAB2XYZb3 = function(L,A,B) {XYZWhiteReference3 * select(pow(funcLAB2XYZa3(L,A,B), 3) - XYZEpsilon, (funcLAB2XYZa3(L,A,B) - 16.0 / 116.0) / 7.787, (funcLAB2XYZa3(L,A,B) - 16.0 / 116.0) / 7.787, pow(funcLAB2XYZa3(L,A,B), 3))}


// input X: between 0 and 100
// input Y: between 0 and 100
// input Z: between 0 and 100
// output R: between 0 and 1
// output G: between 0 and 1
// output B: between 0 and 1
// Note that out-of-range colors are *NOT* corrected. You will have to do this yourself.
#declare funcXYZ2RGBa1 = function(X,Y,Z) {X/100 *  3.2406 + Y/100 * -1.5372 + Z/100 * -0.4986}
#declare funcXYZ2RGBa2 = function(X,Y,Z) {X/100 * -0.9689 + Y/100 *  1.8758 + Z/100 *  0.0415}
#declare funcXYZ2RGBa3 = function(X,Y,Z) {X/100 *  0.0557 + Y/100 * -0.2040 + Z/100 *  1.0570}
#declare funcXYZ2RGBb1 = function(X,Y,Z) {select(funcXYZ2RGBa1(X,Y,Z) - 0.0031308, 12.92 * funcXYZ2RGBa1(X,Y,Z), 12.92 * funcXYZ2RGBa1(X,Y,Z), 1.055 * pow(funcXYZ2RGBa1(X,Y,Z), 1 / 2.4) - 0.055)}
#declare funcXYZ2RGBb2 = function(X,Y,Z) {select(funcXYZ2RGBa2(X,Y,Z) - 0.0031308, 12.92 * funcXYZ2RGBa2(X,Y,Z), 12.92 * funcXYZ2RGBa2(X,Y,Z), 1.055 * pow(funcXYZ2RGBa2(X,Y,Z), 1 / 2.4) - 0.055)}
#declare funcXYZ2RGBb3 = function(X,Y,Z) {select(funcXYZ2RGBa3(X,Y,Z) - 0.0031308, 12.92 * funcXYZ2RGBa3(X,Y,Z), 12.92 * funcXYZ2RGBa3(X,Y,Z), 1.055 * pow(funcXYZ2RGBa3(X,Y,Z), 1 / 2.4) - 0.055)}


// input R: between 0 and 1
// input G: between 0 and 1
// input B: between 0 and 1
// output R: between 0 and 1
// output G: between 0 and 1
// output B: between 0 and 1
// only uses these when converting individual values, not inside isosurfaces or function patterns
#declare correctRGBa1 = function(R,G,B) {selectRGBa(R,selectRGBa(G,selectRGBa(B,R)))}
#declare correctRGBa2 = function(R,G,B) {selectRGBa(R,selectRGBa(G,selectRGBa(B,G)))}
#declare correctRGBa3 = function(R,G,B) {selectRGBa(R,selectRGBa(G,selectRGBa(B,B)))}
#declare correctRGBb1 = function(R,G,B) {selectRGBb(R,selectRGBb(G,selectRGBb(B,correctRGBa1(R,G,B))))}
#declare correctRGBb2 = function(R,G,B) {selectRGBb(R,selectRGBb(G,selectRGBb(B,correctRGBa2(R,G,B))))}
#declare correctRGBb3 = function(R,G,B) {selectRGBb(R,selectRGBb(G,selectRGBb(B,correctRGBa3(R,G,B))))}


// input L = between 0 and 100
// input A = between -128 and +128
// input B = between -128 and +128
// output R: between 0 and 1
// output G: between 0 and 1
// output B: between 0 and 1
#declare convertLAB2RGBa1 = function(L,A,B) {funcXYZ2RGBb1(funcLAB2XYZb1(L,A,B),funcLAB2XYZb2(L,A,B),funcLAB2XYZb3(L,A,B))}
#declare convertLAB2RGBa2 = function(L,A,B) {funcXYZ2RGBb2(funcLAB2XYZb1(L,A,B),funcLAB2XYZb2(L,A,B),funcLAB2XYZb3(L,A,B))}
#declare convertLAB2RGBa3 = function(L,A,B) {funcXYZ2RGBb3(funcLAB2XYZb1(L,A,B),funcLAB2XYZb2(L,A,B),funcLAB2XYZb3(L,A,B))}


// input L = between 0 and 100
// input U = between -128 and +128
// input V = between -128 and +128
// output R: between 0 and 1
// output G: between 0 and 1
// output B: between 0 and 1
#declare convertLUV2RGBa1 = function(L,U,V) {funcXYZ2RGBb1(funcLUV2XYZc1(L,U,V),funcLUV2XYZc2(L,U,V),funcLUV2XYZc3(L,U,V))}
#declare convertLUV2RGBa2 = function(L,U,V) {funcXYZ2RGBb2(funcLUV2XYZc1(L,U,V),funcLUV2XYZc2(L,U,V),funcLUV2XYZc3(L,U,V))}
#declare convertLUV2RGBa3 = function(L,U,V) {funcXYZ2RGBb3(funcLUV2XYZc1(L,U,V),funcLUV2XYZc2(L,U,V),funcLUV2XYZc3(L,U,V))}


// testing
#local VLCH = <25, 50, 270>;
#local VLAB = <funcLCH2LABa1(VLCH.x,VLCH.y,VLCH.z), funcLCH2LABa2(VLCH.x,VLCH.y,VLCH.z), funcLCH2LABa3(VLCH.x,VLCH.y,VLCH.z)>;
#local VXYZ = <funcLAB2XYZb1(VLAB.x,VLAB.y,VLAB.z), funcLAB2XYZb2(VLAB.x,VLAB.y,VLAB.z), funcLAB2XYZb3(VLAB.x,VLAB.y,VLAB.z)>;
#local VRGB = <funcXYZ2RGBb1(VXYZ.x,VXYZ.y,VXYZ.z), funcXYZ2RGBb2(VXYZ.x,VXYZ.y,VXYZ.z), funcXYZ2RGBb3(VXYZ.x,VXYZ.y,VXYZ.z)>;
#local VCOR = < correctRGBb1(VRGB.x,VRGB.y,VRGB.z),  correctRGBb2(VRGB.x,VRGB.y,VRGB.z),  correctRGBb3(VRGB.x,VRGB.y,VRGB.z)>;
#debug concat("\nVLCH = ", vstr(3, VLCH, ", ", 0, 5), "\n")
#debug concat("\nVLAB = ", vstr(3, VLAB, ", ", 0, 5), "\n")
#debug concat("\nVXYZ = ", vstr(3, VXYZ, ", ", 0, 5), "\n")
#debug concat("\nVRGB = ", vstr(3, VRGB, ", ", 0, 5), "\n")
#debug concat("\nVCOR = ", vstr(3, VCOR, ", ", 0, 5), "\n")

//#local VCON = <convertLCH2RGBa1(50, 50, 90),convertLCH2RGBa2(50, 50, 90),convertLCH2RGBa3(50, 50, 90)>;
//#local VCON = <convertLCH2RGBb1(50, 50, 180),convertLCH2RGBb2(50, 50, 180),convertLCH2RGBb3(50, 50, 180)>;
//#debug concat("\nVCON = ", vstr(3, VCON, ", ", 0, 5), "\n")

// input L = between 0 and 100
// input A = between -128 and +128
// input B = between -128 and +128


#declare fD = function(C) {abs(C-0.5)-0.5}
#declare fDist = function(Dr,Dg,Db) {max(Dr,Dg,Db)}
#declare fOutR = function(L,U,V) {convertLUV2RGBa1(L,U,V)}
#declare fOutG = function(L,U,V) {convertLUV2RGBa2(L,U,V)}
#declare fOutB = function(L,U,V) {convertLUV2RGBa3(L,U,V)}
#declare fInpL = function(x,y,z) {y*100}
#declare fInpA = function(x,y,z) {x*256-128}
#declare fInpB = function(x,y,z) {z*256-128}

#declare fFinal = function(L,U,V)
{
	fDist
	(
		fD(fOutR(L,U,V)),
		fD(fOutG(L,U,V)),
		fD(fOutB(L,U,V))
	)
}

#declare pigmentR = pigment
{
	function {linearRGBc(fOutR(fInpL(x,y,z),fInpA(x,y,z),fInpB(x,y,z)))}
	color_map
	{
		[0 color rgb <0,0,0>]
		[1 color rgb <3,0,0>]
	}
}
#declare pigmentG = pigment
{
	function {linearRGBc(fOutG(fInpL(x,y,z),fInpA(x,y,z),fInpB(x,y,z)))}
	color_map
	{
		[0 color rgb <0,0,0>]
		[1 color rgb <0,3,0>]
	}
}
#declare pigmentB = pigment
{
	function {linearRGBc(fOutB(fInpL(x,y,z),fInpA(x,y,z),fInpB(x,y,z)))}
	color_map
	{
		[0 color rgb <0,0,0>]
		[1 color rgb <0,0,3>]
	}
}
#declare pigmentRGB = pigment
{
	average
	pigment_map
	{
		[1 pigmentR]
		[1 pigmentG]
		[1 pigmentB]
	}
}
#declare pigmentCyl = pigment
{
	cylindrical
	color_map
	{
		[ 0/20 color rgbt 10/10]
		[ 1/20 color rgbt 10/10]
		[ 1/20 color rgbt  9/10]
		[20/20 color rgbt  9/10]
	}
	frequency 5
	scale 1.0001
}
#declare pigmentGradY = pigment
{
	gradient y
	pigment_map
	{
		[ 0/20 pigmentCyl]
		[ 1/20 pigmentCyl]
		[ 1/20 color rgbt  9/10]
		[20/20 color rgbt  9/10]
	}
	frequency 5
	scale 1.0001
}


#declare isoShape = isosurface
{
	function {fFinal(fInpL(x,y,z),fInpA(x,y,z),fInpB(x,y,z))}
	threshold	0
	accuracy	0.001
	contained_by
	{
		box {<0,0,0>,<1,1,1>}
	}
	max_gradient	100	// was 20000
}

#declare cubShape = box
{
	0, 1
//	hollow
}

polygon
{
	4, <0, 0,>, <1, 0,>, <1, 1,>, <0, 1,>
	pigment {image_map {png "axis_a_b_normal.png"}}
//	translate	<-0.5, -0.5>
	scale		440/320
	translate	-y * ((440-320)/2 * 440/320)/440
	translate	-x * ((440-320)/2 * 440/320)/440
	rotate		+x * 90
	scale		ScaleAmount
	scale		AxesScale
	translate	-0.000001
}

polygon
{
	4, <0, 0,>, <1, 0,>, <1, 1,>, <0, 1,>
	pigment {image_map {png "axis_L_normal.png"}}
//	translate	<-0.5, -0.5>
	scale		440/320
	translate	-y * ((440-320)/2 * 440/320)/440
	translate	-x * ((440-320)/2 * 440/320)/440
//	rotate		+x * 90
	scale		ScaleAmount
	scale		AxesScale
	translate	-0.000001
}

object
{
	isoShape
	pigment {pigmentRGB}
	finish {emission 0.3 diffuse 1}
	scale ScaleAmount
}

object
{
	cubShape
//	pigment {pigmentGradY}
	pigment {color rgbt  7/8}
	scale ScaleAmount
}