using System;
using System.Linq;
using System.Collections.Generic;

namespace FullScreenShaders.Structure
{
    static class GaussianKernel
    {
        static readonly Dictionary<int, double> idealStandardDeviation = new Dictionary<int,double>
        {
            { 5, 7.013915463849E-01 }, 
            { 7, 8.09171316279E-01 },
            { 9, 9.0372907227E-01 }, 
            { 11, 9.890249035E-01 }, 
            { 13, 1.067359295 }, 
            { 15, 1.1402108 },
            { 17, 1.2086 }
        };

        public static double IdealStandardDeviation(int taps)
        {
            if (!idealStandardDeviation.ContainsKey(taps))
                throw new ArgumentException("No ideal standard deviation for this tap count : " + taps, "taps");

            return idealStandardDeviation[taps];
        }

        public static float[] GetWeights(int taps)
        {
            return GetWeights(taps, IdealStandardDeviation(taps));
        }
        public static float[] GetWeights(int taps, double standardDeviation)
        {
            if (taps < 1 || taps % 2 == 0)
                throw new ArgumentException("Invalid tap count : " + taps, "taps");

            return GetWeightsInternal(taps, standardDeviation).Cast<float>().ToArray();
        }

        public static void GetFastWeights(int taps, out float[] weights, out float[] offsets)
        {
            GetFastWeights(taps, IdealStandardDeviation(taps), out weights, out offsets);
        }
        public static void GetFastWeights(int taps, double standardDeviation, out float[] weights, out float[] offsets)
        {
            if (taps < 1 || taps % 2 == 0)
                throw new ArgumentException("Invalid tap count : " + taps, "taps");

            double[] tempWeights = GetWeightsInternal(taps, standardDeviation);

            bool hasEndTap = tempWeights.Length % 2 == 0;
            int joinedTaps = (tempWeights.Length - 1) / 2,
                totalTaps = joinedTaps + 1 + (hasEndTap ? 1 : 0);

            offsets = new float[totalTaps];
            weights = new float[totalTaps];

            weights[0] = (float)tempWeights[0];
            offsets[0] = 0;

            for (int i = 0; i < joinedTaps; i++)
            {
                double sum = tempWeights[i * 2 + 1] + tempWeights[i * 2 + 2];
                weights[i + 1] = (float)sum;
                offsets[i + 1] = (float)(0.5 - tempWeights[i * 2 + 1] / sum);
            }

            if (hasEndTap)
                weights[weights.Length - 1] = (float)tempWeights[tempWeights.Length - 1];
        }

        static double[] GetWeightsInternal(int taps, double standardDeviation)
        {
            int halfTaps = (taps - 1) / 2 + 1;

            var weights = new double[halfTaps];
            for (int i = 0; i < halfTaps; i++)
                weights[i] = Gaussian(i, standardDeviation);

            double augmentationFactor = 1 + LostLight(taps, standardDeviation);
            for (int i = 0; i < halfTaps; i++)
                weights[i] *= augmentationFactor;

            return weights;
        }

        public static double LostLight(int taps, double standardDeviation)
        {
            int halfTaps = (taps - 1) / 2 + 1;

            double sum = Gaussian(0, standardDeviation);
            for (int i = 1; i < halfTaps; i++)
                sum += Gaussian(i, standardDeviation) * 2;

            return -1 + 1 / sum;
        }

        public static double BoxLikeness(int taps, double standardDeviation)
        {
            int halfTaps = (taps - 1) / 2 + 1;

            var weights = new double[halfTaps];
            for (int i = 0; i < halfTaps; i++)
                weights[i] = Gaussian(i, standardDeviation);

            const double RelativeEpsilon = 1E-14;
            double invTaps = 1.0 / taps;
            IEnumerable<double> boxWeights = Enumerable.Repeat(invTaps + RelativeEpsilon, 1)
                .Concat(Enumerable.Repeat(invTaps - RelativeEpsilon / (halfTaps - 1), halfTaps - 1));

            double correlation = 1 - weights.Correlation(boxWeights);

            const double minAlpha = -1.08975982634075E-03, minA = 1.56809355440769E-02, minB = -3.46280631112407E-02, minC = 6.75067153072151E-02;
            const double maxAlpha = 1.78780869266509E-03, maxA = -3.68630036049787E-02, maxB = 2.59171745554970E-01, maxC = -2.49537606767719E-01;

            double minCorrected = correlation - (minAlpha * Math.Pow(halfTaps, 3) + minA * Math.Pow(halfTaps, 2) + minB * halfTaps + minC);
            double maxCorrected = minCorrected / (maxAlpha * Math.Pow(halfTaps, 3) + maxA * Math.Pow(halfTaps, 2) + maxB * halfTaps + maxC);

            return Math.Min(Math.Max(maxCorrected, 0), 1);
        }

        static double Gaussian(int distance, double standardDeviation)
        {
            return 1 / (Math.Sqrt(2 * Math.PI) * standardDeviation) *
                Math.Exp(-Math.Pow(distance, 2) / (2 * Math.Pow(standardDeviation, 2)));
        }
    }
}
