#include <stdio.h>
#include <stdlib.h>
#include <math.h>
  
#include "shadeop.h"
  
// File initialization
FILE     *f = NULL;
  
//------------------------------------------------------------------------------
//                        Post process of shadeops
//------------------------------------------------------------------------------
  
SHADEOP_CLEANUP(hatchingCurves_cleanup) {
if(f != NULL) {
    fprintf(f, "\nTransformEnd\n");
    fclose(f);
    }
f = NULL;
} 
  
//------------------------------------------------------------------------------
//                        Table for actual functions
//------------------------------------------------------------------------------
  
SHADEOP_TABLE(hatchingCurves) =
{
    { "void hatchingCurves_simple (float, point, normal, vector, float, float)", "", "hatchingCurves_cleanup" }, 
    { "void hatchingCurves_adv (float, point, normal, vector, float, float, float, float, float, float, float, string)", "", "hatchingCurves_cleanup" }, 
    {""}
};
  
//------------------------------------------------------------------------------
//                        Functions used in Shadeops
//------------------------------------------------------------------------------
  
// (Not used)
float dot_product(float *a,float *b,int size) //size: size of the array.
{
    float dp = 0.0f;
    int i;
    for (i=0;i<size;i++)
        dp += a[i] * b[i];
    return dp;
}
  
//------------------------------------------------------------------------------
  
void cross_product(float *a, float *b, float *cp)
{
    cp[0] = a[1]*b[2] - a[2]*b[1];
    cp[1] = a[2]*b[0] - a[0]*b[2];
    cp[2] = a[0]*b[1] - a[1]*b[0];
}
  
//------------------------------------------------------------------------------
  
float frand(void)
{
    return (float)rand()/(float)RAND_MAX;
}
  
//------------------------------------------------------------------------------
  
float rand_length(float min, float max)
{
    float f_rand = 0.0;
    float result = 0.0;
    f_rand = frand();
    result = f_rand*(max - min) + min;
    result = result > 0 ? result : -result;
    return result;
}
  
//------------------------------------------------------------------------------
  
void rotateAlongNormal(float *pnt, float *center, float *n, float rotation)
{
    float cp[3] = {0,0,0};
    float pnt_temp[3] = {pnt[0], pnt[1], pnt[2]};
    cross_product(n, pnt_temp, cp);
    pnt[0] += sin(rotation)*cp[0];
    pnt[1] += sin(rotation)*cp[1];
    pnt[2] += sin(rotation)*cp[2];
    
    pnt[0] += (cos(rotation)-1)*pnt_temp[0];
    pnt[1] += (cos(rotation)-1)*pnt_temp[1];
    pnt[2] += (cos(rotation)-1)*pnt_temp[2];    
}
  
//------------------------------------------------------------------------------
  
// (Not used)
float pointsDist(float *p1, float *p2)
{
    float result = sqrt((p1[0]-p2[0])*(p1[0]-p2[0]) + (p1[1]-p2[1])*(p1[1]-p2[1]) + (p1[2]-p2[2])*(p1[2]-p2[2]));
    return result;
}
  
//------------------------------------------------------------------------------
  
void jitPoint(float *p, float jit)
{
    
    p[0] = p[0] + (frand() * 2.0 - 1.0) * jit;
    p[1] = p[1] + (frand() * 2.0 - 1.0) * jit;
    p[2] = p[2] + (frand() * 2.0 - 1.0) * jit;
}
  
//------------------------------------------------------------------------------
//                                Shadeops
//------------------------------------------------------------------------------
  
SHADEOP (hatchingCurves_simple)
{
  
float rate = *(float*)argv[1];
float *pnt = (float*)argv[2];
float *nf = (float*)argv[3];
float *i = (float*)argv[4];
float minLength = *(float*)argv[5];
float maxLength = *(float*)argv[6];
  
float cp[3] = {0,0,0};
float length = 0.0;
float p1[3] = {0,0,0};
float p2[3] = {0,0,0};
  
  
// Do we have access to an open file?
if(f == NULL) 
    {
    f = fopen("./hatching_curves.rib", "w");
    fprintf(f, "\nTransformBegin\n");
    }
// Write the xyz position of a point
  
if(frand() < rate)
{
  
    cross_product(nf, i, cp);
    length = rand_length(minLength, maxLength);
  
    cp[0] *= length/2.0;
    cp[1] *= length/2.0;
    cp[2] *= length/2.0;
  
    p1[0] = pnt[0] + cp[0];
    p1[1] = pnt[1] + cp[1];
    p1[2] = pnt[2] + cp[2];
    
    p2[0] = pnt[0] - cp[0];
    p2[1] = pnt[1] - cp[1];
    p2[2] = pnt[2] - cp[2];
    
    fprintf(f, "Curves \"linear\" [2] \"nonperiodic\" \"P\" [");
    fprintf(f, "%f %f %f %f %f %f", p1[0], p1[1], p1[2], p2[0], p2[1], p2[2]);
    fprintf(f, "] \"constantwidth\" [0.005]\n");
}
  
return 0;
}
  
//------------------------------------------------------------------------------
  
SHADEOP (hatchingCurves_adv)
{
  
float rate = *(float*)argv[1];
float *pnt = (float*)argv[2];
float *nf = (float*)argv[3];
float *i = (float*)argv[4];
float minLength = *(float*)argv[5];
float maxLength = *(float*)argv[6];
  
float rand = *(float*)argv[7]; // randomness of rotation
float jit = *(float*)argv[8];
float complexity = *(float*)argv[9];
float density = *(float*)argv[10];
float width = *(float*)argv[11];
STRING_DESC    *filepath = (STRING_DESC *)argv[12];
  
float temp = 0.0;
float cp[3] = {0,0,0};
float length = 0.0;
float pStart[3] = {0,0,0};
float pEnd[3] = {0,0,0};
float pStartOrg[3] = {0,0,0};
float pEndOrg[3] = {0,0,0};
int intComplexity = (int)complexity;
float unitLength = 0.0;
int cvNum = 0;
int k = 0;
float unitVec[3] = {0,0,0};
float eachCv[3] = {0,0,0};
float rotation = 0.0;
  
float PI = 3.14159265;
  
  
// Makes sure maxLength is higher than minLength
if(minLength > maxLength)
{
    temp = minLength;
    minLength = maxLength;
    maxLength = temp;
}
  
// Do we have access to an open file?
if(f == NULL) 
    {
    f = fopen(filepath->s, "w");
    fprintf(f, "TransformBegin\n");
    }
    
// Decides to create a curve or not
if(frand() < rate*density)
{
  
    cross_product(nf, i, cp);
    length = rand_length(minLength, maxLength);
    cp[0] *= length/2.0;
    cp[1] *= length/2.0;
    cp[2] *= length/2.0;
  
    pStart[0] = pnt[0] + cp[0];
    pStart[1] = pnt[1] + cp[1];
    pStart[2] = pnt[2] + cp[2];
    
    pEnd[0] = pnt[0] - cp[0];
    pEnd[1] = pnt[1] - cp[1];
    pEnd[2] = pnt[2] - cp[2];
    
// How many units can the curve have?
    if(intComplexity <= 0)
    {
        unitLength = maxLength;
        cvNum = 0;
    }
    else
    {
        unitLength = maxLength / (intComplexity + 1);
        cvNum = (int)(length/unitLength);
    }
  
// Randomizes rotation for start and end points
  
    //--------------- Now not working crrectly -------------------------
    rotation = (frand() * 2 -1) * rand * PI;
    rotateAlongNormal(pStart, pnt, nf, rotation);
    rotateAlongNormal(pEnd, pnt, nf, rotation);
  
// Store the start and end point positions before jittering    
    pStartOrg[0] = pStart[0];
    pStartOrg[1] = pStart[1];
    pStartOrg[2] = pStart[2];
    
    pEndOrg[0] = pEnd[0];
    pEndOrg[1] = pEnd[1];
    pEndOrg[2] = pEnd[2];
    
// jiting for start and end points
    jitPoint(pStart, jit);
    jitPoint(pEnd, jit);
    
// Writes the begining of curves statement
    if(cvNum <= 2)
        fprintf(f, "Curves \"linear\" [%d] \"nonperiodic\" \"P\" [%f %f %f ", (2 + cvNum), pStart[0], pStart[1], pStart[2]);
    else
        fprintf(f, "Curves \"linear\" [%d] \"nonperiodic\" \"P\" [%f %f %f ", (2 + cvNum), pStart[0], pStart[1], pStart[2]);
    
// Writes the intermediate CVs
    unitVec[0] = (pStartOrg[0] - pEndOrg[0])/(float)(cvNum + 1);
    unitVec[1] = (pStartOrg[1] - pEndOrg[1])/(float)(cvNum + 1);
    unitVec[2] = (pStartOrg[2] - pEndOrg[2])/(float)(cvNum + 1);
    
    eachCv[0] = pStartOrg[0];
    eachCv[1] = pStartOrg[1];
    eachCv[2] = pStartOrg[2];
    for(k = 0; k < cvNum; k++)
    {
        float tempCv[3] = {0,0,0};
        
        eachCv[0] += unitVec[0];
        eachCv[1] += unitVec[1];
        eachCv[2] += unitVec[2];
        
        tempCv[0] = eachCv[0];
        tempCv[1] = eachCv[1];
        tempCv[2] = eachCv[2];
        
        jitPoint(tempCv, jit);
        fprintf(f, "%f %f %f ", tempCv[0], tempCv[1], tempCv[2]);
    }
  
// Writes the end of curve statement 
    fprintf(f, "%f %f %f] \"constantwidth\" [%f]\n", pEnd[0], pEnd[1], pEnd[2], width);
}
  
return 0;
}