public static double GradientCoherentNoise3D (double x, double y, double z, int seed,
         NoiseQuality noiseQuality)
   {
      
      // Create a unit-length cube aligned along an integer boundary.  This cube
      // surrounds the input point.
      int x0 = (x > 0.0? (int)x: (int)x - 1);
      int x1 = x0 + 1;
      int y0 = (y > 0.0? (int)y: (int)y - 1);
      int y1 = y0 + 1;
      int z0 = (z > 0.0? (int)z: (int)z - 1);
      int z1 = z0 + 1;

      // Map the difference between the coordinates of the input value and the
      // coordinates of the cube's outer-lower-left vertex onto an S-curve.
      double xs = 0, ys = 0, zs = 0;
      switch (noiseQuality)
      {
         case QUALITY_FAST:
            xs = (x - (double)x0);
            ys = (y - (double)y0);
            zs = (z - (double)z0);
            break;
         case QUALITY_STD:
            xs = Interp.SCurve3 (x - (double)x0);
            ys = Interp.SCurve3 (y - (double)y0);
            zs = Interp.SCurve3 (z - (double)z0);
            break;
         case QUALITY_BEST:
            xs = Interp.SCurve5 (x - (double)x0);
            ys = Interp.SCurve5 (y - (double)y0);
            zs = Interp.SCurve5 (z - (double)z0);
            break;
      }

      // Now calculate the noise values at each vertex of the cube.  To generate
      // the coherent-noise value at the input point, interpolate these eight
      // noise values using the S-curve value as the interpolant (trilinear
      // interpolation.)
      double n0, n1, ix0, ix1, iy0, iy1;
      n0   = GradientNoise3D (x, y, z, x0, y0, z0, seed);
      n1   = GradientNoise3D (x, y, z, x1, y0, z0, seed);
      ix0  = Interp.linearInterp (n0, n1, xs);
      n0   = GradientNoise3D (x, y, z, x0, y1, z0, seed);
      n1   = GradientNoise3D (x, y, z, x1, y1, z0, seed);
      ix1  = Interp.linearInterp (n0, n1, xs);
      iy0  = Interp.linearInterp (ix0, ix1, ys);
      n0   = GradientNoise3D (x, y, z, x0, y0, z1, seed);
      n1   = GradientNoise3D (x, y, z, x1, y0, z1, seed);
      ix0  = Interp.linearInterp (n0, n1, xs);
      n0   = GradientNoise3D (x, y, z, x0, y1, z1, seed);
      n1   = GradientNoise3D (x, y, z, x1, y1, z1, seed);
      ix1  = Interp.linearInterp (n0, n1, xs);
      iy1  = Interp.linearInterp (ix0, ix1, ys);

      return Interp.linearInterp (iy0, iy1, zs);
   }