//This is the source fragment for the project CMPBRDF, published at EGSR 2016:
// http://dcgi.felk.cvut.cz/home/havran/CMPBRDF/
//
// @Article{havran16perceptually,
//  Title                    = {{Perceptually Motivated {BRDF} Comparison using Single Image}},
//  Author                   = {Havran, Vlastimil and Filip, Jiri and Myszkowski, Karol},
//  Journal                  = {Computer Graphics Forum},
//  Year                     = {2016},
//  Doi                      = {10.1111/cgf.12944},
//  ISSN                     = {1467-8659},
//  Publisher                = {The Eurographics Association and John Wiley \& Sons Ltd.}
//}

  // Parameters of a currently generated surfaces, the values in the paper appendix:
  float scale;
  float nl; // waves per circle
  float nc; // waves per radius
  float R; // radius of sphere

  // cropping between x to range <0, 1>
  float crop(float x) const { return (x<0.f) ? 0.f : (x > 1.0f ? 1.0f : x);}

  // ramp function of x between e0 and e1, output is 0 to 1
  float cx(float x, float e0, float e1) const { return crop((x-e0)/(e1-e0));}

  // smooth transition function between e0 and e1 on x-axis, output is 0 to 1
  float tx(float x, float e0, float e1) const {
    float y = cx(x,e0,e1);
    return y*y*y*(y*(y*6.0f - 15.0) + 10.0f);
  }

  // 3D function defining the surfaces used in the paper. Using polar coordinates
  // R,phi and generating Z-coordinate. X and Y coordinate are computed
  // from polar ones as: X= R*cos(phi), Y=R*sin(phi)
  virtual float surfaceFunctionPolar(float r, float phi) {
    if (surfaceType == 2) { // surface 2
      float rads = 0.09f; // sphere radius - fixed size
      float rad2 = 1.0f * rads; // transition starts
      float rad3 = 1.2f * rads; // transition ends
      // kind of spherical surface - z-coordinate
      float rr = rads*0.4f + ( (r < rads) ? sqrt(rads*rads-r*r) : ( (r < 2.0f*rads) ? -sqrt(rads*rads - sqr(2.0f*rads-r)) : -rads ) );
      float r2 = (r < rad2) ? 0.f : ((r-rad2)*nl+0.25f)*M_PI; // radial waves with phase offset
      float Z = scale*tx(r, rad2, rad3)*sin(nc*2.0*phi)*cos(r2) +
        rr*(1.0-tx(r, rad2, rad3)) + sqrt(R*R-r*r)-R;
      return Z;
    }
    if (surfaceType == 3) { // surface 3
      float rads = 0.09f; // sphere radius - fixed size
      float rr = rads*0.25f + ( (r < rads) ? sqrt(rads*rads-r*r) : ( (r < 2.0f*rads) ? -sqrt(rads*rads - sqr(2.0f*rads-r)) : -rads ));
      if (r < rads) r = 0;     
      float Z = scale*10.0f*r*r*sin(nc*phi + sin(nl*M_PI*r))*sqrt(R*R-r*r) + rr;
      return Z;
    }
  }