MakeInterpolationTable:
PROC [colorOperator: ColorOperator, transform: ColorTransform, maxIn: PixelProc]
RETURNS [InterpolationTable] ~ {
colorSpace: ColorSpace ~ transform.domain;
dim: NAT ~ ColorSpaceDimension[colorSpace];
samplesPerPixelOut: NAT ~ transform.rangeMax.dim;
pixelEncoding: ImagerColor.PixelEncoding ~ GetPixelEncoding[colorOperator];
result: InterpolationTable ¬ NEW[InterpolationTableRep];
vScratch: ColorPoint ¬ MakeColorPoint[dim];
CornerIndex: TYPE ~ [0..8);
d: ARRAY CornerIndex OF ColorPoint;
BuildMatrix:
PROC [ix1, ix2: CornerIndex]
RETURNS [
REF MatrixRep] ~ {
ERROR;
Body
};
FOR i:
NAT
IN [0..8)
DO
d[i] ¬ MakeColorPoint[samplesPerPixelOut];
ENDLOOP;
FOR i:
NAT
IN [0..colorOperator.samplesPerPixelIn)
DO
tableMin: REAL ¬ Real.LargestNumber;
tableMax: REAL ¬ -Real.LargestNumber;
sMax: NAT ~ maxIn[i];
encode: REF EncoderArray ¬ NEW[EncoderArray[sMax+1]];
IF pixelEncoding =
NIL
THEN {
tableMin ¬ 0;
tableMax ¬ sMax;
}
ELSE {
FOR j:
NAT
IN [0..sMax]
DO
s: REAL ~ ApplyPixelEncoding[pixelEncoding, i, j];
IF s < tableMin THEN tableMin ¬ s;
IF s > tableMax THEN tableMax ¬ s;
ENDLOOP;
};
IF tableMin < tableMax
THEN {
FOR j:
NAT
IN [0..sMax]
DO
s: REAL ~ ApplyPixelEncoding[pixelEncoding, i, j];
encode[j] ¬ RealInline.MCRound[((s-tableMin)/(tableMax-tableMin)) * encodedMax];
ENDLOOP;
};
encode.min ¬ tableMin;
encode.max ¬ tableMax;
result.encode[i] ¬ encode;
ENDLOOP;
FOR xu: [0..16)
IN [0..16)
DO
FOR yu: [0..16)
IN [0..16)
DO
FOR zu: [0..16)
IN [0..16)
DO
cc: ColorCube ~ NEW[ColorCubeRep];
FOR xf: [0..1]
IN [0..1]
DO
FOR yf: [0..1]
IN [0..1]
DO
FOR zf: [0..1]
IN [0..1]
DO
pt: [0..8) ~ xf*4 + yf*2 + zf;
v: ColorPoint ~ vScratch;
pixelIn: TupleProc ~ {
s: [0..encodedMax] ~
SELECT i
FROM
0 => xu*encodedUnit+xf*(encodedUnit-1),
1 => yu*encodedUnit+yf*(encodedUnit-1),
ENDCASE => zu*encodedUnit+zf*(encodedUnit-1);
e: REF EncoderArray ~ result.encode[i];
RETURN [(REAL[s]/REAL[encodedMax])*(e.max-e.min) + e.min]
};
class: ColorOperatorClass ~ colorOperator.class;
class.apply[colorOperator, pixelIn, colorSpace, v];
transform.proc[transform, v, d[pt]];
ENDLOOP;
ENDLOOP;
ENDLOOP;
cc.matrix[0] ¬ BuildMatrix[4,6];
cc.matrix[1] ¬ BuildMatrix[4,5];
cc.matrix[2] ¬ BuildMatrix[1,5];
cc.matrix[3] ¬ BuildMatrix[2,6];
cc.matrix[4] ¬ BuildMatrix[2,3];
cc.matrix[5] ¬ BuildMatrix[1,3];
cc.base ¬ d[0];
result.gamutSamples[xu][yu][zu] ¬ cc;
d[0] ¬ MakeColorPoint[samplesPerPixelOut];
ENDLOOP;
ENDLOOP;
ENDLOOP;
RETURN [result]
};
ApplyInterpolation:
PROC [interpolationTable: InterpolationTable, colorPoint: ColorPoint]
RETURNS [result: ColorPoint] ~ {
x: [0..encodedMax] ~ interpolationTable.encode[0][RealInline.MCRound[colorPoint[0]]];
y: [0..encodedMax] ~ interpolationTable.encode[1][RealInline.MCRound[colorPoint[1]]];
z: [0..encodedMax] ~ interpolationTable.encode[2][RealInline.MCRound[colorPoint[2]]];
colorCube: ColorCube ~ interpolationTable.gamutSamples[x/encodedUnit][y/encodedUnit][z/encodedUnit];
xr: [0..encodedUnit) ~ x MOD encodedUnit;
yr: [0..encodedUnit) ~ y MOD encodedUnit;
zr: [0..encodedUnit) ~ z MOD encodedUnit;
The parallipiped is split into 6 tets by subdividing it along the planes x=y, x=z, and y=z.
Tetrahedron 0: x > y, x > z, y > z. Vertices are points [0,4,6,7].
Tetrahedron 1: x > y, x > z, y <= z. Vertices are points [0,4,5,7].
Tetrahedron 2: x > y, x <= z, y <= z. Vertices are points [0,1,5,7].
Tetrahedron 3: x <= y, x > z, y > z. Vertices are points [0,2,6,7].
Tetrahedron 4: x <= y, x <= z, y > z. Vertices are points [0,2,3,7].
Tetrahedron 5: x <= y, x <= z, y <= z. Vertices are points [0,1,3,7].
A cell can be subdivided into 5 tets, not 6: the beauty of the 6 tet subdivision is that the cutting planes are so simple that determining which tet a point is inside is very quick.
index: [0..6) ~
IF xr > yr
THEN (
IF xr > zr
THEN (IF yr > zr THEN 0 ELSE 1)
ELSE 2
)
ELSE (
IF xr > zr
THEN 3
ELSE (IF yr > zr THEN 4 ELSE 5)
);
matrix: REF MatrixRep ~ colorCube.matrix[index];
result ¬ MakeColorPoint[colorCube.base.dim];
FOR j:
NAT
IN [0..result.dim)
DO
v: REF ARRAY [0..3) OF REAL ~ matrix[j];
result[j] ¬ colorCube.base[j] + xr*v[0] + yr*v[1] + zr*v[2];
ENDLOOP;
};