/////////////////////////////////////////////////////////////////////////////
// noises.h - various noise-based patterns
//
// Much of this is adapted from Apodaca & Gritz, "Advanced RenderMan",
// Morgan-Kaufmann 1999.
//
/////////////////////////////////////////////////////////////////////////////
// Copyright 2004 NVIDIA Corporation.  All Rights Reserved.
// 
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
// 
// * Redistributions of source code must retain the above copyright
//   notice, this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright
//   notice, this list of conditions and the following disclaimer in the
//   documentation and/or other materials provided with the distribution.
// * Neither the name of NVIDIA nor the names of its contributors
//   may be used to endorse or promote products derived from this software
//   without specific prior written permission.
// 
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// (This is the Modified BSD License)
/////////////////////////////////////////////////////////////////////////////


#ifndef NOISES_H
#define NOISES_H 1


#ifndef PATTERNS_H
#include "patterns.h"    /* Needed for filteredabs */
#endif



#ifndef FILTLOW
#  define FILTLOW 0.2
#endif
#ifndef FILTHIGH
#  define FILTHIGH 0.75
#endif

// If we know the filter size, we can crudely antialias snoise by fading
// to its average value at approximately the Nyquist limit.
float filteredsnoise (float p, float width) {
    return snoise(p) * (1-smoothstep (FILTLOW,FILTHIGH,width));
}
float filteredsnoise (point p, float width) {
    return snoise(p) * (1-smoothstep (FILTLOW,FILTHIGH,width));
}
float filteredsnoise (vector p, float width) {
    return snoise(p) * (1-smoothstep (FILTLOW,FILTHIGH,width));
}

vector filteredvsnoise (float p, float width) {
    return ((vector)snoise(p)) * (1-smoothstep (FILTLOW,FILTHIGH,width));
}
vector filteredvsnoise (point p, float width) {
    return ((vector)snoise(p)) * (1-smoothstep (FILTLOW,FILTHIGH,width));
}
vector filteredvsnoise (vector p, float width) {
    return ((vector)snoise(p)) * (1-smoothstep (FILTLOW,FILTHIGH,width));
}

float filteredpsnoise (float p, float per, float width) {
    return psnoise(p, per) * (1-smoothstep (FILTLOW,FILTHIGH,width));
}
float filteredpsnoise (point p, point per, float width) {
    return psnoise(p, per) * (1-smoothstep (FILTLOW,FILTHIGH,width));
}




// Antialiased fractional Brownian motion
// Inputs: 
//    p, filtwidth   position and approximate inter-pixel spacing
//    octaves        max # of octaves to calculate
//    lacunarity     frequency spacing between successive octaves
//    gain           scaling factor between successive octaves
float
fBm (point p, float filtwidth,
     float octaves, float lacunarity, float gain)
{
    float amp = 1;
    point pp = p;
    float sum = 0, fw = filtwidth;
    float i;

    for (i = 0;  i < octaves;  i += 1) {
	sum += amp * filteredsnoise (pp, fw);
	amp *= gain;  fw *= lacunarity;
        pp = pp * lacunarity + 1;  // multiply and offset
    }
    return sum;
}



// Simplified version for easy use (careful -- it takes derivatives)
float fBm (point p) { return fBm (p, sqrt(area(p)), 4.0, 2.0, 0.5); }




// fBm that takes a float rather than a position
float
fBm (float p, float filtwidth,
     float octaves, float lacunarity, float gain)
{
    float amp = 1;
    float pp = p;
    float sum = 0, fw = filtwidth;
    float i;

    for (i = 0;  i < octaves;  i += 1) {
	sum += amp * filteredsnoise (pp, fw);
	amp *= gain;  fw *= lacunarity;
        pp = pp * lacunarity + 1;  // multiply and offset
    }
    return sum;
}



// Simplified version for easy use (careful -- it takes derivatives)
float fBm (float p) { return fBm (p, sqrt(area(p)), 4.0, 2.0, 0.5); }




// A vector-valued antialiased fBm.
vector
vfBm (point p, float filtwidth,
      float octaves, float lacunarity, float gain)
{
    float amp = 1;
    point pp = p;
    vector sum = 0;
    float fw = filtwidth;
    float i;

    for (i = 0;  i < octaves;  i += 1) {
	sum += amp * filteredvsnoise (pp, fw);
	amp *= gain;  fw *= lacunarity;
        pp = pp * lacunarity + 1;  // multiply and offset
    }
    return sum;
}


// Simplified version for easy use (careful -- it takes derivatives)
vector vfBm (point p) { return vfBm (p, sqrt(area(p)), 4.0, 2.0, 0.5); }



// The stuff that Ken Musgrave calls "VLNoise"
#define VLNoise(Pt,scale) (snoise(((vector)snoise(Pt))*scale+Pt))
#define filteredVLNoise(Pt,fwidth,scale) \
            (filteredsnoise(filteredvsnoise(Pt,fwidth)*scale+Pt,fwidth))


float
VLfBm (point p, float filtwidth,
       float octaves, float lacunarity,
       float gain, float scale)
{
    float amp = 1;
    point pp = p;
    float sum = 0;
    float fw = filtwidth;
    float i;

    for (i = 0;  i < octaves;  i += 1) {
	sum += amp * filteredVLNoise (pp, fw, scale);
	amp *= gain;  fw *= lacunarity;
        pp = pp * lacunarity + 1;  // multiply and offset
    }
    return sum;
}


// Typical use of vfBm:
float VLfBm(point p) { return VLfBm (p, sqrt(area(p)), 4.0, 2.0, 0.5, 1.0); }



// antialiased periodic fBm, float version
// Inputs: 
//    p              position
//    period         period for position (should be integer value)
//    filtwidth      approximate inter-pixel spacing
//    octaves        max # of octaves to calculate
//    lacunarity     frequency spacing between successive octaves
//    gain           scaling factor between successive octaves
float pfBm (float p, float period, float filtwidth,
            float octaves, float lacunarity, float gain)
{
    float amp = 1;
    float pp = p;
    float pper = period;
    float sum = 0; float fw = filtwidth;
    float i;

    for (i = 0;  i < octaves;  i += 1) {
	sum += amp * filteredpsnoise (pp, pper, fw);
	amp *= gain;
        pp *= lacunarity;  pper *= lacunarity;  fw *= lacunarity;
    }
    return sum;
}

float pfBm (float p, float period) {
    return pfBm (p, period, filterwidth(p), 4, 2, 0.5);
}




// antialiased periodic fBm, point version
float pfBm (point p, point period, float filtwidth,
            float octaves, float lacunarity, float gain)
{
    float amp = 1;
    point pp = p;
    point pper = period;
    float sum = 0; float fw = filtwidth;
    float i;

    for (i = 0;  i < octaves;  i += 1) {
	sum += amp * filteredpsnoise (pp, pper, fw);
	amp *= gain;  
        pp *= lacunarity;  pper *= lacunarity;  fw *= lacunarity;
    }
    return sum;
}

float pfBm (point p, point period) {
    return pfBm (p, period, sqrt(area(p)), 4.0, 2.0, 0.5);
}




// Antialiased turbulence.  Watch out -- the abs() call introduces infinite
// frequency content, which makes our antialiasing efforts much trickier!
float turbulence (point p, float filtwidth, float octaves,
                  float lacunarity, float gain)
{
    float amp = 1;
    point pp = p;
    float sum = 0;
    float fw = filtwidth;
    float i;

    for (i = 0;  i < octaves;  i += 1) {
	float n = filteredsnoise (pp, fw);
	sum += amp * filteredabs (n, filterwidth(n));
	amp *= gain;  fw *= lacunarity;
        pp = pp * lacunarity + 1;  // multiply and offset
    }
    return sum;
}


// Typical use of turbulence:
#define turbulence_default(p)  turbulence (p, filterwidth(p), 4, 2, 0.5)



/////////////////////////////////////////////////////////////////////////////
// Voronoi cell noise (a.k.a. Worley noise) functions
//
// These functions assume that space is filled with "features" (points
// of interest).  There are interesting patterns we can make by
// figuring out which feature we are closest to, or to what extent
// we're on the boundary between two features.  Several varieties of
// these computations are below, categorized by the dimension of their
// domains, and the number of close features they are interested in.
//
// All these functions have similar inputs:
//   P      - position to test (for 3-D varieties; 2-D varieties use ss,tt)
//   jitter - how much to jitter the cell center positions (1 is typical,
//             smaller values make a more regular pattern, larger values
//             make a more jagged pattern; use jitter >1 at your risk!).
// And outputs:
//   f_n    - distance to the nth nearest feature (f1 is closest, f2 is
//            the distance to the 2nd closest, etc.)
//   pos_n  - the position of the nth nearest feature.  For 2-D varieties,
//            these are instead spos_n and tpos_n.
/////////////////////////////////////////////////////////////////////////////

// Voronoi cell noise (a.k.a. Worley noise) -- 3-D, 1-feature version.
void
voronoi (point P, float jitter,
         output float f1, output point pos1)
{
    point thiscell = point (floor(P[0]), floor(P[1]), floor(P[2])) + 0.5;
    f1 = 1000;
    float i, j, k;
    for (i = -1;  i <= 1;  i += 1) {
        for (j = -1;  j <= 1;  j += 1) {
            for (k = -1;  k <= 1;  k += 1) {
		point testcell = thiscell + vector(i,j,k);
                point pos = testcell + 
		            jitter * ((vector)cellnoise (testcell) - 0.5);
		vector offset = pos - P;
                float dist = dot(offset,offset); // actually dist^2
                if (dist < f1) {
                    f1 = dist;  pos1 = pos;
                }
            }
	}
    }
    f1 = sqrt(f1);
}


/* Voronoi cell noise (a.k.a. Worley noise) -- 3-D, 2-feature version. */
void
voronoi (point P, float jitter,
         output float f1, output point pos1,
         output float f2, output point pos2)
{
    point thiscell = point (floor(P[0]), floor(P[1]), floor(P[2])) + 0.5;
    f1 = f2 = 1000;
    float i, j, k;
    for (i = -1;  i <= 1;  i += 1) {
        for (j = -1;  j <= 1;  j += 1) {
            for (k = -1;  k <= 1;  k += 1) {
		point testcell = thiscell + vector(i,j,k);
                point pos = testcell + 
		            jitter * ((vector)cellnoise(testcell) - 0.5);
		vector offset = pos - P;
                float dist = dot(offset,offset); // actually dist^2
                if (dist < f1) {
                    f2 = f1;  pos2 = pos1;
                    f1 = dist;  pos1 = pos;
                } else if (dist < f2) {
                    f2 = dist;  pos2 = pos;
		}
            }
	}
    }
    f1 = sqrt(f1);  f2 = sqrt(f2);
}


// Voronoi cell noise (a.k.a. Worley noise) -- 2-D, 1-feature version.
void
voronoi (float ss, float tt, float jitter,
         output float f1, output float spos1, output float tpos1)
{
    float sthiscell = floor(ss)+0.5, tthiscell = floor(tt)+0.5;
    f1 = 1000;
    float i, j;
    for (i = -1;  i <= 1;  i += 1) {
	float stestcell = sthiscell + i;
        for (j = -1;  j <= 1;  j += 1) {
	    float ttestcell = tthiscell + j;
	    float spos = stestcell +
		     jitter * (cellnoise(stestcell, ttestcell) - 0.5);
	    float tpos = ttestcell +
		 jitter * (cellnoise(stestcell+23, ttestcell-87) - 0.5);
	    float soffset = spos - ss;
	    float toffset = tpos - tt;
	    float dist = soffset*soffset + toffset*toffset;
	    if (dist < f1) {
		f1 = dist;
		spos1 = spos;  tpos1 = tpos;
	    }
	}
    }
    f1 = sqrt(f1);
}


// Voronoi cell noise (a.k.a. Worley noise) -- 2-D, 2-feature version.
void
voronoi (float ss, float tt, float jitter,
         output float f1, output float spos1, output float tpos1,
         output float f2, output float spos2, output float tpos2)
{
    float sthiscell = floor(ss)+0.5, tthiscell = floor(tt)+0.5;
    f1 = f2 = 1000;
    float i, j;
    for (i = -1;  i <= 1;  i += 1) {
	float stestcell = sthiscell + i;
        for (j = -1;  j <= 1;  j += 1) {
	    float ttestcell = tthiscell + j;
	    float spos = stestcell +
		   jitter * (cellnoise(stestcell, ttestcell) - 0.5);
	    float tpos = ttestcell +
		   jitter * (cellnoise(stestcell+23, ttestcell-87) - 0.5);
	    float soffset = spos - ss;
	    float toffset = tpos - tt;
	    float dist = soffset*soffset + toffset*toffset;
	    if (dist < f1) {
		f2 = f1;  spos2 = spos1;  tpos2 = tpos1;
		f1 = dist;  spos1 = spos;  tpos1 = tpos;
	    } else if (dist < f2) {
		f2 = dist;
		spos2 = spos;  tpos2 = tpos;
	    }
	}
    }
    f1 = sqrt(f1);  f2 = sqrt(f2);
}


#endif /* NOISES_H */
