/////////////////////////////////////////////////////////////////////////////
// patterns.h - lots of handy pattern functions for building complex
//              shaders.
//
// 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 PATTERNS_H
#define PATTERNS_H 1

// 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 filteredabs (float x, float dx)
{
    float integral (float t) {
	return sign(t) * 0.5 * t*t;
    }

    float x0 = x - 0.5*dx;
    float x1 = x0 + dx;
    float result;
    // If the interval doesn't cross 0, then the box filter is just abs(x).
    // If dx is small (vs x), the integral is numerically unstable;
    // in these cases the interval doesn't cross 0.
    // If dx is 0, the integral is NaN.
    if (x0 * x1 >= 0 || dx == 0.0)
        result = abs(x);
    else
        result = (integral(x1) - integral(x0)) / dx;
    return result;
}



// Antialiased smoothstep(e0,e1,x).  
// Compute the box filter of smoothstep(e0,e1,t) from x-dx/2 to x+dx/2.
// Strategy: divide domain into 3 regions: t < e0, e0 <= t <= e1,
// and t > e1.  Region 1 has integral 0.  Region 2 is computed by
// analytically integrating smoothstep, which is -2t^3+3t^2.  Region 3
// is trivially 1.
//
float filteredsmoothstep (float e0, float e1, float x, float dx)
{
    float integral (float t) {
	return t*t*t * (1 - 0.5*t);
    }

    // Compute x0, x1 bounding region of integration, and normalize so
    // that e0==0, e1==1
    float edgediff = e1 - e0;
    float x0 = (x-e0)/edgediff;
    float fw = dx / edgediff;
    x0 -= 0.5*fw;
    float x1 = x0 + fw;

    // Region 1 always contributes nothing
    float integ = 0;
    // Region 2 - compute integral in region between 0 and 1
    if (x0 < 1 && x1 > 0)
	integ += integral(min(x1,1)) - integral(max(x0,0));
    // Region 3 - is 1.0
    if (x1 > 1)
	integ += x1-max(1,x0);
    return integ / fw;
}



// A 1-D pulse pattern:  return 1 if edge0 <= x <= edge1, otherwise 0
float pulse (float edge0, float edge1, float x)
{
    return step(edge0,x) - step(edge1,x);
}



float filteredpulse (float edge0, float edge1, float x, float dx)
{
    float x0 = x - dx/2;
    float x1 = x0 + dx;
    return max (0, (min(x1,edge1)-max(x0,edge0)) / dx);
}



// A pulse train: a signal that repeats with a given period, and is
// 0 when 0 <= mod(x,period) < edge, and 1 when mod(x,period) > edge.
//
float pulsetrain (float edge, float period, float x)
{
    return pulse (edge, period, mod(x,period));
}



// Filtered pulse train: it's not as simple as just returning the mod
// of filteredpulse -- you have to take into account that the filter may
// cover multiple pulses in the train.
// Strategy: consider the function that is the integral of the pulse
// train from 0 to x. Just subtract!
//
float filteredpulsetrain (float edge, float period, float x, float dx)
{
    // First, normalize so period == 1 and our domain of interest is > 0
    float w = dx/period;
    float x0 = x/period - w/2;
    float x1 = x0+w;
    float nedge = edge / period;   // normalized edge value

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

    // Now we want to integrate the normalized pulsetrain over [x0,x1]
    return (integral(x1) - integral(x0)) / w;
}



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


float filteredsmoothpulse (float e0, float e1, float e2, float e3,
                           float x, float dx)
{
    return filteredsmoothstep(e0,e1,x,dx) - filteredsmoothstep(e2,e3,x,dx);
}



// A pulse train of smoothsteps: a signal that repeats with a given
// period, and is 0 when 0 <= mod(x/period,1) < edge, and 1 when
// mod(x/period,1) > edge.  
//
float smoothpulsetrain (float e0, float e1, float e2, float e3,
                        float period, float x)
{
    return smoothpulse (e0, e1, e2, e3, mod(x,period));
}



// (Box-)Filtered sine wave.
// Strategy: integral of sin(x)dx is cos(x)+C
// FIXME Do better than box filter.  For large dx, this should return 0.
//
float filteredsin (float x, float dx)
{
    // For small dx, integral is numerically unstable.
    if (abs(dx) < 0.01 /*arbitrary*/)
        return sin(x);
    float x0 = x-dx/2;
    float x1 = x0+dx;
    return (cos(x0) - cos(x1))/dx;
}




// varyEach takes a computed color, then tweaks each indexed item
// separately to add some variation.  Hue, saturation, and lightness
// are all independently controlled.  Hue adds, but saturation and
// lightness multiply.
//
color varyEach (color Cin, float index,
                float varyhue, float varysat, float varylum)
{
    // Convert to "hsl" space, it's more convenient
    color Chsl = transformc ("hsl", Cin);
    float h = Chsl[0], s = Chsl[1], l = Chsl[2];
    // Modify Chsl by adding Cvary scaled by our separate h,s,l controls
    h += varyhue * (cellnoise(index+3)-0.5);
    s *= 1 - varysat * (cellnoise(index-14)-0.5);
    l *= 1 - varylum * (cellnoise(index+37)-0.5);
    Chsl = color (mod(h,1), clamp(s,0,1), clamp(l,0,1));
    // Clamp hsl and transform back to rgb space
    return transformc ("hsl", "rgb", clamp(Chsl, (color)0, (color)1));
}




// Given 2-D texture coordinates s,t and their filter widths ds, dt,
// and the width and height of the grooves between tiles (assuming that
// tile spacing is 1.0), figure out which (integer indexed) tile we are
// on and what coordinates (on [0,1]) within our individual tile we are
// shading.
//
float tilepattern (float s, float t, float ds, float dt,
                   float groovewidth, float grooveheight,
                   output float swhichtile, output float twhichtile,
                   output float stile, output float ttile)
{
    swhichtile = floor (s);
    twhichtile = floor (t);
    stile = s - swhichtile;
    ttile = t - twhichtile;
    return filteredpulsetrain (groovewidth, 1, s+groovewidth/2, ds)
             * filteredpulsetrain (grooveheight, 1, t+grooveheight/2, dt);
}



// Given 2-D texture coordinates s,t and their filter widths ds, dt,
// figure out which (integer indexed) tile we are on and what
// coordinates (on [0,1]) within our individual tile we are shading.
// Return 0 or 1 (and smoothly in between at edges) depending on which
// tile in an alternating checker pattern we're on.  Tiles are assumed
// to be 1x1 units large (so the checker pattern repeats every 2 units).
//
float checkerpattern (float s, float t, float ds, float dt,
                      output float swhichtile, output float twhichtile,
                      output float stile, output float ttile)
{
    swhichtile = floor (s);
    twhichtile = floor (t);
    stile = s - swhichtile;
    ttile = t - twhichtile;
    float x = filteredpulsetrain (1, 2, s, ds);
    float y = filteredpulsetrain (1, 2, t, dt);
    return x*y + (1-x)*(1-y);
}




#if 0
// wrapcubetexturef -
// One way to wrap a texture around an object that has no clean st
// space is to project it from all 6 sides of a cube, smoothly
// blending at the seams.  Given a point p and normal n, this function
// returns such a blended lookup from a 1-channel texture.
float
wrapcubetexturef (string texturename, point p, normal n, float blur)
{
    float sum = 0, sumw = 0;

    float w = abs (n[0]);
    w *= smoothstep (.5, .65, w);
    if (gridany(w)) {
        float s = 0.5*p[2]+0.5;
        float t = 0.5*p[1]+0.5;
        sum += w * texture (texturename, s, t,
                            "wrap", "periodic", "blur", blur);
        sumw += w;
    }

    w = abs (n[1]);
    w *= smoothstep (.5, .65, w);
    if (gridany(w)) {
        float s = 0.5*p[0]+0.5;
        float t = 0.5*p[2]+0.5;
        sum += w * texture (texturename, s, t,
                            "wrap", "periodic", "blur", blur);
        sumw += w;
    }

    w = abs (n[2]);
    w *= smoothstep (.5, .65, w);
    if (gridany(w)) {
        float s = 0.5*p[0]+0.5;
        float t = 0.5*p[1]+0.5;
        sum += w * texture (texturename, s, t,
                            "wrap", "periodic", "blur", blur);
        sumw += w;
    }

    if (sumw != 0)
	sum /= sumw;

    return sum;
}


// If blur is not specified, use 0 blur
float
wrapcubetexturef (string texturename, point p, normal n)
{
    return wrapcubetexturef (texturename, p, n, 0);
}
#endif


// basic brick tiling pattern --
//   inputs:
//      s, t                    positions on a 2-D surface
//      ds, dt                  filter widths of s and t
//      tilewidth, tileheight   dimensions of each tile
//      rowstagger              how much does each row stagger relative to
//                                   the previous row
//      rowstaggervary          how much should rowstagger randomly vary
//      jaggedfreq, jaggedamp   adds noise to the edge between the tiles
//   outputs:
//      swhichtile, twhichtile  index which tile the sample is in
//      stile, ttile            position within this tile (0-1)
//
void brickpattern (float s, float t, float ds, float dt,
                   float tilewidth, float tileheight,
                   float rowstagger, float rowstaggervary,
                   float jaggedfreq, float jaggedamp,
                   output float swhichtile, output float twhichtile,
                   output float stile, output float ttile)
{
    float scoord = s, tcoord = t;

    if (jaggedamp != 0.0) {
	// Make the shapes of the bricks vary just a bit
	point PP = noise (s*jaggedfreq/tilewidth, t*jaggedfreq/tileheight);
	scoord += jaggedamp * PP[0];
	tcoord += jaggedamp * PP[1];
    }

    stile = scoord / tilewidth;
    ttile = tcoord / tileheight;
    twhichtile = floor (ttile);   // which brick row?

    // Shift the columns randomly by row
    stile += mod (rowstagger * twhichtile, 1);
    stile += rowstaggervary * (noise (twhichtile+0.5) - 0.5);

    swhichtile = floor (stile);
    stile -= swhichtile;
    ttile -= twhichtile;
}

#endif /* PATTERNS_H */
