gronod/argyllcms
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 11 Wiki Activity
Files
v3.5.0-ICCery
argyllcms/icc/icc_util.c
T
gronod 45bdf51201 Initial commit of v3.5.0
2026-08-20 20:28:40 +01:00

4667 lines
115 KiB
C
Raw Permalink Blame History

/*
* International Color Consortium Format Library (icclib)
* Utility functions used by icclib and other code.
*
* Author: Graeme W. Gill
* Date: 2002/04/22
* Version: 3.0.0
*
* Copyright 1997 - 2022 Graeme W. Gill
*
* This material is licensed with an "MIT" free use license:-
* see the License4.txt file in this directory for licensing details.
*/
/*
* The distinction between "core" and "utility" functions is somewhat
* arbitrary. We do this just to make filesize a bit more managable.
* This file is #included in icc.c.
*/
/* This file is #included in icc.c */
/*
* TTBD:
*
*/
#ifndef M_PI
# define M_PI 3.14159265358979323846
#endif
/* Round a number to the same quantization as a S15Fixed16 */
static double round_S15Fixed16Number(double d) {
d = floor(d * 65536.0 + 0.5); /* Beware! (int)(d + 0.5) doesn't work for -ve nummbets ! */
d = d/65536.0;
return d;
}
/* Macro version */
#define RND_S15FIXED16(xxx) ((xxx) > 0.0 ? (int)((xxx) * 65536.0 + 0.5)/65536.0 \
: (int)((xxx) * 65536.0 - 0.5)/65536.0)
/* Clamp a 3 vector to be +ve */
void icmClamp3(double out[3], double in[3]) {
int i;
for (i = 0; i < 3; i++)
out[i] = in[i] < 0.0 ? 0.0 : in[i];
}
/* Invert (negate) a 3 vector */
void icmInv3(double out[3], double in[3]) {
int i;
for (i = 0; i < 3; i++)
out[i] = -in[i];
}
/* Add two 3 vectors */
void icmAdd3(double out[3], double in1[3], double in2[3]) {
out[0] = in1[0] + in2[0];
out[1] = in1[1] + in2[1];
out[2] = in1[2] + in2[2];
}
/* Subtract two 3 vectors, out = in1 - in2 */
void icmSub3(double out[3], double in1[3], double in2[3]) {
out[0] = in1[0] - in2[0];
out[1] = in1[1] - in2[1];
out[2] = in1[2] - in2[2];
}
/* Divide two 3 vectors, out = in1/in2 */
void icmDiv3(double out[3], double in1[3], double in2[3]) {
out[0] = in1[0]/in2[0];
out[1] = in1[1]/in2[1];
out[2] = in1[2]/in2[2];
}
/* Multiply two 3 vectors, out = in1 * in2 */
void icmMul3(double out[3], double in1[3], double in2[3]) {
out[0] = in1[0] * in2[0];
out[1] = in1[1] * in2[1];
out[2] = in1[2] * in2[2];
}
/* Take (signed) values to power */
void icmPow3(double out[3], double in[3], double p) {
int i;
for (i = 0; i < 3; i++) {
if (in[i] < 0.0)
out[i] = -pow(-in[i], p);
else
out[i] = pow(in[i], p);
}
}
/* Square values */
void icmSqr3(double out[3], double in[3]) {
int i;
for (i = 0; i < 3; i++)
out[i] = in[i] * in[i];
}
/* Square root of values */
void icmSqrt3(double out[3], double in[3]) {
int i;
for (i = 0; i < 3; i++)
out[i] = sqrt(in[i]);
}
/* Take absolute of a 3 vector */
void icmAbs3(double out[3], double in[3]) {
out[0] = fabs(in[0]);
out[1] = fabs(in[1]);
out[2] = fabs(in[2]);
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Dump a 3x3 matrix to a stream */
void icmDump3x3(FILE *fp, char *id, char *pfx, double a[3][3]) {
int i, j;
fprintf(fp, "%s%s[%d][%d]\n",pfx,id,3,3);
for (j = 0; j < 3; j++) {
fprintf(fp, "%s ",pfx);
for (i = 0; i < 3; i++)
fprintf(fp, "%f%s",a[j][i], i < (3-1) ? ", " : "");
fprintf(fp, "\n");
}
}
/* Set a 3x3 matrix to a value */
void icmSetVal3x3(double mat[3][3], double val) {
int i, j;
for (j = 0; j < 3; j++) {
for (i = 0; i < 3; i++) {
mat[j][i] = val;
}
}
}
/* Set a 3x3 matrix to unity */
void icmSetUnity3x3(double mat[3][3]) {
int i, j;
for (j = 0; j < 3; j++) {
for (i = 0; i < 3; i++) {
if (i == j)
mat[j][i] = 1.0;
else
mat[j][i] = 0.0;
}
}
}
/* Copy a 3x3 transform matrix */
void icmCpy3x3(double dst[3][3], double src[3][3]) {
int i, j;
for (j = 0; j < 3; j++) {
for (i = 0; i < 3; i++)
dst[j][i] = src[j][i];
}
}
/* Add one 3x3 to another */
/* dst = src1 + src2 */
void icmAdd3x3(double dst[3][3], double src1[3][3], double src2[3][3]) {
int i, j;
for (j = 0; j < 3; j++) {
for (i = 0; i < 3; i++)
dst[j][i] = src1[j][i] + src2[j][i];
}
}
/* Scale each element of a 3x3 transform matrix */
void icmScale3x3(double dst[3][3], double src[3][3], double scale) {
int i, j;
for (j = 0; j < 3; j++) {
for (i = 0; i < 3; i++)
dst[j][i] = src[j][i] * scale;
}
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/*
mat in out
[ ] [] []
[ ] * [] => []
[ ] [] []
So indexes are mat[out][in]
*/
/* Multiply 3 array by 3x3 transform matrix */
void icmMulBy3x3(double out[3], double mat[3][3], double in[3]) {
double tt[3];
tt[0] = mat[0][0] * in[0] + mat[0][1] * in[1] + mat[0][2] * in[2];
tt[1] = mat[1][0] * in[0] + mat[1][1] * in[1] + mat[1][2] * in[2];
tt[2] = mat[2][0] * in[0] + mat[2][1] * in[1] + mat[2][2] * in[2];
out[0] = tt[0];
out[1] = tt[1];
out[2] = tt[2];
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Tensor product. Multiply two 3 vectors to form a 3x3 matrix */
/* src1[] forms the colums, and src2[] forms the rows in the result */
void icmTensMul3(double dst[3][3], double src1[3], double src2[3]) {
int i, j;
for (j = 0; j < 3; j++) {
for (i = 0; i < 3; i++)
dst[j][i] = src1[j] * src2[i];
}
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Multiply one 3x3 with another */
/* dst = src * dst */
void icmMul3x3(double dst[3][3], double src[3][3]) {
int i, j, k;
double td[3][3]; /* Temporary dest */
for (j = 0; j < 3; j++) {
for (i = 0; i < 3; i++) {
double tt = 0.0;
for (k = 0; k < 3; k++)
tt += src[j][k] * dst[k][i];
td[j][i] = tt;
}
}
/* Copy result out */
for (j = 0; j < 3; j++)
for (i = 0; i < 3; i++)
dst[j][i] = td[j][i];
}
/* Multiply one 3x3 with another #2 */
/* dst = src1 * src2 */
void icmMul3x3_2(double dst[3][3], double src1[3][3], double src2[3][3]) {
int i, j, k;
double td[3][3]; /* Temporary dest */
for (j = 0; j < 3; j++) {
for (i = 0; i < 3; i++) {
double tt = 0.0;
for (k = 0; k < 3; k++)
tt += src1[j][k] * src2[k][i];
td[j][i] = tt;
}
}
/* Copy result out */
for (j = 0; j < 3; j++)
for (i = 0; i < 3; i++)
dst[j][i] = td[j][i];
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Matrix Inversion
by Richard Carling
from "Graphics Gems", Academic Press, 1990
*/
/*
* adjoint( original_matrix, inverse_matrix )
*
* calculate the adjoint of a 3x3 matrix
*
* Let a denote the minor determinant of matrix A obtained by
* ij
*
* deleting the ith row and jth column from A.
*
* i+j
* Let b = (-1) a
* ij ji
*
* The matrix B = (b ) is the adjoint of A
* ij
*/
#define det2x2(a, b, c, d) (a * d - b * c)
static void adjoint(
double out[3][3],
double in[3][3]
) {
double a1, a2, a3, b1, b2, b3, c1, c2, c3;
/* assign to individual variable names to aid */
/* selecting correct values */
a1 = in[0][0]; b1 = in[0][1]; c1 = in[0][2];
a2 = in[1][0]; b2 = in[1][1]; c2 = in[1][2];
a3 = in[2][0]; b3 = in[2][1]; c3 = in[2][2];
/* row column labeling reversed since we transpose rows & columns */
out[0][0] = det2x2(b2, b3, c2, c3);
out[1][0] = - det2x2(a2, a3, c2, c3);
out[2][0] = det2x2(a2, a3, b2, b3);
out[0][1] = - det2x2(b1, b3, c1, c3);
out[1][1] = det2x2(a1, a3, c1, c3);
out[2][1] = - det2x2(a1, a3, b1, b3);
out[0][2] = det2x2(b1, b2, c1, c2);
out[1][2] = - det2x2(a1, a2, c1, c2);
out[2][2] = det2x2(a1, a2, b1, b2);
}
/*
* double = icmDet3x3( a1, a2, a3, b1, b2, b3, c1, c2, c3 )
*
* calculate the determinant of a 3x3 matrix
* in the form
*
* | a1, b1, c1 |
* | a2, b2, c2 |
* | a3, b3, c3 |
*/
double icmDet3x3(double in[3][3]) {
double a1, a2, a3, b1, b2, b3, c1, c2, c3;
double ans;
a1 = in[0][0]; b1 = in[0][1]; c1 = in[0][2];
a2 = in[1][0]; b2 = in[1][1]; c2 = in[1][2];
a3 = in[2][0]; b3 = in[2][1]; c3 = in[2][2];
ans = a1 * det2x2(b2, b3, c2, c3)
- b1 * det2x2(a2, a3, c2, c3)
+ c1 * det2x2(a2, a3, b2, b3);
return ans;
}
/*
* inverse( original_matrix, inverse_matrix )
*
* calculate the inverse of a 3x3 matrix
*
* -1
* A = ___1__ adjoint A
* det A
*/
/* Return non-zero if not invertable */
int icmInverse3x3(
double out[3][3],
double in[3][3]
) {
int i, j;
double det;
/* calculate the 3x3 determinant
* if the determinant is zero,
* then the inverse matrix is not unique.
*/
det = icmDet3x3(in);
if (fabs(det) < ICM_SMALL_NUMBER)
return 1;
/* calculate the adjoint matrix */
adjoint(out, in);
/* scale the adjoint matrix to get the inverse */
for (i = 0; i < 3; i++)
for(j = 0; j < 3; j++)
out[i][j] /= det;
return 0;
}
/* Set a 2x2 matrix to unity */
void icmSetUnity2x2(double mat[2][2]) {
int i, j;
for (j = 0; j < 2; j++) {
for (i = 0; i < 2; i++) {
if (i == j)
mat[j][i] = 1.0;
else
mat[j][i] = 0.0;
}
}
}
/* Invert a 2x2 transform matrix. Return 1 if error. */
int icmInverse2x2(double out[2][2], double in[2][2]) {
double det = det2x2(in[0][0], in[0][1], in[1][0], in[1][1]);
if (fabs(det) < ICM_SMALL_NUMBER)
return 1;
det = 1.0/det;
out[0][0] = det * in[1][1];
out[0][1] = det * -in[0][1];
out[1][0] = det * -in[1][0];
out[1][1] = det * in[0][0];
return 0;
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Transpose a 3x3 matrix */
void icmTranspose3x3(double out[3][3], double in[3][3]) {
int i, j;
if (out != in) {
for (i = 0; i < 3; i++)
for (j = 0; j < 3; j++)
out[i][j] = in[j][i];
} else {
double tt[3][3];
for (i = 0; i < 3; i++)
for (j = 0; j < 3; j++)
tt[i][j] = in[j][i];
for (i = 0; i < 3; i++)
for (j = 0; j < 3; j++)
out[i][j] = tt[i][j];
}
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Compute the dot product of two 3 vectors */
double icmDot3(double in1[3], double in2[3]) {
return in1[0] * in2[0] + in1[1] * in2[1] + in1[2] * in2[2];
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Compute the cross product of two 3D vectors, out = in1 x in2 */
void icmCross3(double out[3], double in1[3], double in2[3]) {
double tt[3];
tt[0] = (in1[1] * in2[2]) - (in1[2] * in2[1]);
tt[1] = (in1[2] * in2[0]) - (in1[0] * in2[2]);
tt[2] = (in1[0] * in2[1]) - (in1[1] * in2[0]);
out[0] = tt[0];
out[1] = tt[1];
out[2] = tt[2];
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Compute the norm (length) squared of a 3 vector */
double icmNorm3sq(double in[3]) {
return in[0] * in[0] + in[1] * in[1] + in[2] * in[2];
}
/* Compute the norm (length) of a 3 vector */
double icmNorm3(double in[3]) {
return sqrt(in[0] * in[0] + in[1] * in[1] + in[2] * in[2]);
}
/* Scale a 3 vector by the given ratio */
void icmScale3(double out[3], double in[3], double rat) {
out[0] = in[0] * rat;
out[1] = in[1] * rat;
out[2] = in[2] * rat;
}
/* Scale a 3 vector by the given ratio and add it */
void icmScaleAdd3(double out[3], double in2[3], double in1[3], double rat) {
out[0] = in2[0] + in1[0] * rat;
out[1] = in2[1] + in1[1] * rat;
out[2] = in2[2] + in1[2] * rat;
}
/* Compute a blend between in0 and in1 */
void icmBlend3(double out[3], double in0[3], double in1[3], double bf) {
out[0] = (1.0 - bf) * in0[0] + bf * in1[0];
out[1] = (1.0 - bf) * in0[1] + bf * in1[1];
out[2] = (1.0 - bf) * in0[2] + bf * in1[2];
}
/* Clip a vector to the range 0.0 .. 1.0 */
void icmClip3(double out[3], double in[3]) {
int j;
for (j = 0; j < 3; j++) {
out[j] = in[j];
if (out[j] < 0.0)
out[j] = 0.0;
else if (out[j] > 1.0)
out[j] = 1.0;
}
}
/* Clip a vector to the range 0.0 .. 1.0 */
/* and retun nz if clipping occured */
int icmClip3sig(double out[3], double in[3]) {
int j;
int clip = 0;
for (j = 0; j < 3; j++) {
out[j] = in[j];
if (out[j] < 0.0) {
out[j] = 0.0;
clip = 1;
} else if (out[j] > 1.0) {
out[j] = 1.0;
clip = 1;
}
}
return clip;
}
/* Clip a vector to the range 0.0 .. 1.0 */
/* and return any clipping margine over the limits. */
double icmClip3marg(double out[3], double in[3]) {
int j;
double tt, marg = 0.0;
for (j = 0; j < 3; j++) {
out[j] = in[j];
if (out[j] < 0.0) {
tt = 0.0 - out[j];
out[j] = 0.0;
if (tt > marg)
marg = tt;
} else if (out[j] > 1.0) {
tt = out[j] - 1.0;
out[j] = 1.0;
if (tt > marg)
marg = tt;
}
}
return marg;
}
/* Normalise a 3 vector to the given length. Return nz if not normalisable */
int icmNormalize3(double out[3], double in[3], double len) {
double tt = sqrt(in[0] * in[0] + in[1] * in[1] + in[2] * in[2]);
if (tt < ICM_SMALL_NUMBER)
return 1;
tt = len/tt;
out[0] = in[0] * tt;
out[1] = in[1] * tt;
out[2] = in[2] * tt;
return 0;
}
/* Compute the norm (length) squared of a vector define by two points */
double icmNorm33sq(double in1[3], double in0[3]) {
int j;
double rv;
for (rv = 0.0, j = 0; j < 3; j++) {
double tt = in1[j] - in0[j];
rv += tt * tt;
}
return rv;
}
/* Compute the norm (length) of a vector define by two points */
double icmNorm33(double in1[3], double in0[3]) {
int j;
double rv;
for (rv = 0.0, j = 0; j < 3; j++) {
double tt = in1[j] - in0[j];
rv += tt * tt;
}
return sqrt(rv);
}
/* Scale a two point vector by the given ratio. in0[] is the origin */
void icmScale33(double out[3], double in1[3], double in0[3], double rat) {
out[0] = in0[0] + (in1[0] - in0[0]) * rat;
out[1] = in0[1] + (in1[1] - in0[1]) * rat;
out[2] = in0[2] + (in1[2] - in0[2]) * rat;
}
/* Normalise a vector from 0->1 to the given length. */
/* The new location of in1[] is returned in out[]. */
/* Return nz if not normalisable */
int icmNormalize33(double out[3], double in1[3], double in0[3], double len) {
int j;
double vl;
for (vl = 0.0, j = 0; j < 3; j++) {
double tt = in1[j] - in0[j];
vl += tt * tt;
}
vl = sqrt(vl);
if (vl < ICM_SMALL_NUMBER)
return 1;
vl = len/vl;
for (j = 0; j < 3; j++) {
out[j] = in0[j] + (in1[j] - in0[j]) * vl;
}
return 0;
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Given two 3D points, create a matrix that rotates */
/* and scales one onto the other about the origin 0,0,0. */
/* The maths is from page 52 of Tomas Moller and Eric Haines "Real-Time Rendering". */
/* s is source vector, t is target vector. */
/* Usage of icmRotMat: */
/* t[0] == mat[0][0] * s[0] + mat[0][1] * s[1] + mat[0][2] * s[2]; */
/* t[1] == mat[1][0] * s[0] + mat[1][1] * s[1] + mat[1][2] * s[2]; */
/* t[2] == mat[2][0] * s[0] + mat[2][1] * s[1] + mat[2][2] * s[2]; */
/* i.e. use icmMulBy3x3 */
void icmRotMat(double m[3][3], double s[3], double t[3]) {
double sl, tl, sn[3], tn[3];
double v[3]; /* Cross product */
double e; /* Dot product */
double h; /* 1-e/Cross product dot squared */
/* Normalise input vectors */
/* The rotation will be about 0,0,0 */
sl = sqrt(s[0] * s[0] + s[1] * s[1] + s[2] * s[2]);
tl = sqrt(t[0] * t[0] + t[1] * t[1] + t[2] * t[2]);
if (sl < 1e-12 || tl < 1e-12) { /* Hmm. Do nothing */
m[0][0] = 1.0;
m[0][1] = 0.0;
m[0][2] = 0.0;
m[1][0] = 0.0;
m[1][1] = 1.0;
m[1][2] = 0.0;
m[2][0] = 0.0;
m[2][1] = 0.0;
m[2][2] = 1.0;
return;
}
sn[0] = s[0]/sl; sn[1] = s[1]/sl; sn[2] = s[2]/sl;
tn[0] = t[0]/tl; tn[1] = t[1]/tl; tn[2] = t[2]/tl;
/* Compute the cross product */
v[0] = (sn[1] * tn[2]) - (sn[2] * tn[1]);
v[1] = (sn[2] * tn[0]) - (sn[0] * tn[2]);
v[2] = (sn[0] * tn[1]) - (sn[1] * tn[0]);
/* Compute the dot product */
e = sn[0] * tn[0] + sn[1] * tn[1] + sn[2] * tn[2];
/* Cross product dot squared */
h = v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
/* If the two input vectors are close to being parallel, */
/* then h will be close to zero. */
if (fabs(h) < 1e-12) {
/* Make sure scale is the correct sign */
if (s[0] * t[0] + s[1] * t[1] + s[2] * t[2] < 0.0)
tl = -tl;
m[0][0] = tl/sl;
m[0][1] = 0.0;
m[0][2] = 0.0;
m[1][0] = 0.0;
m[1][1] = tl/sl;
m[1][2] = 0.0;
m[2][0] = 0.0;
m[2][1] = 0.0;
m[2][2] = tl/sl;
} else {
/* 1-e/Cross product dot squared */
h = (1.0 - e) / h;
m[0][0] = tl/sl * (e + h * v[0] * v[0]);
m[0][1] = tl/sl * (h * v[0] * v[1] - v[2]);
m[0][2] = tl/sl * (h * v[0] * v[2] + v[1]);
m[1][0] = tl/sl * (h * v[0] * v[1] + v[2]);
m[1][1] = tl/sl * (e + h * v[1] * v[1]);
m[1][2] = tl/sl * (h * v[1] * v[2] - v[0]);
m[2][0] = tl/sl * (h * v[0] * v[2] - v[1]);
m[2][1] = tl/sl * (h * v[1] * v[2] + v[0]);
m[2][2] = tl/sl * (e + h * v[2] * v[2]);
}
#ifdef NEVER /* Check result */
{
double tt[3];
icmMulBy3x3(tt, m, s);
if (icmLabDEsq(t, tt) > 1e-4) {
printf("icmRotMat error t, is %f %f %f\n",tt[0],tt[1],tt[2]);
printf(" should be %f %f %f\n",t[0],t[1],t[2]);
}
}
#endif /* NEVER */
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/*
mat in out
[ ] [] []
[ ] * [] => []
[ ] [] []
[ ] [] []
*/
/* Multiply 4 array by 4x4 transform matrix */
void icmMulBy4x4(double out[4], double mat[4][4], double in[4]) {
double tt[4];
tt[0] = mat[0][0] * in[0] + mat[0][1] * in[1] + mat[0][2] * in[2] + mat[0][3] * in[3];
tt[1] = mat[1][0] * in[0] + mat[1][1] * in[1] + mat[1][2] * in[2] + mat[1][3] * in[3];
tt[2] = mat[2][0] * in[0] + mat[2][1] * in[1] + mat[2][2] * in[2] + mat[2][3] * in[3];
tt[3] = mat[3][0] * in[0] + mat[3][1] * in[1] + mat[3][2] * in[2] + mat[3][3] * in[3];
out[0] = tt[0];
out[1] = tt[1];
out[2] = tt[2];
out[3] = tt[3];
}
/* Transpose a 4x4 matrix */
void icmTranspose4x4(double out[4][4], double in[4][4]) {
int i, j;
if (out != in) {
for (i = 0; i < 4; i++)
for (j = 0; j < 4; j++)
out[i][j] = in[j][i];
} else {
double tt[4][4];
for (i = 0; i < 4; i++)
for (j = 0; j < 4; j++)
tt[i][j] = in[j][i];
for (i = 0; i < 4; i++)
for (j = 0; j < 4; j++)
out[i][j] = tt[i][j];
}
}
/* Clip a vector to the range 0.0 .. 1.0 */
/* and return any clipping margine */
double icmClip4marg(double out[4], double in[4]) {
int j;
double tt, marg = 0.0;
for (j = 0; j < 4; j++) {
out[j] = in[j];
if (out[j] < 0.0) {
tt = 0.0 - out[j];
out[j] = 0.0;
if (tt > marg)
marg = tt;
} else if (out[j] > 1.0) {
tt = out[j] - 1.0;
out[j] = 1.0;
if (tt > marg)
marg = tt;
}
}
return marg;
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Copy a 3x4 transform matrix */
void icmCpy3x4(double dst[3][4], double src[3][4]) {
int i, j;
for (j = 0; j < 3; j++) {
for (i = 0; i < 4; i++)
dst[j][i] = src[j][i];
}
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Multiply 3 array by 3x4 transform matrix */
void icmMul3By3x4(double out[3], double mat[3][4], double in[3]) {
double tt[3];
tt[0] = mat[0][0] * in[0] + mat[0][1] * in[1] + mat[0][2] * in[2] + mat[0][3];
tt[1] = mat[1][0] * in[0] + mat[1][1] * in[1] + mat[1][2] * in[2] + mat[1][3];
tt[2] = mat[2][0] * in[0] + mat[2][1] * in[1] + mat[2][2] * in[2] + mat[2][3];
out[0] = tt[0];
out[1] = tt[1];
out[2] = tt[2];
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Given two 3D vectors, create a matrix that translates, */
/* rotates and scales one onto the other. */
/* The maths is from page 52 of Tomas Moller and Eric Haines */
/* "Real-Time Rendering". */
/* s0 -> s1 is source vector, t0 -> t1 is target vector. */
/* Usage of icmRotMat: */
/* t[0] = mat[0][0] * s[0] + mat[0][1] * s[1] + mat[0][2] * s[2] + mat[0][3]; */
/* t[1] = mat[1][0] * s[0] + mat[1][1] * s[1] + mat[1][2] * s[2] + mat[1][3]; */
/* t[2] = mat[2][0] * s[0] + mat[2][1] * s[1] + mat[2][2] * s[2] + mat[2][3]; */
/* i.e. use icmMul3By3x4 */
void icmVecRotMat(double m[3][4], double s1[3], double s0[3], double t1[3], double t0[3]) {
int i, j;
double ss[3], tt[3], rr[3][3];
/* Create the rotation matrix: */
for (i = 0; i < 3; i++) {
ss[i] = s1[i] - s0[i];
tt[i] = t1[i] - t0[i];
}
icmRotMat(rr, ss, tt);
/* Create rotated version of s0: */
icmMulBy3x3(ss, rr, s0);
/* Create 4x4 matrix */
for (j = 0; j < 3; j++) {
for (i = 0; i < 4; i++) {
if (i < 3 && j < 3)
m[j][i] = rr[j][i];
else if (i == 3 && j < 3)
m[j][i] = t0[j] - ss[j];
else if (i == j)
m[j][i] = 1.0;
else
m[j][i] = 0.0;
}
}
#ifdef NEVER /* Check result */
{
double tt0[3], tt1[3];
icmMul3By3x4(tt0, m, s0);
if (icmLabDEsq(t0, tt0) > 1e-4) {
printf("icmVecRotMat error t0, is %f %f %f\n",tt0[0],tt0[1],tt0[2]);
printf(" should be %f %f %f\n",t0[0],t0[1],t0[2]);
}
icmMul3By3x4(tt1, m, s1);
if (icmLabDEsq(t1, tt1) > 1e-4) {
printf("icmVecRotMat error t1, is %f %f %f\n",tt1[0],tt1[1],tt1[2]);
printf(" should be %f %f %f\n",t1[0],t1[1],t1[2]);
}
}
#endif /* NEVER */
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Compute the 3D intersection of a vector and a plane */
/* return nz if there is no intersection */
int icmVecPlaneIsect(
double isect[3], /* return intersection point */
double pl_const, /* Plane equation constant */
double pl_norm[3], /* Plane normal vector */
double ve_1[3], /* Point on line */
double ve_0[3] /* Second point on line */
) {
double den; /* denominator */
double rv; /* Vector parameter value */
double vvec[3]; /* Vector vector */
double ival[3]; /* Intersection value */
/* Compute vector between the two points */
vvec[0] = ve_1[0] - ve_0[0];
vvec[1] = ve_1[1] - ve_0[1];
vvec[2] = ve_1[2] - ve_0[2];
/* Compute the denominator */
den = pl_norm[0] * vvec[0] + pl_norm[1] * vvec[1] + pl_norm[2] * vvec[2];
/* Too small to intersect ? */
if (fabs(den) < 1e-12) {
return 1;
}
/* Compute the parameterized intersction point */
rv = -(pl_norm[0] * ve_0[0] + pl_norm[1] * ve_0[1] + pl_norm[2] * ve_0[2] + pl_const)/den;
/* Compute the actual intersection point */
isect[0] = ve_0[0] + rv * vvec[0];
isect[1] = ve_0[1] + rv * vvec[1];
isect[2] = ve_0[2] + rv * vvec[2];
return 0;
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Compute the closest point on a line to a point. */
/* Return closest point and parameter value if not NULL. */
/* Return nz if the line length is zero */
int icmLinePointClosest(double cp[3], double *pa,
double la0[3], double la1[3], double pp[3]) {
double va[3], vp[3];
double val; /* Vector length squared */
double a; /* Parameter value */
icmSub3(va, la1, la0); /* Line vector */
val = icmNorm3sq(va); /* Vector length squared */
if (val < 1e-12)
return 1;
icmSub3(vp, pp, la0); /* Point vector to line base */
a = icmDot3(vp, va) / val; /* Normalised dist of point projected onto line */
if (cp != NULL)
icmBlend3(cp, la0, la1, a);
if (pa != NULL)
*pa = a;
return 0;
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Compute the closest points between two lines a and b. */
/* Return closest points and parameter values if not NULL. */
/* Return nz if they are paralel. */
/* The maths is from page 338 of Tomas Moller and Eric Haines "Real-Time Rendering". */
int icmLineLineClosest(double ca[3], double cb[3], double *pa, double * pb,
double la0[3], double la1[3], double lb0[3], double lb1[3]) {
double va[3], vb[3];
double vvab[3], vvabns; /* Cross product of va and vb and its norm squared */
double vba0[3]; /* lb0 - la0 */
double tt[3];
double a, b; /* Parameter values */
icmSub3(va, la1, la0);
icmSub3(vb, lb1, lb0);
icmCross3(vvab, va, vb);
vvabns = icmNorm3sq(vvab);
if (vvabns < 1e-12)
return 1;
icmSub3(vba0, lb0, la0);
icmCross3(tt, vba0, vb);
a = icmDot3(tt, vvab) / vvabns;
icmCross3(tt, vba0, va);
b = icmDot3(tt, vvab) / vvabns;
if (pa != NULL)
*pa = a;
if (pb != NULL)
*pb = b;
if (ca != NULL) {
ca[0] = la0[0] + a * va[0];
ca[1] = la0[1] + a * va[1];
ca[2] = la0[2] + a * va[2];
}
if (cb != NULL) {
cb[0] = lb0[0] + b * vb[0];
cb[1] = lb0[1] + b * vb[1];
cb[2] = lb0[2] + b * vb[2];
}
#ifdef NEVER /* Verify */
{
double vab[3]; /* Vector from ca to cb */
vab[0] = lb0[0] + b * vb[0] - la0[0] - a * va[0];
vab[1] = lb0[1] + b * vb[1] - la0[1] - a * va[1];
vab[2] = lb0[2] + b * vb[2] - la0[2] - a * va[2];
if (icmDot3(va, vab) > 1e-6
|| icmDot3(vb, vab) > 1e-6)
warning("icmLineLineClosest verify failed\n");
}
#endif
return 0;
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Given 3 3D points, compute a plane equation. */
/* The normal will be right handed given the order of the points */
/* The plane equation will be the 3 normal components and the constant. */
/* Return nz if any points are cooincident or co-linear */
int icmPlaneEqn3(double eq[4], double p0[3], double p1[3], double p2[3]) {
double ll, v1[3], v2[3];
/* Compute vectors along edges */
v2[0] = p1[0] - p0[0];
v2[1] = p1[1] - p0[1];
v2[2] = p1[2] - p0[2];
v1[0] = p2[0] - p0[0];
v1[1] = p2[1] - p0[1];
v1[2] = p2[2] - p0[2];
/* Compute cross products v1 x v2, which will be the normal */
eq[0] = v1[1] * v2[2] - v1[2] * v2[1];
eq[1] = v1[2] * v2[0] - v1[0] * v2[2];
eq[2] = v1[0] * v2[1] - v1[1] * v2[0];
/* Normalise the equation */
ll = sqrt(eq[0] * eq[0] + eq[1] * eq[1] + eq[2] * eq[2]);
if (ll < 1e-10) {
return 1;
}
eq[0] /= ll;
eq[1] /= ll;
eq[2] /= ll;
/* Compute the plane equation constant */
eq[3] = - (eq[0] * p0[0])
- (eq[1] * p0[1])
- (eq[2] * p0[2]);
return 0;
}
/* Given a 3D point and a plane equation, return the signed */
/* distance from the plane */
double icmPlaneDist3(double eq[4], double p[3]) {
double rv;
rv = eq[0] * p[0]
+ eq[1] * p[1]
+ eq[2] * p[2]
+ eq[3];
return rv;
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Compute the norm (length) of a vector define by two points */
double icmNorm22(double in1[2], double in0[2]) {
int j;
double rv;
for (rv = 0.0, j = 0; j < 2; j++) {
double tt = in1[j] - in0[j];
rv += tt * tt;
}
return sqrt(rv);
}
/* Compute the norm (length) squared of of a vector defined by two points */
double icmNorm22sq(double in1[2], double in0[2]) {
int j;
double rv;
for (rv = 0.0, j = 0; j < 2; j++) {
double tt = in1[j] - in0[j];
rv += tt * tt;
}
return rv;
}
/* Compute the dot product of two 2 vectors */
double icmDot2(double in1[2], double in2[2]) {
return in1[0] * in2[0] + in1[1] * in2[1];
}
/* Compute the dot product of two 2 vectors defined by 4 points */
/* 1->2 and 3->4 */
double icmDot22(double in1[2], double in2[2], double in3[2], double in4[2]) {
return (in2[0] - in1[0]) * (in4[0] - in3[0])
+ (in2[1] - in1[1]) * (in4[1] - in3[1]);
}
/* Normalise a 2 vector to the given length. Return nz if not normalisable */
int icmNormalize2(double out[2], double in[2], double len) {
double tt = sqrt(in[0] * in[0] + in[1] * in[1]);
if (tt < ICM_SMALL_NUMBER)
return 1;
tt = len/tt;
out[0] = in[0] * tt;
out[1] = in[1] * tt;
return 0;
}
/* Return an orthogonal vector */
/* (90 degree rotate = swap and negate X ) */
void icmOrthog2(double out[2], double in[2]) {
double tt;
tt = in[0];
out[0] = -in[1];
out[1] = tt;
}
/* Given 2 2D points, compute a plane equation (implicit line equation). */
/* The normal will be right handed given the order of the points */
/* The plane equation will be the 2 normal components and the constant. */
/* Return nz if any points are cooincident or co-linear */
int icmPlaneEqn2(double eq[3], double p0[2], double p1[2]) {
double ll, v1[3];
/* Compute vectors along edge */
v1[0] = p1[0] - p0[0];
v1[1] = p1[1] - p0[1];
/* Normal to vector */
eq[0] = v1[1];
eq[1] = -v1[0];
/* Normalise the equation */
ll = sqrt(eq[0] * eq[0] + eq[1] * eq[1]);
if (ll < 1e-10) {
return 1;
}
eq[0] /= ll;
eq[1] /= ll;
/* Compute the plane equation constant */
eq[2] = - (eq[0] * p0[0])
- (eq[1] * p0[1]);
return 0;
}
/* Given a 2D point and a plane equation (implicit line), return the signed */
/* distance from the plane. The distance will be +ve if the point */
/* is to the right of the plane formed by two points in order */
double icmPlaneDist2(double eq[3], double p[2]) {
double rv;
rv = eq[0] * p[0]
+ eq[1] * p[1]
+ eq[2];
return rv;
}
/* Return the closest point on an implicit line to a point. */
/* Also return the absolute distance */
double icmImpLinePointClosest2(double cp[2], double eq[3], double pp[2]) {
double q; /* Closest distance to line */
q = eq[0] * pp[0]
+ eq[1] * pp[1]
+ eq[2];
cp[0] = pp[0] - q * eq[0];
cp[1] = pp[1] - q * eq[1];
return fabs(q);
}
/* Return the point of intersection of two implicit lines . */
/* Return nz if there is no intersection (lines are parallel) */
int icmImpLineIntersect2(double res[2], double eq1[3], double eq2[3]) {
double num;
num = eq1[0] * eq2[1] - eq2[0] * eq1[1];
if (fabs(num) < 1e-10)
return 1;
res[0] = (eq2[2] * eq1[1] - eq1[2] * eq2[1])/num;
res[1] = (eq1[2] * eq2[0] - eq2[2] * eq1[0])/num;
return 0;
}
/* Compute the closest point on a line to a point. */
/* Return closest point and parameter value if not NULL. */
/* Return nz if the line length is zero */
int icmLinePointClosest2(double cp[2], double *pa,
double la0[2], double la1[2], double pp[2]) {
double va[2], vp[2];
double val; /* Vector length squared */
double a; /* Parameter value */
va[0] = la1[0] - la0[0]; /* Line vector */
va[1] = la1[1] - la0[1];
val = va[0] * va[0] + va[1] * va[1];
if (val < 1e-12)
return 1;
vp[0] = pp[0] - la0[0]; /* Point vector to line base */
vp[1] = pp[1] - la0[1];
a = (vp[0] * va[0] + vp[1] * va[1]) / val; /* Normalised dist of point projected onto line */
if (cp != NULL) {
cp[0] = (1.0 - a) * la0[0] + a * la1[0];
cp[1] = (1.0 - a) * la0[1] + a * la1[1];
}
if (pa != NULL)
*pa = a;
return 0;
}
/* Given two infinite 2D lines define by 4 points, compute the intersection. */
/* Return nz if there is no intersection (lines are parallel) */
int icmLineIntersect2(double res[2], double p1[2], double p2[2], double p3[2], double p4[2]) {
/* Compute by determinants */
double x1y2_y1x2 = p1[0] * p2[1] - p1[1] * p2[0];
double x3y4_y3x4 = p3[0] * p4[1] - p3[1] * p4[0];
double x1_x2 = p1[0] - p2[0];
double y1_y2 = p1[1] - p2[1];
double x3_x4 = p3[0] - p4[0];
double y3_y4 = p3[1] - p4[1];
double num; /* Numerator */
num = x1_x2 * y3_y4 - y1_y2 * x3_x4;
if (fabs(num) < 1e-10)
return 1;
res[0] = (x1y2_y1x2 * x3_x4 - x1_x2 * x3y4_y3x4)/num;
res[1] = (x1y2_y1x2 * y3_y4 - y1_y2 * x3y4_y3x4)/num;
return 0;
}
/* Given two finite 2D lines define by 4 points, compute their paramaterized intersection. */
/* res[] and/or aprm[] may be NULL. Param is prop. from p1 -> p2, p3 -> p4 */
/* Return 2 if there is no intersection (lines are parallel) */
/* Return 1 lines do not cross within their length */
int icmParmLineIntersect2(double res[2], double aprm[2], double p1[2], double p2[2], double p3[2], double p4[2]) {
double _prm[2];
double *prm = aprm != NULL ? aprm : _prm;
double x21 = p2[0] - p1[0];
double y21 = p2[1] - p1[1];
double x31 = p3[0] - p1[0];
double y31 = p3[1] - p1[1];
double x43 = p4[0] - p3[0];
double y43 = p4[1] - p3[1];
double num; /* Numerator */
num = x43 * y21 - x21 * y43;
if (fabs(num) < 1e-10)
return 2;
prm[0] = (x43 * y31 - x31 * y43)/num; /* Parameter of 1->2 */
prm[1] = (x21 * y31 - x31 * y21)/num; /* Parameter of 3->4 */
if (res != NULL) {
res[0] = x21 * prm[0] + p1[0];
res[1] = y21 * prm[0] + p1[1];
}
if (prm[0] < -1e-10 || prm[0] > (1.0 + 1e-10)
|| prm[1] < -1e-10 || prm[1] > (1.0 + 1e-10))
return 1;
return 0;
}
/* Compute a blend between in0 and in1 for bl 0..1 */
void icmBlend2(double out[2], double in0[2], double in1[2], double bf) {
out[0] = (1.0 - bf) * in0[0] + bf * in1[0];
out[1] = (1.0 - bf) * in0[1] + bf * in1[1];
}
/* Scale a 2 vector by the given ratio */
void icmScale2(double out[2], double in[2], double rat) {
out[0] = in[0] * rat;
out[1] = in[1] * rat;
}
/* Scale a 2 vector by the given ratio and add it */
void icmScaleAdd2(double out[3], double in2[2], double in1[2], double rat) {
out[0] = in2[0] + in1[0] * rat;
out[1] = in2[1] + in1[1] * rat;
}
/* Convert orientation 0 = right, 1 = down, 2 = left, 3 = right */
/* into rotation matrix */
static void icmOr2mat(double mat[2][2], int or) {
if (or == 0) {
mat[0][0] = 1.0;
mat[0][1] = 0.0;
mat[1][0] = 0.0;
mat[1][1] = 1.0;
} else if (or == 1) {
mat[0][0] = 0.0;
mat[0][1] = 1.0;
mat[1][0] = -1.0;
mat[1][1] = 0.0;
} else if (or == 2) {
mat[0][0] = -1.0;
mat[0][1] = 0.0;
mat[1][0] = 0.0;
mat[1][1] = -1.0;
} else {
mat[0][0] = 0.0;
mat[0][1] = -1.0;
mat[1][0] = 1.0;
mat[1][1] = 0.0;
}
}
/* Convert an angle in radians into rotation matrix (0 = right) */
void icmRad2mat(double mat[2][2], double rad) {
double sinv = sin(rad);
double cosv = cos(rad);
mat[0][0] = cosv;
mat[0][1] = -sinv;
mat[1][0] = sinv;
mat[1][1] = cosv;
}
/* Convert an angle in degrees (0 = right) into rotation matrix */
void icmDeg2mat(double mat[2][2], double a) {
double rad = a * 3.1415926535897932384626433832795/180.0;
icmRad2mat(mat, rad);
}
/* Convert a vector into a rotation matrix */
void icmVec2mat(double mat[2][2], double dx, double dy) {
double rad = atan2(dy, dx);
icmRad2mat(mat, rad);
}
/* Multiply 2 array by 2x2 transform matrix */
void icmMulBy2x2(double out[2], double mat[2][2], double in[2]) {
double tt[2];
tt[0] = mat[0][0] * in[0] + mat[0][1] * in[1];
tt[1] = mat[1][0] * in[0] + mat[1][1] * in[1];
out[0] = tt[0];
out[1] = tt[1];
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Set a vector to a single value */
void icmSetN(double *dst, double src, int len) {
int j;
for (j = 0; j < len; j++)
dst[j] = src;
}
/* Copy a vector */
#ifdef NEVER /* Using macro */
void icmCpyN(double *dst, double *src, int len) {
int j;
if (dst != src)
for (j = 0; j < len; j++)
dst[j] = src[j];
}
#endif /* NEVER */
/* Clip a vector to the range 0.0 .. 1.0 */
void icmClipN(double out[3], double in[3], unsigned int len) {
int j;
for (j = 0; j < len; j++) {
out[j] = in[j];
if (out[j] < 0.0)
out[j] = 0.0;
else if (out[j] > 1.0)
out[j] = 1.0;
}
}
/* Clip a vector to the range 0.0 .. 1.0 */
/* and return any clipping margine over the limit */
double icmClipNmarg(double *out, double *in, unsigned int len) {
int j;
double tt, marg = 0.0;
for (j = 0; j < len; j++) {
out[j] = in[j];
if (out[j] < 0.0) {
tt = 0.0 - out[j];
out[j] = 0.0;
if (tt > marg)
marg = tt;
} else if (out[j] > 1.0) {
tt = out[j] - 1.0;
out[j] = 1.0;
if (tt > marg)
marg = tt;
}
}
return marg;
}
/* Return the magnitude (norm) of the difference between two vectors */
double icmDiffN(double *s1, double *s2, int len) {
int i;
double rv = 0.0;
for (i = 0; i < len; i++) {
double tt = s1[i] - s2[i];
rv += tt * tt;
}
return sqrt(rv);
}
/* ----------------------------------------------- */
/* Pseudo - Hilbert count sequencer */
/* This multi-dimensional count sequence is a distributed */
/* Gray code sequence, with direction reversal on every */
/* alternate power of 2 scale. */
/* It is intended to aid cache coherence in multi-dimensional */
/* regular sampling. It approximates the Hilbert curve sequence. */
/* Initialise, returns total usable count */
unsigned
psh_initN(
psh *p, /* Pointer to structure to initialise */
int di, /* Dimensionality */
unsigned int res[], /* Size per coordinate */
int co[] /* Coordinates to initialise (May be NULL) */
) {
int e;
memset(p, 0, sizeof(psh));
p->di = di;
for (e = 0; e < di; e++)
p->res[e] = res[e];
p->xbits = 0;
/* Compute bits per axis and total bits */
for (p->tbits = e = 0; e < di; e++) {
for (p->bits[e] = 0; (1u << p->bits[e]) < p->res[e]; p->bits[e]++)
;
p->tbits += p->bits[e];
if (p->bits[e] > p->xbits)
p->xbits = p->bits[e];
}
if (p->tbits > 32) /* Hmm */
return 0;
/* Compute the total count mask */
p->tmask = ((((unsigned)1) << p->tbits)-1);
/* Compute usable count */
p->count = 1;
for (e = 0; e < di; e++)
p->count *= p->res[e];
p->ix = 0; /* Initial binary index */
if (co != NULL) {
for (e = 0; e < di; e++)
co[e] = 0;
}
return p->count;
}
/* Backwards compatible same res per axis init */
unsigned
psh_init(
psh *p, /* Pointer to structure to initialise */
int di, /* Dimensionality */
unsigned int res, /* Size */
int co[] /* Coordinates to initialise (May be NULL) */
) {
int e;
unsigned int mres[MAX_CHAN];
for (e = 0; e < di; e++)
mres[e] = res;
return psh_initN(p, di, mres, co);
}
/* Diagnostic init - psh returns the value in co[] and then terminates */
unsigned
psh_init_diag(
psh *p, /* Pointer to structure to initialise */
int di, /* Dimensionality */
int co[] /* Coordinates to start on */
) {
int e;
memset(p, 0, sizeof(psh));
p->di = di;
p->diag = 1;
return 1;
}
/* Reset the counter */
void
psh_reset(
psh *p /* Pointer to structure */
) {
p->ix = 0;
}
/* Increment pseudo-hilbert coordinates */
/* Return non-zero if count rolls over to 0 */
int
psh_inc(
psh *p, /* Pointer to structure */
int co[] /* Coordinates to return */
) {
int di = p->di;
int e;
if (p->diag)
return 1;
do {
unsigned int b;
int gix; /* Gray code index */
p->ix = (p->ix + 1) & p->tmask;
gix = p->ix ^ (p->ix >> 1); /* Convert to gray code index */
for (e = 0; e < di; e++)
co[e] = 0;
for (b = 0; b < p->xbits; b++) { /* Distribute bits */
if (b & 1) {
for (e = di-1; e >= 0; e--) { /* In reverse order */
if (b < p->bits[e]) {
co[e] |= (gix & 1) << b;
gix >>= 1;
}
}
} else {
for (e = 0; e < di; e++) { /* In normal order */
if (b < p->bits[e]) {
co[e] |= (gix & 1) << b;
gix >>= 1;
}
}
}
}
/* Convert from Gray to binary coordinates */
for (e = 0; e < di; e++) {
unsigned int sh, tv;
for(sh = 1, tv = co[e];; sh <<= 1) {
unsigned ptv = tv;
tv ^= (tv >> sh);
if (ptv <= 1 || sh == 16)
break;
}
if (tv >= p->res[e]) /* Dumbo filter - increment again if outside cube range */
break;
co[e] = tv;
}
} while (e < di);
return (p->ix == 0);
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* CIE Y (range 0 .. 1) to perceptual CIE 1976 L* (range 0 .. 100) */
double
icmY2L(double val) {
if (val > 0.008856451586)
val = pow(val,1.0/3.0);
else
val = 7.787036979 * val + 16.0/116.0;
val = (116.0 * val - 16.0);
return val;
}
/* Perceptual CIE 1976 L* (range 0 .. 100) to CIE Y (range 0 .. 1) */
double
icmL2Y(double val) {
val = (val + 16.0)/116.0;
if (val > 24.0/116.0)
val = pow(val,3.0);
else
val = (val - 16.0/116.0)/7.787036979;
return val;
}
/* CIE XYZ (range 0 .. 1) to perceptual CIE 1976 L*a*b* */
void
icmXYZ2Lab(icmXYZNumber *w, double *out, double *in) {
double X = in[0], Y = in[1], Z = in[2];
double x,y,z,fx,fy,fz;
x = X/w->X;
y = Y/w->Y;
z = Z/w->Z;
if (x > 0.008856451586)
fx = pow(x,1.0/3.0);
else
fx = 7.787036979 * x + 16.0/116.0;
if (y > 0.008856451586)
fy = pow(y,1.0/3.0);
else
fy = 7.787036979 * y + 16.0/116.0;
if (z > 0.008856451586)
fz = pow(z,1.0/3.0);
else
fz = 7.787036979 * z + 16.0/116.0;
out[0] = 116.0 * fy - 16.0;
out[1] = 500.0 * (fx - fy);
out[2] = 200.0 * (fy - fz);
}
/* Perceptual CIE 1976 L*a*b* to CIE XYZ (range 0 .. 1) */
void
icmLab2XYZ(icmXYZNumber *w, double *out, double *in) {
double L = in[0], a = in[1], b = in[2];
double x,y,z,fx,fy,fz;
fy = (L + 16.0)/116.0;
fx = a/500.0 + fy;
fz = fy - b/200.0;
if (fy > 24.0/116.0)
y = pow(fy,3.0);
else
y = (fy - 16.0/116.0)/7.787036979;
if (fx > 24.0/116.0)
x = pow(fx,3.0);
else
x = (fx - 16.0/116.0)/7.787036979;
if (fz > 24.0/116.0)
z = pow(fz,3.0);
else
z = (fz - 16.0/116.0)/7.787036979;
out[0] = x * w->X;
out[1] = y * w->Y;
out[2] = z * w->Z;
}
/*
* This is a modern update to L*a*b*, based on IPT space.
*
* Differences to L*a*b* and IPT:
* Using inverse CIE 2012 2 degree LMS to XYZ matrix instead of Hunt-Pointer-Estevez.
* Von Kries chromatic adapation in LMS cone space.
* Using L* compression rather than IPT pure 0.43 power.
* Tweaked LMS' to IPT matrix to account for change in XYZ to LMS matrix.
* Output scaled to L*a*b* type ranges, to maintain 1 JND scale.
* (Watch out - L* value is not a non-linear Y value though!
* - Interesting that Dolby force L to be just dependent on Y
* by making L = 0.5 L | 0.5 M in ICtCp space).
*/
/* CIE XYZ to perceptual Lpt */
void
icmXYZ2Lpt(icmXYZNumber *w, double *out, double *in) {
double wxyz[3];
double wlms[3];
double lms[3];
double xyz2lms[3][3] = {
{ 0.2052445519046028, 0.8334486497310412, -0.0386932016356441 },
{ -0.4972221301804286, 1.4034846060306130, 0.0937375241498157 },
{ 0.0000000000000000, 0.0000000000000000, 1.0000000000000000 }
};
double lms2ipt[3][3] = {
{ 0.6585034777870502, 0.1424555300344579, 0.1990409921784920 },
{ 5.6413505933276049, -6.1697985811414187, 0.5284479878138138 },
{ 1.6370552576322106, 0.0192823194340315, -1.6563375770662419 }
};
int j;
/* White point in Cone space */
wxyz[0] = w->X;
wxyz[1] = w->Y;
wxyz[2] = w->Z;
icmMulBy3x3(wlms, xyz2lms, wxyz);
/* Incoming XYZ to Cone space */
icmMulBy3x3(lms, xyz2lms, in);
for (j = 0; j < 3; j++) {
/* Von Kries chromatic adapation */
lms[j] /= wlms[j];
/* Non-linearity */
if (lms[j] > 0.008856451586)
lms[j] = pow(lms[j],1.0/3.0);
else
lms[j] = 7.787036979 * lms[j] + 16.0/116.0;
lms[j] = 116.0 * lms[j] - 16.0;
}
/* IPT */
icmMulBy3x3(out, lms2ipt, lms);
}
void
icmLpt2XYZ(icmXYZNumber *w, double *out, double *in) {
double wxyz[3];
double wlms[3];
double lms[3];
double xyz2lms[3][3] = {
{ 0.2052445519046028, 0.8334486497310412, -0.0386932016356441 },
{ -0.4972221301804286, 1.4034846060306130, 0.0937375241498157 },
{ 0.0000000000000000, 0.0000000000000000, 1.0000000000000000 }
};
double ipt2lms[3][3] = {
{ 1.0000000000000000, 0.0234881527511557, 0.1276631419615779 },
{ 1.0000000000000000, -0.1387534648407132, 0.0759005921388901 },
{ 1.0000000000000000, 0.0215994105411036, -0.4766811148374502 }
};
double lms2xyz[3][3] = {
{ 1.9979376130193824, -1.1864600428553205, 0.1885224298359384 },
{ 0.7078230795296872, 0.2921769204703129, -0.0000000000000000 },
{ 0.0000000000000000, 0.0000000000000000, 1.0000000000000000 }
};
int j;
wxyz[0] = w->X;
wxyz[1] = w->Y;
wxyz[2] = w->Z;
icmMulBy3x3(wlms, xyz2lms, wxyz);
icmMulBy3x3(lms, ipt2lms, in);
for (j = 0; j < 3; j++) {
lms[j] = (lms[j] + 16.0)/116.0;
if (lms[j] > 24.0/116.0)
lms[j] = pow(lms[j], 3.0);
else
lms[j] = (lms[j] - 16.0/116.0)/7.787036979;
lms[j] *= wlms[j];
}
icmMulBy3x3(out, lms2xyz, lms);
}
/* LCh to Lab (general to polar, works with Lpt, Luv too) */
void icmLCh2Lab(double *out, double *in) {
double C, h;
C = in[1];
h = M_PI/180.0 * in[2];
out[0] = in[0];
out[1] = C * cos(h);
out[2] = C * sin(h);
}
/* Lab to LCh (general to polar, works with Lpt, Luv too) */
void icmLab2LCh(double *out, double *in) {
double C, h;
C = sqrt(in[1] * in[1] + in[2] * in[2]);
h = (180.0/M_PI) * atan2(in[2], in[1]);
h = (h < 0.0) ? h + 360.0 : h;
out[0] = in[0];
out[1] = C;
out[2] = h;
}
/* CIE XYZ to perceptual CIE 1976 L*u*v* */
extern ICCLIB_API void icmXYZ2Luv(icmXYZNumber *w, double *out, double *in) {
double X = in[0], Y = in[1], Z = in[2];
double un, vn, u, v, fl, fu, fv;
un = (4.0 * w->X) / (w->X + 15.0 * w->Y + 3.0 * w->Z);
vn = (9.0 * w->Y) / (w->X + 15.0 * w->Y + 3.0 * w->Z);
u = (4.0 * X) / (X + 15.0 * Y + 3.0 * Z);
v = (9.0 * Y) / (X + 15.0 * Y + 3.0 * Z);
Y /= w->Y;
if (Y > 0.008856451586)
fl = pow(Y,1.0/3.0);
else
fl = 7.787036979 * Y + 16.0/116.0;
fu = u - un;
fv = v - vn;
out[0] = 116.0 * fl - 16.0;
out[1] = 13.0 * out[0] * fu;
out[2] = 13.0 * out[0] * fv;
}
/* Perceptual CIE 1976 L*u*v* to CIE XYZ */
extern ICCLIB_API void icmLuv2XYZ(icmXYZNumber *w, double *out, double *in) {
double un, vn, u, v, fl, fu, fv, sum, X, Y, Z;
fl = (in[0] + 16.0)/116.0;
fu = in[1] / (13.0 * in[0]);
fv = in[2] / (13.0 * in[0]);
un = (4.0 * w->X) / (w->X + 15.0 * w->Y + 3.0 * w->Z);
vn = (9.0 * w->Y) / (w->X + 15.0 * w->Y + 3.0 * w->Z);
u = fu + un;
v = fv + vn;
if (fl > 24.0/116.0)
Y = pow(fl,3.0);
else
Y = (fl - 16.0/116.0)/7.787036979;
Y *= w->Y;
sum = (9.0 * Y)/v;
X = (u * sum)/4.0;
Z = (sum - X - 15.0 * Y)/3.0;
out[0] = X;
out[1] = Y;
out[2] = Z;
}
/* XYZ to Yxy */
extern ICCLIB_API void icmXYZ2Yxy(double *out, double *in) {
double sum = in[0] + in[1] + in[2];
double Y, x, y;
if (sum < 1e-9) {
Y = 0.0;
y = 1.0/3.0;
x = 1.0/3.0;
} else {
Y = in[1];
x = in[0]/sum;
y = in[1]/sum;
}
out[0] = Y;
out[1] = x;
out[2] = y;
}
/* XYZ to xy */
extern ICCLIB_API void icmXYZ2xy(double *out, double *in) {
double sum = in[0] + in[1] + in[2];
double x, y;
if (sum < 1e-9) {
y = 1.0/3.0;
x = 1.0/3.0;
} else {
x = in[0]/sum;
y = in[1]/sum;
}
out[0] = x;
out[1] = y;
}
/* Yxy to XYZ */
extern ICCLIB_API void icmYxy2XYZ(double *out, double *in) {
double Y = in[0];
double x = in[1];
double y = in[2];
double z = 1.0 - x - y;
double sum;
if (y < 1e-9) {
out[0] = out[1] = out[2] = 0.0;
} else {
sum = Y/y;
out[0] = x * sum;
out[1] = Y;
out[2] = z * sum;
}
}
/* Y & xy to XYZ */
extern ICCLIB_API void icmY_xy2XYZ(double *out, double *xy, double Y) {
double x = xy[0];
double y = xy[1];
double z = 1.0 - x - y;
double sum;
if (y < 1e-9) {
out[0] = out[1] = out[2] = 0.0;
} else {
sum = Y/y;
out[0] = x * sum;
out[1] = Y;
out[2] = z * sum;
}
}
/* CIE XYZ to perceptual CIE 1976 UCS diagram Yu'v'*/
/* (Yu'v' is a better linear chromaticity space than Yxy) */
extern ICCLIB_API void icmXYZ21976UCS(double *out, double *in) {
double X = in[0], Y = in[1], Z = in[2];
double den, u, v;
den = (X + 15.0 * Y + 3.0 * Z);
if (den < 1e-9) {
Y = 0.0;
u = 4.0/19.0;
v = 9.0/19.0;
} else {
u = (4.0 * X) / den;
v = (9.0 * Y) / den;
}
out[0] = Y;
out[1] = u;
out[2] = v;
}
/* Perceptual CIE 1976 UCS diagram Yu'v' to CIE XYZ */
extern ICCLIB_API void icm1976UCS2XYZ(double *out, double *in) {
double u, v, fl, fu, fv, sum, X, Y, Z;
Y = in[0];
u = in[1];
v = in[2];
if (v < 1e-9) {
X = Y = Z = 0.0;
} else {
X = ((9.0 * u * Y)/(4.0 * v));
Z = -(((20.0 * v + 3.0 * u - 12.0) * Y)/(4.0 * v));
}
out[0] = X;
out[1] = Y;
out[2] = Z;
}
/* CIE 1976 UCS diagram Y & u'v' to XYZ */
extern ICCLIB_API void icm1976UCSY_uv2XYZ(double *out, double *uv, double Y) {
double u, v, fl, fu, fv, sum, X, Z;
u = uv[0];
v = uv[1];
if (v < 1e-9) {
X = Y = Z = 0.0;
} else {
X = ((9.0 * u * Y)/(4.0 * v));
Z = -(((20.0 * v + 3.0 * u - 12.0) * Y)/(4.0 * v));
}
out[0] = X;
out[1] = Y;
out[2] = Z;
}
/* CIE XYZ to perceptual CIE 1976 UCS diagram u'v'*/
/* (u'v' is a better linear chromaticity space than xy) */
extern ICCLIB_API void icmXYZ21976UCSuv(double *out, double *in) {
double X = in[0], Y = in[1], Z = in[2];
double den, u, v;
den = (X + 15.0 * Y + 3.0 * Z);
if (den < 1e-9) {
u = 4.0/19.0;
v = 9.0/19.0;
} else {
u = (4.0 * X) / den;
v = (9.0 * Y) / den;
}
out[0] = u;
out[1] = v;
}
/* CIE XYZ to perceptual CIE 1960 UCS */
/* (This was obsoleted by the 1976UCS, but is still used */
/* in computing color temperatures.) */
extern ICCLIB_API void icmXYZ21960UCS(double *out, double *in) {
double X = in[0], Y = in[1], Z = in[2];
double den, u, v;
den = (X + 15.0 * Y + 3.0 * Z);
if (den < 1e-9) {
Y = 0.0;
u = 4.0/19.0;
v = 6.0/19.0;
} else {
u = (4.0 * X) / den;
v = (6.0 * Y) / den;
}
out[0] = Y;
out[1] = u;
out[2] = v;
}
/* Perceptual CIE 1960 UCS to CIE XYZ */
extern ICCLIB_API void icm1960UCS2XYZ(double *out, double *in) {
double u, v, fl, fu, fv, sum, X, Y, Z;
Y = in[0];
u = in[1];
v = in[2];
if (v < 1e-9) {
X = Y = Z = 0.0;
} else {
X = ((3.0 * u * Y)/(2.0 * v));
Z = -(((10.0 * v + u - 4.0) * Y)/(2.0 * v));
}
out[0] = X;
out[1] = Y;
out[2] = Z;
}
/* CIE XYZ to perceptual CIE 1964 WUV (U*V*W*) */
/* (This is obsolete but still used in computing CRI) */
extern ICCLIB_API void icmXYZ21964WUV(icmXYZNumber *w, double *out, double *in) {
double W, U, V;
double wucs[3];
double iucs[3];
icmXYZ2Ary(wucs, *w);
icmXYZ21960UCS(wucs, wucs);
icmXYZ21960UCS(iucs, in);
W = 25.0 * pow(iucs[0] * 100.0/wucs[0], 1.0/3.0) - 17.0;
U = 13.0 * W * (iucs[1] - wucs[1]);
V = 13.0 * W * (iucs[2] - wucs[2]);
out[0] = W;
out[1] = U;
out[2] = V;
}
/* Perceptual CIE 1964 WUV (U*V*W*) to CIE XYZ */
extern ICCLIB_API void icm1964WUV2XYZ(icmXYZNumber *w, double *out, double *in) {
double W, U, V;
double wucs[3];
double iucs[3];
icmXYZ2Ary(wucs, *w);
icmXYZ21960UCS(wucs, wucs);
W = in[0];
U = in[1];
V = in[2];
iucs[0] = pow((W + 17.0)/25.0, 3.0) * wucs[0]/100.0;
iucs[1] = U / (13.0 * W) + wucs[1];
iucs[2] = V / (13.0 * W) + wucs[2];
icm1960UCS2XYZ(out, iucs);
}
/* CIE CIE1960 UCS to perceptual CIE 1964 WUV (U*V*W*) */
/* (This is used in computing CRI) */
extern ICCLIB_API void icm1960UCS21964WUV(icmXYZNumber *w, double *out, double *in) {
double W, U, V;
double wucs[3];
icmXYZ2Ary(wucs, *w);
icmXYZ21960UCS(wucs, wucs);
W = 25.0 * pow(in[0] * 100.0/wucs[0], 1.0/3.0) - 17.0;
U = 13.0 * W * (in[1] - wucs[1]);
V = 13.0 * W * (in[2] - wucs[2]);
out[0] = W;
out[1] = U;
out[2] = V;
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* NOTE :- that these values are for the 1931 standard observer. */
/* Since they are an arbitrary 4 decimal place accuracy, we round */
/* them to be exactly the same as ICC header encoded values, */
/* to avoid any slight discrepancy with PCS white from profiles. */
/* available D50 Illuminant */
icmXYZNumber icmD50 = { /* Profile illuminant - D50 */
RND_S15FIXED16(0.9642),
RND_S15FIXED16(1.0000),
RND_S15FIXED16(0.8249)
};
icmXYZNumber icmD50_100 = { /* Profile illuminant - D50, scaled to 100 */
RND_S15FIXED16(0.9642) * 100.0,
RND_S15FIXED16(1.0000) * 100.0,
RND_S15FIXED16(0.8249) * 100.0
};
double icmD50_ary3[3] = { /* Profile illuminant - D50 */
RND_S15FIXED16(0.9642),
RND_S15FIXED16(1.0000),
RND_S15FIXED16(0.8249)
};
double icmD50_100_ary3[3] = { /* Profile illuminant - D50, scaled to 100 */
RND_S15FIXED16(0.9642) * 100.0,
RND_S15FIXED16(1.0000) * 100.0,
RND_S15FIXED16(0.8249) * 100.0
};
/* available D65 Illuminant */
icmXYZNumber icmD65 = { /* Profile illuminant - D65 */
RND_S15FIXED16(0.9505),
RND_S15FIXED16(1.0000),
RND_S15FIXED16(1.0890)
};
icmXYZNumber icmD65_100 = { /* Profile illuminant - D65, scaled to 100 */
RND_S15FIXED16(0.9505) * 100.0,
RND_S15FIXED16(1.0000) * 100.0,
RND_S15FIXED16(1.0890) * 100.0
};
double icmD65_ary3[3] = { /* Profile illuminant - D65 */
RND_S15FIXED16(0.9505),
RND_S15FIXED16(1.0000),
RND_S15FIXED16(1.0890)
};
double icmD65_100_ary3[3] = { /* Profile illuminant - D65, scaled to 100 */
RND_S15FIXED16(0.9505) * 100.0,
RND_S15FIXED16(1.0000) * 100.0,
RND_S15FIXED16(1.0890) * 100.0
};
/* Default black point */
icmXYZNumber icmBlack = {
0.0000, 0.0000, 0.0000
};
/* The Standard ("wrong Von-Kries") chromatic transform matrix */
double icmWrongVonKries[3][3] = {
{ 1.0000, 0.0000, 0.0000 },
{ 0.0000, 1.0000, 0.0000 },
{ 0.0000, 0.0000, 1.0000 }
};
/* The Bradford chromatic transform matrix */
double icmBradford[3][3] = {
{ RND_S15FIXED16( 0.8951), RND_S15FIXED16( 0.2664), RND_S15FIXED16(-0.1614) },
{ RND_S15FIXED16(-0.7502), RND_S15FIXED16( 1.7135), RND_S15FIXED16( 0.0367) },
{ RND_S15FIXED16( 0.0389), RND_S15FIXED16(-0.0685), RND_S15FIXED16( 1.0296) }
};
/* Return the normal Delta E given two Lab values */
double icmLabDE(double *Lab0, double *Lab1) {
double rv = 0.0, tt;
tt = Lab0[0] - Lab1[0];
rv += tt * tt;
tt = Lab0[1] - Lab1[1];
rv += tt * tt;
tt = Lab0[2] - Lab1[2];
rv += tt * tt;
return sqrt(rv);
}
/* Return the normal Delta E squared, given two Lab values */
double icmLabDEsq(double *Lab0, double *Lab1) {
double rv = 0.0, tt;
tt = Lab0[0] - Lab1[0];
rv += tt * tt;
tt = Lab0[1] - Lab1[1];
rv += tt * tt;
tt = Lab0[2] - Lab1[2];
rv += tt * tt;
return rv;
}
/* Return the normal Delta E squared given two XYZ values */
extern ICCLIB_API double icmXYZLabDEsq(icmXYZNumber *w, double *in0, double *in1) {
double lab0[3], lab1[3], rv;
icmXYZ2Lab(w, lab0, in0);
icmXYZ2Lab(w, lab1, in1);
rv = icmLabDEsq(lab0, lab1);
return rv;
}
/* Return the normal Delta E given two XYZ values */
extern ICCLIB_API double icmXYZLabDE(icmXYZNumber *w, double *in0, double *in1) {
double lab0[3], lab1[3], rv;
icmXYZ2Lab(w, lab0, in0);
icmXYZ2Lab(w, lab1, in1);
rv = icmLabDE(lab0, lab1);
return rv;
}
/* Return the normal Delta E squared given two XYZ values */
extern ICCLIB_API double icmXYZLptDEsq(icmXYZNumber *w, double *in0, double *in1) {
double lab0[3], lab1[3], rv;
icmXYZ2Lpt(w, lab0, in0);
icmXYZ2Lpt(w, lab1, in1);
rv = icmLabDEsq(lab0, lab1);
return rv;
}
/* Return the normal Delta E given two XYZ values */
extern ICCLIB_API double icmXYZLptDE(icmXYZNumber *w, double *in0, double *in1) {
double lab0[3], lab1[3], rv;
icmXYZ2Lpt(w, lab0, in0);
icmXYZ2Lpt(w, lab1, in1);
rv = icmLabDE(lab0, lab1);
return rv;
}
/* (Note that CIE94 can give odd results for very large delta E's, */
/* when one of the two points is near the neutral axis: */
/* ie DE(A,B + del) != DE(A,B) + DE(del) */
#ifdef NEVER
{
double w1[3] = { 99.996101, -0.058417, -0.010557 };
double c1[3] = { 60.267956, 98.845863, -61.163277 };
double w2[3] = { 100.014977, -0.138339, 0.089744 };
double c2[3] = { 60.294464, 98.117104, -60.843159 };
printf("DE 1 = %f, 2 = %f\n", icmLabDE(c1, w1), icmLabDE(c2, w2));
printf("DE94 1 = %f, 2 = %f\n", icmCIE94(c1, w1), icmCIE94(c2, w2));
}
#endif
/* Return the CIE94 Delta E color difference measure, squared */
double icmCIE94sq(double Lab0[3], double Lab1[3]) {
double desq, dhsq;
double dlsq, dcsq;
double c12;
{
double dl, da, db;
dl = Lab0[0] - Lab1[0];
dlsq = dl * dl; /* dl squared */
da = Lab0[1] - Lab1[1];
db = Lab0[2] - Lab1[2];
/* Compute normal Lab delta E squared */
desq = dlsq + da * da + db * db;
}
{
double c1, c2, dc;
/* Compute chromanance for the two colors */
c1 = sqrt(Lab0[1] * Lab0[1] + Lab0[2] * Lab0[2]);
c2 = sqrt(Lab1[1] * Lab1[1] + Lab1[2] * Lab1[2]);
c12 = sqrt(c1 * c2); /* Symetric chromanance */
/* delta chromanance squared */
dc = c1 - c2;
dcsq = dc * dc;
}
/* Compute delta hue squared */
if ((dhsq = desq - dlsq - dcsq) < 0.0)
dhsq = 0.0;
{
double sc, sh;
/* Weighting factors for delta chromanance & delta hue */
sc = 1.0 + 0.045 * c12;
sh = 1.0 + 0.015 * c12;
return dlsq + dcsq/(sc * sc) + dhsq/(sh * sh);
}
}
/* Return the CIE94 Delta E color difference measure */
double icmCIE94(double Lab0[3], double Lab1[3]) {
return sqrt(icmCIE94sq(Lab0, Lab1));
}
/* Return the CIE94 Delta E color difference measure for two XYZ values */
extern ICCLIB_API double icmXYZCIE94(icmXYZNumber *w, double *in0, double *in1) {
double lab0[3], lab1[3];
icmXYZ2Lab(w, lab0, in0);
icmXYZ2Lab(w, lab1, in1);
return sqrt(icmCIE94sq(lab0, lab1));
}
/* From the paper "The CIEDE2000 Color-Difference Formula: Implementation Notes, */
/* Supplementary Test Data, and Mathematical Observations", by */
/* Gaurav Sharma, Wencheng Wu and Edul N. Dalal, */
/* Color Res. Appl., vol. 30, no. 1, pp. 21-30, Feb. 2005. */
/* Return the CIEDE2000 Delta E color difference measure squared, for two Lab values */
double icmCIE2Ksq(double *Lab0, double *Lab1) {
double C1, C2;
double h1, h2;
double dL, dC, dH;
double dsq;
/* The trucated value of PI is needed to ensure that the */
/* test cases pass, as one of them lies on the edge of */
/* a mathematical discontinuity. The precision is still */
/* enough for any practical use. */
#define RAD2DEG(xx) (180.0/M_PI * (xx))
#define DEG2RAD(xx) (M_PI/180.0 * (xx))
/* Compute Cromanance and Hue angles */
{
double C1ab, C2ab;
double Cab, Cab7, G;
double a1, a2;
C1ab = sqrt(Lab0[1] * Lab0[1] + Lab0[2] * Lab0[2]);
C2ab = sqrt(Lab1[1] * Lab1[1] + Lab1[2] * Lab1[2]);
Cab = 0.5 * (C1ab + C2ab);
Cab7 = pow(Cab,7.0);
G = 0.5 * (1.0 - sqrt(Cab7/(Cab7 + 6103515625.0)));
a1 = (1.0 + G) * Lab0[1];
a2 = (1.0 + G) * Lab1[1];
C1 = sqrt(a1 * a1 + Lab0[2] * Lab0[2]);
C2 = sqrt(a2 * a2 + Lab1[2] * Lab1[2]);
if (C1 < 1e-9)
h1 = 0.0;
else {
h1 = RAD2DEG(atan2(Lab0[2], a1));
if (h1 < 0.0)
h1 += 360.0;
}
if (C2 < 1e-9)
h2 = 0.0;
else {
h2 = RAD2DEG(atan2(Lab1[2], a2));
if (h2 < 0.0)
h2 += 360.0;
}
}
/* Compute delta L, C and H */
{
double dh;
dL = Lab1[0] - Lab0[0];
dC = C2 - C1;
if (C1 < 1e-9 || C2 < 1e-9) {
dh = 0.0;
} else {
dh = h2 - h1;
if (dh > 180.0)
dh -= 360.0;
else if (dh < -180.0)
dh += 360.0;
}
dH = 2.0 * sqrt(C1 * C2) * sin(DEG2RAD(0.5 * dh));
}
{
double L, C, h, T;
double hh, ddeg;
double C7, RC, L50sq, SL, SC, SH, RT;
double dLsq, dCsq, dHsq, RCH;
L = 0.5 * (Lab0[0] + Lab1[0]);
C = 0.5 * (C1 + C2);
if (C1 < 1e-9 || C2 < 1e-9) {
h = h1 + h2;
} else {
h = h1 + h2;
if (fabs(h1 - h2) > 180.0) {
if (h < 360.0)
h += 360.0;
else if (h >= 360.0)
h -= 360.0;
}
h *= 0.5;
}
T = 1.0 - 0.17 * cos(DEG2RAD(h-30.0)) + 0.24 * cos(DEG2RAD(2.0 * h))
+ 0.32 * cos(DEG2RAD(3.0 * h + 6.0)) - 0.2 * cos(DEG2RAD(4.0 * h - 63.0));
L50sq = (L - 50.0) * (L - 50.0);
SL = 1.0 + (0.015 * L50sq)/sqrt(20.0 + L50sq);
SC = 1.0 + 0.045 * C;
SH = 1.0 + 0.015 * C * T;
dLsq = dL/SL;
dCsq = dC/SC;
dHsq = dH/SH;
hh = (h - 275.0)/25.0;
ddeg = 30.0 * exp(-hh * hh);
C7 = pow(C, 7.0);
RC = 2.0 * sqrt(C7/(C7 + 6103515625.0));
RT = -sin(DEG2RAD(2 * ddeg)) * RC;
RCH = RT * dCsq * dHsq;
dLsq *= dLsq;
dCsq *= dCsq;
dHsq *= dHsq;
dsq = dLsq + dCsq + dHsq + RCH;
}
return dsq;
#undef RAD2DEG
#undef DEG2RAD
}
/* Return the CIE2DE000 Delta E color difference measure for two Lab values */
double icmCIE2K(double *Lab0, double *Lab1) {
return sqrt(icmCIE2Ksq(Lab0, Lab1));
}
/* Return the CIEDE2000 Delta E color difference measure for two XYZ values */
ICCLIB_API double icmXYZCIE2K(icmXYZNumber *w, double *in0, double *in1) {
double lab0[3], lab1[3];
icmXYZ2Lab(w, lab0, in0);
icmXYZ2Lab(w, lab1, in1);
return sqrt(icmCIE2Ksq(lab0, lab1));
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Independent chromatic adaptation transform utility. */
/* Return a 3x3 chromatic adaptation matrix */
/* Use icmMulBy3x3(dst, mat, src) */
/* NOTE that to transform primaries they */
/* must be mat[XYZ][RGB] format! */
void icmChromAdaptMatrix(
int flags, /* Use Bradford flag, Transform given matrix flag */
icmXYZNumber d_wp, /* Destination white point */
icmXYZNumber s_wp, /* Source white point */
double mat[3][3] /* Destination matrix */
) {
double dst[3], src[3]; /* Source & destination white points */
double vkmat[3][3]; /* Von Kries matrix */
static int inited = 0; /* Compute inverse bradford once */
static double ibradford[3][3]; /* Inverse Bradford */
/* Set initial matrix to unity if creating from scratch */
if (!(flags & ICM_CAM_MULMATRIX)) {
icmSetUnity3x3(mat);
}
icmXYZ2Ary(src, s_wp);
icmXYZ2Ary(dst, d_wp);
if (flags & ICM_CAM_BRADFORD) {
icmMulBy3x3(src, icmBradford, src);
icmMulBy3x3(dst, icmBradford, dst);
}
/* Setup the Von Kries white point adaptation matrix */
vkmat[0][0] = dst[0]/src[0];
vkmat[1][1] = dst[1]/src[1];
vkmat[2][2] = dst[2]/src[2];
vkmat[0][1] = vkmat[0][2] = 0.0;
vkmat[1][0] = vkmat[1][2] = 0.0;
vkmat[2][0] = vkmat[2][1] = 0.0;
/* Transform to Bradford space if requested */
if (flags & ICM_CAM_BRADFORD) {
icmMul3x3(mat, icmBradford);
}
/* Apply chromatic adaptation */
icmMul3x3(mat, vkmat);
/* Transform from Bradford space */
if (flags & ICM_CAM_BRADFORD) {
if (inited == 0) {
icmInverse3x3(ibradford, icmBradford);
inited = 1;
}
icmMul3x3(mat, ibradford);
}
/* We're done */
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* RGB XYZ primaries device to RGB->XYZ transform matrix. */
/* We assume that the device is perfectly additive, but that */
/* there may be a scale factor applied to the channels to */
/* match the white point at RGB = 1. */
/* Use icmMulBy3x3(dst, mat, src) */
/* Return non-zero if matrix would be singular */
int icmRGBXYZprim2matrix(
double red[3], /* Red colorant */
double green[3], /* Green colorant */
double blue[3], /* Blue colorant */
double white[3], /* White point */
double mat[3][3] /* Destination matrix[RGB][XYZ] */
) {
double tmat[3][3];
double t[3];
/* Assemble the colorants into a matrix */
tmat[0][0] = red[0];
tmat[0][1] = green[0];
tmat[0][2] = blue[0];
tmat[1][0] = red[1];
tmat[1][1] = green[1];
tmat[1][2] = blue[1];
tmat[2][0] = red[2];
tmat[2][1] = green[2];
tmat[2][2] = blue[2];
/* Compute the inverse */
if (icmInverse3x3(mat, tmat))
return 1;
/* Compute scale vector that maps colorants to white point */
t[0] = mat[0][0] * white[0]
+ mat[0][1] * white[1]
+ mat[0][2] * white[2];
t[1] = mat[1][0] * white[0]
+ mat[1][1] * white[1]
+ mat[1][2] * white[2];
t[2] = mat[2][0] * white[0]
+ mat[2][1] * white[1]
+ mat[2][2] * white[2];
/* Now formulate the transform matrix */
mat[0][0] = red[0] * t[0];
mat[0][1] = green[0] * t[1];
mat[0][2] = blue[0] * t[2];
mat[1][0] = red[1] * t[0];
mat[1][1] = green[1] * t[1];
mat[1][2] = blue[1] * t[2];
mat[2][0] = red[2] * t[0];
mat[2][1] = green[2] * t[1];
mat[2][2] = blue[2] * t[2];
return 0;
}
/* RGB Yxy primaries to device to RGB->XYZ transform matrix */
/* Return non-zero if matrix would be singular */
/* Use icmMulBy3x3(dst, mat, src) */
int icmRGBYxyprim2matrix(
double red[3], /* Red colorant */
double green[3], /* Green colorant */
double blue[3], /* Blue colorant */
double white[3], /* White point */
double mat[3][3], /* Return matrix[RGB][XYZ] */
double wXYZ[3] /* Return white XYZ */
) {
double r[3], g[3], b[3];
icmYxy2XYZ(r, red);
icmYxy2XYZ(g, green);
icmYxy2XYZ(b, blue);
icmYxy2XYZ(wXYZ, white);
return icmRGBXYZprim2matrix(r, g, b, wXYZ, mat);
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Pre-round a 3x3 matrix to ensure that the product of */
/* the matrix and the input value is the same as */
/* the quantized matrix product. This is used to improve accuracy */
/* of 'chad' tag in computing absolute white point. */
void icmQuantize3x3S15Fixed16(
double targ[3], /* Target of product */
double mat[3][3], /* matrix[][] to be quantized */
double in[3] /* Input of product - must not be 0.0! */
) {
int i, j;
double sum[3]; /* == target */
double tmp[3]; /* Uncorrected sum */
// printf("In = %.8f %.8f %.8f\n",in[0], in[1], in[2]);
// printf("Target = %.8f %.8f %.8f\n",targ[0], targ[1], targ[2]);
for (j = 0; j < 3; j++)
sum[j] = targ[j];
/* Pre-quantize the matrix, and then ensure that the */
/* sum of the product of the quantized values is the quantized */
/* sum by assigning the rounding error to the largest component. */
for (i = 0; i < 3; i++) {
int bix = 0;
double bval = -1e9;
/* locate the largest and quantize each matrix component */
for (j = 0; j < 3; j++) {
if (fabs(mat[i][j]) > bval) { /* Locate largest */
bix = j;
bval = fabs(mat[i][j]);
}
mat[i][j] = round_S15Fixed16Number(mat[i][j]);
}
/* Check the product of the uncorrected quantized values */
tmp[i] = 0.0;
for (j = 0; j < 3; j++)
tmp[i] += mat[i][j] * in[j];
/* Compute the value the largest has to be */
/* to ensure that sum of the quantized mat[][] times in[] is */
/* equal to the quantized sum. */
for (j = 0; j < 3; j++) {
if (j == bix)
continue;
sum[i] -= mat[i][j] * in[j];
}
mat[i][bix] = round_S15Fixed16Number(sum[i]/in[i]);
/* Check the product of the corrected quantized values */
sum[i] = 0.0;
for (j = 0; j < 3; j++)
sum[i] += mat[i][j] * in[j];
}
// printf("Q Sum = %.8f %.8f %.8f\n",tmp[0], tmp[1], tmp[2]);
// printf("Q cor Sum = %.8f %.8f %.8f\n",sum[0], sum[1], sum[2]);
}
/* Pre-round RGB device primary values to ensure that */
/* the sum of the quantized primaries is the same as */
/* the quantized sum. */
/* [Note matrix is transposed compared to quantize3x3S15Fixed16() ] */
void quantizeRGBprimsS15Fixed16(
double mat[3][3] /* matrix[RGB][XYZ] */
) {
int i, j;
double sum[3];
// printf("D50 = %f %f %f\n",icmD50.X, icmD50.Y, icmD50.Z);
/* Compute target sum of primary XYZ */
for (i = 0; i < 3; i++) {
sum[i] = 0.0;
for (j = 0; j < 3; j++)
sum[i] += mat[j][i];
}
// printf("Sum = %f %f %f\n",sum[0], sum[1], sum[2]);
/* Pre-quantize the primary XYZ's, and then ensure that the */
/* sum of the quantized values is the quantized sum by assigning */
/* the rounding error to the largest component. */
for (i = 0; i < 3; i++) {
int bix = 0;
double bval = -1e9;
/* locate the largest and quantize each component */
for (j = 0; j < 3; j++) {
if (fabs(mat[j][i]) > bval) { /* Locate largest */
bix = j;
bval = fabs(mat[j][i]);
}
mat[j][i] = round_S15Fixed16Number(mat[j][i]);
}
/* Compute the value the largest has to be */
/* to ensure that sum of the quantized values is */
/* equal to the quantized sum */
for (j = 0; j < 3; j++) {
if (j == bix)
continue;
sum[i] -= mat[j][i];
}
mat[bix][i] = round_S15Fixed16Number(sum[i]);
/* Check the sum of the quantized values */
// sum[i] = 0.0;
// for (j = 0; j < 3; j++)
// sum[i] += mat[j][i];
}
// printf("Q cor Sum = %f %f %f\n",sum[0], sum[1], sum[2]);
}
/* - - - - - - - - - - - - - - - - - - - - - - - - */
/* Some PCS utility functions */
/* Clip Lab, while maintaining hue angle. */
/* Return nz if clipping occured */
int icmClipLab(double out[3], double in[3]) {
double C;
out[0] = in[0];
out[1] = in[1];
out[2] = in[2];
if (out[0] >= 0.0 && out[0] <= 100.0
&& out[1] >= -128.0 && out[1] <= 127.0
&& out[2] >= -128.0 && out[2] <= 127.0)
return 0;
/* Clip L */
if (out[0] < 0.0)
out[0] = 0.0;
else if (out[0] > 100.0)
out[0] = 100.0;
C = out[1];
if (fabs(out[2]) > fabs(C))
C = out[2];
if (C < -128.0 || C > 127.0) {
if (C < 0.0)
C = -128.0/C;
else
C = 127.0/C;
out[1] *= C;
out[2] *= C;
}
return 1;
}
/* Clip XYZ, while maintaining hue angle */
/* Return nz if clipping occured */
int icmClipXYZ(double out[3], double in[3]) {
out[0] = in[0];
out[1] = in[1];
out[2] = in[2];
if (out[0] >= 0.0 && out[0] <= 1.9999
&& out[1] >= 0.0 && out[1] <= 1.9999
&& out[2] >= 0.0 && out[2] <= 1.9999)
return 0;
/* Clip Y and scale X and Z similarly */
if (out[1] > 1.9999) {
out[0] *= 1.9999/out[1];
out[2] *= 1.9999/out[1];
out[1] = 1.9999;
} else if (out[1] < 0.0) {
out[0] = 0.0;
out[1] = 0.0;
out[2] = 0.0;
}
if (out[0] < 0.0 || out[0] > 1.9999
|| out[2] < 0.0 || out[2] > 1.9999) {
double D50[3] = { 0.9642, 1.0000, 0.8249 };
double bb = 0.0;
/* Scale the D50 so that it has the same Y value as our color */
D50[0] *= out[1];
D50[1] *= out[1];
D50[2] *= out[1];
/* Figure out what blend factor with Y scaled D50, brings our */
/* color X and Z values into range */
if (out[0] < 0.0) {
double b;
b = (0.0 - out[0])/(D50[0] - out[0]);
if (b > bb)
bb = b;
} else if (out[0] > 1.9999) {
double b;
b = (1.9999 - out[0])/(D50[0] - out[0]);
if (b > bb)
bb = b;
}
if (out[2] < 0.0) {
double b;
b = (0.0 - out[2])/(D50[2] - out[2]);
if (b > bb)
bb = b;
} else if (out[2] > 1.9999) {
double b;
b = (1.9999 - out[2])/(D50[2] - out[2]);
if (b > bb)
bb = b;
}
/* Do the desaturation */
out[0] = D50[0] * bb + (1.0 - bb) * out[0];
out[2] = D50[2] * bb + (1.0 - bb) * out[2];
}
return 1;
}
/* --------------------------------------------------------------- */
/* Some video specific functions */
/* Should add ST.2048 log functions */
/* Convert Lut table index/value to YPbPr */
/* (Same as Lut_Lut2YPbPr() ) */
void icmLut2YPbPr(double *out, double *in) {
out[0] = in[0]; /* Y */
out[1] = in[1] - 0.5; /* Cb */
out[2] = in[2] - 0.5; /* Cr */
}
/* Convert YPbPr to Lut table index/value */
/* (Same as Lut_YPbPr2Lut() ) */
void icmYPbPr2Lut(double *out, double *in) {
out[0] = in[0]; /* Y */
out[1] = in[1] + 0.5; /* Cb */
out[2] = in[2] + 0.5; /* Cr */
}
/* Convert Rec601 RGB' into YPbPr, or "full range YCbCr" */
/* where input 0..1, output 0..1, -0.5 .. 0.5, -0.5 .. 0.5 */
/* [From the Rec601 spec. ] */
void icmRec601_RGBd_2_YPbPr(double out[3], double in[3]) {
double tt[3];
tt[0] = 0.299 * in[0] + 0.587 * in[1] + 0.114 * in[2];
tt[1] = -0.299 /1.772 * in[0]
+ -0.587 /1.772 * in[1]
+ (1.0-0.114)/1.772 * in[2];
tt[2] = (1.0-0.299)/1.402 * in[0]
+ -0.587 /1.402 * in[1]
+ -0.114 /1.402 * in[2];
out[0] = tt[0];
out[1] = tt[1];
out[2] = tt[2];
}
/* Convert Rec601 YPbPr to RGB' (== "full range YCbCr") */
/* where input 0..1, -0.5 .. 0.5, -0.5 .. 0.5, output 0.0 .. 1 */
/* [Inverse of above] */
void icmRec601_YPbPr_2_RGBd(double out[3], double in[3]) {
double tt[3];
tt[0] = 1.000000000 * in[0] + 0.000000000 * in[1] + 1.402000000 * in[2];
tt[1] = 1.000000000 * in[0] + -0.344136286 * in[1] + -0.714136286 * in[2];
tt[2] = 1.000000000 * in[0] + 1.772000000 * in[1] + 0.000000000 * in[2];
out[0] = tt[0];
out[1] = tt[1];
out[2] = tt[2];
}
/* Convert Rec709 1150/60/2:1 RGB' into YPbPr, or "full range YCbCr" */
/* where input 0..1, output 0..1, -0.5 .. 0.5, -0.5 .. 0.5 */
/* (This is for digital Rec709 and is very close to analog Rec709 60Hz.) */
/* [From the Rec709 specification] */
void icmRec709_RGBd_2_YPbPr(double out[3], double in[3]) {
double tt[3];
tt[0] = 0.2126 * in[0] + 0.7152 * in[1] + 0.0722 * in[2];
tt[1] = 1.0/1.8556 * -0.2126 * in[0]
+ 1.0/1.8556 * -0.7152 * in[1]
+ 1.0/1.8556 * (1.0-0.0722) * in[2];
tt[2] = 1.0/1.5748 * (1.0-0.2126) * in[0]
+ 1.0/1.5748 * -0.7152 * in[1]
+ 1.0/1.5748 * -0.0722 * in[2];
out[0] = tt[0];
out[1] = tt[1];
out[2] = tt[2];
}
/* Convert Rec709 1150/60/2:1 YPbPr to RGB' (== "full range YCbCr") */
/* where input 0..1, -0.5 .. 0.5, -0.5 .. 0.5, output 0.0 .. 1 */
/* (This is for digital Rec709 and is very close to analog Rec709 60Hz.) */
/* [Inverse of above] */
void icmRec709_YPbPr_2_RGBd(double out[3], double in[3]) {
double tt[3];
tt[0] = 1.000000000 * in[0] + 0.000000000 * in[1] + 1.574800000 * in[2];
tt[1] = 1.000000000 * in[0] + -0.187324273 * in[1] + -0.468124273 * in[2];
tt[2] = 1.000000000 * in[0] + 1.855600000 * in[1] + 0.000000000 * in[2];
out[0] = tt[0];
out[1] = tt[1];
out[2] = tt[2];
}
/* Convert Rec709 1250/50/2:1 RGB' into YPbPr, or "full range YCbCr" */
/* where input 0..1, output 0..1, -0.5 .. 0.5, -0.5 .. 0.5 */
/* (This is for analog Rec709 50Hz) */
/* [From the Rec709 specification] */
void icmRec709_50_RGBd_2_YPbPr(double out[3], double in[3]) {
double tt[3];
tt[0] = 0.299 * in[0] + 0.587 * in[1] + 0.114 * in[2];
tt[1] = 0.564 * -0.299 * in[0]
+ 0.564 * -0.587 * in[1]
+ 0.564 * (1.0-0.114) * in[2];
tt[2] = 0.713 * (1.0-0.299) * in[0]
+ 0.713 * -0.587 * in[1]
+ 0.713 * -0.114 * in[2];
out[0] = tt[0];
out[1] = tt[1];
out[2] = tt[2];
}
/* Convert Rec709 1250/50/2:1 YPbPr to RGB' (== "full range YCbCr") */
/* where input 0..1, -0.5 .. 0.5, -0.5 .. 0.5, output 0.0 .. 1 */
/* (This is for analog Rec709 50Hz) */
/* [Inverse of above] */
void icmRec709_50_YPbPr_2_RGBd(double out[3], double in[3]) {
double tt[3];
tt[0] = 1.000000000 * in[0] + 0.000000000 * in[1] + 1.402524544 * in[2];
tt[1] = 1.000000000 * in[0] + -0.344340136 * in[1] + -0.714403473 * in[2];
tt[2] = 1.000000000 * in[0] + 1.773049645 * in[1] + 0.000000000 * in[2];
out[0] = tt[0];
out[1] = tt[1];
out[2] = tt[2];
}
/* Convert Rec2020 RGB' into Non-constant luminance YPbPr, or "full range YCbCr" */
/* where input 0..1, output 0..1, -0.5 .. 0.5, -0.5 .. 0.5 */
/* [From the Rec2020 specification] */
void icmRec2020_NCL_RGBd_2_YPbPr(double out[3], double in[3]) {
double tt[3];
tt[0] = 0.2627 * in[0] + 0.6780 * in[1] + 0.0593 * in[2];
tt[1] = 1/1.8814 * -0.2627 * in[0]
+ 1/1.8814 * -0.6780 * in[1]
+ 1/1.8814 * (1.0-0.0593) * in[2];
tt[2] = 1/1.4746 * (1.0-0.2627) * in[0]
+ 1/1.4746 * -0.6780 * in[1]
+ 1/1.4746 * -0.0593 * in[2];
out[0] = tt[0];
out[1] = tt[1];
out[2] = tt[2];
}
/* Convert Rec2020 Non-constant luminance YPbPr into RGB' (== "full range YCbCr") */
/* where input 0..1, -0.5 .. 0.5, -0.5 .. 0.5, output 0.0 .. 1 */
/* [Inverse of above] */
void icmRec2020_NCL_YPbPr_2_RGBd(double out[3], double in[3]) {
double tt[3];
tt[0] = 1.000000000 * in[0] + 0.000000000 * in[1] + 1.474600000 * in[2];
tt[1] = 1.000000000 * in[0] + -0.164553127 * in[1] + -0.571353127 * in[2];
tt[2] = 1.000000000 * in[0] + 1.881400000 * in[1] + 0.000000000 * in[2];
out[0] = tt[0];
out[1] = tt[1];
out[2] = tt[2];
}
/* Convert Rec2020 RGB' into Constant luminance YPbPr, or "full range YCbCr" */
/* where input 0..1, output 0..1, -0.5 .. 0.5, -0.5 .. 0.5 */
/* [From the Rec2020 specification] */
void icmRec2020_CL_RGBd_2_YPbPr(double out[3], double in[3]) {
int i;
double tt[3];
/* Convert RGB' to RGB */
for (i = 0; i < 3; i++) {
if (in[i] < (4.5 * 0.0181))
tt[i] = in[i]/4.5;
else
tt[i] = pow((in[i] + 0.0993)/1.0993, 1.0/0.45);
}
/* Y value */
tt[0] = 0.2627 * tt[0] + 0.6780 * tt[1] + 0.0593 * tt[2];
/* Y' value */
if (tt[0] < 0.0181)
tt[0] = tt[0] * 4.5;
else
tt[0] = 1.0993 * pow(tt[0], 0.45) - 0.0993;
tt[1] = in[2] - tt[0];
if (tt[1] <= 0.0)
tt[1] /= 1.9404;
else
tt[1] /= 1.5816;
tt[2] = in[0] - tt[0];
if (tt[2] <= 0.0)
tt[2] /= 1.7184;
else
tt[2] /= 0.9936;
out[0] = tt[0];
out[1] = tt[1];
out[2] = tt[2];
}
/* Convert Rec2020 Constant luminance YPbPr into RGB' (== "full range YCbCr") */
/* where input 0..1, -0.5 .. 0.5, -0.5 .. 0.5, output 0.0 .. 1 */
/* [Inverse of above] */
void icmRec2020_CL_YPbPr_2_RGBd(double out[3], double in[3]) {
int i;
double tin[3], tt[3];
/* Y' */
tin[0] = in[0];
/* B' - Y' */
if (in[1] <= 0.0)
tin[1] = 1.9404 * in[1];
else
tin[1] = 1.5816 * in[1];
/* R' - Y' */
if (in[2] <= 0.0)
tin[2] = 1.7184 * in[2];
else
tin[2] = 0.9936 * in[2];
/* R' */
tt[0] = tin[2] + tin[0];
/* Y' */
tt[1] = tin[0];
/* B' */
tt[2] = tin[1] + tin[0];
/* Convert RYB' to RYB */
for (i = 0; i < 3; i++) {
if (tt[i] < (4.5 * 0.0181))
tin[i] = tt[i]/4.5;
else
tin[i] = pow((tt[i] + 0.0993)/1.0993, 1.0/0.45);
}
/* G */
tt[1] = (tin[1] - 0.2627 * tin[0] - 0.0593 * tin[2])/0.6780;
/* G' */
if (tt[1] < 0.0181)
tt[1] = tt[1] * 4.5;
else
tt[1] = 1.0993 * pow(tt[1], 0.45) - 0.0993;
out[0] = tt[0];
out[1] = tt[1];
out[2] = tt[2];
}
/* Convert Rec601/Rec709/Rec2020 YPbPr to YCbCr Video range. */
/* input 0..1, -0.5 .. 0.5, -0.5 .. 0.5, */
/* output 16/255 .. 235/255, 16/255 .. 240/255, 16/255 .. 240/255 */
void icmRecXXX_YPbPr_2_YCbCr(double out[3], double in[3]) {
out[0] = ((235.0 - 16.0) * in[0] + 16.0)/255.0;
out[1] = ((128.0 - 16.0) * 2.0 * in[1] + 128.0)/255.0;
out[2] = ((128.0 - 16.0) * 2.0 * in[2] + 128.0)/255.0;
}
/* Convert Rec601/Rec709/Rec2020 Video YCbCr to YPbPr range. */
/* input 16/255 .. 235/255, 16/255 .. 240/255, 16/255 .. 240/255 */
/* output 0..1, -0.5 .. 0.5, -0.5 .. 0.5, */
void icmRecXXX_YCbCr_2_YPbPr(double out[3], double in[3]) {
out[0] = (255.0 * in[0] - 16.0)/(235.0 - 16.0);
out[1] = (255.0 * in[1] - 128.0)/(2.0 * (128.0 - 16.0));
out[2] = (255.0 * in[2] - 128.0)/(2.0 * (128.0 - 16.0));
}
/* Convert full range RGB to Video range 16..235 RGB */
void icmRGB_2_VidRGB(double out[3], double in[3]) {
out[0] = ((235.0 - 16.0) * in[0] + 16.0)/255.0;
out[1] = ((235.0 - 16.0) * in[1] + 16.0)/255.0;
out[2] = ((235.0 - 16.0) * in[2] + 16.0)/255.0;
}
/* Convert Video range 16..235 RGB to full range RGB */
/* Return nz if outside RGB range */
void icmVidRGB_2_RGB(double out[3], double in[3]) {
out[0] = (255.0 * in[0] - 16.0)/(235.0 - 16.0);
out[1] = (255.0 * in[1] - 16.0)/(235.0 - 16.0);
out[2] = (255.0 * in[2] - 16.0)/(235.0 - 16.0);
}
/* =============================================================== */
/* PS 3.14-2009, Digital Imaging and Communications in Medicine */
/* (DICOM) Part 14: Grayscale Standard Display Function */
/* JND index value 1..1023 to L 0.05 .. 3993.404 cd/m^2 */
static double icmDICOM_fwd_nl(double jnd) {
double a = -1.3011877;
double b = -2.5840191e-2;
double c = 8.0242636e-2;
double d = -1.0320229e-1;
double e = 1.3646699e-1;
double f = 2.8745620e-2;
double g = -2.5468404e-2;
double h = -3.1978977e-3;
double k = 1.2992634e-4;
double m = 1.3635334e-3;
double jj, num, den, rv;
jj = jnd = log(jnd);
num = a;
den = 1.0;
num += c * jj;
den += b * jj;
jj *= jnd;
num += e * jj;
den += d * jj;
jj *= jnd;
num += g * jj;
den += f * jj;
jj *= jnd;
num += m * jj;
den += h * jj;
jj *= jnd;
den += k * jj;
rv = pow(10.0, num/den);
return rv;
}
/* JND index value 1..1023 to L 0.05 .. 3993.404 cd/m^2 */
double icmDICOM_fwd(double jnd) {
if (jnd < 0.5)
jnd = 0.5;
if (jnd > 1024.0)
jnd = 1024.0;
return icmDICOM_fwd_nl(jnd);
}
/* L 0.05 .. 3993.404 cd/m^2 to JND index value 1..1023 */
/* This is not super accurate - typically to 0.03 .. 0.1 jne. */
static double icmDICOM_bwd_apx(double L) {
double A = 71.498068;
double B = 94.593053;
double C = 41.912053;
double D = 9.8247004;
double E = 0.28175407;
double F = -1.1878455;
double G = -0.18014349;
double H = 0.14710899;
double I = -0.017046845;
double rv, LL;
if (L < 0.049982) { /* == jnd 0.5 */
return 0.5;
}
if (L > 4019.354716) /* == jnd 1024 */
L = 4019.354716;
LL = L = log10(L);
rv = A;
rv += B * LL;
LL *= L;
rv += C * LL;
LL *= L;
rv += D * LL;
LL *= L;
rv += E * LL;
LL *= L;
rv += F * LL;
LL *= L;
rv += G * LL;
LL *= L;
rv += H * LL;
LL *= L;
rv += I * LL;
return rv;
}
/* L 0.05 .. 3993.404 cd/m^2 to JND index value 1..1023 */
/* Polish the aproximate solution twice using Newton's itteration */
double icmDICOM_bwd(double L) {
double rv, Lc, prv, pLc, de;
int i;
if (L < 0.045848) /* == jnd 0.5 */
L = 0.045848;
if (L > 4019.354716) /* == jnd 1024 */
L = 4019.354716;
/* Approx solution */
rv = icmDICOM_bwd_apx(L);
/* Compute aprox derivative */
Lc = icmDICOM_fwd_nl(rv);
prv = rv + 0.01;
pLc = icmDICOM_fwd_nl(prv);
do {
de = (rv - prv)/(Lc - pLc);
prv = rv;
rv -= (Lc - L) * de;
pLc = Lc;
Lc = icmDICOM_fwd_nl(rv);
} while (fabs(Lc - L) > 1e-8);
return rv;
}
/* ---------------------------------------------------------- */
/* Convert an angle in radians into chromatic RGB values */
/* in a simple geometric fashion, with 0 = Red. */
void icmRad2RGB(double rgb[3], double ang) {
double th1 = 1.0/3.0 * 2.0 * M_PI;
double th2 = 2.0/3.0 * 2.0 * M_PI;
double bl;
while (ang < 0.0)
ang += 2.0 * M_PI;
while (ang >= (2.0 * M_PI))
ang -= 2.0 * M_PI;
if (ang < th1) {
bl = ang/th1;
rgb[0] = (1.0 - bl);
rgb[1] = bl;
rgb[2] = 0.0;
} else if (ang < th2) {
bl = (ang - th1)/th1;
rgb[0] = 0.0;
rgb[1] = (1.0 - bl);
rgb[2] = bl;
} else {
bl = (ang - th2)/th1;
rgb[0] = bl;
rgb[1] = 0.0;
rgb[2] = (1.0 - bl);
}
}
/* ---------------------------------------------------------- */
/* Print an int vector to a string. */
/* Returned static buffer is re-used every 5 calls. */
char *icmPiv(int di, int *p) {
# define BUFSZ (MAX_CHAN * 8 * 16)
char *fmt = "%d";
static char buf[5][BUFSZ];
static int ix = 0;
int brem = BUFSZ;
int e;
char *bp;
if (p == NULL)
return "(null)";
if (++ix >= 5)
ix = 0;
bp = buf[ix];
for (e = 0; e < di && brem > 10; e++) {
int tt;
if (e > 0)
*bp++ = ' ', brem--;
tt = snprintf(bp, brem, fmt, p[e]);
if (tt < 0 || tt >= brem)
break; /* Run out of room... */
bp += tt;
brem -= tt;
}
return buf[ix];
# undef BUFSZ
}
/* Print a double color vector to a string. */
/* Returned static buffer is re-used every 5 calls. */
char *icmPdv(int di, double *p) {
return icmPdvf(di, NULL, p);
}
/* Print a double color vector to a string with format. */
/* Returned static buffer is re-used every 5 calls. */
char *icmPdvf(int di, char *fmt, double *p) {
# define NBUF 10
# define BUFSZ (MAX_CHAN * 16)
static char buf[NBUF][BUFSZ];
static int ix = 0;
int brem = BUFSZ;
int e;
char *bp;
if (p == NULL)
return "(null)";
if (fmt == NULL)
fmt = "%.8f";
if (++ix >= NBUF)
ix = 0;
bp = buf[ix];
for (e = 0; e < di && brem > 10; e++) {
int tt;
if (e > 0)
*bp++ = ' ', brem--;
tt = snprintf(bp, brem, fmt, p[e]);
if (tt < 0 || tt >= brem)
break; /* Run out of room... */
bp += tt;
brem -= tt;
}
return buf[ix];
# undef BUFSZ
# undef NBUF
}
/* Print a float color vector to a string. */
/* Returned static buffer is re-used every 5 calls. */
char *icmPfv(int di, float *p) {
# define BUFSZ (MAX_CHAN * 8 * 16)
char *fmt = "%.8f";
static char buf[5][BUFSZ];
static int ix = 0;
int brem = BUFSZ;
int e;
char *bp;
if (p == NULL)
return "(null)";
if (++ix >= 5)
ix = 0;
bp = buf[ix];
for (e = 0; e < di && brem > 10; e++) {
int tt;
if (e > 0)
*bp++ = ' ', brem--;
tt = snprintf(bp, brem, fmt, p[e]);
if (tt < 0 || tt >= brem)
break; /* Run out of room... */
bp += tt;
brem -= tt;
}
return buf[ix];
# undef BUFSZ
}
/* Print an XYZ */
/* Returned static buffer is re-used every 5 calls. */
char *icmPXYZ(icmXYZNumber *p) {
double xyz[3];
icmXYZ2Ary(xyz, *p);
return icmPdv(3, xyz);
}
/* Print an 0..1 range XYZ as a D50 Lab string */
/* Returned static buffer is re-used every 5 calls. */
char *icmPLab(double *p) {
double lab[3];
icmXYZ2Lab(&icmD50, lab, p);
return icmPdv(3, lab);
}
/* ---------------------------------------------------------- */
/* Unicode conversions */
/*
* Copyright 2001-2004 Unicode, Inc.
*
* Disclaimer
*
* This source code is provided as is by Unicode, Inc. No claims are
* made as to fitness for any particular purpose. No warranties of any
* kind are expressed or implied. The recipient agrees to determine
* applicability of information provided. If this file has been
* purchased on magnetic or optical media from Unicode, Inc., the
* sole remedy for any claim will be exchange of defective media
* within 90 days of receipt.
*
* Limitations on Rights to Redistribute This Code
*
* Unicode, Inc. hereby grants the right to freely use the information
* supplied in this file in the creation of products supporting the
* Unicode Standard, and to make copies of this file in any form
* for internal or external distribution as long as this notice
* remains attached.
*/
/* ---------------------------------------------------------------------
Conversions between UTF32, UTF-16, and UTF-8. Source code file.
Author: Mark E. Davis, 1994.
Rev History: Rick McGowan, fixes & updates May 2001.
Sept 2001: fixed const & error conditions per
mods suggested by S. Parent & A. Lillich.
June 2002: Tim Dodd added detection and handling of incomplete
source sequences, enhanced error detection, added casts
to eliminate compiler warnings.
July 2003: slight mods to back out aggressive FFFE detection.
Jan 2004: updated switches in from-UTF8 conversions.
Oct 2004: updated to use UNI_MAX_LEGAL_UTF32 in UTF-32 conversions.
Dec 2022 Modified for use in icclib by Graeme W. Gill.
------------------------------------------------------------------------ */
/* Some fundamental constants */
#define UNI_REPLACEMENT_CHAR (icmUTF32)0x0000FFFD
#define UNI_BOM_CHAR (icmUTF32)0x0000FEFF
#define UNI_MAX_BMP (icmUTF32)0x0000FFFF
#define UNI_MAX_UTF16 (icmUTF32)0x0010FFFF
#define UNI_MAX_UTF32 (icmUTF32)0x7FFFFFFF
#define UNI_MAX_LEGAL_UTF32 (icmUTF32)0x0010FFFF
static const int halfShift = 10; /* used for shifting by 10 bits */
static const icmUTF32 halfBase = 0x0010000;
static const icmUTF32 halfMask = 0x3FF;
#define UNI_SUR_HIGH_START (icmUTF32)0xD800
#define UNI_SUR_HIGH_END (icmUTF32)0xDBFF
#define UNI_SUR_LOW_START (icmUTF32)0xDC00
#define UNI_SUR_LOW_END (icmUTF32)0xDFFF
/* --------------------------------------------------------------------- */
/*
* Index into the table below with the first byte of a UTF-8 sequence to
* get the number of trailing bytes that are supposed to follow it.
* Note that *legal* UTF-8 values can't have 4 or 5-bytes. The table is
* left as-is for anyone who may want to do such conversion, which was
* allowed in earlier algorithms.
*/
static const char trailingBytesForUTF8[256] = {
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2, 3,3,3,3,3,3,3,3,4,4,4,4,5,5,5,5
};
/* Mask of bits that the extra bytes is needed for */
static const icmUTF32 msbitsTrailingBytesForUTF8[6] = {
0x0000007f, /* Bits 0..6 */
0x00000780, /* Bits 7..10 */
0x0000f800, /* Bits 11..15 */
0x001f0000, /* Bits 16..20 */
0x03e00000, /* Bits 21..25 */
0x7c000000 /* Bits 26..30 */
};
/* Mask to leave just payload bits given a UTF-8 byte */
static const ORD8 paylloadMaskForUTF8[256] = {
0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,
0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,
0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,
0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,
0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,
0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,
0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,
0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,0x7f,
0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,
0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,
0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,
0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,0x3f,
0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,
0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,0x1f,
0x0f,0x0f,0x0f,0x0f,0x0f,0x0f,0x0f,0x0f,0x0f,0x0f,0x0f,0x0f,0x0f,0x0f,0x0f,0x0f,
0x07,0x07,0x07,0x07,0x07,0x07,0x07,0x07,0x03,0x03,0x03,0x03,0x01,0x01,0x01,0x01
};
/*
* Magic values subtracted from a buffer value during UTF8 conversion.
* This table contains as many values as there might be trailing bytes
* in a UTF-8 sequence.
*/
static const icmUTF32 offsetsFromUTF8[6] = {
0x00000000UL, 0x00003080UL, 0x000E2080UL,
0x03C82080UL, 0xFA082080UL, 0x82082080UL
};
/*
* Once the bits are split out into bytes of UTF-8, this is a mask OR-ed
* into the first byte, depending on how many bytes follow. There are
* as many entries in this table as there are UTF-8 sequence types.
* (I.e., one byte sequence, two byte... etc.). Remember that sequencs
* for *legal* UTF-8 will be 4 or fewer bytes total.
*/
static const icmUTF8 firstByteMark[7] = { 0x00, 0x00, 0xC0, 0xE0, 0xF0, 0xF8, 0xFC };
/* --------------------------------------------------------------------- */
/* Convert ICC buffer of UTF16 to UTF8 using icmSn serialisation. */
/* Input is expected to be nul terminated on the last character. */
/* Checks and repairs UTF16 string to be legal UTF-8. */
/* Output will be nul terminated. */
/* len is number of UTF16 bytes including nul. */
/* Return resulting size in bytes. */
/* out pointer can be NULL. */
/* pillegal can be NULL. */
/* If nonul is nz then UTF16 is not expected to have a nul terminator */
size_t icmUTF16SntoUTF8(icmUTFerr *pillegal, icmUTF8 *out, icmFBuf *bin, size_t len, int nonul) {
int fix = 1; /* We always fix */
icmUTFerr illegal = icmUTF_ok;
size_t ilen = len;
icmUTF8 *iout = out;
//printf("\n~1 icmUTF16SntoUTF8: len %d bytes\n",len);
if (len & 1)
illegal |= icmUTF_odd_bytes;
for (;;) {
unsigned short bytesToWrite = 0;
icmUTF32 ch;
if (len < 2) {
if (!nonul) {
illegal |= icmUTF_no_nul; /* No nul terminator before end */
//printf(" ~1 icmUTF_no_nul\n");
}
break;
}
icmSn_ui_UInt16(bin, &ch);
len -= 2;
//printf(" ~1 ch = 0x%x remaining len = %d bytes\n",ch, len);
if (ch == 0) {
if (len >= 2) { /* nul terminator prior to end */
illegal |= icmUTF_prem_nul;
//printf(" ~1 icmUTF_prem_nul\n");
}
if (nonul) {
illegal |= icmUTF_unex_nul;
}
break;
}
/* ICC UTF-16 should not have a Byte Order Mark, */
/* and resulting UTF-8 should not reproduce it. */
if (ch == UNI_BOM_CHAR && (ilen == (len+2))) {
illegal |= icmUTF_unn_bom;
//printf(" ~1 icmUTF_unn_bom\n");
continue;
}
/* If we have a surrogate pair, decode them into UTF32 */
if (ch >= UNI_SUR_HIGH_START && ch <= UNI_SUR_HIGH_END) {
icmUTF32 ch2; /* Next UTF16 */
//printf(" ~1 found high surrogate\n");
if (len < 2) {
illegal |= icmUTF_bad_surr;
if (fix) ch = UNI_REPLACEMENT_CHAR;
//printf(" ~1 icmUTF_bad_surr 1\n");
break;
}
icmSn_ui_UInt16(bin, &ch2);
//printf(" ~1 next ch = 0x%x\n",ch2);
/* If the next input is the low pair, convert to UTF32. */
if (ch2 >= UNI_SUR_LOW_START && ch2 <= UNI_SUR_LOW_END) {
//printf(" ~1 found low surrogate\n");
ch = ((ch - UNI_SUR_HIGH_START) << halfShift)
+ (ch2 - UNI_SUR_LOW_START) + halfBase;
len -= 2;
/* A stranded surrogate high */
} else {
icmSn_rseek(bin, -2); /* Put back the character */
illegal |= icmUTF_bad_surr;
if (fix) ch = UNI_REPLACEMENT_CHAR;
//printf(" ~1 icmUTF_bad_surr 2\n");
}
/* An unexpected surrogate low */
} else if (ch >= UNI_SUR_LOW_START && ch <= UNI_SUR_LOW_END) {
//printf(" ~1 found low surrogate\n");
illegal |= icmUTF_bad_surr;
if (fix) ch = UNI_REPLACEMENT_CHAR;
//printf(" ~1 icmUTF_bad_surr 3\n");
}
/* We have a UTF32 value. Convert to UTF8 */
if (ch == 0) {
illegal |= icmUTF_emb_nul; /* Embedded nul */
ch = UNI_REPLACEMENT_CHAR;
//printf(" ~1 icmUTF_emb_nul\n");
}
/* Figure out how many bytes the result will require */
if (ch < (icmUTF32)0x80) { bytesToWrite = 1;
} else if (ch < (icmUTF32)0x800) { bytesToWrite = 2;
} else if (ch < (icmUTF32)0x10000) { bytesToWrite = 3;
} else if (ch < (icmUTF32)0x110000) { bytesToWrite = 4;
} else { bytesToWrite = 3;
/* Shouldn't happen ? */
illegal |= icmUTF_bad_surr;
ch = UNI_REPLACEMENT_CHAR;
//printf(" ~1 icmUTF_bad_surr 4\n");
}
if (iout != NULL) {
const icmUTF32 byteMask = 0xBF;
const icmUTF32 byteMark = 0x80;
out += bytesToWrite;
switch (bytesToWrite) { /* note: everything falls through. */
case 4:
*--out = (icmUTF8)((ch | byteMark) & byteMask); ch >>= 6;
case 3:
*--out = (icmUTF8)((ch | byteMark) & byteMask); ch >>= 6;
case 2:
*--out = (icmUTF8)((ch | byteMark) & byteMask); ch >>= 6;
case 1:
*--out = (icmUTF8)(ch | firstByteMark[bytesToWrite]);
}
}
out += bytesToWrite;
}
if (iout != NULL)
*out = 0;
out++;
if (pillegal != NULL)
*pillegal = illegal;
return (out - iout) * sizeof(icmUTF8);
}
/* General convert UTF16 to UTF8. */
/* Repairs UTF16 string to be legal UTF-8. */
/* Input is expected to be nul terminated. */
/* Output will be nul terminated. */
/* Return resulting size in bytes including nul. */
/* out pointer can be NULL to size buffer. */
/* pillegal can be NULL. */
size_t icmUTF16toUTF8(icmUTFerr *pillegal, icmUTF8 *out, icmUTF16 *in) {
int fix = 1; /* Always fix */
icmUTFerr illegal = icmUTF_ok;
ORD16 *iin = in;
icmUTF8 *iout = out;
for (;;) {
unsigned short bytesToWrite = 0;
icmUTF32 ch;
ch = *in++;;
if (ch == 0) {
break;
}
/* Once in machine representation, a Byte Order Mark is not needed, */
/* and resulting UTF-8 should not reproduce it. */
if (ch == UNI_BOM_CHAR && iin == (in-2)) {
illegal |= icmUTF_unn_bom;
continue;
}
/* If we have a surrogate pair, decode them int UTF32 */
if (ch >= UNI_SUR_HIGH_START && ch <= UNI_SUR_HIGH_END) {
icmUTF32 ch2 = 0; /* Next UTF16 */
ch2 = *in;
/* If the next input is the low pair, convert to UTF32. */
if (ch2 >= UNI_SUR_LOW_START && ch2 <= UNI_SUR_LOW_END) {
ch = ((ch - UNI_SUR_HIGH_START) << halfShift)
+ (ch2 - UNI_SUR_LOW_START) + halfBase;
in++;
/* A stranded surrogate high */
} else {
illegal |= icmUTF_bad_surr;
if (fix) ch = UNI_REPLACEMENT_CHAR;
}
/* An unexpected surrogate low */
} else if (ch >= UNI_SUR_LOW_START && ch <= UNI_SUR_LOW_END) {
illegal |= icmUTF_bad_surr;
if (fix) ch = UNI_REPLACEMENT_CHAR;
}
/* We have a UTF32 value. Convert to UTF8 */
if (ch == 0) {
illegal |= icmUTF_emb_nul; /* Embedded nul */
ch = UNI_REPLACEMENT_CHAR;
}
/* Figure out how many bytes the result will require */
if (ch < (icmUTF32)0x80) { bytesToWrite = 1;
} else if (ch < (icmUTF32)0x800) { bytesToWrite = 2;
} else if (ch < (icmUTF32)0x10000) { bytesToWrite = 3;
} else if (ch < (icmUTF32)0x110000) { bytesToWrite = 4;
} else { bytesToWrite = 3;
/* Shouldn't happen ? */
illegal |= icmUTF_bad_surr;
ch = UNI_REPLACEMENT_CHAR;
}
if (iout != NULL) {
const icmUTF32 byteMask = 0xBF;
const icmUTF32 byteMark = 0x80;
out += bytesToWrite;
switch (bytesToWrite) { /* note: everything falls through. */
case 4:
*--out = (icmUTF8)((ch | byteMark) & byteMask); ch >>= 6;
case 3:
*--out = (icmUTF8)((ch | byteMark) & byteMask); ch >>= 6;
case 2:
*--out = (icmUTF8)((ch | byteMark) & byteMask); ch >>= 6;
case 1:
*--out = (icmUTF8)(ch | firstByteMark[bytesToWrite]);
}
}
out += bytesToWrite;
}
if (iout != NULL)
*out = 0; /* nul terminator */
out++;
if (pillegal != NULL)
*pillegal = illegal;
return (out - iout) * sizeof(icmUTF8);
}
/* --------------------------------------------------------------------- */
/* Convert UTF-8 to ICC buffer of UTF-16 using icmSn serialisation. */
/* Input must be nul terminated, and so will output. */
/* len is number of UTF8 characters including nul. */
/* Return resulting size in bytes including nul. */
/* bout may be NULL to size buffer. */
/* If nonul is nz then UTF16 is not expected to have a nul terminator */
size_t icmUTF8toUTF16Sn(icmUTFerr *pillegal, icmFBuf *bout, icmUTF8 *in, size_t len, int nonul) {
icmUTFerr illegal = icmUTF_ok;
icmUTF16 val;
size_t osize = 0;
//printf("\n~1 icmUTF8toUTF16Sn:\n");
if (in == NULL) {
//printf("~1 in == NULL, setting len = 0\n");
len = 0;
}
for (;;) {
unsigned short extraBytesToRead;
icmUTF32 ch;
//printf("~1 len %zu\n",len);
if (len == 0) {
illegal |= icmUTF_no_nul; /* No nul terminator before end */
break;
}
ch = *in++;
len--;
//printf(" ~1 raw ch = 0x%x\n",ch);
if (ch == 0) {
break; /* We're done */
}
if ((ch & 0xc0) == 0x80) {
illegal |= icmUTF_unex_cont; /* Unexpected continuation byte */
ch = UNI_REPLACEMENT_CHAR;
} else {
int i;
extraBytesToRead = trailingBytesForUTF8[ch];
ch &= paylloadMaskForUTF8[ch]; /* Just the payload bits */
if (extraBytesToRead > 3)
illegal |= icmUTF_tmany_cont; /* Too many continuation bytes */
/* (we'll decode them anyway) */
/* Read trailing bytes while checking for nul terminator and bad encoding */
for (i = 0; i < extraBytesToRead; i++) {
icmUTF32 tch;
if (len == 0) {
illegal |= icmUTF_no_nul; /* No nul terminator before end */
break;
}
tch = *in++;
if (tch == 0) {
illegal |= icmUTF_short_cont; /* Cut short by end of string */
ch = UNI_REPLACEMENT_CHAR;
break;
}
/* If byte hasn't got continuation code */
if ((tch & 0xc0) != 0x80) {
illegal |= icmUTF_short_cont; /* Cut short by next non-continuation */
ch = UNI_REPLACEMENT_CHAR;
break;
}
tch &= 0x3f; /* Just payload bits */
ch = (ch << 6) + tch;
}
/* See if the encoding is over long */
if (extraBytesToRead > 0
&& (ch & msbitsTrailingBytesForUTF8[extraBytesToRead]) == 0) {
illegal |= icmUTF_overlong; /* overlong (doesn't affect decoding) */
}
}
if (len == 0
&& (illegal & icmUTF_no_nul)) /* No nul terminator before end */
break;
/* We now have UTF32 values, so convert to UTF16 */
//printf(" ~1 cooked ch = 0x%x\n",ch);
if (ch == 0) {
illegal |= icmUTF_emb_nul; /* Embedded nul */
ch = UNI_REPLACEMENT_CHAR;
}
/* Is a character <= 0xFFFF */
if (ch <= UNI_MAX_BMP) {
/* UTF-16 surrogate values are illegal in UTF-32 */
if (ch >= UNI_SUR_HIGH_START && ch <= UNI_SUR_LOW_END) {
illegal |= icmUTF_inv_cdpnt; /* Bad code point - surrogate */
if (bout != NULL) {
val = UNI_REPLACEMENT_CHAR;
icmSn_us_UInt16(bout, &val);
//printf(" ~1 output 0x%x\n",val);
}
osize += 2;
/* Normal UTF-16 character */
} else {
if (bout != NULL) {
val = (icmUTF16)ch;
icmSn_us_UInt16(bout, &val);
//printf(" ~1 output 0x%x\n",val);
}
osize += 2;
}
/* Outside unicode range */
} else if (ch > UNI_MAX_UTF16) {
illegal |= icmUTF_ovr_cdpnt; /* Bad code point - beyond unicode */
if (bout != NULL) {
val = UNI_REPLACEMENT_CHAR;
icmSn_us_UInt16(bout, &val);
//printf(" ~1 output 0x%x\n",val);
}
osize += 2;
/* Is a character in range 0x10000 - 0x10FFFF, */
/* so output a surrogate pair */
} else {
ch -= halfBase;
if (bout != NULL) {
val = (icmUTF16)((ch >> halfShift) + UNI_SUR_HIGH_START);
icmSn_us_UInt16(bout, &val);
//printf(" ~1 output 0x%x\n",val);
val = (icmUTF16)((ch & halfMask) + UNI_SUR_LOW_START);
icmSn_us_UInt16(bout, &val);
//printf(" ~1 output 0x%x\n",val);
}
osize += 4;
}
}
if (!nonul) {
if (bout != NULL) {
val = 0;
icmSn_us_UInt16(bout, &val);
//printf(" ~1 output 0x%x\n",val);
}
osize += 2;
}
if (pillegal != NULL)
*pillegal = illegal;
//printf(" ~1 returning osize %zu with illegal = 0x%X\n",osize * sizeof(ORD8),illegal);
return osize * sizeof(ORD8);
}
/* General convert UTF-8 to UTF-16. */
/* Input must be nul terminated, and so will output. */
/* Return resulting size in bytes including nul. */
/* out can be NULL. */
size_t icmUTF8toUTF16(icmUTFerr *pillegal, icmUTF16 *out, icmUTF8 *in) {
icmUTFerr illegal = icmUTF_ok;
icmUTF16 *iout = out;
for (;;) {
unsigned short extraBytesToRead;
icmUTF32 ch;
ch = *in++;
if (ch == 0)
break; /* We're done */
if ((ch & 0xc0) == 0x80) {
illegal |= icmUTF_unex_cont; /* Unexpected continuation byte */
ch = UNI_REPLACEMENT_CHAR;
} else {
int i;
extraBytesToRead = trailingBytesForUTF8[ch];
ch &= paylloadMaskForUTF8[ch]; /* Just the payload bits */
if (extraBytesToRead > 3)
illegal |= icmUTF_tmany_cont; /* Too many extra bytes */
/* (we'll decode them anyway) */
/* Read trailing bytes while checking for nul terminator and bad encoding */
for (i = 0; i < extraBytesToRead; i++) {
icmUTF32 tch;
tch = *in++;
if (tch == 0) {
illegal |= icmUTF_short_cont; /* Cut short by end of string */
ch = UNI_REPLACEMENT_CHAR;
break;
}
/* If byte hasn't got continuation code */
if ((tch & 0xc0) != 0x80) {
illegal |= icmUTF_short_cont; /* Cut short by next non-continuation */
ch = UNI_REPLACEMENT_CHAR;
break;
}
tch &= 0x3f; /* Just payload bits */
ch = (ch << 6) + tch;
}
/* See if the encoding is over long */
if (extraBytesToRead > 0
&& (ch & msbitsTrailingBytesForUTF8[extraBytesToRead]) == 0) {
illegal |= icmUTF_overlong; /* overlong (doesn't affect decoding) */
}
}
/* We now have UTF32 values, so convert to UTF16 */
if (ch == 0) {
illegal |= icmUTF_emb_nul; /* Embedded nul */
ch = UNI_REPLACEMENT_CHAR;
}
/* Is a character <= 0xFFFF */
if (ch <= UNI_MAX_BMP) {
/* UTF-16 surrogate values are illegal in UTF-32 */
if (ch >= UNI_SUR_HIGH_START && ch <= UNI_SUR_LOW_END) {
illegal |= icmUTF_inv_cdpnt; /* Bad code point - surrogate */
if (iout != NULL)
*out = UNI_REPLACEMENT_CHAR;
out++;
/* Normal UTF-16 character */
} else {
if (iout != NULL)
*out = (icmUTF16)ch;
out++;
}
/* Outside unicode range */
} else if (ch > UNI_MAX_UTF16) {
illegal |= icmUTF_ovr_cdpnt; /* Bad code point - beyond unicode */
if (iout != NULL)
*out = UNI_REPLACEMENT_CHAR;
out++;
/* Is a character in range 0x10000 - 0x10FFFF, */
/* so output a surrogate pair */
} else {
ch -= halfBase;
if (iout != NULL) {
out[0] = (icmUTF16)((ch >> halfShift) + UNI_SUR_HIGH_START);
out[1] = (icmUTF16)((ch & halfMask) + UNI_SUR_LOW_START);
}
out += 2;
}
}
if (iout != NULL)
*out = 0;
out++;
if (pillegal != NULL)
*pillegal = illegal;
return (out - iout) * sizeof(icmUTF16);
}
/* --------------------------------------------------------------------- */
/* Convert UTF-8 to printable ASCII + HTML escapes. */
/* Input must be nul terminated, and so will output. */
/* Return resulting size in bytes (including terminating nul). */
/* *pillegal will have icmUTF_non_ascii set if any non-ASCII found. */
/* All output pointers may be NULL */
size_t icmUTF8toHTMLESC(icmUTFerr *pillegal, char *out, icmUTF8 *in) {
icmUTFerr illegal = icmUTF_ok;
int u8 = 0;
char *iout = out;
for (;;) {
unsigned short extraBytesToRead;
icmUTF32 ch;
ch = *in++;
if (ch == 0)
break; /* We're done */
if ((ch & 0xc0) == 0x80) {
illegal |= icmUTF_unex_cont; /* Unexpected continuation byte */
ch = UNI_REPLACEMENT_CHAR;
} else {
int i;
extraBytesToRead = trailingBytesForUTF8[ch];
ch &= paylloadMaskForUTF8[ch]; /* Just the payload bits */
if (extraBytesToRead > 3)
illegal |= icmUTF_tmany_cont; /* Too many extra bytes */
/* (we'll decode them anyway) */
/* Read trailing bytes while checking for nul terminator and bad encoding */
for (i = 0; i < extraBytesToRead; i++) {
icmUTF32 tch;
tch = *in++;
if (tch == 0) {
illegal |= icmUTF_short_cont; /* Cut short by end of string */
ch = UNI_REPLACEMENT_CHAR;
break;
}
/* If byte hasn't got continuation code */
if ((tch & 0xc0) != 0x80) {
illegal |= icmUTF_short_cont; /* Cut short by next non-continuation */
ch = UNI_REPLACEMENT_CHAR;
break;
}
tch &= 0x3f; /* Just payload bits */
ch = (ch << 6) + tch;
}
/* See if the encoding is over long */
if (extraBytesToRead > 0
&& (ch & msbitsTrailingBytesForUTF8[extraBytesToRead]) == 0) {
illegal |= icmUTF_overlong; /* overlong (doesn't affect decoding) */
}
}
/* We now have UTF32 values */
if (ch > 0x7f)
illegal |= icmUTF_non_ascii; /* non-ASCII found */
/* HTML special characters */
if (ch == '&') {
if (iout != NULL)
memmove(out, "&amp", 4);
out += 4;
} else if (ch == '<') {
if (iout != NULL)
memmove(out, "&lt", 3);
out += 3;
} else if (ch == '>') {
if (iout != NULL)
memmove(out, "&gt", 3);
out += 3;
} else if (ch == '"') {
if (iout != NULL)
memmove(out, "&quot", 5);
out += 5;
} else if (ch == '\'') {
if (iout != NULL)
memmove(out, "&#39", 4);
out += 4;
/* Printable ASCII */
} else if (ch >= 0x20 && ch < 0x80) {
if (iout != NULL) *out = ch;
out++;
/* HTML Excape form "&#DDDD;" */
} else {
char buf[17], *bp;
if (iout != NULL) *out = '&';
out++;
if (iout != NULL) *out = '#';
out++;
/* Do a simple utoa */
bp = &buf[17];
*--bp = '\000';
if (ch == 0)
*--bp = '0';
for (; ch; ch /= 10)
*--bp = '0' + ch % 10;
/* Copy it out */
for ( ;*bp != '\000'; bp++) {
if (iout != NULL) *out = *bp;
out++;
}
if (iout != NULL) *out = ';';
out++;
}
}
if (iout != NULL)
*out = 0;
out++;
if (pillegal != NULL)
*pillegal = illegal;
return (out - iout) * sizeof(char);
}
/* --------------------------------------------------------------------- */
/* Convert UTF-8 to ASCII. Out of range characters are replaced by '?'. */
/* Input must be nul terminated, and so will output. */
/* Return resulting size in bytes (including terminating nul). */
/* *pillegal will have icmUTF_non_ascii set if any non-ASCII found. */
/* All output pointers may be NULL */
size_t icmUTF8toASCII(icmUTFerr *pillegal, char *out, icmUTF8 *in) {
icmUTFerr illegal = icmUTF_ok;
int u8 = 0;
char *iout = out;
char replacement_char = '?';
for (;;) {
unsigned short extraBytesToRead;
icmUTF32 ch;
ch = *in++;
if (ch == 0)
break; /* We're done */
if ((ch & 0xc0) == 0x80) {
illegal |= icmUTF_unex_cont; /* Unexpected continuation byte */
ch = replacement_char;
} else {
int i;
extraBytesToRead = trailingBytesForUTF8[ch];
ch &= paylloadMaskForUTF8[ch]; /* Just the payload bits */
if (extraBytesToRead > 3)
illegal |= icmUTF_tmany_cont; /* Too many extra bytes */
/* (we'll decode them anyway) */
/* Read trailing bytes while checking for nul terminator and bad encoding */
for (i = 0; i < extraBytesToRead; i++) {
icmUTF32 tch;
tch = *in++;
if (tch == 0) {
illegal |= icmUTF_short_cont; /* Cut short by end of string */
ch = replacement_char;
break;
}
/* If byte hasn't got continuation code */
if ((tch & 0xc0) != 0x80) {
illegal |= icmUTF_short_cont; /* Cut short by next non-continuation */
ch = replacement_char;
break;
}
tch &= 0x3f; /* Just payload bits */
ch = (ch << 6) + tch;
}
/* See if the encoding is over long */
if (extraBytesToRead > 0
&& (ch & msbitsTrailingBytesForUTF8[extraBytesToRead]) == 0) {
illegal |= icmUTF_overlong; /* overlong (doesn't affect decoding) */
}
}
/* We now have UTF32 values */
if (ch > 0x7f) {
ch = replacement_char;
illegal |= icmUTF_non_ascii; /* non-ASCII found */
}
if (iout != NULL) *out = ch;
out++;
}
if (iout != NULL)
*out = 0;
out++;
if (pillegal != NULL)
*pillegal = illegal;
return (out - iout) * sizeof(char);
}
/* Convert ICC buffer of ASCIIZ to UTF8 using icmSn serialisation. */
/* Input is expected to be nul terminated on the last character. */
/* ilen is the expected number of ASCIIZ input bytes including nul. */
/* If fxlen is > 0, then the input is from a fixed length buffer of fxlen. */
/* If fxlen is < 0, then the input has a maximum length of -fxlen */
/* (If fxlen != 0 then usually ilen = abs(fxlen)) */
/* Output will be nul terminated. */
/* Return resulting size in bytes. */
/* out pointer can be NULL. */
/* pillegal can be NULL. */
size_t icmASCIIZSntoUTF8(icmUTFerr *pillegal, icmUTF8 *out, icmFBuf *bin, size_t ilen, int fxlen) {
icmUTFerr illegal = icmUTF_ok;
char replacement_char = '?';
icmUTF8 *iout = out;
size_t remin = ilen; /* fixed/max length buffer bytes left */
icmUTF32 ch;
//printf("~1 icmASCIIZSntoUTF8 got ilen %zu fxlen %zu\n",ilen,fxlen);
if (fxlen > 0) {
remin = fxlen;
} else if (fxlen < 0) {
remin = -fxlen;
}
if (fxlen != 0 && ilen > remin) {
illegal |= icmUTF_ascii_toolong; /* Can't be longer than fixed/max length */
ilen = remin;
}
for (;;) {
if (ilen < 1) {
illegal |= icmUTF_no_nul; /* No nul terminator before end */
break;
}
icmSn_ui_UInt8(bin, &ch);
ilen --;
if (fxlen != 0)
remin--;
if (ch == 0) { /* Got to the end */
if (fxlen == 0 && ilen > 0)
illegal |= icmUTF_prem_nul;
break;
}
/* Make sure this is legal ASCII */
if (ch > 0x7f) {
ch = replacement_char;
illegal |= icmUTF_non_ascii; /* non-ASCII found */
}
/* Since this is ASCII, we don't need any conversion to UTF-8 */
if (iout != NULL)
*out = (icmUTF8)ch;
out++;
}
/* Use up the rest of a fixed buffer */
while (fxlen > 0 && remin > 0) {
icmSn_ui_UInt8(bin, &ch);
remin--;
}
if (iout != NULL) /* Add nul terminator */
*out = 0;
out++;
if (pillegal != NULL)
*pillegal = illegal;
//printf("~1 returning %zu\n",(out - iout) * sizeof(icmUTF8));
return (out - iout) * sizeof(icmUTF8);
}
/* Convert UTF-8 to ICC buffer of ASCIIZ using icmSn serialisation. */
/* Input is expected to be nul terminated, and so will output. */
/* ilen is number of UTF8 characters including nul of the input. */
/* If fxlen is > 0, then the output is to a fixed length buffer of that length, */
/* and will be padded with 0 up to that length. */
/* If fxlen < 0, then the otput has a maxumum length of -fxlen. */
/* Return resulting size in bytes. */
/* bout may be NULL. */
size_t icmUTF8toASCIIZSn(icmUTFerr *pillegal, icmFBuf *bout, icmUTF8 *in, size_t ilen, int fxlen) {
icmUTFerr illegal = icmUTF_ok;
icmUTF32 replacement_char = '?';
ORD8 val;
size_t osize = 0;
size_t remout; /* fixed/max length buffer bytes left */
icmUTF32 ch;
if (fxlen > 0) {
remout = fxlen;
} else if (fxlen < 0) {
remout = -fxlen;
}
if (in == NULL)
ilen = 0;
for (;;) {
unsigned short extraBytesToRead;
if (ilen == 0) {
illegal |= icmUTF_no_nul; /* No nul terminator before end */
break;
}
ch = *in++;
ilen--;
if (ch == 0) {
if (ilen > 0)
illegal |= icmUTF_prem_nul;
break; /* We're done */
}
if ((ch & 0xc0) == 0x80) {
illegal |= icmUTF_unex_cont; /* Unexpected continuation byte */
ch = replacement_char;
} else {
int i;
extraBytesToRead = trailingBytesForUTF8[ch];
ch &= paylloadMaskForUTF8[ch]; /* Just the payload bits */
if (extraBytesToRead > 3)
illegal |= icmUTF_tmany_cont; /* Too many continuation bytes */
/* (we'll decode them anyway) */
/* Read trailing bytes while checking for nul terminator and bad encoding */
for (i = 0; i < extraBytesToRead; i++) {
icmUTF32 tch;
if (ilen == 0) {
illegal |= icmUTF_no_nul; /* No nul terminator before end */
break;
}
tch = *in++;
if (tch == 0) {
illegal |= icmUTF_short_cont; /* Cut short by end of string */
ch = replacement_char;
break;
}
/* If byte hasn't got continuation code */
if ((tch & 0xc0) != 0x80) {
illegal |= icmUTF_short_cont; /* Cut short by next non-continuation */
ch = replacement_char;
break;
}
tch &= 0x3f; /* Just payload bits */
ch = (ch << 6) + tch;
}
/* See if the encoding is over long */
if (extraBytesToRead > 0
&& (ch & msbitsTrailingBytesForUTF8[extraBytesToRead]) == 0) {
illegal |= icmUTF_overlong; /* overlong (doesn't affect decoding) */
}
}
if (ilen == 0
&& (illegal & icmUTF_no_nul)) /* No nul terminator before end */
break;
/* We now have UTF32 values, so convert to UTF16 */
if (ch > 0x7f) {
ch = replacement_char;
illegal |= icmUTF_non_ascii; /* non-ASCII found */
}
if (fxlen != 0 && remout <= 1) { /* Won't be room for final nul */
illegal |= icmUTF_ascii_toolong; /* Can't be longer than fixed/max buffer */
break;
}
if (bout != NULL) {
val = (ORD8) ch;
icmSn_uc_UInt8(bout, &val);
}
osize++;
if (fxlen != 0)
remout--;
}
/* add nul terminator */
if (bout != NULL) {
val = 0;
icmSn_uc_UInt8(bout, &val);
}
osize++;
if (fxlen != 0)
remout--;
/* if fixed length, pad */
while (fxlen > 0 && remout > 0) {
ch = 0;
if (bout != NULL) {
icmSn_ui_UInt8(bout, &ch);
}
remout--;
}
if (pillegal != NULL)
*pillegal = illegal;
return osize * sizeof(ORD8);
}
/* --------------------------------------------------------------------- */
/* ICC ScriptCode support */
/* Convert ICC buffer of ScriptCode to UTF8 using icmSn serialisation. */
/* Input is expected to be nul terminated on the last character. */
/* len is number of ASCIIZ bytes including nul. */
/* len may be 0, in which case there is no ScriptCode */
/* Output will be nul terminated. */
/* Return resulting size in bytes. */
/* out pointer can be NULL. */
/* pillegal can be NULL. */
size_t icmSntoScriptCode(icmUTFerr *pillegal, icmUTF8 *out, icmFBuf *bin, size_t len) {
icmUTFerr illegal = icmUTF_ok;
char replacement_char = '?';
size_t ilen = len;
icmUTF8 *iout = out;
icmUTF32 ch;
size_t remin = 67;
if (len > remin) {
illegal |= icmUTF_sc_toolong; /* Can't be longer than 67 */
len = remin;
}
for (;;) {
if (len < 1) {
if (ilen > 0)
illegal |= icmUTF_no_nul; /* No nul terminator before end */
break;
}
icmSn_ui_UInt8(bin, &ch);
len--;
remin--;
if (ch == 0) {
if (len > 0)
illegal |= icmUTF_no_nul;
break;
}
if (iout != NULL)
*out = (icmUTF8)ch;
out++;
}
/* Use up the rest of the fixed ScriptCode buffer length */
while (remin > 0) {
/* Note this won't pick up an un-padded script code when it is */
/* simply concatenated within a larger tag such as ProfileSequenceDescType */
if (bin->get_space(bin) == 0) { /* Wasn't padded */
illegal |= icmUTF_sc_tooshort;
break;
}
icmSn_ui_UInt8(bin, &ch); /* Could check if ch == 0 ? */
remin--;
}
if (ilen > 0) {
/* Add nul */
if (iout != NULL)
*out = 0;
out++;
}
if (pillegal != NULL)
*pillegal = illegal;
return (out - iout) * sizeof(icmUTF8);
}
/* Convert ScriptCode to ICC buffer of ScriptCode using icmSn serialisation. */
/* Input must be nul terminated, and so will output. */
/* len is number of UTF8 characters including nul. */
/* Input may be NULL and len == 0 in which case there is no ScriptCode */
/* Buffer output will be nul terminated and always 67 bytes */
/* If buffer is NULL or len is 0, the output will be 67 bytes of 0 */
/* Return resulting string size in bytes. (file may be larger with padding) */
/* bout may be NULL. */
size_t icmScriptCodetoSn(icmUTFerr *pillegal, icmFBuf *bout, ORD8 *in, size_t len) {
icmUTFerr illegal = icmUTF_ok;
icmUTF32 replacement_char = '?';
ORD8 ch;
size_t osize = 0;
size_t remout = 67; /* Remaining output space */
if (in == NULL)
len = 0;
for (;;) {
if (len == 0) {
if (in != NULL)
illegal |= icmUTF_no_nul; /* No nul terminator before end */
break;
}
ch = *in++;
len--;
if (ch == 0)
break; /* We're done */
if (remout <= 1) { /* No room for final nul */
illegal |= icmUTF_sc_toolong; /* Can't be longer than 67 */
break;
}
if (bout != NULL) {
icmSn_uc_UInt8(bout, &ch);
}
osize++;
remout--;
}
if (in != NULL) {
/* nul terminator */
if (bout != NULL) {
ch = 0;
icmSn_uc_UInt8(bout, &ch);
}
osize++;
remout--;
}
/* pad until remout== 0 */
while (remout > 0) {
if (bout != NULL) {
ch = 0;
icmSn_uc_UInt8(bout, &ch);
}
remout--;
}
if (pillegal != NULL)
*pillegal = illegal;
return osize * sizeof(ORD8);
}
/* --------------------------------------------------------------------- */
/* Convert icmUTFerr to a comma separated list of error descriptions */
struct {
icmUTFerr mask;
char *string;
} UTFerrStrings[] = {
{ icmUTF_emb_nul, "embedded nul" },
{ icmUTF_no_nul, "no nul terminator" },
{ icmUTF_unex_nul, "unexpected nul terminator" },
{ icmUTF_prem_nul, "premature nul terminator" },
{ icmUTF_bad_surr, "bad surrogate" },
{ icmUTF_unn_bom, "unnecessary BOM" },
{ icmUTF_unex_cont, "unexpected continuation byte" },
{ icmUTF_tmany_cont, "to many continuation bytes" },
{ icmUTF_short_cont, "not enough continuation bytes" },
{ icmUTF_overlong, "overlong encoding" },
{ icmUTF_inv_cdpnt, "invalid codepoint" },
{ icmUTF_ovr_cdpnt, "over-range codepoint" },
{ icmUTF_non_ascii, "non-ASCII characters" },
{ icmUTF_ascii_toolong,"ASCII longer than fixed sized buffer" },
{ icmUTF_sc_tooshort, "ScriptCode occupies less than 67 bytes" },
{ icmUTF_sc_toolong, "ScriptCode longer than 67 bytes" },
{ 0, NULL },
};
char *icmUTFerr2str(icmUTFerr err) {
static char buf[400];
char *bp = buf;
int ne = 0;
int i;
if (err == icmUTF_ok)
return "ok";
for (i = 0; UTFerrStrings[i].string != NULL; i++) {
if (err & UTFerrStrings[i].mask) {
if (ne)
bp += sprintf(bp, ", ");
bp += sprintf(bp, "%s",UTFerrStrings[i].string);
ne = 1;
}
}
return buf;
}
Reference in New Issue View Git Blame Copy Permalink