float	ycbcrCoeffs[3] = { .299, .587, .114 };

/* default coding range is CCIR Rec 601-1 with no headroom/footroom */

unsigned long refBlackWhite[6] = { 0, 255, 128, 255, 128, 255 };

#define	LumaRed		ycbcrCoeffs[0]

#define	LumaGreen	ycbcrCoeffs[1]

#define	LumaBlue	ycbcrCoeffs[2]

long	eRtotal = 0;

long	eGtotal = 0;

long	eBtotal = 0;

long	preveRtotal = 0;

long	preveGtotal = 0;

long	preveBtotal = 0;

unsigned long AbseRtotal = 0;

unsigned long AbseGtotal = 0;

unsigned long AbseBtotal = 0;

unsigned long eCodes = 0;

unsigned long preveCodes = 0;

unsigned long eBits = 0;

unsigned long preveBits = 0;

static	void setupLumaTables();

static int abs(int v) { return (v < 0 ? -v : v); }

static double pct(int v,double range) { return (v*100. / range); }

static void check(int R, int G, int B);

float	D1, D2;

float	D3, D4;

float	D5, D6;

int

main(int argc, char** argv)

{

    int R, G, B;

    if (argc > 1) {

	refBlackWhite[0] = 16;

	refBlackWhite[1] = 235;

	refBlackWhite[2] = 128;

	refBlackWhite[3] = 240;

	refBlackWhite[4] = 128;

	refBlackWhite[5] = 240;

    }

    D3 = 2 - 2*LumaRed;

    D4 = 2 - 2*LumaBlue;

    D1 = 1. / D3;

    D2 = 1. / D4;

    D5 = D3*LumaRed / LumaGreen;

    D6 = D4*LumaBlue / LumaGreen;

    setupLumaTables();

    for (R = 0; R < 256; R++) {

	for (G = 0; G < 256; G++)

	    for (B = 0; B < 256; B++)

		check(R, G, B);

	printf("[%3u] c %u/%u b %u/%u (R %u/%d/%u G %u/%d/%u B %u/%d/%u)\n"

	    , R

	    , eCodes - preveCodes, eCodes

	    , eBits - preveBits, eBits

	    , abs(AbseRtotal - preveRtotal), eRtotal , AbseRtotal

	    , abs(AbseGtotal - preveGtotal), eGtotal , AbseGtotal

	    , abs(AbseBtotal - preveBtotal), eBtotal , AbseBtotal

	);

	preveRtotal = AbseRtotal;

	preveGtotal = AbseGtotal;

	preveBtotal = AbseBtotal;

	preveCodes = eCodes;

	preveBits = eBits;

    }

    printf("%u total codes\n", 256*256*256);

    printf("total error: %u codes %u bits (R %d/%u G %d/%u B %d/%u)\n"

	, eCodes

	, eBits

	, eRtotal , AbseRtotal

	, eGtotal , AbseGtotal

	, eBtotal , AbseBtotal

    );

    return (0);

}

float	*lumaRed;

float	*lumaGreen;

float	*lumaBlue;

static float*

setupLuma(float c)

{

    float *v = (float *)_TIFFmalloc(256 * sizeof (float));

    int i;

    for (i = 0; i < 256; i++)

	v[i] = c * i;

    return (v);

}

static void

setupLumaTables(void)

{

    lumaRed = setupLuma(LumaRed);

    lumaGreen = setupLuma(LumaGreen);

    lumaBlue = setupLuma(LumaBlue);

}

static unsigned

V2Code(float f, unsigned long RB, unsigned long RW, int CR)

{

    unsigned int c = (unsigned int)((((f)*(RW-RB)/CR)+RB)+.5);

    return (c > 255 ? 255 : c);

}

#define	Code2V(c, RB, RW, CR)	((((c)-(int)RB)*(float)CR)/(float)(RW-RB))

#define	CLAMP(f,min,max) \

    (int)((f)+.5 < (min) ? (min) : (f)+.5 > (max) ? (max) : (f)+.5)

static

void

check(int R, int G, int B)

{

    float Y, Cb, Cr;

    int iY, iCb, iCr;

    float rY, rCb, rCr;

    float rR, rG, rB;

    int eR, eG, eB;

    Y = lumaRed[R] + lumaGreen[G] + lumaBlue[B];

    Cb = (B - Y)*D2;

    Cr = (R - Y)*D1;

    iY = V2Code(Y, refBlackWhite[0], refBlackWhite[1], 255);

    iCb = V2Code(Cb, refBlackWhite[2], refBlackWhite[3], 127);

    iCr = V2Code(Cr, refBlackWhite[4], refBlackWhite[5], 127);

    rCb = Code2V(iCb, refBlackWhite[2], refBlackWhite[3], 127);

    rCr = Code2V(iCr, refBlackWhite[4], refBlackWhite[5], 127);

    rY = Code2V(iY, refBlackWhite[0], refBlackWhite[1], 255);

    rR = rY + rCr*D3;

    rB = rY + rCb*D4;

    rG = rY - rCb*D6 - rCr*D5;

    eR = R - CLAMP(rR,0,255);

    eG = G - CLAMP(rG,0,255);

    eB = B - CLAMP(rB,0,255);

    if (abs(eR) > 1 || abs(eG) > 1 || abs(eB) > 1) {

	printf("R %u G %u B %u", R, G, B);

	printf(" Y %g Cb %g Cr %g", Y, Cb, Cr);

	printf(" iY %u iCb %u iCr %u", iY, iCb, iCr);

	printf("\n -> Y %g Cb %g Cr %g", rY, rCb, rCr);

	printf(" R %g (%u) G %g (%u) B %g (%u) E=[%d %d %d])\n"

	    , rR, CLAMP(rR,0,255)

	    , rG, CLAMP(rG,0,255)

	    , rB, CLAMP(rB,0,255)

	    , eR, eG, eB

	);

    }

    eRtotal += eR;

    eGtotal += eG;

    eBtotal += eB;

    AbseRtotal += abs(eR);

    AbseGtotal += abs(eG);

    AbseBtotal += abs(eB);

    if (eR | eG | eB)

	eCodes++;

    eBits += abs(eR) + abs(eG) + abs(eB);

}

/*

 * Local Variables:

 * mode: c

 * c-basic-offset: 8

 * fill-column: 78

 * End:

 */