/* Copyright (C) 2005-2011 Fabio Riccardi */
package com.lightcrafts.utils;
import Jama.Matrix;
import java.awt.image.DataBuffer;
import java.awt.color.ICC_ProfileRGB;
import com.lightcrafts.jai.JAIContext;
import com.lightcrafts.utils.Color.BlackBody;
public class ColorScience {
static final float[][] rgbXYZ;
static final float[] wtptXYZ;
static final float whitePointTemperature;
public static final float Wr;
public static final float Wg;
public static final float Wb;
public static final float W[];
static final float[][] Cxy;
public static final float[][] XYZToRGBMat;
public static final float[][] RGBToXYZMat;
public static float[] XYZ2xy(float[] XYZ) {
return new float[] {
XYZ[0] / (XYZ[0] + XYZ[1] + XYZ[2]),
XYZ[1] / (XYZ[0] + XYZ[1] + XYZ[2])
};
}
public static float[] xy2XYZ(float[] xy) {
return new float[] {
xy[0]/xy[1],
1,
(1 - xy[0] - xy[1]) / xy[1]
};
}
public static float[] XYZ2xyY(float[] XYZ) {
return new float[] {
XYZ[0] / (XYZ[0] + XYZ[1] + XYZ[2]),
XYZ[1] / (XYZ[0] + XYZ[1] + XYZ[2]),
XYZ[1]
};
}
public static float[] xyY2XYZ(float[] xyY) {
return new float[] {
xyY[2] * xyY[0]/xyY[1],
xyY[2],
xyY[2] * (1 - xyY[0] - xyY[1]) / xyY[1]
};
}
static {
ICC_ProfileParameters pp = new ICC_ProfileParameters((ICC_ProfileRGB) JAIContext.linearProfile);
rgbXYZ = pp.rgbXYZ;
wtptXYZ = pp.wtptXYZ;
Cxy = pp.Cxy;
whitePointTemperature = pp.whitePointTemperature;
W = pp.W;
Wr = pp.W[0];
Wg = pp.W[1];
Wb = pp.W[2];
RGBToXYZMat = pp.RGBToXYZMat;
XYZToRGBMat = pp.XYZToRGBMat;
}
// TODO: we have to finish cleaning up this.
public static class ICC_ProfileParameters {
public float rgbXYZ[][] = new float[3][3];
public float wtptXYZ[];
public float Cxy[][];
public float whitePointTemperature;
public float W[];
public float RGBToXYZMat[][];
public float XYZToRGBMat[][];
public ICC_ProfileParameters(ICC_ProfileRGB profile) {
// Extract the rgbXYZ from the current linear color profile
float[][] rgbXYZt = (profile).getMatrix();
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
rgbXYZ[i][j] = rgbXYZt[j][i];
// Same for the white point
wtptXYZ = profile.getMediaWhitePoint();
// compute the xy coordinates of the workspace primaries form XYZ
float Cr = rgbXYZ[0][0] + rgbXYZ[0][1] + rgbXYZ[0][2];
float Cg = rgbXYZ[1][0] + rgbXYZ[1][1] + rgbXYZ[1][2];
float Cb = rgbXYZ[2][0] + rgbXYZ[2][1] + rgbXYZ[2][2];
Cxy = new float[][]{
{rgbXYZ[0][0] / Cr, rgbXYZ[0][1] / Cr},
{rgbXYZ[1][0] / Cg, rgbXYZ[1][1] / Cg},
{rgbXYZ[2][0] / Cb, rgbXYZ[2][1] / Cb}
};
float Zr[] = new float[] {Cxy[0][0], Cxy[0][1], 1 - Cxy[0][0] - Cxy[0][1]};
float Zg[] = new float[] {Cxy[1][0], Cxy[1][1], 1 - Cxy[1][0] - Cxy[1][1]};
float Zb[] = new float[] {Cxy[2][0], Cxy[2][1], 1 - Cxy[2][0] - Cxy[2][1]};
// wtptXYZ -> wtptxy
final float[] wtptxy = XYZ2xy(wtptXYZ);
whitePointTemperature = CCTX(wtptxy[0]);
W = W(whitePointTemperature, Cxy);
RGBToXYZMat = new Matrix (new float[][] {
mul(W[0] / Zr[1], Zr),
mul(W[1] / Zg[1], Zg),
mul(W[2] / Zb[1], Zb)
}).transpose().getArrayFloat();
XYZToRGBMat = new Matrix(RGBToXYZMat).inverse().getArrayFloat();
}
}
// Return x,y coordinates of the CIE D illuminant specified by t, good between 4000K and 25000K
/* static float[] D(float t) {
double x;
if (t <= 7000)
x = -4.6070E9 * Math.pow(t, -3.) + 2.9678E6 * Math.pow(t, -2.) + 0.09911E3 * Math.pow(t, -1.) + 0.244063;
else
x = -2.0064E9 * Math.pow(t, -3.) + 1.9018E6 * Math.pow(t, -2.) + 0.24748E3 * Math.pow(t, -1.) + 0.237040;
double y = -3.000 * Math.pow(x, 2.) + 2.870 * x - 0.275;
return new float[]{(float) x, (float) y};
} */
// Return x,y coordinates of the CIE D illuminant specified by t, good between 1000K and 40000K
static float[] D(float t) {
return BlackBody.xy(t);
}
public static float CCTX(float x) {
return BlackBody.t(x);
}
public static float[][] RGBtoZYX() {
double mdata[][] = new double[3][3];
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
mdata[i][j] = rgbXYZ[i][j];
Matrix XYZRGB = new Matrix(mdata);
float[][] S = new Matrix(new double[][] {{wtptXYZ[0], wtptXYZ[1], wtptXYZ[2]}}).times(XYZRGB.inverse()).getArrayFloat();
return new float[][] {
{S[0][0] * rgbXYZ[0][0], S[0][0] * rgbXYZ[0][1], S[0][0] * rgbXYZ[0][2]},
{S[0][1] * rgbXYZ[1][0], S[0][1] * rgbXYZ[1][1], S[0][1] * rgbXYZ[1][2]},
{S[0][2] * rgbXYZ[2][0], S[0][2] * rgbXYZ[2][1], S[0][2] * rgbXYZ[2][2]}
};
}
static float[] mul(float c, float[] v) {
float[] r = new float[v.length];
for (int i = 0; i < v.length; i++)
r[i] = c * v[i];
return r;
}
public static float[] W(float T) {
return W(T, Cxy);
}
public static float[] W(float T, float Cxy[][]) {
float DT[] = D(T);
// Compute z-coordinates
float R[] = {Cxy[0][0], Cxy[0][1], 1 - Cxy[0][0] - Cxy[0][1]};
float G[] = {Cxy[1][0], Cxy[1][1], 1 - Cxy[1][0] - Cxy[1][1]};
float B[] = {Cxy[2][0], Cxy[2][1], 1 - Cxy[2][0] - Cxy[2][1]};
float W[] = {DT[0], DT[1], 1 - DT[0] - DT[1]};
// Compute luminance weights for the primaries using the white point
Matrix RGB = vec(R, G, B);
Matrix WGB = vec(W, G, B);
Matrix RWB = vec(R, W, B);
Matrix RGW = vec(R, G, W);
float rgbDet = (float) RGB.det();
return new float[] {(R[1] * (float) WGB.det()) / (W[1] * rgbDet),
(G[1] * (float) RWB.det()) / (W[1] * rgbDet),
(B[1] * (float) RGW.det()) / (W[1] * rgbDet)};
}
/*
8.2 - HSI, HSL, HSV, and related color spaces
The representation of the colors in the RGB and CMY(K) color spaces are designed
for specific devices. But for a human observer, they are not useful definitions.
For user interfaces a more intuitive color space, designed for the way we actually
think about color is to be preferred. Such a color space is HSI; Hue, Saturation and
Intensity, which can be thought of as a RGB cube tipped up onto one corner. The line
from RGB=min to RGB=max becomes verticle and is the intensity axis. The position of a
point on the circumference of a circle around this axis is the hue and the saturation
is the radius from the central intensity axis to the color.
Green
/\
/ \ ^
/V=1 x \ \ Hue (angle, so that Hue(Red)=0,
Blue -------------- Red Hue(Green)=120, and Hue(blue)=240 deg)
\ | /
\ |-> Saturation (distance from the central axis)
\ | /
\ | /
\ | /
\ |/
V=0 x (Intensity=0 at the top of the apex and =1 at the base of the cone)
The big disadvantage of this model is the conversion which is mainly because the hue is
expressed as an angle. The transforms are given below:
Hue = (Alpha-arctan((Red-intensity)*(3^0.5)/(Green-Blue)))/(2*PI)
with { Alpha=PI/2 if Green>Blue
{ Aplha=3*PI/2 if Green<Blue
{ Hue=1 if Green=Blue
Saturation = (Red^2+Green^2+Blue^2-Red*Green-Red*Blue-Blue*Green)^0.5
Intensity = (Red+Green+Blue)/3
Note that you have to compute Intensity *before* Hue. If not, you must assume that:
Hue = (Alpha-arctan((2*Red-Green-Blue)/((Green-Blue)*(3^0.5))))/(2*PI).
I assume that H, S, L, R, G, and B are within the range of [0;1]. Another point of view of
this cone is to project the coordinates onto the base. The 2D projection is:
Red: (1;0)
Green: (cos(120 deg);sin(120 deg)) = (-0.5; 0.866)
Blue: (cos(240 deg);sin(240 deg)) = (-0.5;-0.866)
Now you need intermediate coordinates:
x = Red-0.5*(Green+Blue)
y = 0.866*(Green-Blue)
Finally, you have:
Hue = arctan2(x,y)/(2*PI) ; Just one formula, always in the correct quadrant
Saturation = (x^2+y^2)^0.5
Intensity = (Red+Green+Blue)/3
The intermediate coordinates YST { I, x, y } are a cool substitute for HSI in most calculations
RGBtoYST = new float[][]{{Wr, Wg, Wb, 0},
{1, -.5, -.5, 0},
{0, .5 * Math.sqrt(3), -.5 * Math.sqrt(3), 0},
{0, 0, 0, 1}};
YSTtoRGB = new float[][] {{1, (Wg+Wb), (Wb - Wg) / Math.sqrt(3), 0},
{1, (-Wr), (2*Wb + Wr) / Math.sqrt(3), 0},
{1, (-Wr), -(2*Wg + Wr) / Math.sqrt(3), 0}};
The code below uses a normalized version of the intermediate coordinates YST to make sure it fits in the
dynamic range of a positive integer representation
*/
public static abstract class LinearTransform {
static double[][] scaleTransform(double[][] t, int dataType) {
switch (dataType) {
case DataBuffer.TYPE_BYTE:
t[1][3] = 0xFF / 2.;
t[2][3] = 0xFF / 2.;
break;
case DataBuffer.TYPE_SHORT:
t[1][3] = Short.MAX_VALUE / 2.;
t[2][3] = Short.MAX_VALUE / 2.;
break;
case DataBuffer.TYPE_USHORT:
t[1][3] = 0xFFFF / 2.;
t[2][3] = 0xFFFF / 2.;
break;
case DataBuffer.TYPE_INT:
t[1][3] = Integer.MAX_VALUE / 2.;
t[2][3] = Integer.MAX_VALUE / 2.;
break;
case DataBuffer.TYPE_FLOAT:
case DataBuffer.TYPE_DOUBLE:
t[1][3] = 0.5;
t[2][3] = 0.5;
break;
}
return t;
}
abstract double [][] getTransform();
public double[][] fromRGB(int dataType) {
double[][] t = scaleTransform(getTransform(), dataType);
return strip(t);
}
public double[][] toRGB(int dataType) {
double[][] t = scaleTransform(getTransform(), dataType);
return strip(new Matrix(t).inverse().getArray());
}
}
public static class XYZ extends LinearTransform {
double [][] getTransform() {
return new double[][]{
{Wr, Wg, Wb, 0},
{0.5, -0.25, -0.25, 0.5},
{0, 0.5, -0.5, 0.5},
{0, 0, 0, 1}
};
}
}
public static class YST extends LinearTransform {
double [][] getTransform() {
return new double[][]{
{Wr, Wg, Wb, 0},
{0.5, -0.25, -0.25, 0.5},
{0, 0.5, -0.5, 0.5},
{0, 0, 0, 1}
};
}
}
public static class LLab extends LinearTransform {
double [][] getTransform() {
return new double[][]{
{Wr, Wg, Wb, 0},
{0.5, -0.5, 0, 0.5},
{0, 0.5, -0.5, 0.5},
{0, 0, 0, 1}
};
}
}
public static class YCC extends LinearTransform {
double [][] getTransform() {
return new double[][]{
{Wr, Wg, Wb, 0},
{-Wr / (2 - 2 * Wb) + .5, -Wg / (2 - 2 * Wb) + .5, (1 - Wb) / (2 - 2 * Wb) + .5, 0},
{(1 - Wr) / (2 - 2 * Wr) + .5, -Wg / (2 - 2 * Wr) + .5, -Wb / (2 - 2 * Wr) + .5, 0},
{0, 0, 0, 1}
};
}
}
static Matrix vec(float[] A, float[] B, float[] C) {
Matrix ABC = new Matrix(3, 3);
for (int i = 0; i < 3; i++)
ABC.set(i, 0, A[i]);
for (int i = 0; i < 3; i++)
ABC.set(i, 1, B[i]);
for (int i = 0; i < 3; i++)
ABC.set(i, 2, C[i]);
return ABC;
}
public static double[][] strip(double[][] x) {
double r[][] = new double[3][];
r[0] = x[0];
r[1] = x[1];
r[2] = x[2];
return r;
}
// Bradford cone response matrix, seems to deliver more consistent results
static float[][] Bradford = {
{ 0.8951f, -0.7502f, 0.0389f},
{ 0.2664f, 1.7135f, -0.0685f},
{-0.1614f, 0.0367f, 1.0296f}
};
static float[][] VonKries = {
{0.40024f, -0.22630f, 0.00000f},
{0.70760f, 1.16532f, 0.00000f},
{-0.08081f, 0.04570f, 0.91822f},
};
// CAT02 matrix, sometimes gives more "neutral" results (no cyan cast)
static float[][] CAT02 = {
{ 0.7328f, -0.7036f, 0.0030f},
{ 0.4296f, 1.6975f, 0.0136f},
{-0.1624f, 0.0061f, 0.9834f}
};
static float[][] Sharp = {
{ 1.2694f, -0.8364f, 0.0297f},
{-0.0988f, 1.8006f, -0.0315f},
{-0.1706f, 0.0357f, 1.0018f}
};
static float[][] CMCCAT2000 = {
{ 0.7982f, -0.5918f, 0.0008f},
{ 0.3389f, 1.5512f, 0.0239f},
{-0.1371f, 0.0406f, 0.9753f}
};
static float[][] XYZScaling = {
{1, 0, 0},
{0, 1, 0},
{0, 0, 1}
};
public enum CAMethod {
Bradford, VonKries, CAT02, Sharp, CMCCAT2000, XYZScaling, Mixed
}
static float mixer(float t) {
return (float) (Math.atan((t - 5000)/100) / Math.PI + 0.5);
}
static float[][] matrix(float t) {
float matrix[][] = new float[3][3];
float m = mixer(t);
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
matrix[i][j] = XYZScaling[i][j] * (1-m) + Bradford[i][j] * m;
return matrix;
}
public static float[][] chromaticAdaptation(float source, float target) {
return chromaticAdaptation(source, target, CAMethod.Bradford);
}
public static float[][] chromaticAdaptation(float source, float target, CAMethod caMethod) {
final float[][] method;
switch (caMethod) {
case Bradford:
method = Bradford;
break;
case VonKries:
method = VonKries;
break;
case CAT02:
method = CAT02;
break;
case Sharp:
method = Sharp;
break;
case CMCCAT2000:
method = CMCCAT2000;
break;
case Mixed:
method = matrix(target);
break;
case XYZScaling:
default:
method = XYZScaling;
break;
}
Matrix B = new Matrix(method);
// source illuminant tristimulus in cone response domain
float[] sXYZ = xy2XYZ(D(source));
sXYZ = new Matrix(new float[][] {{sXYZ[0], sXYZ[1], sXYZ[2]}}).times(B).getArrayFloat()[0];
// target illuminant tristimulus in cone response domain
float[] tXYZ = xy2XYZ(D(target));
tXYZ = new Matrix(new float[][] {{tXYZ[0], tXYZ[1], tXYZ[2]}}).times(B).getArrayFloat()[0];
// scaling matrix for the colors
float[][] diag = new float[][]{
{sXYZ[0] / tXYZ[0], 0, 0},
{0, sXYZ[1] / tXYZ[1], 0},
{0, 0, sXYZ[2] / tXYZ[2]}
};
// total tansform
return B.times(new Matrix(diag)).times(B.inverse()).getArrayFloat();
}
public static double saturation(double r, double g, double b) {
double min = Math.min(r, Math.min(g, b));
double max = Math.max(r, Math.max(g, b));
return max != 0 ? 1 - min / max : 0;
}
private static Matrix RGBtoZYX = new Matrix(RGBtoZYX()).transpose();
private static Matrix XYZtoRGB = RGBtoZYX.inverse();
public static float[] neutralTemperature(float rgb[], float refT, CAMethod caMethod) {
float sat = Float.MAX_VALUE;
float minT = 0;
double wbr = 0, wbg = 0, wbb = 0;
float tint = 0;
Matrix color = new Matrix(new double[][]{{rgb[0]}, {rgb[1]}, {rgb[2]}});
for (int t = 1000; t < 40000; t+= (0.001 * t)) {
Matrix B = new Matrix(chromaticAdaptation(t, refT, caMethod));
Matrix combo = XYZtoRGB.times(B.times(RGBtoZYX));
Matrix adapdedColor = combo.times(color);
double r = clip(adapdedColor.get(0, 0));
double g = clip(adapdedColor.get(1, 0));
double b = clip(adapdedColor.get(2, 0));
float tSat = (float) saturation(r, g, b);
if (tSat < sat) {
sat = tSat;
minT = t;
wbr = r;
wbg = g;
wbb = b;
}
}
if (wbr != 0 || wbg != 0 || wbb != 0) {
tint = (float) (- (wbg - (wbr + wbb) / 2));
}
return new float[] {minT, tint};
}
static double clip(double x) {
return Math.min(Math.max(0, x), 1);
}
public static float findTemperature(float rgb[], float refT, CAMethod caMethod) {
float minDiff = Float.MAX_VALUE;
float minT = 0;
float[] xyzRef = JAIContext.linearColorSpace.toCIEXYZ(rgb);
float[] labRef = JAIContext.labColorSpace.fromCIEXYZ(xyzRef);
Matrix gray = new Matrix(new double[][]{{0.18}, {0.18}, {0.18}});
for (int t = 1000; t < 40000; t+= (0.001 * t)) {
Matrix B = new Matrix(chromaticAdaptation(t, refT, caMethod));
Matrix combo = XYZtoRGB.times(B.times(RGBtoZYX));
gray = combo.times(gray);
double r = clip(gray.get(0, 0));
double g = clip(gray.get(1, 0));
double b = clip(gray.get(2, 0));
float[] xyzGray = JAIContext.linearColorSpace.toCIEXYZ(new float[] {(float) r, (float) g, (float) b});
float[] labGray = JAIContext.labColorSpace.fromCIEXYZ(xyzGray);
float diff = 0;
for (int i = 1; i < 3; i++) {
float di = labGray[i] - labRef[i];
diff += di * di;
}
diff = (float) Math.sqrt(diff);
if (diff < minDiff) {
minDiff = diff;
minT = t;
/* wbr = r / 256;
wbg = g / 256;
wbb = b / 256; */
}
}
/* if (wbr != 0 || wbg != 0 || wbb != 0) {
tint = (float) (- (wbg - (wbr + wbb) / 2));
} */
return minT;
}
public static void main( final String[] args ) {
for (int i = 2000; i < 25000; i += 500)
System.out.println("m(" + i + ") : " + mixer(i));
float D65[] = D(whitePointTemperature);
System.out.println("D65: " + D65[0] + ", " + D65[1] + ", " + (1 - D65[0] - D65[1]));
System.out.println("xr: " + Cxy[0][0] + ", yr: " + Cxy[0][1]);
System.out.println("xg: " + Cxy[1][0] + ", yg: " + Cxy[1][1]);
System.out.println("xb: " + Cxy[2][0] + ", yb: " + Cxy[2][1]);
System.out.println("W: " + Wr + ", " + Wg + ", " + Wb);
/* Matrix ww = new Matrix(RGBtoZYX()).transpose();
ww.print(8, 6);
Matrix ca = new Matrix(chromaticAdaptation(7500, 5000));
float[] result = new Matrix(new float[][] {{0.2, 0.2, 0.2}}).times(ca).getArray()[0];
for (int j = 0; j < 3; j++)
System.out.print(" " + result[j]);
System.out.println(); */
System.out.println("RGBToXYZ");
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++)
System.out.print(" " + RGBToXYZMat[i][j]);
System.out.println();
}
float[][] rgb2xyz = RGBtoZYX();
System.out.println("rgb2xyz");
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++)
System.out.print(" " + rgb2xyz[i][j]);
System.out.println();
}
System.out.println("XYZToRGB");
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++)
System.out.print(" " + XYZToRGBMat[i][j]);
System.out.println();
}
}
}