surface occlgather(float samples = 256)
{
      normal Ns = faceforward(normalize(N),I);

      float hits = 0;
      gather("illuminance", P, Ns, PI/2, samples, "distribution", "cosine") {
        hits += 1;
      }

      float occlusion = hits / samples;

      Ci = (1.0 - occlusion) * Cs;
      Oi = 1;
}