IDT_INV_REC709.dctl

// Inverse Rec.709 to ACES AP0 DCTL

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

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

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;

#define AP0_2_XYZ_MAT		RGBtoXYZ( AP0)
#define XYZ_2_AP0_MAT		XYZtoRGB( AP0)
#define AP1_2_XYZ_MAT		RGBtoXYZ( AP1)
#define XYZ_2_AP1_MAT		XYZtoRGB( AP1)
#define AP0_2_AP1_MAT		mult_f33_f33( AP0_2_XYZ_MAT, XYZ_2_AP1_MAT)
#define AP1_2_AP0_MAT		mult_f33_f33( 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 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		_exp10f(_log10f(0.02f))

__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__ 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} };

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

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

__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__ 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__ 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__ mat3 XYZtoRGB( Chromaticities N) {
mat3 M = invert_f33(RGBtoXYZ(N));
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 bt1886_f( float V, float gamma, float Lw, float Lb) {
float a = _powf( _powf( Lw, 1.0f/gamma) - _powf( Lb, 1.0f/gamma), gamma);
float b = _powf( Lb, 1.0f/gamma) / ( _powf( Lw, 1.0f/gamma) - _powf( Lb, 1.0f/gamma));
float L = a * _powf( _fmaxf( V + b, 0.0f), gamma);
return L;
}

__DEVICE__ float3 bt1886_f_f3( float3 V, float gamma, float Lw, float Lb) {
float3 L;
L.x = bt1886_f( V.x, gamma, Lw, Lb);L.y = bt1886_f( V.y, gamma, Lw, Lb);L.z = bt1886_f( V.z, gamma, Lw, Lb);
return L;
}

__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 * _exp2f(-15.0f), 0.0001f}, 
{0.18f, 4.8f}, {0.18f * _exp2f(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 _exp10f(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 _exp10f( 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 * _exp2f(-6.5f) ),  0.02f},{segmented_spline_c5_fwd( 0.18f ), 4.8f},
{segmented_spline_c5_fwd( 0.18f * _exp2f(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 _exp10f( logx);
}

__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__ float3 InvODT_Rec709( float3 outputCV) {
const Chromaticities DISPLAY_PRI = REC709_PRI;
const mat3 DISPLAY_PRI_2_XYZ_MAT = RGBtoXYZ(DISPLAY_PRI);
const float DISPGAMMA = 2.4f; 
float L_W = 1.0f;
float L_B = 0.0f;
float3 linearCV = bt1886_f_f3( outputCV, DISPGAMMA, L_W, L_B);
float3 XYZ = mult_f3_f33( linearCV, DISPLAY_PRI_2_XYZ_MAT);
XYZ = mult_f3_f33( XYZ, invert_f33( D60_2_D65_CAT));
linearCV = mult_f3_f33( XYZ, XYZ_2_AP1_MAT);
linearCV = mult_f3_f33( linearCV, invert_f33( ODT_SAT_MAT));
float3 rgbPre = linCV_2_Y_f3( linearCV, CINEMA_WHITE, CINEMA_BLACK);
float3 rgbPost;
rgbPost = segmented_spline_c9_rev_f3( rgbPre);
float3 oces = mult_f3_f33( rgbPost, AP1_2_AP0_MAT);
return oces;
}

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

__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_Rec709(aces);
aces = inverseRRT(aces);
return aces;
}