/* */
/* Little cms - profiler construction set */
/* Copyright (C) 1998-2001 Marti Maria <marti@littlecms.com> */
/* */
/* THIS SOFTWARE IS PROVIDED "AS-IS" AND WITHOUT WARRANTY OF ANY KIND, */
/* EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY */
/* WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. */
/* */
/* IN NO EVENT SHALL MARTI MARIA BE LIABLE FOR ANY SPECIAL, INCIDENTAL, */
/* INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
/* OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, */
/* WHETHER OR NOT ADVISED OF THE POSSIBILITY OF DAMAGE, AND ON ANY THEORY OF */
/* LIABILITY, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE */
/* OF THIS SOFTWARE. */
/* */
/* This file is free software; you can redistribute it and/or modify it */
/* under the terms of the GNU General Public License as published by */
/* the Free Software Foundation; either version 2 of the License, or */
/* (at your option) any later version. */
/* */
/* This program is distributed in the hope that it will be useful, but */
/* WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU */
/* General Public License for more details. */
/* */
/* You should have received a copy of the GNU General Public License */
/* along with this program; if not, write to the Free Software */
/* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */
/* */
/* As a special exception to the GNU General Public License, if you */
/* distribute this file as part of a program that contains a */
/* configuration script generated by Autoconf, you may include it under */
/* the same distribution terms that you use for the rest of that program. */
/* */
/* Version 1.09a */
#include "lcmsprf.h"
LPGAMMATABLE cdecl cmsxEstimateGamma(LPSAMPLEDCURVE X, LPSAMPLEDCURVE Y, int nResultingPoints);
void cdecl cmsxCompleteLabOfPatches(LPMEASUREMENT m, SETOFPATCHES Valids, int Medium);
void cdecl cmsxComputeLinearizationTables(LPMEASUREMENT m,
int ColorSpace,
LPGAMMATABLE Lin[3],
int nResultingPoints,
int Medium);
void cdecl cmsxApplyLinearizationTable(double In[3],
LPGAMMATABLE Gamma[3],
double Out[3]);
void cdecl cmsxApplyLinearizationGamma(WORD In[3], LPGAMMATABLE Gamma[3], WORD Out[3]);
/* ------------------------------------------------------------- Implementation */
#define EPSILON 0.00005
#define LEVENBERG_MARQUARDT_ITERATE_MAX 150
/* In order to track linearization tables, we use following procedure */
/* */
/* We first assume R', G' and B' does exhibit a non-linear behaviour */
/* that can be separated for each channel as Yr(R'), Yg(G'), Yb(B') */
/* This is the shaper step */
/* */
/* R = Lr(R') */
/* G = Lg(G') */
/* B = Lb(B') (0.0) */
/* */
/* After this step, RGB is converted to XYZ by a matrix multiplication */
/* */
/* |X| |R| */
/* |Y| = [M]�|G| */
/* |Z| |B| (1.0) */
/* */
/* In order to extract Lr,Lg,Lb tables, we are interested only on Y part */
/* */
/* Y = (m1 * R + m2 * G + m3 * B) (1.1) */
/* */
/* The total intensity for maximum RGB = (1, 1, 1) should be 1, */
/* */
/* 1 = m1 * 1 + m2 * 1 + m3 * 1, so */
/* */
/* m1 + m2 + m3 = 1.0 (1.2) */
/* */
/* We now impose that for neutral (gray) patches, RGB components must be equal */
/* */
/* R = G = B = Gray */
/* */
/* So, substituting in (1.1): */
/* */
/* Y = (m1 + m2 + m3) Gray */
/* */
/* and for (1.2), (m1+m2+m3) = 1, so */
/* */
/* Y = Gray = Lr(R') = Lg(G') = Lb(B') */
/* */
/* That is, after prelinearization, RGB of gray patches should give */
/* same values for R, G and B. And this value is Y. */
/* */
/* */
static
LPSAMPLEDCURVE NormalizeTo(LPSAMPLEDCURVE X, double N, BOOL lAddEndPoint)
{
int i, nItems;
LPSAMPLEDCURVE XNorm;
nItems = X ->nItems;
if (lAddEndPoint) nItems++;
XNorm = cmsAllocSampledCurve(nItems);
for (i=0; i < X ->nItems; i++) {
XNorm ->Values[i] = X ->Values[i] / N;
}
if (lAddEndPoint)
XNorm -> Values[X ->nItems] = 1.0;
return XNorm;
}
/* */
/* ------------------------------------------------------------------------------ */
/* */
/* Our Monitor model. We assume gamma has a general expression of */
/* */
/* Fn(x) = (Gain * x + offset) ^ gamma | for x >= 0 */
/* Fn(x) = 0 | for x < 0 */
/* */
/* First partial derivatives are */
/* */
/* dFn/dGamma = Fn * ln(Base) */
/* dFn/dGain = gamma * x * ((Gain * x + Offset) ^ (gamma -1)) */
/* dFn/dOffset = gamma * ((Gain * x + Offset) ^ (gamma -1)) */
/* */
static
void GammaGainOffsetFn(double x, double *a, double *y, double *dyda, int na)
{
double Gamma,Gain,Offset;
double Base;
Gamma = a[0];
Gain = a[1];
Offset = a[2];
Base = Gain * x + Offset;
if (Base < 0) {
Base = 0.0;
*y = 0.0;
dyda[0] = 0.0;
dyda[1] = 0.0;
dyda[2] = 0.0;
} else {
/* The function itself */
*y = pow(Base, Gamma);
/* dyda[0] is partial derivative across Gamma */
dyda[0] = *y * log(Base);
/* dyda[1] is partial derivative across gain */
dyda[1] = (x * Gamma) * pow(Base, Gamma-1.0);
/* dyda[2] is partial derivative across offset */
dyda[2] = Gamma * pow(Base, Gamma-1.0);
}
}
/* Fit curve to our gamma-gain-offset model. */
static
BOOL OneTry(LPSAMPLEDCURVE XNorm, LPSAMPLEDCURVE YNorm, double a[])
{
LCMSHANDLE h;
double ChiSq, OldChiSq;
int i;
BOOL Status = true;
/* initial guesses */
a[0] = 3.0; /* gamma */
a[1] = 4.0; /* gain */
a[2] = 6.0; /* offset */
a[3] = 0.0; /* Thereshold */
a[4] = 0.0; /* Black */
/* Significance = 0.02 gives good results */
h = cmsxLevenbergMarquardtInit(XNorm, YNorm, 0.02, a, 3, GammaGainOffsetFn);
if (h == NULL) return false;
OldChiSq = cmsxLevenbergMarquardtChiSq(h);
for(i = 0; i < LEVENBERG_MARQUARDT_ITERATE_MAX; i++) {
if (!cmsxLevenbergMarquardtIterate(h)) {
Status = false;
break;
}
ChiSq = cmsxLevenbergMarquardtChiSq(h);
if(OldChiSq != ChiSq && (OldChiSq - ChiSq) < EPSILON)
break;
OldChiSq = ChiSq;
}
cmsxLevenbergMarquardtFree(h);
return Status;
}
/* Tries to fit gamma as per IEC 61966-2.1 using Levenberg-Marquardt method */
/* */
/* Y = (aX + b)^Gamma | X >= d */
/* Y = cX | X < d */
LPGAMMATABLE cmsxEstimateGamma(LPSAMPLEDCURVE X, LPSAMPLEDCURVE Y, int nResultingPoints)
{
double a[5];
LPSAMPLEDCURVE XNorm, YNorm;
double e, Max;
/* Coarse approximation, to find maximum. */
/* We have only a portion of curve. It is likely */
/* maximum will not fall on exactly 100. */
if (!OneTry(X, Y, a))
return 0;
/* Got parameters. Compute maximum. */
e = a[1]* 255.0 + a[2];
if (e < 0) return 0;
Max = pow(e, a[0]);
/* Normalize values to maximum */
XNorm = NormalizeTo(X, 255.0, false);
YNorm = NormalizeTo(Y, Max, false);
/* Do the final fitting */
if (!OneTry(XNorm, YNorm, a))
return 0;
/* Type 3 = IEC 61966-2.1 (sRGB) */
/* Y = (aX + b)^Gamma | X >= d */
/* Y = cX | X < d */
return cmsBuildParametricGamma(nResultingPoints, 3, a);
}
/* A dumb bubble sort */
static
void Bubble(LPSAMPLEDCURVE C, LPSAMPLEDCURVE L)
{
#define SWAP(a, b) { tmp = (a); (a) = (b); (b) = tmp; }
BOOL lSwapped;
int i, nItems;
double tmp;
nItems = C -> nItems;
do {
lSwapped = false;
for (i= 0; i < nItems - 1; i++) {
if (C->Values[i] > C->Values[i+1]) {
SWAP(C->Values[i], C->Values[i+1]);
SWAP(L->Values[i], L->Values[i+1]);
lSwapped = true;
}
}
} while (lSwapped);
#undef SWAP
}
/* Check for monotonicity. Force it if is not the case. */
static
void CheckForMonotonicSampledCurve(LPSAMPLEDCURVE t)
{
int n = t ->nItems;
int i;
double last;
last = t ->Values[n-1];
for (i = n-2; i >= 0; --i) {
if (t ->Values[i] > last)
t ->Values[i] = last;
else
last = t ->Values[i];
}
}
/* The main gamma inferer. Tries first by gamma-gain-offset, */
/* if not proper reverts to curve guessing. */
static
LPGAMMATABLE BuildGammaTable(LPSAMPLEDCURVE C, LPSAMPLEDCURVE L, int nResultingPoints)
{
LPSAMPLEDCURVE Cw, Lw, Cn, Ln;
LPSAMPLEDCURVE out;
LPGAMMATABLE Result;
double Lmax, Lend, Cmax;
/* Try to see if it can be fitted */
Result = cmsxEstimateGamma(C, L, nResultingPoints);
if (Result)
return Result;
/* No... build curve from scratch. Since we have not */
/* endpoints, a coarse linear extrapolation should be */
/* applied in order to get the expected maximum. */
Cw = cmsDupSampledCurve(C);
Lw = cmsDupSampledCurve(L);
Bubble(Cw, Lw);
/* Get endpoint */
Lmax = Lw->Values[Lw ->nItems - 1];
Cmax = Cw->Values[Cw ->nItems - 1];
/* Linearly extrapolate */
Lend = (255 * Lmax) / Cmax;
Ln = NormalizeTo(Lw, Lend, true);
Cn = NormalizeTo(Cw, 255.0, true);
cmsFreeSampledCurve(Cw);
cmsFreeSampledCurve(Lw);
/* Add endpoint */
out = cmsJoinSampledCurves(Cn, Ln, nResultingPoints);
cmsFreeSampledCurve(Cn);
cmsFreeSampledCurve(Ln);
CheckForMonotonicSampledCurve(out);
cmsSmoothSampledCurve(out, nResultingPoints*4.);
cmsClampSampledCurve(out, 0, 1.0);
Result = cmsConvertSampledCurveToGamma(out, 1.0);
cmsFreeSampledCurve(out);
return Result;
}
void cmsxCompleteLabOfPatches(LPMEASUREMENT m, SETOFPATCHES Valids, int Medium)
{
LPPATCH White;
cmsCIEXYZ WhiteXYZ;
int i;
if (Medium == MEDIUM_REFLECTIVE_D50)
{
WhiteXYZ.X = D50X * 100.;
WhiteXYZ.Y = D50Y * 100.;
WhiteXYZ.Z = D50Z * 100.;
}
else {
White = cmsxPCollFindWhite(m, Valids, NULL);
if (!White) return;
WhiteXYZ = White ->XYZ;
}
/* For all patches with XYZ and without Lab, add Lab values. */
/* Transmissive profiles does need to locate its own white */
/* point for device gray. Reflective does use D50 */
for (i=0; i < m -> nPatches; i++) {
if (Valids[i]) {
LPPATCH p = m -> Patches + i;
if ((p ->dwFlags & PATCH_HAS_XYZ) &&
(!(p ->dwFlags & PATCH_HAS_Lab) || (Medium == MEDIUM_TRANSMISSIVE))) {
cmsXYZ2Lab(&WhiteXYZ, &p->Lab, &p->XYZ);
p -> dwFlags |= PATCH_HAS_Lab;
}
}
}
}
/* Compute linearization tables, trying to fit in a pure */
/* exponential gamma. If gamma cannot be accurately infered, */
/* then does build a smooth, monotonic curve that does the job. */
void cmsxComputeLinearizationTables(LPMEASUREMENT m,
int ColorSpace,
LPGAMMATABLE Lin[3],
int nResultingPoints,
int Medium)
{
LPSAMPLEDCURVE R, G, B, L;
LPGAMMATABLE gr, gg, gb;
SETOFPATCHES Neutrals;
int nGrays;
int i;
/* We need Lab for grays. */
cmsxCompleteLabOfPatches(m, m->Allowed, Medium);
/* Add neutrals, normalize to max */
Neutrals = cmsxPCollBuildSet(m, false);
cmsxPCollPatchesNearNeutral(m, m ->Allowed, 15, Neutrals);
nGrays = cmsxPCollCountSet(m, Neutrals);
R = cmsAllocSampledCurve(nGrays);
G = cmsAllocSampledCurve(nGrays);
B = cmsAllocSampledCurve(nGrays);
L = cmsAllocSampledCurve(nGrays);
nGrays = 0;
/* Collect patches */
for (i=0; i < m -> nPatches; i++) {
if (Neutrals[i]) {
LPPATCH gr = m -> Patches + i;
R -> Values[nGrays] = gr -> Colorant.RGB[0];
G -> Values[nGrays] = gr -> Colorant.RGB[1];
B -> Values[nGrays] = gr -> Colorant.RGB[2];
L -> Values[nGrays] = gr -> XYZ.Y;
nGrays++;
}
}
gr = BuildGammaTable(R, L, nResultingPoints);
gg = BuildGammaTable(G, L, nResultingPoints);
gb = BuildGammaTable(B, L, nResultingPoints);
cmsFreeSampledCurve(R);
cmsFreeSampledCurve(G);
cmsFreeSampledCurve(B);
cmsFreeSampledCurve(L);
if (ColorSpace == PT_Lab) {
LPGAMMATABLE Gamma3 = cmsBuildGamma(nResultingPoints, 3.0);
Lin[0] = cmsJoinGammaEx(gr, Gamma3, nResultingPoints);
Lin[1] = cmsJoinGammaEx(gg, Gamma3, nResultingPoints);
Lin[2] = cmsJoinGammaEx(gb, Gamma3, nResultingPoints);
cmsFreeGamma(gr); cmsFreeGamma(gg); cmsFreeGamma(gb);
cmsFreeGamma(Gamma3);
}
else {
LPGAMMATABLE Gamma1 = cmsBuildGamma(nResultingPoints, 1.0);
Lin[0] = cmsJoinGammaEx(gr, Gamma1, nResultingPoints);
Lin[1] = cmsJoinGammaEx(gg, Gamma1, nResultingPoints);
Lin[2] = cmsJoinGammaEx(gb, Gamma1, nResultingPoints);
cmsFreeGamma(gr); cmsFreeGamma(gg); cmsFreeGamma(gb);
cmsFreeGamma(Gamma1);
}
}
/* Apply linearization. WORD encoded version */
void cmsxApplyLinearizationGamma(WORD In[3], LPGAMMATABLE Gamma[3], WORD Out[3])
{
L16PARAMS Lut16;
cmsCalcL16Params(Gamma[0] -> nEntries, &Lut16);
Out[0] = cmsLinearInterpLUT16(In[0], Gamma[0] -> GammaTable, &Lut16);
Out[1] = cmsLinearInterpLUT16(In[1], Gamma[1] -> GammaTable, &Lut16);
Out[2] = cmsLinearInterpLUT16(In[2], Gamma[2] -> GammaTable, &Lut16);
}
/* Apply linearization. double version */
void cmsxApplyLinearizationTable(double In[3], LPGAMMATABLE Gamma[3], double Out[3])
{
WORD rw, gw, bw;
double rd, gd, bd;
L16PARAMS Lut16;
cmsCalcL16Params(Gamma[0] -> nEntries, &Lut16);
rw = (WORD) floor(_cmsxSaturate255To65535(In[0]) + .5);
gw = (WORD) floor(_cmsxSaturate255To65535(In[1]) + .5);
bw = (WORD) floor(_cmsxSaturate255To65535(In[2]) + .5);
rd = cmsLinearInterpLUT16(rw , Gamma[0] -> GammaTable, &Lut16);
gd = cmsLinearInterpLUT16(gw, Gamma[1] -> GammaTable, &Lut16);
bd = cmsLinearInterpLUT16(bw, Gamma[2] -> GammaTable, &Lut16);
Out[0] = _cmsxSaturate65535To255(rd); /* back to 0..255 */
Out[1] = _cmsxSaturate65535To255(gd);
Out[2] = _cmsxSaturate65535To255(bd);
}