/*
* Copyright (c) 1998 Thomas E. Burge. All rights reserved.
*
* FILE: showk.sl
*
* DESCRIPTION:
* Calculates different types of curvature. Right now this only works
* on PRMan, but I'll publish it here in case someone finds it of
* interest. It needs more testing, but the numbers do seem to match
* Alias's numbers.
*
* TO-DO:
* The way a curvature value is matched to a color needs to be improved
* and I need to double check the expressions with the terms Puu.NN,
* Puv.NN and Pvv.NN.
*
* Do note that the second derivatives seem to cause some renderers
* problems as the values do not appear too valid at times.
*/
/*#define DEBUG*/
/* Note that the value 1/Calck(dPdu, Du(dPdu)) is the radius
* of curvature in the u direction of a surface. So
* a cylinder (with an unscaled cirular cross section)
* would have a radius at P equal to the radius of
* curvature at P.
* Q1 is a first derivative and Q2 is the matching second
* derivative Q1'.
*/
float Calck( vector Q1; vector Q2 )
{
return length((Q1^Q2)/pow(length(Q1),3));
}
/* Solves the quadratic equation and finds two roots. */
void CalcQuadEquation( float a; float b; float c;
output float r1; output float r2 )
{
float q;
q = b*b - 4*a*c;
if ( q < 0 )
{
r1 = 0;
r2 = 0;
}
else
{
q = -0.5*( b + sign(b) * sqrt(q));
r1 = q / a;
r2 = c / q;
}
}
/* Solves a quadratic equation for principle curvatures k1 and k2. */
#define CALC_K1K2() \
CalcQuadEquation( ((dPdu.dPdu)*(dPdv.dPdv) - pow(dPdu.dPdv, 2)), \
-((dPdu.dPdu)*(Pvv.NN)+(dPdv.dPdv)*(Puu.NN) \
- 2*(dPdu.dPdv)*(Puv.NN)), \
((Puu.NN)*(Pvv.NN) - pow(Puv.NN, 2)), \
k1,k2)
/* The way the color ramp is handled is not very good. You need to
* find the min and max values of the type of curvature you are
* viewing, otherwise one color gets used for the entire surface.
*/
color CalcColor( color c1; color c2; color c3;
float min; float max;
float index )
{
float i;
color c;
i = clamp( index, min, max );
i = (i-min)/(max-min);
c = spline( i, c1, 0.8*c1, c2, 0.8*c3, c3 );
return c;
}
surface showk( float Kd=0.5;
float Ka=0.1;
float roughness=.1;
/* gaussian, mean, min, max, ku, or kv */
string CurvatureType = "gaussian";
color c1 = color(0,1,0);
color c2 = color(1,0,0);
color c3 = color(0,0,1);
float min = -10.0;
float max = 10.0;
)
{
normal NN;
normal Nf;
vector V;
vector Puu, Puv, Pvv;
float k, k1, k2;
color Ct;
float G, M;
NN = normalize(N);
Nf = faceforward( NN, I );
V = -normalize(I);
if ( CurvatureType == "gaussian")
{
Puu = Du(dPdu);
Puv = Du(dPdv);
Pvv = Dv(dPdv);
k = ((Puu.NN)*(Pvv.NN) - pow(Puv.NN, 2))
/ ((dPdu.dPdu)*(dPdv.dPdv) - pow(dPdu.dPdv, 2));
#ifdef DEBUG
/* The Gaussian curvature is the product of the principle
* curvatures at point P. Note that the Gaussian curvature
* can be calculated as above or from the principle
* curvatures referred to below as k1 and k2. The values
* can be calculated by solving a quadratic equation or
* a simpler equation when dPdu.dPdv==Puv.NN==0.
*/
CALC_K1K2();
if ( abs(dPdu.dPdv)<0.000001 && abs(Puv.NN)<0.000001 )
printf(
"dPdu.dPdv:%g Puv.NN:%g k:%g=k1*k2=%g (k1=%g=%g)(k2=%g=%g)\n",
dPdu.dPdv, Puv.NN, k, k1*k2,
k1,(Puu.NN)/(dPdu.dPdu),
k2,(Pvv.NN)/(dPdv.dPdv) );
else
printf(" dPdu.dPdv:%g Puv.NN:%g k:%g = k1*k2 = %g\n",
dPdu.dPdv, Puv.NN, k, k1*k2 );
#endif
Ct = CalcColor(c1,c2,c3, min,max, k);
}
else if ( CurvatureType == "mean")
{
Puu = Du(dPdu);
Puv = Du(dPdv);
Pvv = Dv(dPdv);
k = ((dPdu.dPdu)*(Pvv.NN)+(dPdv.dPdv)*(Puu.NN)
- 2*(dPdu.dPdv)*(Puv.NN))
/ (2*((dPdu.dPdu)*(dPdv.dPdv) - pow(dPdu.dPdv, 2)));
#ifdef DEBUG
/* The Mean curvature is the average of the principle
* curvatures at point P. Note that the Mean curvature
* can be calculated as above or from the principle
* curvatures referred to below as k1 and k2. The values
* can be calculated by solving a quadratic equation or
* a simpler equation when dPdu.dPdv==Puv.NN==0.
*/
CALC_K1K2();
if ( abs(dPdu.dPdv)<0.000001 && abs(Puv.NN)<0.000001 )
printf(
"dPdu.dPdv:%g Puv.NN:%g k:%g=(k1+k2)/2:%g (k1=%g=%g)(k2=%g=%g)\n",
dPdu.dPdv, Puv.NN, k, (k1+k2)/2,
k1,(Puu.NN)/(dPdu.dPdu),
k2,(Pvv.NN)/(dPdv.dPdv) );
else
printf(" dPdu.dPdv:%g Puv.NN:%g k:%g = (k1+k2)/2 = %g\n",
dPdu.dPdv, Puv.NN, k, (k1+k2)/2 );
#endif
Ct = CalcColor(c1,c2,c3, min,max, k);
}
else if ( CurvatureType == "ku")
{
k = Calck(dPdu, Du(dPdu));
Ct = CalcColor(c1,c2,c3, min,max, k);
}
else if ( CurvatureType == "kv")
{
k = Calck(dPdv, Du(dPdv));
Ct = CalcColor(c1,c2,c3, min,max, k);
}
else if ( CurvatureType == "min")
{
Puu = Du(dPdu);
Puv = Du(dPdv);
Pvv = Dv(dPdv);
CALC_K1K2();
k = min(k1,k2);
#ifdef DEBUG
printf( "min k:%g k1:%g k2:%g\n", k, k1, k2 );
#endif
Ct = CalcColor(c1,c2,c3, min,max, k);
}
else if ( CurvatureType == "max")
{
Puu = Du(dPdu);
Puv = Du(dPdv);
Pvv = Dv(dPdv);
CALC_K1K2();
k = max(k1,k2);
#ifdef DEBUG
printf( "max k:%g\n", k );
#endif
Ct = CalcColor(c1,c2,c3, min,max, k);
}
else /* information printed out */
{
/* Print surface information for testing by adding the following
* line to a RIB file such as the one attached below:
*
* Surface "showk" "CurvatureType" ""
*
* Using the above Surface statement in the RIB file below (bumpy.rib),
* the following file creates file listing Gaussian, Mean and
* principle curvatures:
*
* prman bumpy.rib > junk.txt
*
*/
Puu = Du(dPdu);
Puv = Du(dPdv);
Pvv = Dv(dPdv);
G = ((Puu.NN)*(Pvv.NN) - pow(Puv.NN, 2))
/ ((dPdu.dPdu)*(dPdv.dPdv) - pow(dPdu.dPdv, 2));
M = ((dPdu.dPdu)*(Pvv.NN)+(dPdv.dPdv)*(Puu.NN)
- 2*(dPdu.dPdv)*(Puv.NN))
/ (2*((dPdu.dPdu)*(dPdv.dPdv) - pow(dPdu.dPdv, 2)));
CALC_K1K2();
printf("u,v(%g,%g)G:%g M:%g k1:%g k2:%g\n", u,v, G, M, k1, k2);
Ct = 1.0;
}
Oi = Os;
Ci = Os * ( Cs * Ct * (Ka*ambient() + Kd*diffuse(Nf)) );
}
#ifdef SAMPLE_RIB
Display "bumpy.rib.tiff" "framebuffer" "rgba"
Format 300 300 1
Projection "perspective" "fov" [30]
Rotate -53.0903 0 0 1
Rotate 18.2711 1 0 0
Rotate 22.6563 0 1 0
Translate -1.13392 0.97189 2.71655
WorldBegin
LightSource "distantlight" 1 "intensity" [1]
LightSource "ambientlight" 2 "intensity" [.1]
Declare "CurvatureType" "uniform string"
Declare "min" "uniform float"
Declare "max" "uniform float"
Surface "showk" "CurvatureType" "gaussian" "min" -600 "max" 90
Surface "showk" "CurvatureType" "mean" "min" -10 "max" 10
Surface "showk" "CurvatureType" "mean" "min" -10 "max" 10
Surface "showk" "CurvatureType" "max" "min" -23 "max" 23
NuPatch 7 4 [0 0 0 0 0.25 0.5 0.751669 1 1 1 1] 0 1 7 4 [-1 -1 -1 -1 -0.746696 -0.5 -0.253304 0 0 0 0] -1 0 "Pw" [-0.5 -0.5 0 1
-0.5 -0.416667 0 1 -0.5 -0.25 0 1 -0.5 0.00055639 0 1
-0.5 0.250556 0 1 -0.5 0.417223 0 1 -0.5 0.5 0 1
-0.415565 -0.5 0 1 -0.415565 -0.416667 0.16867 1 -0.415565 -0.25 0.227614 1
-0.415565 0.00055639 -0.0832115 1 -0.415565 0.250556 -0.355132 1 -0.415565 0.417223 -0.217987 1
-0.415565 0.5 0 1 -0.248899 -0.5 0 1 -0.248899 -0.416667 0.223213 1
-0.248899 -0.25 0.272004 1 -0.248899 0.00055639 -0.20521 1 -0.248899 0.250556 -0.586827 1
-0.248899 0.417223 -0.346515 1 -0.248899 0.5 0 1 0 -0.5 0 1
0 -0.416667 -0.0812803 1 0 -0.25 -0.203201 1 0 0.00055639 -0.264297 1
0 0.250556 -0.202929 1 0 0.417223 -0.0807376 1 0 0.5 0 1
0.248899 -0.5 0 1 0.248899 -0.416667 -0.348844 1 0.248899 -0.25 -0.586082 1
0.248899 0.00055639 -0.2033 1 0.248899 0.250556 0.273169 1 0.248899 0.417223 0.221723 1
0.248899 0.5 0 1 0.415565 -0.5 0 1 0.415565 -0.416667 -0.219452 1
0.415565 -0.25 -0.354569 1 0.415565 0.00055639 -0.0819158 1 0.415565 0.250556 0.228346 1
0.415565 0.417223 0.167544 1 0.415565 0.5 0 1 0.5 -0.5 0 1
0.5 -0.416667 0 1 0.5 -0.25 0 1 0.5 0.00055639 0 1
0.5 0.250556 0 1 0.5 0.417223 0 1 0.5 0.5 0 1]
WorldEnd
#endif
[Affine Toolkit]
[RIB Utilities]
[Bitmap Utilities]
[Handy Little Utilities]
[Libraries]
[Using the Libraries]