/*
  Metric
  Copyright (C) 2006 Yangli Hector Yee

  This program 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/

#include "Metric.h"
#include "CompareArgs.h"
#include "RGBAImage.h"
#include "LPyramid.h"
#include <math.h>
#include "pdiff.h"

#ifndef M_PI
#define M_PI 3.14159265f
#endif

/*
 * Given the adaptation luminance, this function returns the
 * threshold of visibility in cd per m^2
 * TVI means Threshold vs Intensity function
 * This version comes from Ward Larson Siggraph 1997
 */

float tvi(float adaptation_luminance)
{
    /* returns the threshold luminance given the adaptation luminance
       units are candelas per meter squared
    */

    float log_a, r, result;
    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;
    } else {
	r = log_a - 1.255f;
    }

    result = powf(10.0f , r);

    return result;

}

/* computes the contrast sensitivity function (Barten SPIE 1989)
 * given the cycles per degree (cpd) and luminance (lum)
 */
float csf(float cpd, float lum)
{
    float a, b, result;

    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;
}

/*
 * Visual Masking Function
 * from Daly 1993
 */
float mask(float contrast)
{
    float a, b, result;
    a = powf(392.498f * contrast,  0.7f);
    b = powf(0.0153f * a, 4.0f);
    result = powf(1.0f + b, 0.25f);

    return result;
}

/* convert Adobe RGB (1998) with reference white D65 to XYZ */
void AdobeRGBToXYZ(float r, float g, float b, float &x, float &y, float &z)
{
    /* matrix is from http://www.brucelindbloom.com/ */
    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;
}

void XYZToLAB(float x, float y, float z, float &L, float &A, float &B)
{
    static float xw = -1;
    static float yw;
    static float zw;
    /* reference white */
    if (xw < 0) {
	AdobeRGBToXYZ(1, 1, 1, xw, yw, zw);
    }
    const float epsilon  = 216.0f / 24389.0f;
    const float kappa = 24389.0f / 27.0f;
    float f[3];
    float r[3];
    r[0] = x / xw;
    r[1] = y / yw;
    r[2] = z / zw;
    for (int i = 0; i < 3; i++) {
	if (r[i] > epsilon) {
	    f[i] = powf(r[i], 1.0f / 3.0f);
	} else {
	    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]);
}

int Yee_Compare_Images(RGBAImage *image_a,
		       RGBAImage *image_b,
		       float gamma,
		       float luminance,
		       float field_of_view,
		       bool verbose);

bool Yee_Compare(CompareArgs &args)
{
    if ((args.ImgA->Get_Width() != args.ImgB->Get_Width()) ||
	(args.ImgA->Get_Height() != args.ImgB->Get_Height())) {
	args.ErrorStr = "Image dimensions do not match\n";
	return false;
    }

    unsigned int i, dim, pixels_failed;
    dim = args.ImgA->Get_Width() * args.ImgA->Get_Height();
    bool identical = true;
    for (i = 0; i < dim; i++) {
	if (args.ImgA->Get(i) != args.ImgB->Get(i)) {
	    identical = false;
	    break;
	}
    }
    if (identical) {
	args.ErrorStr = "Images are binary identical\n";
	return true;
    }

    pixels_failed = Yee_Compare_Images (args.ImgA, args.ImgB,
					args.Gamma, args.Luminance,
					args.FieldOfView, args.Verbose);

    if (pixels_failed < args.ThresholdPixels) {
	args.ErrorStr = "Images are perceptually indistinguishable\n";
	return true;
    }

    char different[100];
    sprintf(different, "%d pixels are different\n", pixels_failed);

    args.ErrorStr = "Images are visibly different\n";
    args.ErrorStr += different;

    if (args.ImgDiff) {
#if IMAGE_DIFF_CODE_ENABLED
	if (args.ImgDiff->WritePPM()) {
	    args.ErrorStr += "Wrote difference image to ";
	    args.ErrorStr+= args.ImgDiff->Get_Name();
	    args.ErrorStr += "\n";
	} else {
	    args.ErrorStr += "Could not write difference image to ";
	    args.ErrorStr+= args.ImgDiff->Get_Name();
	    args.ErrorStr += "\n";
	}
#endif
	args.ErrorStr += "Generation of image \"difference\" is currently disabled\n";
    }
    return false;
}

int
pdiff_compare (cairo_surface_t *surface_a,
	       cairo_surface_t *surface_b,
	       double gamma,
	       double luminance,
	       double field_of_view)
{
    RGBAImage *image_a, *image_b;

    image_a = new RGBACairoImage (surface_a);
    image_b = new RGBACairoImage (surface_b);

    return Yee_Compare_Images (image_a, image_b,
			       gamma, luminance,
			       field_of_view, false);
}

int Yee_Compare_Images(RGBAImage *image_a,
		       RGBAImage *image_b,
		       float gamma,
		       float luminance,
		       float field_of_view,
		       bool verbose)
{
    unsigned int i, dim;
    dim = image_a->Get_Width() * image_a->Get_Height();

    /* assuming colorspaces are in Adobe RGB (1998) convert to XYZ */
    float *aX = new float[dim];
    float *aY = new float[dim];
    float *aZ = new float[dim];
    float *bX = new float[dim];
    float *bY = new float[dim];
    float *bZ = new float[dim];
    float *aLum = new float[dim];
    float *bLum = new float[dim];

    float *aA = new float[dim];
    float *bA = new float[dim];
    float *aB = new float[dim];
    float *bB = new float[dim];

    if (verbose) printf("Converting RGB to XYZ\n");

    unsigned int x, y, w, h;
    w = image_a->Get_Width();
    h = image_a->Get_Height();
    for (y = 0; y < h; y++) {
	for (x = 0; x < w; x++) {
	    float r, g, b, l;
	    i = x + y * w;
	    r = powf(image_a->Get_Red(i) / 255.0f, gamma);
	    g = powf(image_a->Get_Green(i) / 255.0f, gamma);
	    b = powf(image_a->Get_Blue(i) / 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(image_b->Get_Red(i) / 255.0f, gamma);
	    g = powf(image_b->Get_Green(i) / 255.0f, gamma);
	    b = powf(image_b->Get_Blue(i) / 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;
	}
    }

    if (verbose) printf("Constructing Laplacian Pyramids\n");

    LPyramid *la = new LPyramid(aLum, w, h);
    LPyramid *lb = new LPyramid(bLum, w, h);

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

    if (verbose) printf("Performing test\n");

    float num_pixels = 1;
    unsigned int 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;
    }

    float cpd[MAX_PYR_LEVELS];
    cpd[0] = 0.5f * pixels_per_degree;
    for (i = 1; i < MAX_PYR_LEVELS; i++) cpd[i] = 0.5f * cpd[i - 1];
    float csf_max = csf(3.248f, 100.0f);

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

    unsigned int pixels_failed = 0;
    for (y = 0; y < h; y++) {
	for (x = 0; x < w; x++) {
	    int index = x + y * w;
	    float contrast[MAX_PYR_LEVELS - 2];
	    float sum_contrast = 0;
	    for (i = 0; i < MAX_PYR_LEVELS - 2; i++) {
		float n1 = fabsf(la->Get_Value(x,y,i) - la->Get_Value(x,y,i + 1));
		float n2 = fabsf(lb->Get_Value(x,y,i) - lb->Get_Value(x,y,i + 1));
		float numerator = (n1 > n2) ? n1 : n2;
		float d1 = fabsf(la->Get_Value(x,y,i+2));
		float d2 = fabsf(lb->Get_Value(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;
	    float F_mask[MAX_PYR_LEVELS - 2];
	    float adapt = la->Get_Value(x,y,adaptation_level) + lb->Get_Value(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));
	    }
	    float 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;
	    float delta = fabsf(la->Get_Value(x,y,0) - lb->Get_Value(x,y,0));
	    bool pass = true;
	    /* pure luminance test */
	    if (delta > factor * tvi(adapt)) {
		pass = false;
	    } else {
		/* CIE delta E test with modifications */
		float color_scale = 1.0f;
		/* ramp down the color test in scotopic regions */
		if (adapt < 10.0f) {
		    color_scale = 1.0f - (10.0f - color_scale) / 10.0f;
		    color_scale = color_scale * color_scale;
		}
		float da = aA[index] - bA[index];
		float db = aB[index] - bB[index];
		da = da * da;
		db = db * db;
		float delta_e = (da + db) * color_scale;
		if (delta_e > factor) {
		    pass = false;
		}
	    }
	    if (!pass)
		pixels_failed++;
#if IMAGE_DIFF_ENABLED
	    if (!pass) {
		if (args.ImgDiff) {
		    args.ImgDiff->Set(255, 0, 0, 255, index);
		}
	    } else {
		if (args.ImgDiff) {
		    args.ImgDiff->Set(0, 0, 0, 255, index);
		}
	    }
#endif
	}
    }

    if (aX) delete[] aX;
    if (aY) delete[] aY;
    if (aZ) delete[] aZ;
    if (bX) delete[] bX;
    if (bY) delete[] bY;
    if (bZ) delete[] bZ;
    if (aLum) delete[] aLum;
    if (bLum) delete[] bLum;
    if (la) delete la;
    if (lb) delete lb;
    if (aA) delete aA;
    if (bA) delete bA;
    if (aB) delete aB;
    if (bB) delete bB;

    return pixels_failed;
}