Inverse_sRGB_to_Clog_P3D60.dctl

// Inverse sRGB to C-Log P3D60

typedef struct {
float2 red; float2 green; float2 blue; float2 white;
} Chromaticities;

typedef struct {
float3 c0, c1, c2;
} mat3;

typedef struct {
float4 c0, c1, c2, c3;
} mat4;

typedef struct {
float x; float y;
} SplineMapPoint;

typedef struct {
float coefsLow[6];float coefsHigh[6];
SplineMapPoint minPoint;SplineMapPoint midPoint;SplineMapPoint maxPoint;
float slopeLow;float slopeHigh;
} SegmentedSplineParams_c5;

typedef struct {
float coefsLow[10];float coefsHigh[10];
SplineMapPoint minPoint;SplineMapPoint midPoint;SplineMapPoint maxPoint;
float slopeLow;float slopeHigh;
} SegmentedSplineParams_c9;

__CONSTANT__ Chromaticities AP0 = { {0.7347f, 0.2653f}, {0.0f, 1.0f}, {0.0001f, -0.077f}, {0.32168f, 0.33767f} };
__CONSTANT__ Chromaticities AP1 = { {0.713f, 0.293f}, {0.165f, 0.83f}, {0.128f, 0.044f}, {0.32168f, 0.33767f} };
__CONSTANT__ Chromaticities REC709_PRI = { {0.64f, 0.33f}, {0.3f, 0.6f}, {0.15f, 0.06f}, {0.3127f, 0.329f} };
__CONSTANT__ Chromaticities P3D65_PRI = { {0.68f, 0.32f}, {0.265f, 0.69f}, {0.15f, 0.06f}, {0.3127f, 0.329f} };
__CONSTANT__ Chromaticities P3D60_PRI = { {0.68f, 0.32f}, {0.265f, 0.69f}, {0.15f, 0.06f}, {0.32168f, 0.33767f} };
__CONSTANT__ float DIM_SURROUND_GAMMA = 0.9811f;
__CONSTANT__ float ODT_SAT_FACTOR = 0.93f;
__CONSTANT__ mat3 MM = {{0.5f,-1.0f,0.5f},{-1.0f,1.0f,0.5f},{0.5f,0.0f,0.0f}};
__CONSTANT__ mat3 CONE_RESP_MAT_BRADFORD = {{0.89510f,-0.75020f,0.03890f},{0.26640f,1.71350f,-0.06850f},{-0.16140f,0.03670f,1.02960f}};

#define AP0_2_XYZ_MAT		RGBtoXYZ_m4( AP0, 1.0f)
#define XYZ_2_AP0_MAT		XYZtoRGB_m4( AP0, 1.0f)
#define AP1_2_XYZ_MAT		RGBtoXYZ_m4( AP1, 1.0f)
#define XYZ_2_AP1_MAT		XYZtoRGB_m4( AP1, 1.0f)
#define AP0_2_AP1_MAT		mult_f44_f44( AP0_2_XYZ_MAT, XYZ_2_AP1_MAT)
#define AP1_2_AP0_MAT		mult_f44_f44( AP1_2_XYZ_MAT, XYZ_2_AP0_MAT)
#define AP1_RGB2Y			make_float3( AP1_2_XYZ_MAT.c0.y, AP1_2_XYZ_MAT.c1.y, AP1_2_XYZ_MAT.c2.y)
#define RGB_2_YAB_MAT		make_mat3( make_float3(1.0f/3.0f, 1.0f/2.0f, 0.0f), make_float3(1.0f/3.0f, -1.0f/4.0f, sqrt3over4), make_float3(1.0f/3.0f, -1.0f/4.0f, -sqrt3over4))
#define REC709_2_XYZ_MAT	RGBtoXYZ( REC709_PRI, 1.0f)   
#define REC709_RGB2Y		make_float3( REC709_2_XYZ_MAT.c0.y, REC709_2_XYZ_MAT.c1.y, REC709_2_XYZ_MAT.c2.y)
#define D60_2_D65_CAT		calculate_cat_matrix( AP0.white, REC709_PRI.white)
#define ODT_SAT_MAT			calc_sat_adjust_matrix( ODT_SAT_FACTOR, AP1_RGB2Y)
#define CINEMA_WHITE		48.0f
#define CINEMA_BLACK		_powf(10.0f, _log10f(0.02f))

__DEVICE__ mat3 make_mat3( float3 A, float3 B, float3 C) {
mat3 D; D.c0 = A; D.c1 = B; D.c2 = C;
return D;
}

__DEVICE__ mat4 make_mat4( float4 A, float4 B, float4 C, float4 D) {
mat4 E; E.c0 = A; E.c1 = B; E.c2 = C; E.c3 = D;
return E;
}

__DEVICE__ float4 make_float4_f3_f( float3 A, float B) {
return make_float4(A.x, A.y, A.z, B);
}

__DEVICE__ mat4 make_mat4_f33_f(mat3 A, float B) {
mat4 C; C.c0 = make_float4_f3_f(A.c0, 0.0f); C.c1 = make_float4_f3_f(A.c1, 0.0f);
C.c2 = make_float4_f3_f(A.c2, 0.0f); C.c3 = make_float4(0.0f, 0.0f, 0.0f, B);
return C;
}

__DEVICE__ float _pow10f(float x) {
return _powf(10.0f, x);
}

__DEVICE__ float dot_f3_f3(float3 x, float3 y) {
return x.x * y.x + x.y * y.y + x.z * y.z;
}

__DEVICE__ float3 mult_f3_f33( float3 X, mat3 A) {
float r[3];
float x[3] = {X.x, X.y, X.z};
float a[3][3] =	{{A.c0.x, A.c0.y, A.c0.z}, {A.c1.x, A.c1.y, A.c1.z}, {A.c2.x, A.c2.y, A.c2.z}};
for( int i = 0; i < 3; ++i){
r[i] = 0.0f;
for( int j = 0; j < 3; ++j){
r[i] = r[i] + x[j] * a[j][i];}}
return make_float3(r[0], r[1], r[2]);
}

__DEVICE__ mat3 mult_f33_f33( mat3 A, mat3 B){
float r[3][3];
float a[3][3] =	{{A.c0.x, A.c0.y, A.c0.z}, {A.c1.x, A.c1.y, A.c1.z}, {A.c2.x, A.c2.y, A.c2.z}};
float b[3][3] =	{{B.c0.x, B.c0.y, B.c0.z}, {B.c1.x, B.c1.y, B.c1.z}, {B.c2.x, B.c2.y, B.c2.z}};
for( int i = 0; i < 3; ++i){
for( int j = 0; j < 3; ++j){
r[i][j] = 0.0f;
for( int k = 0; k < 3; ++k){
r[i][j] = r[i][j] + a[i][k] * b[k][j];}}}
mat3 R = make_mat3(make_float3(r[0][0], r[0][1], r[0][2]), 
make_float3(r[1][0], r[1][1], r[1][2]), make_float3(r[2][0], r[2][1], r[2][2]));
return R;
}

__DEVICE__ float3 mult_f3_f44( float3 X, mat4 A) {
float r[3];
float x[3] = {X.x, X.y, X.z};
float a[4][4] =	{{A.c0.x, A.c0.y, A.c0.z, A.c0.w},
{A.c1.x, A.c1.y, A.c1.z, A.c1.w},
{A.c2.x, A.c2.y, A.c2.z, A.c2.w},
{A.c3.x, A.c3.y, A.c3.z, A.c3.w}};
for( int i = 0; i < 3; ++i){
r[i] = 0.0f;
for( int j = 0; j < 3; ++j){
r[i] = r[i] + x[j] * a[j][i];}
r[i] = r[i] + a[3][i];}
float s = 1.0f / (x[0] * a[0][3] + x[1] * a[1][3] + x[2] * a[2][3] + a[3][3]);
for( int k = 0; k < 3; ++k){
r[k] = r[k] * s;}
return make_float3(r[0], r[1], r[2]);
}

__DEVICE__ mat3 mult_f_f33( float f, mat3 A) {
float r[3][3];
float a[3][3] =	{{A.c0.x, A.c0.y, A.c0.z}, {A.c1.x, A.c1.y, A.c1.z}, {A.c2.x, A.c2.y, A.c2.z}};
for( int i = 0; i < 3; ++i ){
for( int j = 0; j < 3; ++j ){
r[i][j] = f * a[i][j];}}
mat3 R = make_mat3(make_float3(r[0][0], r[0][1], r[0][2]), 
make_float3(r[1][0], r[1][1], r[1][2]), make_float3(r[2][0], r[2][1], r[2][2]));
return R;
}

__DEVICE__ mat4 mult_f44_f44 (mat4 A, mat4 B) {
float r[4][4];
float a[4][4] = {{A.c0.x, A.c0.y, A.c0.z, A.c0.w},{A.c1.x, A.c1.y, A.c1.z, A.c1.w},
{A.c2.x, A.c2.y, A.c2.z, A.c2.w},{A.c3.x, A.c3.y, A.c3.z, A.c3.w}};
float b[4][4] =	{{B.c0.x, B.c0.y, B.c0.z, B.c0.w},{B.c1.x, B.c1.y, B.c1.z, B.c1.w},
{B.c2.x, B.c2.y, B.c2.z, B.c2.w},{B.c3.x, B.c3.y, B.c3.z, B.c3.w}};
for( int i = 0; i < 4; ++i){
for( int j = 0; j < 4; ++j){
r[i][j] = 0.0f;
for( int k = 0; k < 4; ++k){
r[i][j] = r[i][j] + a[i][k] * b[k][j];}}}
mat4 R = make_mat4(make_float4(r[0][0], r[0][1], r[0][2], r[0][3]), 
make_float4(r[1][0], r[1][1], r[1][2], r[1][3]), make_float4(r[2][0], r[2][1], r[2][2], r[2][3]),
make_float4(r[3][0], r[3][1], r[3][2], r[3][3]));
return R;
}

__DEVICE__ float3 MAX( float3 rgb, float max)
{
return make_float3(_fmaxf(rgb.x, max),_fmaxf(rgb.y, max),_fmaxf(rgb.z, max));
}

__DEVICE__ mat3 invert_f33 (mat3 A) {
mat3 R;
float result[3][3];
float a[3][3] =	{{A.c0.x, A.c0.y, A.c0.z},{A.c1.x, A.c1.y, A.c1.z},{A.c2.x, A.c2.y, A.c2.z}};
float det =   a[0][0] * a[1][1] * a[2][2] + a[0][1] * a[1][2] * a[2][0] + a[0][2] * a[1][0] * a[2][1]
- a[2][0] * a[1][1] * a[0][2] - a[2][1] * a[1][2] * a[0][0] - a[2][2] * a[1][0] * a[0][1];
if( det != 0.0 ){
result[0][0] = a[1][1] * a[2][2] - a[1][2] * a[2][1];result[0][1] = a[2][1] * a[0][2] - a[2][2] * a[0][1];
result[0][2] = a[0][1] * a[1][2] - a[0][2] * a[1][1];result[1][0] = a[2][0] * a[1][2] - a[1][0] * a[2][2];
result[1][1] = a[0][0] * a[2][2] - a[2][0] * a[0][2];result[1][2] = a[1][0] * a[0][2] - a[0][0] * a[1][2];
result[2][0] = a[1][0] * a[2][1] - a[2][0] * a[1][1];result[2][1] = a[2][0] * a[0][1] - a[0][0] * a[2][1];
result[2][2] = a[0][0] * a[1][1] - a[1][0] * a[0][1];
R = make_mat3(make_float3(result[0][0], result[0][1], result[0][2]), 
make_float3(result[1][0], result[1][1], result[1][2]), make_float3(result[2][0], result[2][1], result[2][2]));
return mult_f_f33( 1.0f / det, R);}
R = make_mat3(make_float3(1.0f, 0.0f, 0.0f), 
make_float3(0.0f, 1.0f, 0.0f), make_float3(0.0f, 0.0f, 1.0f));
return R;
}

__DEVICE__ mat3 transpose_f33 (mat3 A) {
float r[3][3];
float a[3][3] =	{{A.c0.x, A.c0.y, A.c0.z},{A.c1.x, A.c1.y, A.c1.z},{A.c2.x, A.c2.y, A.c2.z}};
for( int i = 0; i < 3; ++i){
for( int j = 0; j < 3; ++j){
r[i][j] = a[j][i];}}
mat3 R = make_mat3(make_float3(r[0][0], r[0][1], r[0][2]), 
make_float3(r[1][0], r[1][1], r[1][2]), make_float3(r[2][0], r[2][1], r[2][2]));
return R;
}

__DEVICE__ float Y_2_linCV( float Y, float Ymax, float Ymin) {
return (Y - Ymin) / (Ymax - Ymin);
}

__DEVICE__ float linCV_2_Y( float linCV, float Ymax, float Ymin) {
return linCV * (Ymax - Ymin) + Ymin;
}

__DEVICE__ float3 Y_2_linCV_f3( float3 Y, float Ymax, float Ymin) {
float3 linCV;
linCV.x = Y_2_linCV( Y.x, Ymax, Ymin);linCV.y = Y_2_linCV( Y.y, Ymax, Ymin);linCV.z = Y_2_linCV( Y.z, Ymax, Ymin);
return linCV;
}

__DEVICE__ float3 linCV_2_Y_f3( float3 linCV, float Ymax, float Ymin) {
float3 Y;
Y.x = linCV_2_Y( linCV.x, Ymax, Ymin);Y.y = linCV_2_Y( linCV.y, Ymax, Ymin);Y.z = linCV_2_Y( linCV.z, Ymax, Ymin);
return Y;
}

__DEVICE__ mat3 RGBtoXYZ( Chromaticities N) {
mat3 M = make_mat3(make_float3(N.red.x, N.red.y, 1.0f - (N.red.x + N.red.y)),
make_float3(N.green.x, N.green.y, 1.0f - (N.green.x + N.green.y)),
make_float3(N.blue.x, N.blue.y, 1.0f - (N.blue.x + N.blue.y)));
float3 wh = make_float3(N.white.x / N.white.y, 1.0f, (1.0f - (N.white.x + N.white.y)) / N.white.y);
wh = mult_f3_f33(wh, invert_f33(M));
mat3 WH =  make_mat3(make_float3(wh.x, 0.0f, 0.0f), 
make_float3(0.0f, wh.y, 0.0f), make_float3(0.0f, 0.0f, wh.z));
M = mult_f33_f33(WH, M);
return M;
}

__DEVICE__ mat4 RGBtoXYZ_m4( Chromaticities N, float W) {
mat3 A = RGBtoXYZ(N);
mat4 M = make_mat4_f33_f(A, W);
return M;
}

__DEVICE__ mat3 XYZtoRGB( Chromaticities N) {
mat3 M = invert_f33(RGBtoXYZ(N));
return M;
}

__DEVICE__ mat4 XYZtoRGB_m4( Chromaticities N, float W) {
mat3 A = XYZtoRGB(N);
mat4 M = make_mat4_f33_f(A, W);
return M;
}

__DEVICE__ float3 XYZ_2_xyY( float3 XYZ) {  
float3 xyY;
float divisor = (XYZ.x + XYZ.y + XYZ.z);
if (divisor == 0.0f) divisor = 1e-10f;
xyY.x = XYZ.x / divisor;
xyY.y = XYZ.y / divisor;  
xyY.z = XYZ.y;
return xyY;
}

__DEVICE__ float3 xyY_2_XYZ( float3 xyY) {
float3 XYZ;
XYZ.x = xyY.x * xyY.z / _fmaxf( xyY.y, 1e-10f);
XYZ.y = xyY.z;  
XYZ.z = (1.0f - xyY.x - xyY.y) * xyY.z / _fmaxf( xyY.y, 1e-10f);
return XYZ;
}

__DEVICE__ mat3 calculate_cat_matrix ( float2 src_xy, float2 des_xy) {
mat3 coneRespMat = CONE_RESP_MAT_BRADFORD;
const float3 src_xyY = { src_xy.x, src_xy.y, 1.0f };
const float3 des_xyY = { des_xy.x, des_xy.y, 1.0f };
float3 src_XYZ = xyY_2_XYZ( src_xyY );
float3 des_XYZ = xyY_2_XYZ( des_xyY );
float3 src_coneResp = mult_f3_f33( src_XYZ, coneRespMat);
float3 des_coneResp = mult_f3_f33( des_XYZ, coneRespMat);
mat3 vkMat = { { des_coneResp.x / src_coneResp.x, 0.0f, 0.0f },
{ 0.0f, des_coneResp.y / src_coneResp.y, 0.0f },{ 0.0f, 0.0f, des_coneResp.z / src_coneResp.z }};
mat3 cat_matrix = mult_f33_f33( coneRespMat, mult_f33_f33( vkMat, invert_f33( coneRespMat ) ) );
return cat_matrix;
}

__DEVICE__ mat3 calc_sat_adjust_matrix ( float sat, float3 rgb2Y) {
float M[3][3];
M[0][0] = (1.0f - sat) * rgb2Y.x + sat;
M[1][0] = (1.0f - sat) * rgb2Y.x;
M[2][0] = (1.0f - sat) * rgb2Y.x;
M[0][1] = (1.0f - sat) * rgb2Y.y;
M[1][1] = (1.0f - sat) * rgb2Y.y + sat;
M[2][1] = (1.0f - sat) * rgb2Y.y;
M[0][2] = (1.0f - sat) * rgb2Y.z;
M[1][2] = (1.0f - sat) * rgb2Y.z;
M[2][2] = (1.0f - sat) * rgb2Y.z + sat;
mat3 R = make_mat3(make_float3(M[0][0], M[0][1], M[0][2]), 
make_float3(M[1][0], M[1][1], M[1][2]), make_float3(M[2][0], M[2][1], M[2][2]));
R = transpose_f33(R);    
return R;
}

__DEVICE__ float moncurve_f( float x, float gamma, float offs ) {
float y;
const float fs = (( gamma - 1.0f) / offs) * _powf( offs * gamma / ( ( gamma - 1.0f) * ( 1.0f + offs)), gamma);
const float xb = offs / ( gamma - 1.0f);
if ( x >= xb) y = _powf( ( x + offs) / ( 1.0f + offs), gamma);
else
y = x * fs;
return y;
}

__DEVICE__ float3 moncurve_f_f3( float3 x, float gamma, float offs) {
float3 y;
y.x = moncurve_f( x.x, gamma, offs);y.y = moncurve_f( x.y, gamma, offs);y.z = moncurve_f( x.z, gamma, offs);
return y;
}

__DEVICE__ float clog_to_lin( float clog) {
float out;
if(clog < 0.12512248f)
out = -( _powf( 10.0f, ( 0.12512248f - clog ) / 0.45310179f ) - 1.0f ) / 10.1596f;
else
out = ( _powf( 10.0f, ( clog - 0.12512248f ) / 0.45310179f ) - 1.0f ) / 10.1596f;
return out;
}

__DEVICE__ float lin_to_clog( float lin) {
float clog;
clog = lin < clog_to_lin(0.0730597f) ? -(0.529136f * (_log10f(-lin * 10.1596f + 1.0f)) - 0.0730597f) : 0.529136f * _log10f(10.1596f * lin + 1.0f) + 0.0730597f;
return clog;
}

__DEVICE__ float3 lin_to_clog_f3( float3 lin) {
lin.x = lin_to_clog(lin.x);lin.y = lin_to_clog(lin.y);lin.z = lin_to_clog(lin.z);
return lin;
}

__DEVICE__ SegmentedSplineParams_c5 RRT_PARAMS() {
SegmentedSplineParams_c5 A = {{ -4.0f, -4.0f, -3.1573765773f, -0.4852499958f, 1.8477324706f, 1.8477324706f}, 
{ -0.7185482425f, 2.0810307172f, 3.6681241237f, 4.0f, 4.0f, 4.0f}, {0.18f * _powf(2.0f, -15.0f), 0.0001f}, 
{0.18f, 4.8f}, {0.18f * _powf(2.0f, 18.0f), 10000.0f}, 0.0f, 0.0f};
return A;
}

__DEVICE__ float segmented_spline_c5_fwd( float x) {
SegmentedSplineParams_c5 C = RRT_PARAMS();
const int N_KNOTS_LOW = 4;
const int N_KNOTS_HIGH = 4;
float logx = _log10f( _fmaxf(x, 0.0f )); float logy;
if ( logx <= _log10f(C.minPoint.x) ) { 
logy = logx * C.slopeLow + ( _log10f(C.minPoint.y) - C.slopeLow * _log10f(C.minPoint.x) );
} else if (( logx > _log10f(C.minPoint.x) ) && ( logx < _log10f(C.midPoint.x) )) {
float knot_coord = (N_KNOTS_LOW-1) * (logx - _log10f(C.minPoint.x))/(_log10f(C.midPoint.x) - _log10f(C.minPoint.x));
int j = knot_coord;
float t = knot_coord - j;
float3 cf = make_float3( C.coefsLow[ j], C.coefsLow[ j + 1], C.coefsLow[ j + 2]);
float3 monomials = make_float3( t * t, t, 1.0f );
logy = dot_f3_f3( monomials, mult_f3_f33( cf, MM));
} else if (( logx >= _log10f(C.midPoint.x) ) && ( logx < _log10f(C.maxPoint.x) )) {
float knot_coord = (N_KNOTS_HIGH-1) * (logx-log10(C.midPoint.x))/(_log10f(C.maxPoint.x)-log10(C.midPoint.x));
int j = knot_coord;
float t = knot_coord - j;
float3 cf = { C.coefsHigh[ j], C.coefsHigh[ j + 1], C.coefsHigh[ j + 2]}; 
float3 monomials = make_float3( t * t, t, 1.0f );logy = dot_f3_f3( monomials, mult_f3_f33( cf, MM));
} else {
logy = logx * C.slopeHigh + ( _log10f(C.maxPoint.y) - C.slopeHigh * _log10f(C.maxPoint.x) );}
return _pow10f(logy);
}

__DEVICE__ float segmented_spline_c5_rev ( float y) {  
SegmentedSplineParams_c5 C = RRT_PARAMS();
const int N_KNOTS_LOW = 4;
const int N_KNOTS_HIGH = 4;
const float KNOT_INC_LOW = (_log10f(C.midPoint.x) - _log10f(C.minPoint.x)) / (N_KNOTS_LOW - 1.0f);
const float KNOT_INC_HIGH = (_log10f(C.maxPoint.x) - _log10f(C.midPoint.x)) / (N_KNOTS_HIGH - 1.0f);
float KNOT_Y_LOW[ N_KNOTS_LOW];
for (int i = 0; i < N_KNOTS_LOW; i = i+1) {
KNOT_Y_LOW[ i] = ( C.coefsLow[i] + C.coefsLow[i+1]) / 2.0f;};
float KNOT_Y_HIGH[ N_KNOTS_HIGH];
for (int i = 0; i < N_KNOTS_HIGH; i = i+1) {
KNOT_Y_HIGH[ i] = ( C.coefsHigh[i] + C.coefsHigh[i+1]) / 2.0f;};
float logy = _log10f( _fmaxf(y, 1e-10f));float logx;
if (logy <= _log10f(C.minPoint.y)) {
logx = _log10f(C.minPoint.x);
} else if ( (logy > _log10f(C.minPoint.y)) && (logy <= _log10f(C.midPoint.y)) ) {
unsigned int j;float3 cf;
if ( logy > KNOT_Y_LOW[ 0] && logy <= KNOT_Y_LOW[ 1]) {
cf.x = C.coefsLow[0];  cf.y = C.coefsLow[1];  cf.z = C.coefsLow[2];  j = 0;
} else if ( logy > KNOT_Y_LOW[ 1] && logy <= KNOT_Y_LOW[ 2]) {
cf.x = C.coefsLow[1];  cf.y = C.coefsLow[2];  cf.z = C.coefsLow[3];  j = 1;
} else if ( logy > KNOT_Y_LOW[ 2] && logy <= KNOT_Y_LOW[ 3]) {
cf.x = C.coefsLow[2];  cf.y = C.coefsLow[3];  cf.z = C.coefsLow[4];  j = 2;} 
const float3 tmp = mult_f3_f33( cf, MM);
float a = tmp.x;float b = tmp.y;float c = tmp.z;c = c - logy;
const float d = _sqrtf( b * b - 4.0f * a * c);
const float t = ( 2.0f * c) / ( -d - b);
logx = _log10f(C.minPoint.x) + ( t + j) * KNOT_INC_LOW;
} else if ( (logy > _log10f(C.midPoint.y)) && (logy < _log10f(C.maxPoint.y)) ) {
unsigned int j;float3 cf;
if ( logy > KNOT_Y_HIGH[ 0] && logy <= KNOT_Y_HIGH[ 1]) {
cf.x = C.coefsHigh[0];  cf.y = C.coefsHigh[1];  cf.z = C.coefsHigh[2];  j = 0;
} else if ( logy > KNOT_Y_HIGH[ 1] && logy <= KNOT_Y_HIGH[ 2]) {
cf.x = C.coefsHigh[1];  cf.y = C.coefsHigh[2];  cf.z = C.coefsHigh[3];  j = 1;
} else if ( logy > KNOT_Y_HIGH[ 2] && logy <= KNOT_Y_HIGH[ 3]) {
cf.x = C.coefsHigh[2];  cf.y = C.coefsHigh[3];  cf.z = C.coefsHigh[4];  j = 2;} 
const float3 tmp = mult_f3_f33( cf, MM);
float a = tmp.x;float b = tmp.y;float c = tmp.z;c = c - logy;
const float d = _sqrtf( b * b - 4.0f * a * c);
const float t = ( 2.0f * c) / ( -d - b);
logx = _log10f(C.midPoint.x) + ( t + j) * KNOT_INC_HIGH;
} else {
logx = _log10f(C.maxPoint.x);}
return _pow10f( logx);
}

__DEVICE__ SegmentedSplineParams_c9 ODT_48nits() {
SegmentedSplineParams_c9 A =
{{ -1.6989700043f, -1.6989700043f, -1.4779f, -1.2291f, -0.8648f, -0.448f, 0.00518f, 0.4511080334f, 0.9113744414f, 0.9113744414f},
{ 0.5154386965f, 0.8470437783f, 1.1358f, 1.3802f, 1.5197f, 1.5985f, 1.6467f, 1.6746091357f, 1.6878733390f, 1.6878733390f },
{segmented_spline_c5_fwd( 0.18f * _powf(2.0f, -6.5f) ),  0.02f},{segmented_spline_c5_fwd( 0.18f ), 4.8f},
{segmented_spline_c5_fwd( 0.18f * _powf(2.0f, 6.5f) ), 48.0f}, 0.0f, 0.04f};
return A;
};

__DEVICE__ float segmented_spline_c9_rev ( float y) {  
SegmentedSplineParams_c9 C = ODT_48nits();
const int N_KNOTS_LOW = 8;
const int N_KNOTS_HIGH = 8;
const float KNOT_INC_LOW = (_log10f(C.midPoint.x) - _log10f(C.minPoint.x)) / (N_KNOTS_LOW - 1.0f);
const float KNOT_INC_HIGH = (_log10f(C.maxPoint.x) - _log10f(C.midPoint.x)) / (N_KNOTS_HIGH - 1.0f);
float KNOT_Y_LOW[ N_KNOTS_LOW];
for (int i = 0; i < N_KNOTS_LOW; i = i+1) {
KNOT_Y_LOW[ i] = ( C.coefsLow[i] + C.coefsLow[i+1]) / 2.0f;};
float KNOT_Y_HIGH[ N_KNOTS_HIGH];
for (int i = 0; i < N_KNOTS_HIGH; i = i+1) {
KNOT_Y_HIGH[ i] = ( C.coefsHigh[i] + C.coefsHigh[i+1]) / 2.0f;};
float logy = _log10f( _fmaxf( y, 1e-10f));
float logx;
if (logy <= _log10f(C.minPoint.y)) {
logx = _log10f(C.minPoint.x);
} else if ( (logy > _log10f(C.minPoint.y)) && (logy <= _log10f(C.midPoint.y)) ) {
unsigned int j;float3 cf;
if ( logy > KNOT_Y_LOW[ 0] && logy <= KNOT_Y_LOW[ 1]) {
cf.x = C.coefsLow[0];  cf.y = C.coefsLow[1];  cf.z = C.coefsLow[2];  j = 0;
} else if ( logy > KNOT_Y_LOW[ 1] && logy <= KNOT_Y_LOW[ 2]) {
cf.x = C.coefsLow[1];  cf.y = C.coefsLow[2];  cf.z = C.coefsLow[3];  j = 1;
} else if ( logy > KNOT_Y_LOW[ 2] && logy <= KNOT_Y_LOW[ 3]) {
cf.x = C.coefsLow[2];  cf.y = C.coefsLow[3];  cf.z = C.coefsLow[4];  j = 2;
} else if ( logy > KNOT_Y_LOW[ 3] && logy <= KNOT_Y_LOW[ 4]) {
cf.x = C.coefsLow[3];  cf.y = C.coefsLow[4];  cf.z = C.coefsLow[5];  j = 3;
} else if ( logy > KNOT_Y_LOW[ 4] && logy <= KNOT_Y_LOW[ 5]) {
cf.x = C.coefsLow[4];  cf.y = C.coefsLow[5];  cf.z = C.coefsLow[6];  j = 4;
} else if ( logy > KNOT_Y_LOW[ 5] && logy <= KNOT_Y_LOW[ 6]) {
cf.x = C.coefsLow[5];  cf.y = C.coefsLow[6];  cf.z = C.coefsLow[7];  j = 5;
} else if ( logy > KNOT_Y_LOW[ 6] && logy <= KNOT_Y_LOW[ 7]) {
cf.x = C.coefsLow[6];  cf.y = C.coefsLow[7];  cf.z = C.coefsLow[8];  j = 6;}
const float3 tmp = mult_f3_f33( cf, MM);
float a = tmp.x;float b = tmp.y;float c = tmp.z;c = c - logy;
const float d = _sqrtf( b * b - 4.0f * a * c);
const float t = ( 2.0f * c) / ( -d - b);
logx = _log10f(C.minPoint.x) + ( t + j) * KNOT_INC_LOW;
} else if ( (logy > _log10f(C.midPoint.y)) && (logy < _log10f(C.maxPoint.y)) ) {
unsigned int j;float3 cf;
if ( logy > KNOT_Y_HIGH[ 0] && logy <= KNOT_Y_HIGH[ 1]) {
cf.x = C.coefsHigh[0];  cf.y = C.coefsHigh[1];  cf.z = C.coefsHigh[2];  j = 0;
} else if ( logy > KNOT_Y_HIGH[ 1] && logy <= KNOT_Y_HIGH[ 2]) {
cf.x = C.coefsHigh[1];  cf.y = C.coefsHigh[2];  cf.z = C.coefsHigh[3];  j = 1;
} else if ( logy > KNOT_Y_HIGH[ 2] && logy <= KNOT_Y_HIGH[ 3]) {
cf.x = C.coefsHigh[2];  cf.y = C.coefsHigh[3];  cf.z = C.coefsHigh[4];  j = 2;
} else if ( logy > KNOT_Y_HIGH[ 3] && logy <= KNOT_Y_HIGH[ 4]) {
cf.x = C.coefsHigh[3];  cf.y = C.coefsHigh[4];  cf.z = C.coefsHigh[5];  j = 3;
} else if ( logy > KNOT_Y_HIGH[ 4] && logy <= KNOT_Y_HIGH[ 5]) {
cf.x = C.coefsHigh[4];  cf.y = C.coefsHigh[5];  cf.z = C.coefsHigh[6];  j = 4;
} else if ( logy > KNOT_Y_HIGH[ 5] && logy <= KNOT_Y_HIGH[ 6]) {
cf.x = C.coefsHigh[5];  cf.y = C.coefsHigh[6];  cf.z = C.coefsHigh[7];  j = 5;
} else if ( logy > KNOT_Y_HIGH[ 6] && logy <= KNOT_Y_HIGH[ 7]) {
cf.x = C.coefsHigh[6];  cf.y = C.coefsHigh[7];  cf.z = C.coefsHigh[8];  j = 6;}
const float3 tmp = mult_f3_f33( cf, MM);
float a = tmp.x;float b = tmp.y;float c = tmp.z;c = c - logy;
const float d = _sqrtf( b * b - 4.0f * a * c);
const float t = ( 2.0f * c) / ( -d - b);
logx = _log10f(C.midPoint.x) + ( t + j) * KNOT_INC_HIGH;
} else {
logx = _log10f(C.maxPoint.x);}
return _pow10f( logx);
}

__DEVICE__ float3 dimSurround_to_darkSurround( float3 linearCV) {
float3 XYZ = mult_f3_f44( linearCV, AP1_2_XYZ_MAT); 
float3 xyY = XYZ_2_xyY(XYZ);
xyY.z = _fmaxf( xyY.z, 0.0f);
xyY.z = _powf( xyY.z, 1.0f / DIM_SURROUND_GAMMA);
XYZ = xyY_2_XYZ(xyY);
return mult_f3_f44( XYZ, XYZ_2_AP1_MAT);
}

__DEVICE__ float3 segmented_spline_c5_rev_f3( float3 rgb) {
rgb.x = segmented_spline_c5_rev( rgb.x);rgb.y = segmented_spline_c5_rev( rgb.y);rgb.z = segmented_spline_c5_rev( rgb.z);
return rgb;
}

__DEVICE__ float3 segmented_spline_c5_fwd_f3( float3 rgb) {
rgb.x = segmented_spline_c5_fwd( rgb.x);rgb.y = segmented_spline_c5_fwd( rgb.y);rgb.z = segmented_spline_c5_fwd( rgb.z);
return rgb;
}

__DEVICE__ float3 segmented_spline_c9_rev_f3( float3 rgb) {
rgb.x = segmented_spline_c9_rev( rgb.x);rgb.y = segmented_spline_c9_rev( rgb.y);rgb.z = segmented_spline_c9_rev( rgb.z);
return rgb;
}

__DEVICE__ float fullRange_to_smpteRange( float rgbIn) {
float REFBLACK =  64.0f / 1023.0f;float REFWHITE = 940.0f / 1023.0f;
return rgbIn * (REFWHITE - REFBLACK) + REFBLACK;
}

__DEVICE__ float3 fullRange_to_smpteRange_f3( float3 rgbIn){
rgbIn.x = fullRange_to_smpteRange(rgbIn.x);rgbIn.y = fullRange_to_smpteRange(rgbIn.y);rgbIn.z = fullRange_to_smpteRange(rgbIn.z);
return rgbIn;
}

__DEVICE__ float3 InvODT_sRGB( float3 outputCV) {
const Chromaticities DISPLAY_PRI = REC709_PRI;
const mat4 DISPLAY_PRI_2_XYZ_MAT = RGBtoXYZ_m4(DISPLAY_PRI, 1.0f);
const float DISPGAMMA = 2.4f; 
const float OFFSET = 0.055f;
float3 linearCV;
linearCV = moncurve_f_f3( outputCV, DISPGAMMA, OFFSET);
float3 XYZ = mult_f3_f44( linearCV, DISPLAY_PRI_2_XYZ_MAT);
XYZ = mult_f3_f33( XYZ, invert_f33( D60_2_D65_CAT));
linearCV = mult_f3_f44( XYZ, XYZ_2_AP1_MAT);
linearCV = mult_f3_f33( linearCV, invert_f33( ODT_SAT_MAT));
//linearCV = dimSurround_to_darkSurround( linearCV);
float3 rgbPre = linCV_2_Y_f3( linearCV, CINEMA_WHITE, CINEMA_BLACK);
float3 rgbPost;
rgbPost = segmented_spline_c9_rev_f3( rgbPre);
float3 oces = mult_f3_f44( rgbPost, AP1_2_AP0_MAT);
return oces;
}

__DEVICE__ float3 inverseRRT( float3 oces) {
float3 rgbPost;
rgbPost = segmented_spline_c5_rev_f3( oces);
rgbPost = MAX(rgbPost, 0.0f);
return rgbPost;
}

__DEVICE__ float3 aces_to_clog( float3 aces) {
const Chromaticities P3_PRI = P3D60_PRI;
const mat4 XYZ_2_P3_MAT = XYZtoRGB_m4(P3_PRI, 1.0f);
float3 rgb = mult_f3_f44( aces, AP0_2_XYZ_MAT);
//rgb = mult_f3_f33( rgb, D60_2_D65_CAT);
rgb = mult_f3_f44( rgb, XYZ_2_P3_MAT);
rgb = lin_to_clog_f3( rgb);
rgb = fullRange_to_smpteRange_f3(rgb);
return rgb;
}

__DEVICE__ float3 transform(int p_Width, int p_Height, int p_X, int p_Y, float p_R, float p_G, float p_B)
{
float3 aces = make_float3(p_R, p_G, p_B);
aces = InvODT_sRGB(aces);
aces = inverseRRT(aces);
aces = aces_to_clog(aces);
return aces;
}