float worley(point p; 
            string spacename; 
            float  freq)
{
point pp = transform(spacename, p) * freq;
point thiscell = point(floor(xcomp(pp))+0.5,
                       floor(ycomp(pp))+0.5,
                       floor(zcomp(pp))+0.5);
float dist2nearest = 1000;
uniform 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 + vector cellnoise(testcell)-0.5;
           float dist = distance(pos,pp);
           if(dist < dist2nearest)
           dist2nearest = dist;
           }
return dist2nearest;
}
  
surface
preBakeSSS(string     filename = "", 
                    displaychannels = ""; 
                    color     Kdt = 1)
{
color irrad, rad_t;
normal Nn = normalize(N);
float a = area(P, "dicing"); // micropolygon area
  
/* Compute direct illumination (ambient and diffuse) */
irrad = ((float worley(P, "world", 0.6) - 0.5) * 0.4) + 0.9 * (ambient() + diffuse(Nn));
  
/* Compute the radiance diffusely transmitted through the surface.
   Kdt is the surface color (could use a texture for e.g. freckles) */
rad_t = Kdt * irrad;
  
/* Store in point cloud file */
bake3d(filename, displaychannels, P, Nn, "interpolate", 1,
   "_area", a, "_radiance_t", rad_t);
  
Ci = rad_t * Cs * Os;
Oi = Os;
}