33 files changed, 784 insertions, 55 deletions

diff --git a/icc/ClayRGB1998.icm b/icc/ClayRGB1998.icm
index 4c8ecb3..4c8ecb3 100644..100755
--- a/icc/ClayRGB1998.icm
+++ b/icc/ClayRGB1998.icm
diff --git a/icc/EBU3213_PAL.icm b/icc/EBU3213_PAL.icm
index b0e67bc..b0e67bc 100644..100755
--- a/icc/EBU3213_PAL.icm
+++ b/icc/EBU3213_PAL.icm
diff --git a/icc/Jamfile b/icc/Jamfile
index 4bedbba..4bedbba 100644..100755
--- a/icc/Jamfile
+++ b/icc/Jamfile
diff --git a/icc/License4.txt b/icc/License4.txt
index 61be41a..61be41a 100644..100755
--- a/icc/License4.txt
+++ b/icc/License4.txt
diff --git a/icc/Makefile b/icc/Makefile
index 85f19be..85f19be 100644..100755
--- a/icc/Makefile
+++ b/icc/Makefile
diff --git a/icc/Makefile.IBMNT b/icc/Makefile.IBMNT
index 26a2bd2..26a2bd2 100644..100755
--- a/icc/Makefile.IBMNT
+++ b/icc/Makefile.IBMNT
diff --git a/icc/Makefile.OSX b/icc/Makefile.OSX
index fb1c88f..fb1c88f 100644..100755
--- a/icc/Makefile.OSX
+++ b/icc/Makefile.OSX
diff --git a/icc/Makefile.UNIX b/icc/Makefile.UNIX
index 91d6005..91d6005 100644..100755
--- a/icc/Makefile.UNIX
+++ b/icc/Makefile.UNIX
diff --git a/icc/Makefile.WNT b/icc/Makefile.WNT
index 4ca291d..4ca291d 100644..100755
--- a/icc/Makefile.WNT
+++ b/icc/Makefile.WNT
diff --git a/icc/ProPhoto.icm b/icc/ProPhoto.icm
index cb6ee56..cb6ee56 100644..100755
--- a/icc/ProPhoto.icm
+++ b/icc/ProPhoto.icm
diff --git a/icc/ProPhotoLin.icm b/icc/ProPhotoLin.icm
index 28369aa..28369aa 100644..100755
--- a/icc/ProPhotoLin.icm
+++ b/icc/ProPhotoLin.icm
diff --git a/icc/Readme.txt b/icc/Readme.txt
index 2b5fd53..2b5fd53 100644..100755
--- a/icc/Readme.txt
+++ b/icc/Readme.txt
diff --git a/icc/Rec2020.icm b/icc/Rec2020.icm
index 0decaf6..0decaf6 100644..100755
--- a/icc/Rec2020.icm
+++ b/icc/Rec2020.icm
diff --git a/icc/Rec709.icm b/icc/Rec709.icm
index abf0a63..abf0a63 100644..100755
--- a/icc/Rec709.icm
+++ b/icc/Rec709.icm
diff --git a/icc/SMPTE431_P3.icm b/icc/SMPTE431_P3.icm
index 8055849..8055849 100644..100755
--- a/icc/SMPTE431_P3.icm
+++ b/icc/SMPTE431_P3.icm
diff --git a/icc/SMPTE_RP145_NTSC.icm b/icc/SMPTE_RP145_NTSC.icm
index 7fc78fb..7fc78fb 100644..100755
--- a/icc/SMPTE_RP145_NTSC.icm
+++ b/icc/SMPTE_RP145_NTSC.icm
diff --git a/icc/afiles b/icc/afiles
index 3ee789a..3ee789a 100644..100755
--- a/icc/afiles
+++ b/icc/afiles
diff --git a/icc/icc.c b/icc/icc.c
index b239060..de8c539 100644..100755
--- 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;
@@ -14063,6 +14163,199 @@ void icmRotMat(double m[3][3], double s[3], double t[3]) {
 }
 
 /* - - - - - - - - - - - - - - - - - - - - - - - - */
+
+/* 
+
+  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]) {
 	double tt[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;
@@ -19269,6 +19839,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. */
 char *icmPiv(int di, int *p) {
diff --git a/icc/icc.h b/icc/icc.h
index e0bf88e..1e2caf8 100644..100755
--- a/icc/icc.h
+++ b/icc/icc.h
@@ -33,7 +33,7 @@
 #define ICCLIB_VERSION 0x020021
 #define ICCLIB_VERSION_STR "2.21"
 
-#undef ENABLE_V4		/* V4 is not fully implemented */
+#undef ENABLE_V4		/* [und] V4 is not fully implemented, but enables parsing */
 
 /*
  *  Note XYZ scaling to 1.0, not 100.0
@@ -360,6 +360,9 @@ typedef int icmSig;	/* Otherwise un-enumerated 4 byte signature */
 /* A monochrome CIE Y space */
 #define icmSigYData ((icColorSpaceSignature) icmMakeTag('Y',' ',' ',' '))
 
+/* A modern take on Lab */
+#define icmSigLptData ((icColorSpaceSignature) icmMakeTag('L','p','t',' '))
+
 
 /* Pseudo Color Space Signatures - just used within icclib */
 
@@ -1778,13 +1781,23 @@ extern ICCLIB_API unsigned int icmCSSig2chanNames( icColorSpaceSignature sig, ch
 #define icmSet3(d_ary, s_val) ((d_ary)[0] = (s_val), (d_ary)[1] = (s_val), \
                               (d_ary)[2] = (s_val))
 
+/* Copy a 2 vector */
+#define icmCpy2(d_ary, s_ary) ((d_ary)[0] = (s_ary)[0], (d_ary)[1] = (s_ary)[1])
+
 /* Copy a 3 vector */
 #define icmCpy3(d_ary, s_ary) ((d_ary)[0] = (s_ary)[0], (d_ary)[1] = (s_ary)[1], \
                               (d_ary)[2] = (s_ary)[2])
 
+/* Copy a 4 vector */
+#define icmCpy4(d_ary, s_ary) ((d_ary)[0] = (s_ary)[0], (d_ary)[1] = (s_ary)[1], \
+                              (d_ary)[2] = (s_ary)[2], (d_ary)[3] = (s_ary)[3])
+
 /* Clamp a 3 vector to be +ve */
 void icmClamp3(double out[3], double in[3]);
 
+/* Invert a 3 vector */
+void icmInv3(double out[3], double in[3]);
+
 /* Add two 3 vectors */
 void icmAdd3(double out[3], double in1[3], double in2[3]);
 
@@ -1805,6 +1818,9 @@ void icmMul3(double out[3], double in1[3], double in2[3]);
 
 #define ICMMUL3(o, i, j) ((o)[0] = (i)[0] * (j)[0], (o)[1] = (i)[1] * (j)[1], (o)[2] = (i)[2] * (j)[2])
 
+/* Take values to power */
+void icmPow3(double out[3], double in[3], double p);
+
 /* Take absolute of a 3 vector */
 void icmAbs3(double out[3], double in[3]);
 
@@ -1816,6 +1832,8 @@ double icmDot3(double in1[3], double in2[3]);
 /* Compute the cross product of two 3D vectors, out = in1 x in2 */
 void icmCross3(double out[3], double in1[3], double in2[3]);
 
+#define icmNorm2(i) sqrt((i)[0] * (i)[0] + (i)[1] * (i)[1])
+
 /* Compute the norm squared (length squared) of a 3 vector */
 double icmNorm3sq(double in[3]);
 
@@ -1831,15 +1849,29 @@ void icmScale3(double out[3], double in[3], double rat);
 
 #define ICMSCALE3(o, i, j) ((o)[0] = (i)[0] * (j), (o)[1] = (i)[1] * (j), (o)[2] = (i)[2] * (j))
 
+/* Scale a 3 vector by the given ratio and add it */
+void icmScaleAdd3(double out[3], double in2[3], double in1[3], double rat);
+
 /* Compute a blend between in0 and in1 */
 void icmBlend3(double out[3], double in0[3], double in1[3], double bf);
 
 /* Clip a vector to the range 0.0 .. 1.0 */
 void icmClip3(double out[3], double in[3]);
 
+/* Clip a vector to the range 0.0 .. 1.0 */
+/* and return nz if clipping occured */
+int icmClip3sig(double out[3], double in[3]);
+
+/* Clip a vector to the range 0.0 .. 1.0 */
+/* and return any clipping margine */
+double icmClip3marg(double out[3], double in[3]);
+
 /* Normalise a 3 vector to the given length. Return nz if not normalisable */
 int icmNormalize3(double out[3], double in[3], double len);
 
+/* Compute the norm (length) of of a vector defined by two points */
+double icmNorm22(double in1[2], double in0[2]);
+
 /* Compute the norm squared (length squared) of a vector defined by two points */
 double icmNorm33sq(double in1[3], double in0[3]);
 
@@ -1854,6 +1886,8 @@ void icmScale33(double out[3], double in1[3], double in0[3], double rat);
 /* Return nz if not normalisable */
 int icmNormalize33(double out[3], double in1[3], double in0[3], double len);
 
+/* Set a 3x3 matrix to a value */
+void icmSetVal3x3(double mat[3][3], double val);
 
 /* Set a 3x3 matrix to unity */
 void icmSetUnity3x3(double mat[3][3]);
@@ -1861,6 +1895,9 @@ void icmSetUnity3x3(double mat[3][3]);
 /* Copy a 3x3 transform matrix */
 void icmCpy3x3(double out[3][3], double mat[3][3]);
 
+/* Add a 3x3 transform matrix to another */
+void icmAdd3x3(double dst[3][3], double src1[3][3], double src2[3][3]);
+
 /* Scale each element of a 3x3 transform matrix */
 void icmScale3x3(double dst[3][3], double src[3][3], double scale);
 
@@ -1868,10 +1905,6 @@ void icmScale3x3(double dst[3][3], double src[3][3], double scale);
 /* Organization is mat[out][in] */
 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]);
-
 /* 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]);
@@ -1937,6 +1970,44 @@ double icmPlaneDist3(double eq[4], double p[3]);
 
 /* - - - - - - - - - - - - - - - - - - - - - - - */
 
+/* Multiply 4 vector by 4x4 transform matrix */
+/* Organization is mat[out][in] */
+void icmMulBy4x4(double out[4], double mat[4][4], double in[4]);
+
+/* Transpose a 4x4 matrix */
+void icmTranspose4x4(double out[4][4], double in[4][4]);
+
+/* Clip a vector to the range 0.0 .. 1.0 */
+/* and return any clipping margine */
+double icmClip4marg(double out[4], double in[4]);
+
+/* - - - - - - - - - - - - - - - - - - - - - - - */
+
+/* Multiply 5 vector by 5x5 transform matrix */
+/* Organization is mat[out][in] */
+void icmMulBy5x5(double out[5], double mat[5][5], double in[5]);
+
+/* Transpose a 5x5 matrix */
+void icmTranspose5x5(double out[5][5], double in[5][5]);
+
+/* Clip a vector to the range 0.0 .. 1.0 */
+/* and return any clipping margine */
+double icmClip5marg(double out[5], double in[5]);
+
+
+/* Multiply 6 vector by 6x6 transform matrix */
+/* Organization is mat[out][in] */
+void icmMulBy6x6(double out[6], double mat[6][6], double in[6]);
+
+/* Transpose a 6x6 matrix */
+void icmTranspose6x6(double out[6][6], double in[6][6]);
+
+/* Clip a vector to the range 0.0 .. 1.0 */
+/* and return any clipping margine */
+double icmClip6marg(double out[6], double in[6]);
+
+/* - - - - - - - - - - - - - - - - - - - - - - - */
+
 /* Given 2 2D points, compute a plane 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. */
@@ -1947,10 +2018,30 @@ int icmPlaneEqn2(double eq[3], double p0[2], double p1[2]);
 /* distance from the plane */
 double icmPlaneDist2(double eq[3], double p[2]);
 
+/* 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]);
+
+/* 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]);
+
+
+/* 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]);
+
 /* Given two infinite 2D lines define by two pairs of 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]);
 
+/* Given two finite 2D lines define by 4 points, compute their paramaterized intersection. */
+/* aprm may be NULL */
+/* Return nz if there is no intersection (lines are parallel or do not cross in length) */
+int icmParmLineIntersect2(double ares[2], double aprm[2], double p1[2], double p2[2], double p3[2], double p4[2]);
+
 /* Multiply 2 array by 2x2 transform matrix */
 void icmMulBy2x2(double out[2], double mat[2][2], double in[2]);
 
@@ -1983,18 +2074,32 @@ extern ICCLIB_API void icmXYZ2Lab(icmXYZNumber *w, double *out, double *in);
 /* Perceptual Lab to CIE XYZ */
 extern ICCLIB_API void icmLab2XYZ(icmXYZNumber *w, double *out, double *in);
 
-/* LCh to Lab */
+/* CIE XYZ to perceptual Lpt */
+extern ICCLIB_API void icmXYZ2Lpt(icmXYZNumber *w, double *out, double *in);
+
+/* Perceptual Lpt to CIE XYZ */
+extern ICCLIB_API void icmLpt2XYZ(icmXYZNumber *w, double *out, double *in);
+
+/* LCh to Lab (general) */
 extern ICCLIB_API void icmLCh2Lab(double *out, double *in);
 
-/* Lab to LCh */
+/* Lab to LCh (general) */
 extern ICCLIB_API void icmLab2LCh(double *out, double *in);
 
+
 /* XYZ to Yxy */
 extern ICCLIB_API void icmXYZ2Yxy(double *out, double *in);
 
 /* Yxy to XYZ */
 extern ICCLIB_API void icmYxy2XYZ(double *out, double *in);
 
+/* XYZ to xy */
+extern ICCLIB_API void icmXYZ2xy(double *out, double *in);
+
+/* Y & xy to XYZ */
+extern ICCLIB_API void icmY_xy2XYZ(double *out, double *xy, double Y);
+
+
 /* CIE XYZ to perceptual Luv */
 extern ICCLIB_API void icmXYZ2Luv(icmXYZNumber *w, double *out, double *in);
 
@@ -2003,12 +2108,16 @@ 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);
 
 /* Perceptual CIE 1976 UCS diagram Yu'v' to CIE XYZ */
 extern ICCLIB_API void icm1976UCS2XYZ(double *out, double *in);
 
+/* 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);
+
 
 /* CIE XYZ to perceptual CIE 1960 UCS */
 /* (This was obsoleted by the 1976UCS, but is still used */
@@ -2072,9 +2181,18 @@ extern ICCLIB_API double icmLabDE(double *in0, double *in1);
 /* Return the normal Delta E squared, given two Lab values */
 extern ICCLIB_API double icmLabDEsq(double *in0, double *in1);
 
+/* Return the normal Delta E squared given two XYZ values */
+extern ICCLIB_API double icmXYZLabDEsq(icmXYZNumber *w, double *in0, double *in1);
+
 /* Return the normal Delta E given two XYZ values */
 extern ICCLIB_API double icmXYZLabDE(icmXYZNumber *w, double *in0, double *in1);
 
+/* Return the normal Delta E squared given two XYZ values */
+extern ICCLIB_API double icmXYZLptDEsq(icmXYZNumber *w, double *in0, double *in1);
+
+/* Return the normal Delta E given two XYZ values */
+extern ICCLIB_API double icmXYZLptDE(icmXYZNumber *w, double *in0, double *in1);
+
 
 /* Return the CIE94 Delta E color difference measure for two Lab values */
 extern ICCLIB_API double icmCIE94(double *in0, double *in1);
@@ -2249,6 +2367,13 @@ double icmDICOM_bwd(double L);
 
 
 /* ---------------------------------------------------------- */
+/* Some utility functions */
+
+/* Convert an angle in radians into chromatic RGB values */
+/* in a simple geometric fashion with 0 = Red */
+void icmRad2RGB(double rgb[3], double ang);
+
+/* ---------------------------------------------------------- */
 /* Print an int vector to a string. */
 /* Returned static buffer is re-used every 5 calls. */
 char *icmPiv(int di, int *p);
diff --git a/icc/iccV42.h b/icc/iccV42.h
index a84f79d..a84f79d 100644..100755
--- a/icc/iccV42.h
+++ b/icc/iccV42.h
diff --git a/icc/iccdump.c b/icc/iccdump.c
index 8574258..8574258 100644..100755
--- a/icc/iccdump.c
+++ b/icc/iccdump.c
diff --git a/icc/icclu.c b/icc/icclu.c
index 5c732d0..5c732d0 100644..100755
--- a/icc/icclu.c
+++ b/icc/icclu.c
diff --git a/icc/iccrw.c b/icc/iccrw.c
index f33c164..f33c164 100644..100755
--- a/icc/iccrw.c
+++ b/icc/iccrw.c
diff --git a/icc/iccstd.c b/icc/iccstd.c
index 6629250..6629250 100644..100755
--- a/icc/iccstd.c
+++ b/icc/iccstd.c
diff --git a/icc/icctest.c b/icc/icctest.c
index 647ba07..647ba07 100644..100755
--- a/icc/icctest.c
+++ b/icc/icctest.c
diff --git a/icc/lab2lab.icm b/icc/lab2lab.icm
index 39a5911..39a5911 100644..100755
--- a/icc/lab2lab.icm
+++ b/icc/lab2lab.icm
diff --git a/icc/log.txt b/icc/log.txt
index b531f5c..b531f5c 100644..100755
--- a/icc/log.txt
+++ b/icc/log.txt
diff --git a/icc/lutest.c b/icc/lutest.c
index fb05ea9..fb05ea9 100644..100755
--- a/icc/lutest.c
+++ b/icc/lutest.c
diff --git a/icc/mcheck.c b/icc/mcheck.c
index 51368ca..51368ca 100644..100755
--- a/icc/mcheck.c
+++ b/icc/mcheck.c
diff --git a/icc/mkDispProf.c b/icc/mkDispProf.c
index b7c32c6..b7c32c6 100644..100755
--- a/icc/mkDispProf.c
+++ b/icc/mkDispProf.c
diff --git a/icc/sRGB.icm b/icc/sRGB.icm
index 59b4507..59b4507 100644..100755
--- a/icc/sRGB.icm
+++ b/icc/sRGB.icm
diff --git a/icc/testDE2K.c b/icc/testDE2K.c
index 11a4201..11a4201 100644..100755
--- a/icc/testDE2K.c
+++ b/icc/testDE2K.c
diff --git a/icc/todo.txt b/icc/todo.txt
index bdbaa5e..bdbaa5e 100644..100755
--- a/icc/todo.txt
+++ b/icc/todo.txt