tvi (float adaptation_luminance)

    log_a = log10f(adaptation_luminance);

    if (log_a < -3.94f) {

	r = -2.86f;

    } else if (log_a < -1.44f) {

	r = powf(0.405f * log_a + 1.6f , 2.18f) - 2.86f;

    } else if (log_a < -0.0184f) {

	r = log_a - 0.395f;

    } else if (log_a < 1.9f) {

	r = powf(0.249f * log_a + 0.65f, 2.7f) - 0.72f;

	r = log_a - 1.255f;

    result = powf(10.0f , r);

    return result;

csf (float cpd, float lum)

    a = 440.0f * powf((1.0f + 0.7f / lum), -0.2f);

    b = 0.3f * powf((1.0f + 100.0f / lum), 0.15f);

    result = a * cpd * expf(-b * cpd) * sqrtf(1.0f + 0.06f * expf(b * cpd));

    return result;

mask (float contrast)

    a = powf(392.498f * contrast,  0.7f);

    b = powf(0.0153f * a, 4.0f);

    result = powf(1.0f + b, 0.25f);

    return result;

AdobeRGBToXYZ (float r, float g, float b, float *x, float *y, float *z)

    *x = r * 0.576700f + g * 0.185556f + b * 0.188212f;

    *y = r * 0.297361f + g * 0.627355f + b * 0.0752847f;

    *z = r * 0.0270328f + g * 0.0706879f + b * 0.991248f;

}

XYZToLAB (float x, float y, float z, float *L, float *A, float *B)

    const float epsilon  = 216.0f / 24389.0f;

    const float kappa = 24389.0f / 27.0f;

    if (xw < 0) {

	AdobeRGBToXYZ(1, 1, 1, &xw, &yw, &zw);

    r[0] = x / xw;

    r[1] = y / yw;

    r[2] = z / zw;

    for (i = 0; i < 3; i++) {

	if (r[i] > epsilon) {

	    f[i] = powf(r[i], 1.0f / 3.0f);

	    f[i] = (kappa * r[i] + 16.0f) / 116.0f;

    *L = 116.0f * f[1] - 16.0f;

    *A = 500.0f * (f[0] - f[1]);

    *B = 200.0f * (f[1] - f[2]);

}

_get_pixel (const uint32_t *data, int i, cairo_format_t format)

    if (format == CAIRO_FORMAT_ARGB32)

	return data[i];

	return data[i] | 0xff000000;

_get_red (const uint32_t *data, int i, cairo_format_t format)

    pixel = _get_pixel (data, i, format);

    alpha = (pixel & 0xff000000) >> 24;

    if (alpha == 0)

	return 0;

	return (((pixel & 0x00ff0000) >> 16) * 255 + alpha / 2) / alpha;

_get_green (const uint32_t *data, int i, cairo_format_t format)

    pixel = _get_pixel (data, i, format);

    alpha = (pixel & 0xff000000) >> 24;

    if (alpha == 0)

	return 0;

	return (((pixel & 0x0000ff00) >> 8) * 255 + alpha / 2) / alpha;

_get_blue (const uint32_t *data, int i, cairo_format_t format)

    pixel = _get_pixel (data, i, format);

    alpha = (pixel & 0xff000000) >> 24;

    if (alpha == 0)

	return 0;

	return (((pixel & 0x000000ff) >> 0) * 255 + alpha / 2) / alpha;

xmalloc (size_t size)

    buf = malloc (size);

    if (buf == NULL) {

    return buf;

pdiff_compare (cairo_surface_t *surface_a,

    unsigned int dim = (cairo_image_surface_get_width (surface_a)

			* cairo_image_surface_get_height (surface_a));

    w = cairo_image_surface_get_width (surface_a);

    h = cairo_image_surface_get_height (surface_a);

    if (w < 3 || h < 3) /* too small for the Laplacian convolution */

	return -1;

    format_a = cairo_image_surface_get_format (surface_a);

    format_b = cairo_image_surface_get_format (surface_b);

    assert (format_a == CAIRO_FORMAT_RGB24 || format_a == CAIRO_FORMAT_ARGB32);

    assert (format_b == CAIRO_FORMAT_RGB24 || format_b == CAIRO_FORMAT_ARGB32);

    aX = xmalloc (dim * sizeof (float));

    aY = xmalloc (dim * sizeof (float));

    aZ = xmalloc (dim * sizeof (float));

    bX = xmalloc (dim * sizeof (float));

    bY = xmalloc (dim * sizeof (float));

    bZ = xmalloc (dim * sizeof (float));

    aLum = xmalloc (dim * sizeof (float));

    bLum = xmalloc (dim * sizeof (float));

    aA = xmalloc (dim * sizeof (float));

    bA = xmalloc (dim * sizeof (float));

    aB = xmalloc (dim * sizeof (float));

    bB = xmalloc (dim * sizeof (float));

    data_a = (uint32_t *) cairo_image_surface_get_data (surface_a);

    data_b = (uint32_t *) cairo_image_surface_get_data (surface_b);

    for (y = 0; y < h; y++) {

	for (x = 0; x < w; x++) {

	    i = x + y * w;

	    r = powf(_get_red (data_a, i, format_a) / 255.0f, gamma);

	    g = powf(_get_green (data_a, i, format_a) / 255.0f, gamma);

	    b = powf(_get_blue (data_a, i, format_a) / 255.0f, gamma);

	    AdobeRGBToXYZ(r,g,b,&aX[i],&aY[i],&aZ[i]);

	    XYZToLAB(aX[i], aY[i], aZ[i], &l, &aA[i], &aB[i]);

	    r = powf(_get_red (data_b, i, format_b) / 255.0f, gamma);

	    g = powf(_get_green (data_b, i, format_b) / 255.0f, gamma);

	    b = powf(_get_blue (data_b, i, format_b) / 255.0f, gamma);

	    AdobeRGBToXYZ(r,g,b,&bX[i],&bY[i],&bZ[i]);

	    XYZToLAB(bX[i], bY[i], bZ[i], &l, &bA[i], &bB[i]);

	    aLum[i] = aY[i] * luminance;

	    bLum[i] = bY[i] * luminance;

    la = lpyramid_create (aLum, w, h);

    lb = lpyramid_create (bLum, w, h);

    num_one_degree_pixels = (float) (2 * tan(field_of_view * 0.5 * M_PI / 180) * 180 / M_PI);

    pixels_per_degree = w / num_one_degree_pixels;

    num_pixels = 1;

    adaptation_level = 0;

    for (i = 0; i < MAX_PYR_LEVELS; i++) {

	adaptation_level = i;

	if (num_pixels > num_one_degree_pixels) break;

	num_pixels *= 2;

    cpd[0] = 0.5f * pixels_per_degree;

    for (i = 1; i < MAX_PYR_LEVELS; i++) cpd[i] = 0.5f * cpd[i - 1];

    csf_max = csf(3.248f, 100.0f);

    for (i = 0; i < MAX_PYR_LEVELS - 2; i++) F_freq[i] = csf_max / csf( cpd[i], 100.0f);

    pixels_failed = 0;

    for (y = 0; y < h; y++) {

	for (x = 0; x < w; x++) {

	    int index = x + y * w;

	    float sum_contrast = 0;

	    for (i = 0; i < MAX_PYR_LEVELS - 2; i++) {

		float n1 = fabsf(lpyramid_get_value (la,x,y,i) - lpyramid_get_value (la,x,y,i + 1));

		float n2 = fabsf(lpyramid_get_value (lb,x,y,i) - lpyramid_get_value (lb,x,y,i + 1));

		float numerator = (n1 > n2) ? n1 : n2;

		float d1 = fabsf(lpyramid_get_value(la,x,y,i+2));

		float d2 = fabsf(lpyramid_get_value(lb,x,y,i+2));

		float denominator = (d1 > d2) ? d1 : d2;

		if (denominator < 1e-5f) denominator = 1e-5f;

		contrast[i] = numerator / denominator;

		sum_contrast += contrast[i];

	    if (sum_contrast < 1e-5) sum_contrast = 1e-5f;

	    adapt = lpyramid_get_value(la,x,y,adaptation_level) + lpyramid_get_value(lb,x,y,adaptation_level);

	    adapt *= 0.5f;

	    if (adapt < 1e-5) adapt = 1e-5f;

	    for (i = 0; i < MAX_PYR_LEVELS - 2; i++) {

		F_mask[i] = mask(contrast[i] * csf(cpd[i], adapt));

	    factor = 0;

	    for (i = 0; i < MAX_PYR_LEVELS - 2; i++) {

		factor += contrast[i] * F_freq[i] * F_mask[i] / sum_contrast;

	    if (factor < 1) factor = 1;

	    if (factor > 10) factor = 10;

	    delta = fabsf(lpyramid_get_value(la,x,y,0) - lpyramid_get_value(lb,x,y,0));

	    pass = true;

	    if (delta > factor * tvi(adapt)) {

		pass = false;

		float color_scale = 1.0f;

		float da = aA[index] - bA[index];

		float db = aB[index] - bB[index];

		if (adapt < 10.0f) {

		    color_scale = 1.0f - (10.0f - color_scale) / 10.0f;

		    color_scale = color_scale * color_scale;

		da = da * da;

		db = db * db;

		delta_e = (da + db) * color_scale;

		if (delta_e > factor) {

		    pass = false;

	    if (!pass)

		pixels_failed++;

    free (aX);

    free (aY);

    free (aZ);

    free (bX);

    free (bY);

    free (bZ);

    free (aLum);

    free (bLum);

    lpyramid_destroy (la);

    lpyramid_destroy (lb);

    free (aA);

    free (bA);

    free (aB);

    free (bB);

    return pixels_failed;