1float	ycbcrCoeffs[3] = { .299, .587, .114 };
2/* default coding range is CCIR Rec 601-1 with no headroom/footroom */
3unsigned long refBlackWhite[6] = { 0, 255, 128, 255, 128, 255 };
4
5#define	LumaRed		ycbcrCoeffs[0]
6#define	LumaGreen	ycbcrCoeffs[1]
7#define	LumaBlue	ycbcrCoeffs[2]
8
9long	eRtotal = 0;
10long	eGtotal = 0;
11long	eBtotal = 0;
12long	preveRtotal = 0;
13long	preveGtotal = 0;
14long	preveBtotal = 0;
15unsigned long AbseRtotal = 0;
16unsigned long AbseGtotal = 0;
17unsigned long AbseBtotal = 0;
18unsigned long eCodes = 0;
19unsigned long preveCodes = 0;
20unsigned long eBits = 0;
21unsigned long preveBits = 0;
22
23static	void setupLumaTables();
24static int abs(int v) { return (v < 0 ? -v : v); }
25static double pct(int v,double range) { return (v*100. / range); }
26static void check(int R, int G, int B);
27
28float	D1, D2;
29float	D3, D4;
30float	D5, D6;
31
32int
33main(int argc, char** argv)
34{
35    int R, G, B;
36
37    if (argc > 1) {
38	refBlackWhite[0] = 16;
39	refBlackWhite[1] = 235;
40	refBlackWhite[2] = 128;
41	refBlackWhite[3] = 240;
42	refBlackWhite[4] = 128;
43	refBlackWhite[5] = 240;
44    }
45    D3 = 2 - 2*LumaRed;
46    D4 = 2 - 2*LumaBlue;
47    D1 = 1. / D3;
48    D2 = 1. / D4;
49    D5 = D3*LumaRed / LumaGreen;
50    D6 = D4*LumaBlue / LumaGreen;
51    setupLumaTables();
52    for (R = 0; R < 256; R++) {
53	for (G = 0; G < 256; G++)
54	    for (B = 0; B < 256; B++)
55		check(R, G, B);
56	printf("[%3u] c %u/%u b %u/%u (R %u/%d/%u G %u/%d/%u B %u/%d/%u)\n"
57	    , R
58	    , eCodes - preveCodes, eCodes
59	    , eBits - preveBits, eBits
60	    , abs(AbseRtotal - preveRtotal), eRtotal , AbseRtotal
61	    , abs(AbseGtotal - preveGtotal), eGtotal , AbseGtotal
62	    , abs(AbseBtotal - preveBtotal), eBtotal , AbseBtotal
63	);
64	preveRtotal = AbseRtotal;
65	preveGtotal = AbseGtotal;
66	preveBtotal = AbseBtotal;
67	preveCodes = eCodes;
68	preveBits = eBits;
69    }
70    printf("%u total codes\n", 256*256*256);
71    printf("total error: %u codes %u bits (R %d/%u G %d/%u B %d/%u)\n"
72	, eCodes
73	, eBits
74	, eRtotal , AbseRtotal
75	, eGtotal , AbseGtotal
76	, eBtotal , AbseBtotal
77    );
78    return (0);
79}
80
81float	*lumaRed;
82float	*lumaGreen;
83float	*lumaBlue;
84
85static float*
86setupLuma(float c)
87{
88    float *v = (float *)_TIFFmalloc(256 * sizeof (float));
89    int i;
90    for (i = 0; i < 256; i++)
91	v[i] = c * i;
92    return (v);
93}
94
95static void
96setupLumaTables(void)
97{
98    lumaRed = setupLuma(LumaRed);
99    lumaGreen = setupLuma(LumaGreen);
100    lumaBlue = setupLuma(LumaBlue);
101}
102
103static unsigned
104V2Code(float f, unsigned long RB, unsigned long RW, int CR)
105{
106    unsigned int c = (unsigned int)((((f)*(RW-RB)/CR)+RB)+.5);
107    return (c > 255 ? 255 : c);
108}
109
110#define	Code2V(c, RB, RW, CR)	((((c)-(int)RB)*(float)CR)/(float)(RW-RB))
111
112#define	CLAMP(f,min,max) \
113    (int)((f)+.5 < (min) ? (min) : (f)+.5 > (max) ? (max) : (f)+.5)
114
115static
116void
117check(int R, int G, int B)
118{
119    float Y, Cb, Cr;
120    int iY, iCb, iCr;
121    float rY, rCb, rCr;
122    float rR, rG, rB;
123    int eR, eG, eB;
124
125    Y = lumaRed[R] + lumaGreen[G] + lumaBlue[B];
126    Cb = (B - Y)*D2;
127    Cr = (R - Y)*D1;
128    iY = V2Code(Y, refBlackWhite[0], refBlackWhite[1], 255);
129    iCb = V2Code(Cb, refBlackWhite[2], refBlackWhite[3], 127);
130    iCr = V2Code(Cr, refBlackWhite[4], refBlackWhite[5], 127);
131    rCb = Code2V(iCb, refBlackWhite[2], refBlackWhite[3], 127);
132    rCr = Code2V(iCr, refBlackWhite[4], refBlackWhite[5], 127);
133    rY = Code2V(iY, refBlackWhite[0], refBlackWhite[1], 255);
134    rR = rY + rCr*D3;
135    rB = rY + rCb*D4;
136    rG = rY - rCb*D6 - rCr*D5;
137    eR = R - CLAMP(rR,0,255);
138    eG = G - CLAMP(rG,0,255);
139    eB = B - CLAMP(rB,0,255);
140    if (abs(eR) > 1 || abs(eG) > 1 || abs(eB) > 1) {
141	printf("R %u G %u B %u", R, G, B);
142	printf(" Y %g Cb %g Cr %g", Y, Cb, Cr);
143	printf(" iY %u iCb %u iCr %u", iY, iCb, iCr);
144	printf("\n -> Y %g Cb %g Cr %g", rY, rCb, rCr);
145	printf(" R %g (%u) G %g (%u) B %g (%u) E=[%d %d %d])\n"
146	    , rR, CLAMP(rR,0,255)
147	    , rG, CLAMP(rG,0,255)
148	    , rB, CLAMP(rB,0,255)
149	    , eR, eG, eB
150	);
151    }
152    eRtotal += eR;
153    eGtotal += eG;
154    eBtotal += eB;
155    AbseRtotal += abs(eR);
156    AbseGtotal += abs(eG);
157    AbseBtotal += abs(eB);
158    if (eR | eG | eB)
159	eCodes++;
160    eBits += abs(eR) + abs(eG) + abs(eB);
161}
162/*
163 * Local Variables:
164 * mode: c
165 * c-basic-offset: 8
166 * fill-column: 78
167 * End:
168 */
169