835
836
wr_xicc, icmFwd, icmDefaultIntent,
837
838
flags, /* Flags */
838
npat, tpat, 0.0, wpscale, smooth, avgdev,
839
npat, npat, tpat, NULL, 0.0, wpscale, smooth, avgdev,
839
840
NULL, NULL, NULL, iquality)) == NULL)
840
841
error("%d, %s",wr_xicc->errc, wr_xicc->err);
846
847
} else { /* cLUT based profile */
849
icxMatrixModel *mm = NULL;
849
851
xicc *wr_xicc; /* extention object */
850
852
icxLuBase *AtoB; /* AtoB ixcLu */
854
double avgdist; /* Average distance between points */
860
859
if (verb) printf("Creating extrapolation black and white points:\n");
862
avgdist = pow(1.0/(double)npat, 1.0/3.0);
865
else if (avgdist > 0.3)
867
//printf("~1 avgdist = %f\n",avgdist);
869
861
if ((mpat = (cow *)malloc(sizeof(cow) * npat)) == NULL)
870
862
error("Malloc failed - mpat[]");
872
/* Select points from full set */
873
for (range = 0.05;; range *= 1.5) {
877
for (nmpat = j = 0; j < npat; j++) {
881
icmCpy3(mpat[nmpat].p, tpat[j].p);
882
icmCpy3(mpat[nmpat].v, tpat[j].v);
884
/* Locate largest/smallest RGB value */
885
mxp = -1e6, mnp = 1e6;
886
for (k = 0; k < 3; k++) {
887
if (tpat[j].p[k] > mxp)
889
if (tpat[j].p[k] < mnp)
892
mxp -= mnp; /* Spread; 0 for R=G=B */
895
mpat[nmpat].w = 1.0 - mxp/range;
864
/* Weight points from full set to build matrix model */
865
/* to extrapolate the neutral axis */
866
for (nmpat = j = 0; j < npat; j++) {
870
icmCpy3(mpat[nmpat].p, tpat[j].p);
871
icmCpy3(mpat[nmpat].v, tpat[j].v);
873
/* Locate largest/smallest RGB value */
874
mxp = -1e6, mnp = 1e6;
875
for (k = 0; k < 3; k++) {
876
if (tpat[j].p[k] > mxp)
878
if (tpat[j].p[k] < mnp)
881
mxp -= mnp; /* Spread; 0 for R=G=B */
883
mpat[nmpat].w = pow(1.1 - mxp, 2.0);
896
884
//printf("~1 added value %d: %f %f %f -> %f %f %f wt %f\n",j, mpat[nmpat].p[0], mpat[nmpat].p[1], mpat[nmpat].p[2], mpat[nmpat].v[0], mpat[nmpat].v[1], mpat[nmpat].v[2],mpat[nmpat].w);
901
if (nmpat >= 16 || range >= 0.99) {
902
//printf("~1 stopping with %d points at range %f\n",nmpat,range);
906
/* Hmm. Not enough points with that range */
909
if (verb) printf("%d/%d patches for extrapolation gamma/matrix model\n",nmpat,npat);
911
888
/* Create gamma/matrix model to extrapolate with. */
912
889
/* (Use ofset & gain, gamma curve as 0th order with 1 harmonic, */
913
/* and heavily smooth it.) */
890
/* and smooth it.) */
914
891
if ((mm = new_MatrixModel(verb, nmpat, mpat, wantLab,
915
892
/* quality */ -1, /* isLinear */ ptype == prof_matonly,
916
893
/* isGamma */ 0, /* isShTRC */ 0,
917
894
/* shape0gam */ 1, /* clipbw */ 0, /* clipprims */ 0,
918
/* smooth */ 1.0, /* scale */ 0.7)) == NULL) {
895
// /* smooth */ 1.0, /* scale */ 0.7)) == NULL) {
896
/* smooth */ 1.0, /* scale */ 1.0)) == NULL) {
919
897
error("Creating extrapolation matrix model failed - memory ?");
900
#ifdef NEVER /* Plot Lab of model */
940
919
do_plot(xx,y0,y1,y2,XRES);
942
921
#endif /* DEBUG_PLOT */
945
/* Create a black and white patch */
946
for (i = 0; i < 2; i++) {
947
int cix; /* Closest point index */
948
int eix; /* End point index */
949
double cde = 1e60; /* Closest point distance */
951
double corr[3], cwt; /* Correction */
953
tpat[npat + nxpat].p[0] =
954
tpat[npat + nxpat].p[1] =
955
tpat[npat + nxpat].p[2] = (double)i;
957
/* Locate closest point */
958
for (nmpat = j = 0; j < npat; j++) {
962
/* Locate largest/smallest RGB value */
963
mxp = -1e6, mnp = 1e6;
964
for (k = 0; k < 3; k++) {
965
if (tpat[j].p[k] > mxp)
967
if (tpat[j].p[k] < mnp)
970
mxp -= mnp; /* Spread; 0 for R=G=B */
972
tt = icmNorm33(tpat[npat + nxpat].p, tpat[j].p);
926
int pcsy; /* Effective PCS L or Y chanel index */
928
double dwhite[MXDI]; /* Device white */
930
double avgdist; /* Average distance between points */
932
/* Figure out the device white point. */
933
/* Note that this is duplicating code in xicc/xmatrix.c */
937
pcsy = 0; /* L or Lab */
939
pcsy = 1; /* Y of XYZ */
941
for (i = 0; i < npat; i++) {
944
/* Create D50 Lab to allow some chromatic sensitivity */
945
/* in picking the white point */
947
icmCpy3(labv, tpat[i].v);
949
icmXYZ2Lab(&icmD50, labv, tpat[i].v);
951
/* Tilt things towards D50 neutral white patches */
952
yv = labv[0] - 0.3 * sqrt(labv[1] * labv[1] + labv[2] * labv[2]);
954
for (j = 0; j < 3; j++)
955
dwhite[j] = tpat[i].p[j];
961
/* Fix extrapolation matrix to be perfect at white point */
962
mm->force(mm, tpat[wix].v, tpat[wix].p);
964
/* Scale the white point to make one dev value 1.0 */
966
for (j = 0; j < 3; j++) {
967
if (dwhite[j] > mxdw)
970
for (j = 0; j < 3; j++) {
974
avgdist = pow(1.0/(double)npat, 1.0/3.0);
977
else if (avgdist > 0.3)
979
//printf("~1 avgdist = %f\n",avgdist);
981
/* For points with white point device ratio, */
982
/* and points with R=G=B ratio, create extrapolation points. */
983
for (ii = 0; ii < 2; ii++) {
986
dwhite[0] = dwhite[1] = dwhite[2] = 1.0;
988
/* Create a series of black and white patch */
989
for (i = 0; i < 2; i++) {
990
int cix; /* Closest point index */
991
int eix; /* End point index */
992
double cde = 1e60; /* Closest point distance */
994
double corr[3], cwt; /* Correction */
998
icmScale3(tpat[eix].p, dwhite, (double)i);
1000
/* Locate closest point */
1001
for (j = 0; j < npat; j++) {
1005
/* Locate largest/smallest RGB value */
1006
mxp = -1e6, mnp = 1e6;
1007
for (k = 0; k < 3; k++) {
1008
if (tpat[j].p[k] > mxp)
1010
if (tpat[j].p[k] < mnp)
1013
mxp -= mnp; /* Spread; 0 for R=G=B */
1015
tt = icmNorm33(tpat[eix].p, tpat[j].p);
981
1024
//printf("~1 closest %d: de %f, %f %f %f -> %f %f %f\n",cix, cde, tpat[cix].p[0], tpat[cix].p[1], tpat[cix].p[2], tpat[cix].v[0], tpat[cix].v[1], tpat[cix].v[2]);
986
1029
//printf("~1 closest gam/matrix -> %f %f %f\n",val[0],val[1],val[2]);
989
/* Lookup matrix value for our new point */
991
mm->lookup(mm, tpat[eix].v, tpat[eix].p);
1032
/* Lookup matrix value for our new point */
1033
mm->lookup(mm, tpat[eix].v, tpat[eix].p);
992
1034
//printf("~1 got value %d: %f %f %f -> %f %f %f\n",i, tpat[eix].p[0], tpat[eix].p[1], tpat[eix].p[2], tpat[eix].v[0], tpat[eix].v[1], tpat[eix].v[2]);
993
/* Weight the extra point so that it doesn't overpower the */
994
/* nearest real point to it too much. */
996
if (tt > avgdist) /* Distance at which sythetic point has 100% weight */
998
tpat[eix].w = EXTRAP_WEIGHT * tt/avgdist;
1035
/* Weight the extra point so that it doesn't overpower the */
1036
/* nearest real point to it too much. */
1038
if (tt > avgdist) /* Distance at which sythetic point has 100% weight */
1040
tpat[eix].w = 0.5 * EXTRAP_WEIGHT * tt/avgdist;
999
1041
//printf("~1 weight %f\n",tpat[eix].w);
1001
printf("Added synthetic point @ %f %f %f, val %f %f %f, weight %f\n",tpat[eix].p[0], tpat[eix].p[1], tpat[eix].p[2], tpat[eix].v[0], tpat[eix].v[1], tpat[eix].v[2],tpat[eix].w);
1004
/* If there is a lot of space, add a second intemediate point */
1043
printf("Added synthetic point @ %f %f %f, val %f %f %f, weight %f\n",tpat[eix].p[0], tpat[eix].p[1], tpat[eix].p[2], tpat[eix].v[0], tpat[eix].v[1], tpat[eix].v[2],tpat[eix].w);
1046
/* If there is a lot of space, add a second intemediate point */
1005
1047
//printf("~1 cde = %f, avgdist = %f\n",cde,avgdist);
1006
if (cde >= (0.5 * avgdist)) {
1007
int nxps; /* Number of extra points including end point */
1008
nxps = 1 + (int)(cde/(0.5 * avgdist));
1009
if (nxps > EXTRAP_MAXPNTS)
1010
nxps = EXTRAP_MAXPNTS;
1048
if (cde >= (0.5 * avgdist)) {
1049
int nxps; /* Number of extra points including end point */
1050
nxps = 1 + (int)(cde/(0.5 * avgdist));
1051
if (nxps > EXTRAP_MAXPNTS)
1052
nxps = EXTRAP_MAXPNTS;
1012
1054
//printf("~1 nxps = %d\n",nxps);
1013
for (j = 1; j < nxps; j++) {
1016
bl = j/(nxps + 1.0);
1055
for (j = 1; j < nxps; j++) {
1058
bl = j/(nxps + 1.0);
1018
ipos = (1.0 - bl) * tpat[eix].p[0]
1019
+ bl * (tpat[cix].p[0] + tpat[cix].p[1] + tpat[cix].p[1])/3.0;
1020
tpat[npat + nxpat].p[0] =
1021
tpat[npat + nxpat].p[1] =
1022
tpat[npat + nxpat].p[2] = ipos;
1024
/* Lookup matrix value for our new point */
1025
mm->lookup(mm, tpat[npat + nxpat].v, tpat[npat + nxpat].p);
1027
/* Weight the extra point so that it doesn't overpower the */
1028
/* nearest real point to it too much. */
1029
cde = icmNorm33(tpat[cix].p, tpat[npat + nxpat].p);
1031
if (cde > avgdist) /* Distance at which sythetic point has 100% weight */
1033
tpat[npat + nxpat].w = EXTRAP_WEIGHT * cde/avgdist;
1035
printf("Added synthetic point @ %f %f %f, val %f %f %f, weight %f\n",tpat[npat + nxpat].p[0], tpat[npat + nxpat].p[1], tpat[npat + nxpat].p[2], tpat[npat + nxpat].v[0], tpat[npat + nxpat].v[1], tpat[npat + nxpat].v[2],tpat[npat + nxpat].w);
1060
ipos = (1.0 - bl) * tpat[eix].p[0]
1061
+ bl * (tpat[cix].p[0] + tpat[cix].p[1] + tpat[cix].p[1])/3.0;
1062
icmScale3(tpat[npat + nxpat].p, dwhite, ipos);
1064
/* Lookup matrix value for our new point */
1065
mm->lookup(mm, tpat[npat + nxpat].v, tpat[npat + nxpat].p);
1067
/* Weight the extra point so that it doesn't overpower the */
1068
/* nearest real point to it too much. */
1069
cde = icmNorm33(tpat[cix].p, tpat[npat + nxpat].p);
1071
if (cde > avgdist) /* Distance at which sythetic point has 100% weight */
1073
tpat[npat + nxpat].w = 0.5 * EXTRAP_WEIGHT * cde/avgdist;
1075
printf("Added synthetic point @ %f %f %f, val %f %f %f, weight %f\n",tpat[npat + nxpat].p[0], tpat[npat + nxpat].p[1], tpat[npat + nxpat].p[2], tpat[npat + nxpat].v[0], tpat[npat + nxpat].v[1], tpat[npat + nxpat].v[2],tpat[npat + nxpat].w);
1043
1083
/* Wrap with an expanded icc */