/*
 * Noise.cs
 * Authors:     August Zinsser
 * 
 * Copyright August Zinsser 2007
 * This program is distributed under the terms of the GNU General Public License
 */


using System;
using System.Collections.Generic;
using System.Text;
using Microsoft.Xna.Framework;


namespace Pina3D
{
    /// <summary>
    /// Generates a multi-dimensional field of solid noise. Range is the canonical random number [0,1). Note that the origin will always have
    /// a value of 0.5f.
    /// </summary>
    public class PseudoRandomMap
    {
        protected Vector3[] mGradientHash;      // Stores hashcodes for a random 3D vectors
        protected int[] mIndices;               // Used to look into the GradientField to produce a pseudorandom Gradient Field
        protected Vector3? mSize;               // Mod lookups by these numbers to make the field tesselate


        /// <summary>
        /// Create a new arbitrarily dimensional field of seeded solid noise. The first 3 Dimensions
        /// will tesselate seamlessly after the specified number of grid points.
        /// </summary>
        /// <param name="amplitude">Maximum value possible (min is 0)</param>
        /// <param name="seed">Seed for the pseudorandom number generator</param>
        /// <param name="xSize">Number of grid points in the X axis</param>
        /// <param name="ySize">Number of grid points in the Y axis</param>
        /// <param name="zSize">Number of grid points in the Z axis</param>
        public PseudoRandomMap(int seed, int xSize, int ySize, int zSize)
            : this(seed)
        {
            mSize = new Vector3(xSize, ySize, zSize);
        }


        /// <summary>
        /// Create a new arbitrarily dimensional field of seeded solid noise.
        /// </summary>
        /// <param name="seed"></param>
        public PseudoRandomMap(int seed)
        {
            Random rand = new Random(seed);

            // Build a random permutation of indices
            int n = 512;
            mIndices = new int[n];
            List<int> availableIndices = new List<int>(n);
            for (int i = 0; i < n; i++)
                availableIndices.Add(i);
            for (int i = 0; i < n; i++)
            {
                int take = (int)(rand.NextDouble() * availableIndices.Count);
                mIndices[i] = availableIndices[take];
                availableIndices.RemoveAt(take);
            }

            // Construct data that maps to the Gradient Field
            mGradientHash = new Vector3[n];
            for (int i = 0; i < n; i++)
            {
                // Calculate points within the unit sphere and then get the unit vector to those points
                // Throwing away the "corners" of the unit cube ensures that all directions have equal probability.
                float lengthSqrd = 2;
                while (lengthSqrd > 1)
                {
                    mGradientHash[i] = new Vector3((float)rand.NextDouble() * 2f - 1, (float)rand.NextDouble() * 2f - 1, (float)rand.NextDouble() * 2f - 1);
                    lengthSqrd = mGradientHash[i].LengthSquared();
                }
                mGradientHash[i].Normalize();
            }
        }


        /// <summary>
        /// Returns the bicubic-interpolated value at the specified coordinates.
        /// </summary>
        public virtual float LookUp(float x, float y, float z)
        {
            // Wrap the coords if possible (creates seamless maps)
            if (mSize.HasValue)
            {
                x %= mSize.Value.X;
                y %= mSize.Value.Y;
                z %= mSize.Value.Z;
            }

            // Find the 8 nearest grid points
            int x0 = (int)Math.Floor(x);
            int x1 = x0 + 1;
            int y0 = (int)Math.Floor(y);
            int y1 = y0 + 1;
            int z0 = (int)Math.Floor(z);
            int z1 = z0 + 1;

            // Calculate Solid Noise (similar to Perlin's algorithm)
            // For each neighbor, find the dot product of the current point (x,y) and the gradient of said neighbor to offset the local max/mins and
            // get a magnitude for the current point. The final value is the weighted sum (using cubic interpolation) of the 4 nearest neighbors'
            // calculated magnitude.
            float result =
                GetWeightedValue(x - x0, y - y0, z - z0, GetGradientVector(x0, y0, z0)) +
                GetWeightedValue(x - x1, y - y0, z - z0, GetGradientVector(x1, y0, z0)) +
                GetWeightedValue(x - x0, y - y1, z - z0, GetGradientVector(x0, y1, z0)) +
                GetWeightedValue(x - x1, y - y1, z - z0, GetGradientVector(x1, y1, z0)) +
                GetWeightedValue(x - x0, y - y0, z - z1, GetGradientVector(x0, y0, z1)) +
                GetWeightedValue(x - x1, y - y0, z - z1, GetGradientVector(x1, y0, z1)) +
                GetWeightedValue(x - x0, y - y1, z - z1, GetGradientVector(x0, y1, z1)) +
                GetWeightedValue(x - x1, y - y1, z - z1, GetGradientVector(x1, y1, z1));

            // Interpolate (again) the result to increase contrast
            return Logistic(result);
        }


        /// <summary>
        /// Returns the bicubic-interpolated value at the specified coordinates.
        /// </summary>
        public virtual float LookUp(float x, float y)
        {
            return LookUp(x, y, 0);
        }


        /// <summary>
        /// Returns the bicubic-interpolated value at the specified coordinates.
        /// </summary>
        public virtual float LookUp(float x)
        {
            return LookUp(x, 0, 0);
        }


        /*
         * Helper Functions
         */

        private float GetWeightedValue(float u, float v, Vector2 value)
        {
            return CubicWeight(u) * CubicWeight(v) * (Vector2.Dot(value, new Vector2(u, v)));
        }


        private float GetWeightedValue(float u, float v, float w, Vector3 value)
        {
            return CubicWeight(u) * CubicWeight(v) * CubicWeight(w) * (Vector3.Dot(value, new Vector3(u, v, w)));
        }


        /// <summary>
        /// Bicubic Convolution Algorithm. Weights each point for bicubic interpolation.
        /// </summary>
        /// <param name="t"></param>
        /// <returns></returns>
        private float CubicWeight(float t)
        {
            t = (float)Math.Abs(t);
            float a = 0f;   // Edge-cases are not handled in LookUp so moving a away from 0 exaggerates artifacts on the gridlines
            if (t <= 1)
                return (float)((a + 2) * Math.Pow(t, 3) - (a + 3) * Math.Pow(t, 2)) + 1;
            else if (t <= 2)
                return (float)(a * Math.Pow(t, 3) - 5 * a * Math.Pow(t, 2) + 8 * a * t - 4 * a);
            else
                return 0;
        }


        /// <summary>
        /// Logistic Function, effectifly pushing the parameter closer to the endpoints
        /// </summary>
        /// <param name="t"></param>
        /// <returns></returns>
        private float Logistic(float t)
        {
            t *= 6;
            t = 1f / (1f + (float)Math.Pow(MathHelper.E, -1 * t));
            return t;
        }


        private Vector3 GetGradientVector(int x, int y, int z)
        {
            return mGradientHash[GetWrappedIndex(x + GetWrappedIndex(y + GetWrappedIndex(z)))];
        }


        private int GetWrappedIndex(int i)
        {
            int index = i % mIndices.Length;
            if (index < 0)
                index += mIndices.Length;
            return mIndices[index];
        }
    }


    /// <summary>
    /// Combines random maps of decreasing frequencies and amplitude to create fractal-like maps of perlin noise.
    /// </summary>
    public class PerlinNoise : PseudoRandomMap
    {
        float mInverseAmplitude;        // Multiply results by this number to scale sampled spots down to the (-1,1) range
        float mInversePersistence;      // Relates the amplitude to the frequency of the current iteration
        int mMaxIterations;             // How many times subdivision may occur
        float mAmplitudeThreshold;      // Cutoff for adding subdivisions


        /// <summary>
        /// High persistence creates smoother fields, and low persistence creates rockier fields. This value must be > 1, and default
        /// value is 2 (values <= 1 are ignored).
        /// </summary>
        public float Persistence { set { if (value > 1) mInversePersistence = 1 / value; } }


        /// <summary>
        /// Create a new arbitrarily dimensional field of seeded perlin noise.
        /// </summary>
        /// <param name="seed"></param>
        public PerlinNoise(int seed)
            : base(seed)
        {
            mInversePersistence = .5f;
            mMaxIterations = 10;
            mAmplitudeThreshold = .05f;

            // Find the maximum value (the sum of the max of each iteration) and then take the inverse
            mInverseAmplitude = 0;
            float currentAmplitude = 1;
            for (int i = 0; i < mMaxIterations && currentAmplitude >= mAmplitudeThreshold; i++)
            {
                currentAmplitude = (float)Math.Pow(mInversePersistence, i);
                mInverseAmplitude += currentAmplitude;
            }
            mInverseAmplitude = 1f / mInverseAmplitude;
        }


        /// <summary>
        /// Return the bicubic-interpolated value at the specified coordinates.
        /// </summary>
        public override float LookUp(float x, float y, float z)
        {
            float result = 0;
            float amplitude = 1;

            for (int i = 0; i < mMaxIterations && amplitude >= mAmplitudeThreshold; i++)
            {
                float frequency = (float)Math.Pow(2, i);
                amplitude = (float)Math.Pow(mInversePersistence, i);
                result += base.LookUp(x * frequency, y * frequency, z * frequency) * amplitude;
            }

            // Interpolate the result to increase contrast
            return result * mInverseAmplitude;
        }


        /// <summary>
        /// Performs a spline interpolation between -1 and 1 to increase contrast
        /// </summary>
        /// <param name="t"></param>
        /// <returns></returns>
        private float SplineInterpolate(float t, float slope)
        {
            return (2 * (float)Math.Pow(t, 3) - 3 * (float)Math.Pow(t, 2));
        }
    }
}