#define  MINFILTWIDTH 1.0e-6
#include "perlin.h"

#ifndef M_PI
   #define       M_PI 3.14159265358979323846
#endif

/* Antialiased abs().
 * Compute the box filter of abs(t) from x-dx/2 to x+dx/2.
 * Hinges on the realization that the indefinite integral of abs(x) is
 * sign(x) * 1/2 x*x;
 */
float integral (float t)
{
        if (t < 0) {return -(0.5 * t*t);}
        else {return (0.5 * t*t);}
}

float length (float v1[3])
{
        float square;
        square = v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2];
        return sqrt(square);
}

float * mix (float c1[3], float c2[3], float mixvalue)
{
        float *v;

        v=malloc(sizeof(float) * 3);
        v[0] = ((1.0 - mixvalue)*c1[0]) + mixvalue*c2[0];
        v[1] = ((1.0 - mixvalue)*c1[1]) + mixvalue*c2[1];
        v[2] = ((1.0 - mixvalue)*c1[2]) + mixvalue*c2[2];

        return v;
}

float smoothstep (float e1, float e2, float x)
{
        float interval, x2, y;

        if (x < e1) return 0.0;
        if (x > e2) return 1.0;
        interval = fabs(e2- e1);
        x2 = (x - e1)/interval;
        y = (x2*x2)*(3 - (2*x2));
        return y;
}

float noise(float p[3])
{
        float intensity;

        intensity = noise3(p[0], p[1], p[2]);

        return intensity + 0.5;
}

float snoise(float p[3])
{
        float intensity;
        intensity = noise3(p[0], p[1], p[2]);
        return intensity;
}

float *vsnoise(float p[3])
{
        float *VSN = malloc(sizeof(float) *3);
        float x = p[0], y = p[1], z = p[2];

        VSN[0] = noise3(x, y, z);
        VSN[1] = noise3(y, z+17.7, x+23.3);
        VSN[2] = noise3(z, x+41.7, y-13.1);

        return VSN;
}


float filterednoise(float p[3], float width)
{
        return (noise(p))*(1-(smoothstep(0.2, 0.75, width)));
}

float filteredsnoise(float p[3], float width)
{
        return (snoise(p))*(1-(smoothstep(0.2, 0.75, width)));
}

float *filteredvsnoise(float p[3], float width)
{
        float *VN;

        VN = vsnoise(p);
        
        VN[0] *= (1 - (smoothstep(0.2, 0.75, width)));
        VN[1] *= (1 - (smoothstep(0.2, 0.75, width)));
        VN[2] *= (1 - (smoothstep(0.2, 0.75, width)));
        return VN;
}


float filteredabs (float x, float dx)
{
    float x0 = x - 0.5*dx;
    float x1 = x0 + dx;
    return (integral(x1) - integral(x0)) / dx;
}

/* fractional Brownian motion. */

float fBm (float p[3], float filtwidth, int octaves, float lacunarity,
                 float gain)
{
        float amp = 1;
        float fw, sum=0;
        int k, limit = octaves;
        float pp[3];

        pp[0] = p[0];
        pp[1] = p[1];
        pp[2] = p[2];
        fw = filtwidth;

        for (k = 0; k < limit; k++)
        {
                sum += amp * filteredsnoise (pp, fw);
                if (limit > 4) {printf("limit = %08i\n", limit); break;}
                amp *= gain;
                pp[0] = pp[0]*lacunarity;
                pp[1] = pp[1]*lacunarity;
                pp[2] = pp[2]*lacunarity;
                if (limit > 4) {printf("after limit = %08i\n", limit); break;}
                fw *= lacunarity;
        }

    return sum;
}


float *vfBm (float p[3], float filtwidth, float octaves, float lacunarity,
                 float gain) {
        float amp = 1;
        float pp[3];
        float *sum=malloc(sizeof(float) * 3);
        float *ttt;
        float fw, i;

        sum[0] = sum[1] = sum[2] = 0;
        pp[0] = p[0];
        pp[1] = p[1];
        pp[2] = p[2];
        fw = filtwidth;

        for (i = 0;  i < octaves;  i += 1)
        {
                ttt = filteredvsnoise (pp, fw);
                sum[0] += (amp*ttt[0]);
                sum[1] += (amp*ttt[1]);
                sum[2] += (amp*ttt[2]);
		free(ttt);
                amp *= gain;
                pp[0] = pp[0]*lacunarity;
                pp[1] = pp[1]*lacunarity;
                pp[2] = pp[2]*lacunarity;
                fw *= lacunarity;
        }
    return sum;
}

void rotateX(float theta, float *vR)
{
        float x, y, z;

        x = vR[0];
        y = vR[1]*cos(theta) - vR[2]*sin(theta);
        z = vR[2]*cos(theta) + vR[1]*sin(theta);
        vR[0] = x;
        vR[1] = y; 
        vR[2] = z;
}

void rotateY(float theta, float *vR)
{
        float x, y, z;

        y = vR[1];
        z = vR[2]*cos(theta) - vR[0]*sin(theta);
        x = vR[0]*cos(theta) + vR[2]*sin(theta);
        vR[0] = x;
        vR[1] = y;
        vR[2] = z;
}

void rotateZ(float theta, float *vR)
{
        float x, y, z;

        z = vR[2];
        x = vR[0]*cos(theta) - vR[1]*sin(theta);
        y = vR[1]*cos(theta) + vR[0]*sin(theta);

        vR[0] =x;
        vR[1] =y;
        vR[2] =z;
}

float RMturbulence (float p[3], float filtwidth,
                        float octaves, float lacunarity, float gain)
{
        float amp = 1;
        float pp[3];
        float sum = 0;
        float fw =filtwidth;
        float i;
        pp[0] = p[0];
        pp[1] = p[1];
        pp[2] = p[2];

        for (i = 0;  i < octaves;  i += 1)
        {
                float n = filteredsnoise (pp, fw);
                sum += amp * filteredabs (n, fw);
                amp *= gain;
                pp[0] = pp[0]*lacunarity;
                pp[1] = pp[1]*lacunarity;
                pp[2] = pp[2]*lacunarity;
                fw *= lacunarity;
    }

    return sum;
}

float cellnoise(float x)
{
    int seed;

        x = floor(x);
        x=(x+1.2)/2.4;

        srand(1937 +(int)(x*7680));
        seed= rand() % 768;
        srand(seed + (int)(x*x*seed));
        seed= rand() % 768;
        srand(seed + (int)(x/seed));
        seed= rand() % 768;
        srand(seed + (int)(x+seed));

        return (rand() % 512) / 512.0;
}

void normalize(float *v1) {
        float len = length(v1);
        v1[0] = v1[0]/len;
        v1[1] = v1[1]/len;
        v1[2] = v1[2]/len;
}

/* Definite integral of normalized pulsetrain from 0 to t */
float def_integral (float t, float nedge)
{
        return ((1-nedge)*floor(t) + MAX2(0, t-floor(t)-nedge));
}

float smoothpulse (float e0, float e1, float e2, float e3, float x)
{
        return smoothstep(e0, e1, x) - smoothstep(e2, e3, x);
}

float filteredpulsetrain (float edge, float period, float x, float dx)
{
        /* First normalize so period == 1 and our domain of interest
           is > 0 */

        float w, x0, x1, nedge;
        w = dx/period;
        x0 = x/period - w/2;
        x1 = x0 + w;
        nedge = edge/period;   /* normalize edge value */

        /* Now we want to integrate the normalized pulsetrain over [x0,x1]*/

        return (def_integral(x1, nedge) - def_integral(x0, nedge))/w;
}

float area(float v1[3], float v2[3])
{
        float vector[3];

        vector[0] = v1[1]*v2[2] - v1[2]*v2[1];
        vector[1] = v1[2]*v2[0] - v1[0]*v2[2];
        vector[2] = v1[0]*v2[1] - v1[1]*v2[0];
        return length(vector);
}


float filterwidthp(float v1[3], float v2[3])
{
        float size;

        size = sqrt(area(v1, v2));

        if (size > MINFILTWIDTH)
        {       return size;}
        else
        {       return MINFILTWIDTH;}
}