From f6b8e0eae4374f339487a33e3e4fe5462d5816e1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B6rg=20Frings-F=C3=BCrst?= Date: Sat, 25 Nov 2017 10:16:00 +0100 Subject: New upstream version 2.0.0 --- icc/icc.c | 698 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 651 insertions(+), 47 deletions(-) mode change 100644 => 100755 icc/icc.c (limited to 'icc/icc.c') diff --git a/icc/icc.c b/icc/icc.c old mode 100644 new mode 100755 index b239060..de8c539 --- a/icc/icc.c +++ b/icc/icc.c @@ -804,7 +804,7 @@ static void read_PCSNumber(icc *icp, icColorSpaceSignature csig, double pcs[3], pcs[1] = read_DCS16Number(p+2); pcs[2] = read_DCS16Number(p+4); } - switch (csig) { + switch ((int)csig) { case icSigXYZData: Lut_Lut2XYZ(pcs, pcs); break; @@ -840,7 +840,7 @@ static int write_PCSNumber(icc *icp, icColorSpaceSignature csig, double pcs[3], csig = icmSigLabV2Data; } - switch (csig) { + switch ((int)csig) { case icSigXYZData: Lut_XYZ2Lut(v, pcs); break; @@ -876,7 +876,7 @@ static int write_PCSNumber(icc *icp, icColorSpaceSignature csig, double pcs[3], /* Public: */ int read_Primitive(icc *icp, icmPrimType ptype, void *prim, char *p) { - switch(ptype) { + switch (ptype) { case icmUInt8Number: *((unsigned int *)prim) = read_UInt8Number(p); return 0; @@ -944,7 +944,7 @@ int read_Primitive(icc *icp, icmPrimType ptype, void *prim, char *p) { /* Public: */ int write_Primitive(icc *icp, icmPrimType ptype, char *p, void *prim) { - switch(ptype) { + switch (ptype) { case icmUInt8Number: return write_UInt8Number(*((unsigned int *)prim), p); case icmUInt16Number: @@ -1073,7 +1073,7 @@ static int check_null_string16(char *cp, int len) { /* Color Space to number of component conversion */ /* Return 0 on error */ static unsigned int number_ColorSpaceSignature(icColorSpaceSignature sig) { - switch(sig) { + switch ((int)sig) { case icSigXYZData: return 3; case icSigLabData: @@ -1134,6 +1134,9 @@ static unsigned int number_ColorSpaceSignature(icColorSpaceSignature sig) { return 1; case icmSigLData: return 1; + case icmSigLptData: + return 3; + case icmSigL8Data: return 1; case icmSigLV2Data: @@ -1171,7 +1174,7 @@ static int chnames_ColorSpaceSignature( icColorSpaceSignature sig, char *cvals[] /* Pointers to return for each channel */ ) { - switch (sig) { + switch ((int)sig) { case icSigXYZData: cvals[0] = "CIE X"; cvals[1] = "CIE Y"; @@ -1246,6 +1249,12 @@ char *cvals[] /* Pointers to return for each channel */ cvals[0] = "CIE L*"; return 2; + case icmSigLptData: + cvals[0] = "L"; + cvals[1] = "p"; + cvals[2] = "t"; + return 2; + default: break; @@ -1382,7 +1391,7 @@ static char *string_AsciiOrBinaryData(unsigned int flags) { /* public tags and sizes */ static const char *string_TagSignature(icTagSignature sig) { static char buf[80]; - switch(sig) { + switch ((int)sig) { case icSigAToB0Tag: return "AToB0 Multidimentional Transform"; case icSigAToB1Tag: @@ -1487,7 +1496,7 @@ static const char *string_TagSignature(icTagSignature sig) { /* technology signature descriptions */ static const char *string_TechnologySignature(icTechnologySignature sig) { static char buf[80]; - switch(sig) { + switch (sig) { case icSigDigitalCamera: return "Digital Camera"; case icSigFilmScanner: @@ -1541,7 +1550,7 @@ static const char *string_TechnologySignature(icTechnologySignature sig) { /* type signatures */ static const char *string_TypeSignature(icTagTypeSignature sig) { static char buf[80]; - switch(sig) { + switch (sig) { case icSigCurveType: return "Curve"; case icSigDataType: @@ -1599,7 +1608,7 @@ static const char *string_TypeSignature(icTagTypeSignature sig) { /* Color Space Signatures */ static const char *string_ColorSpaceSignature(icColorSpaceSignature sig) { static char buf[80]; - switch(sig) { + switch ((int)sig) { case icSigXYZData: return "XYZ"; case icSigLabData: @@ -1662,6 +1671,9 @@ static const char *string_ColorSpaceSignature(icColorSpaceSignature sig) { return "L"; case icmSigL8Data: return "L"; + case icmSigLptData: + return "Lpt"; + case icmSigLV2Data: return "L"; case icmSigLV4Data: @@ -1692,7 +1704,7 @@ char *ColorSpaceSignature2str(icColorSpaceSignature sig) { /* profileClass enumerations */ static const char *string_ProfileClassSignature(icProfileClassSignature sig) { static char buf[80]; - switch(sig) { + switch (sig) { case icSigInputClass: return "Input"; case icSigDisplayClass: @@ -1716,7 +1728,7 @@ static const char *string_ProfileClassSignature(icProfileClassSignature sig) { /* Platform Signatures */ static const char *string_PlatformSignature(icPlatformSignature sig) { static char buf[80]; - switch(sig) { + switch ((int)sig) { case icSigMacintosh: return "Macintosh"; case icSigMicrosoft: @@ -1738,7 +1750,7 @@ static const char *string_PlatformSignature(icPlatformSignature sig) { /* Measurement Geometry, used in the measurmentType tag */ static const char *string_MeasurementGeometry(icMeasurementGeometry sig) { static char buf[30]; - switch(sig) { + switch (sig) { case icGeometryUnknown: return "Unknown"; case icGeometry045or450: @@ -1754,7 +1766,7 @@ static const char *string_MeasurementGeometry(icMeasurementGeometry sig) { /* Rendering Intents, used in the profile header */ static const char *string_RenderingIntent(icRenderingIntent sig) { static char buf[30]; - switch(sig) { + switch((int)sig) { case icPerceptual: return "Perceptual"; case icRelativeColorimetric: @@ -5130,7 +5142,7 @@ double *in /* Input array[outputChan] */ int rv = 0; double *gp; /* Pointer to grid cube base */ double co[MAX_CHAN]; /* Coordinate offset with the grid cell */ - int si[MAX_CHAN]; /* co[] Sort index, [0] = smalest */ + int si[MAX_CHAN]; /* co[] Sort index, [0] = smallest */ /* We are using a simplex (ie. tetrahedral for 3D input) interpolation. */ /* This method is more appropriate for XYZ/RGB/CMYK input spaces, */ @@ -5388,7 +5400,7 @@ double *in /* Input array[outputChan] */ int rv = 0; double *gp; /* Pointer to grid cube base */ double co[MAX_CHAN]; /* Coordinate offset with the grid cell */ - int si[MAX_CHAN]; /* co[] Sort index, [0] = smalest */ + int si[MAX_CHAN]; /* co[] Sort index, [0] = smallest */ /* We are using a simplex (ie. tetrahedral for 3D input) interpolation. */ /* This method is more appropriate for XYZ/RGB/CMYK input spaces, */ @@ -5825,7 +5837,7 @@ int icmSetMultiLutTables( iv = _iv + maxchan; /* Allow for "index under" and smoothing radius values */ /* Setup input table value min-max */ - if (inmin == NULL || imax == NULL) { + if (inmin == NULL || inmax == NULL) { #ifdef SYMETRICAL_DEFAULT_LAB_RANGE /* Symetrical default range. */ /* We are assuming V2 Lab16 encoding, since this is a lut16type that always uses */ /* this encoding */ @@ -13548,6 +13560,13 @@ void icmClamp3(double out[3], double in[3]) { out[i] = in[i] < 0.0 ? 0.0 : in[i]; } +/* Invert 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]; @@ -13576,6 +13595,18 @@ void icmMul3(double out[3], double in1[3], double in2[3]) { out[2] = in1[2] * in2[2]; } +/* Take 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); + } +} + /* Take absolute of a 3 vector */ void icmAbs3(double out[3], double in[3]) { out[0] = fabs(in[0]); @@ -13585,6 +13616,16 @@ void icmAbs3(double out[3], double in[3]) { /* - - - - - - - - - - - - - - - - - - - - - - - - */ +/* 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; @@ -13609,6 +13650,17 @@ void icmCpy3x3(double dst[3][3], double src[3][3]) { } } +/* 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; @@ -13645,17 +13697,6 @@ void icmMulBy3x3(double out[3], double mat[3][3], double in[3]) { } /* - - - - - - - - - - - - - - - - - - - - - - - - */ -/* 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]; - } -} - /* - - - - - - - - - - - - - - - - - - - - - - - - */ /* Tensor product. Multiply two 3 vectors to form a 3x3 matrix */ /* src1[] forms the colums, and src2[] forms the rows in the result */ @@ -13885,6 +13926,13 @@ void icmScale3(double out[3], double in[3], double 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]; @@ -13904,6 +13952,47 @@ void icmClip3(double out[3], double in[3]) { } } +/* 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 */ +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]); @@ -13917,6 +14006,17 @@ int icmNormalize3(double out[3], double in[3], double len) { return 0; } +/* 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 a vector define by two points */ double icmNorm33sq(double in1[3], double in0[3]) { int j; @@ -14062,6 +14162,199 @@ void icmRotMat(double m[3][3], double s[3], double t[3]) { #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; +} + +/* - - - - - - - - - - - - - - - - - - - - - - - - */ + +/* + + mat in out + +[ ] [] [] +[ ] [] [] +[ ] * [] => [] +[ ] [] [] +[ ] [] [] + + */ + +/* Multiply 5 array by 5x5 transform matrix */ +void icmMulBy5x5(double out[5], double mat[5][5], double in[5]) { + int i, j; + double tt[5]; + + for (i = 0; i < 5; i++) { + tt[i] = 0.0; + for (j = 0; j < 5; j++) + tt[i] += mat[i][j] * in[j]; + } + + for (i = 0; i < 5; i++) + out[i] = tt[i]; +} + +/* Transpose a 5x5 matrix */ +void icmTranspose5x5(double out[5][5], double in[5][5]) { + int i, j; + if (out != in) { + for (i = 0; i < 5; i++) + for (j = 0; j < 5; j++) + out[i][j] = in[j][i]; + } else { + double tt[5][5]; + for (i = 0; i < 5; i++) + for (j = 0; j < 5; j++) + tt[i][j] = in[j][i]; + for (i = 0; i < 5; i++) + for (j = 0; j < 5; j++) + out[i][j] = tt[i][j]; + } +} + +/* Clip a vector to the range 0.0 .. 1.0 */ +/* and return any clipping margine */ +double icmClip5marg(double out[5], double in[5]) { + int j; + double tt, marg = 0.0; + for (j = 0; j < 5; 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; +} + + +/* Multiply 6 array by 6x6 transform matrix */ +void icmMulBy6x6(double out[6], double mat[6][6], double in[6]) { + int i, j; + double tt[6]; + + for (i = 0; i < 6; i++) { + tt[i] = 0.0; + for (j = 0; j < 6; j++) + tt[i] += mat[i][j] * in[j]; + } + + for (i = 0; i < 6; i++) + out[i] = tt[i]; +} + +/* Transpose a 6x6 matrix */ +void icmTranspose6x6(double out[6][6], double in[6][6]) { + int i, j; + if (out != in) { + for (i = 0; i < 6; i++) + for (j = 0; j < 6; j++) + out[i][j] = in[j][i]; + } else { + double tt[6][6]; + for (i = 0; i < 6; i++) + for (j = 0; j < 6; j++) + tt[i][j] = in[j][i]; + for (i = 0; i < 6; i++) + for (j = 0; j < 6; j++) + out[i][j] = tt[i][j]; + } +} + +/* Clip a vector to the range 0.0 .. 1.0 */ +/* and return any clipping margine */ +double icmClip6marg(double out[6], double in[6]) { + int j; + double tt, marg = 0.0; + for (j = 0; j < 6; 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; +} + /* - - - - - - - - - - - - - - - - - - - - - - - - */ /* Multiply 2 array by 2x2 transform matrix */ void icmMulBy2x2(double out[2], double mat[2][2], double in[2]) { @@ -14207,11 +14500,11 @@ double ve_0[3] /* Second point on line */ 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 */ + double val; /* Vector length squared */ double a; /* Parameter value */ icmSub3(va, la1, la0); /* Line vector */ - val = icmNorm3(va); /* Vector length */ + val = icmNorm3sq(va); /* Vector length squared */ if (val < 1e-12) return 1; @@ -14345,7 +14638,7 @@ double icmPlaneDist3(double eq[4], double p[3]) { } /* - - - - - - - - - - - - - - - - - - - - - - - - */ -/* Given 2 2D points, compute a plane equation. */ +/* 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 */ @@ -14375,7 +14668,7 @@ int icmPlaneEqn2(double eq[3], double p0[2], double p1[2]) { return 0; } -/* Given a 2D point and a plane equation, return the signed */ +/* 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]) { @@ -14388,6 +14681,70 @@ double icmPlaneDist2(double eq[3], double p[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]) { @@ -14411,6 +14768,41 @@ int icmLineIntersect2(double res[2], double p1[2], double p2[2], double p3[2], d return 0; } +/* Given two finite 2D lines define by 4 points, compute their paramaterized intersection. */ +/* aprm may be NULL */ +/* 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; +} + /* - - - - - - - - - - - - - - - - - - - - - - - - */ /* CIE Y (range 0 .. 1) to perceptual CIE 1976 L* (range 0 .. 100) */ @@ -14497,26 +14889,122 @@ icmLab2XYZ(icmXYZNumber *w, double *out, double *in) { 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 2degree LMS to XYZ matrix instead of Hunt-Pointer-Estevez + * Von Kries chromatic adapation in LMS 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!). + */ + +/* 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; -/* LCh to Lab */ + /* 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 = 3.14159265359/180.0 * in[2]; + 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 Luv too) */ +/* 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/3.14159265359) * atan2(in[2], in[1]); + h = (180.0/M_PI) * atan2(in[2], in[1]); h = (h < 0.0) ? h + 360.0 : h; out[0] = in[0]; @@ -14531,8 +15019,8 @@ extern ICCLIB_API void icmXYZ2Yxy(double *out, double *in) { if (sum < 1e-9) { Y = 0.0; - y = 0.333333333; - x = 0.333333333; + y = 1.0/3.0; + x = 1.0/3.0; } else { Y = in[1]; x = in[0]/sum; @@ -14543,6 +15031,22 @@ extern ICCLIB_API void icmXYZ2Yxy(double *out, double *in) { 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]; @@ -14560,6 +15064,21 @@ extern ICCLIB_API void icmYxy2XYZ(double *out, double *in) { } } +/* 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 L*u*v* */ extern ICCLIB_API void icmXYZ2Luv(icmXYZNumber *w, double *out, double *in) { @@ -14616,7 +15135,7 @@ extern ICCLIB_API void icmLuv2XYZ(icmXYZNumber *w, double *out, double *in) { } /* CIE XYZ to perceptual CIE 1976 UCS diagram Yu'v'*/ -/* (Yu'v' is a better chromaticity space than Yxy) */ +/* (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; @@ -14657,6 +15176,27 @@ extern ICCLIB_API void icm1976UCS2XYZ(double *out, double *in) { 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.) */ @@ -14862,6 +15402,16 @@ double icmLabDEsq(double *Lab0, double *Lab1) { 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; @@ -14872,6 +15422,26 @@ extern ICCLIB_API double icmXYZLabDE(icmXYZNumber *w, double *in0, double *in1) 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) */ @@ -14962,8 +15532,8 @@ double icmCIE2Ksq(double *Lab0, double *Lab1) { /* 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/3.14159265358979 * (xx)) -#define DEG2RAD(xx) (3.14159265358979/180.0 * (xx)) +#define RAD2DEG(xx) (180.0/M_PI * (xx)) +#define DEG2RAD(xx) (M_PI/180.0 * (xx)) /* Compute Cromanance and Hue angles */ { @@ -18487,7 +19057,7 @@ static icmLuBase* icc_get_luobj ( if (intent == icmDefaultIntent) intent = icPerceptual; /* Make this the default */ - switch (intent) { + switch ((int)intent) { case icAbsoluteColorimetric: ttag = icSigAToB1Tag; fbtag = icSigAToB0Tag; @@ -18588,7 +19158,7 @@ static icmLuBase* icc_get_luobj ( if (intent == icmDefaultIntent) intent = icPerceptual; /* Make this the default */ - switch (intent) { + switch ((int)intent) { case icAbsoluteColorimetric: ttag = icSigBToA1Tag; fbtag = icSigBToA0Tag; @@ -18708,7 +19278,7 @@ static icmLuBase* icc_get_luobj ( if (intent == icmDefaultIntent) intent = icPerceptual; /* Make this the default */ - switch (intent) { + switch ((int)intent) { case icRelativeColorimetric: case icAbsoluteColorimetric: ttag = icSigAToB1Tag; @@ -18778,7 +19348,7 @@ static icmLuBase* icc_get_luobj ( if (intent == icmDefaultIntent) intent = icPerceptual; /* Make this the default */ - switch (intent) { + switch ((int)intent) { case icRelativeColorimetric: case icAbsoluteColorimetric: ttag = icSigBToA1Tag; @@ -18852,7 +19422,7 @@ static icmLuBase* icc_get_luobj ( return NULL; } #else /* Be more forgiving */ - switch (intent) { + switch ((int)intent) { case icAbsoluteColorimetric: case icmAbsolutePerceptual: /* Special icclib intent */ case icmAbsoluteSaturation: /* Special icclib intent */ @@ -18879,7 +19449,7 @@ static icmLuBase* icc_get_luobj ( case icmPreview: /* PCS to PCS */ - switch (intent) { + switch ((int)intent) { case icRelativeColorimetric: ttag = icSigPreview1Tag; break; @@ -19268,6 +19838,40 @@ icmAlloc *al /* Memory allocator */ } +/* ---------------------------------------------------------- */ + +/* 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. */ -- cgit v1.2.3