Perlin noise interpolation

Question

I've been writing a Perlin noise generator in F# using this and this, and have been successful up until the interpolation part of the algorithm. This is the working code so far (you probably don't need to read the next 2 code chunks as they're just here for context, so don't be scared away):

// PerlinNoiseProvider3D.fs
open System
open System.Collections.Generic
open Axiom.Math
open Ops

// Instances return a pseudorandom noise value between -1 and 1. scale defines how far apart the grid points are spaced.
type PerlinNoiseProvider3D(scale) =
  let randGen = new Random()
  // Each point in the grid has a random gradient
  let grid = new Dictionary<Vector3, Vector3>()

  // I stole this handy algorithm from this SE question:
  // http://math.stackexchange.com/questions/44689/how-to-find-a-random-axis-or-unit-vector-in-3d
  let randomAngle() =
    // Random value between 0 and 2π
    let θ = randGen.NextDouble() * Math.PI * 2.0
    // Random value between -1 and 1
    let z = randGen.NextDouble() * 2.0 - 1.0
    // Compute the resulting point
    (sqrt(1.0 - (z**2.0)) * cos(θ)) @@ (sqrt(1.0 - (z**2.0)) * sin(θ)) @@ z

  // Returns the gradient (just a vector pointing in a particular direction) at the given GRID point (not world coordinates)
  let getGradientAt v =
    try
      grid.[v]
    with
      | :? KeyNotFoundException -> grid.[v] <- randomAngle(); grid.[v]

  // Calculates the influence of the gradient g at gp on p, which is the dot product of gp and a vector facing from gp to p. Note that gp is _not_ in grid coordinates!!
  let influence (p: Vector3) (gp: Vector3) (g: Vector3) =
    // Find the vector going from the gp to p 
    let f = p - gp
    // Dot product on the grandient and the vector facing p from gp
    float <| (g.x * f.x) + (g.y * f.y) + (g.z * f.y)

  member this.GetValue (point: Vector3) =
    let v: Vector3 = point / scale
    // There's gotta be a shorter way to do this...
    let extract (i: Vector3) (min: Vector3) (max: Vector3) =
      let x =
        match int i.x with
        | 0 -> min.x
        | 1 -> max.x
        | _ -> failwith "Bad x index (shouldn't happen!)"
      let y =
        match int i.y with
        | 0 -> min.y
        | 1 -> max.y
        | _ -> failwith "Bad y index (shouldn't happen!)"
      let z =
        match int i.z with
        | 0 -> min.z
        | 1 -> max.z
        | _ -> failwith "Bad z index (shouldn't happen!)"
      x @@ y @@ z
    let min = (floor <| float v.x) @@ (floor <| float v.y) @@ (floor <| float v.z)
    let max = (ceil <| float v.x) @@ (ceil <| float v.y) @@ (ceil <| float v.z)
    // Relative 3D array of neighboring points (fst in tuple) and data (snd in tuples)
    let neighbors = Array3D.init 2 2 2 (fun xi yi zi -> extract (xi @@ yi @@ zi) min max, 0 @@ 0 @@ 0)
    let pInfluence = influence v
    let influences = neighbors |> Array3D.map (fun (p, g) -> pInfluence p (getGradientAt p))

and

// Ops.fs
module Ops
  let average (numbers: float list) =
    let rec average' (numbers': float list) =
      if numbers'.Length > 1 then
        average' [
          for i in 0..2..(numbers'.Length - 1) ->
            let a = numbers'.[i]
            try
              (a + numbers'.[i + 1]) / 2.0
            with
            | :? System.ArgumentException -> a
        ]
      else
        numbers'
    (average' numbers).[0]

  // Temporary 3D array average
  let average3D numbers =
    // Simply flatten the list and use the average function defined earlier
    average [
      for x in 0..(Array3D.length1 numbers - 1) do
        for y in 0..(Array3D.length2 numbers - 1) do
          for z in 0..(Array3D.length3 numbers - 1) ->
            numbers.[x, y, z]
    ]

I know that's a lot to take in, but none of the above code given is broken; it works as expected to produce incomplete images like this (it's just missing some):

This is what I've come up with for the interpolation (just the other code + this, removing the average3D call as the result):

// PerlinNoiseProvider3D.fs
// ...
let interp((a: float, b), axis) = a + (axis * (b - a))
let Sx, Sy, Sz = S <| float v.x, S <| float v.y, S <| float v.z
let zE1, zE2, zE3, zE4 =
  (influences.[0, 0, 0], influences.[0, 0, 1]),
  (influences.[0, 1, 0], influences.[0, 1, 1]),
  (influences.[1, 0, 0], influences.[1, 0, 1]),
  (influences.[1, 1, 0], influences.[1, 1, 1])
// Interpolate the points along the z-axis
let xE1, xE2 =
  (interp(zE1, Sz), interp(zE2, Sz)),
  (interp(zE3, Sz), interp(zE4, Sz))
// Interpolate the edges along the x-axis
let yE = (interp(xE1, Sx), interp(xE2, Sx))
// And finally, interpolate the y edge to yield our final value
interp(yE, Sy)

and

// Ops.fs
// ..
// Ease curve
let ease p = (3.0 * (p ** 2.0)) - (2.0 * (p ** 3.0))
let S x = ease (x - (floor x))

Which gives this result, which looks even worse:

I thought that you're basically just supposed to interpolate each of the axes -- in the case of 3D, all the 4 z-axes, say, then the resulting 2 x-axes, and finally the resulting y axis to get the final value. And it doesn't seem to work. I must be misunderstanding something here; it just doesn't work! Maybe my interpolation function is wrong. Maybe my application of it is wrong. Any help appreciated. Even an explanation of precisely what how you're supposed to do this -- most other sources for the final step say something like, "then interpolate these dot products together to get the final value."

PS: I'm using structures from the game library Axiom (a C# port of Ogre), mainly Vector3, so I defined an @@ operator that is used throughout this code and that is used to create Vector3's easily like so: 1 @@ 2 @@ 3

Answer 1

From here.

11.1.3 Normal vector interpolation This technique uses a similar interpolation method, except that the quantity interpolated across the polygon is the surface normal vector, instead of the computed intensity itself. (This method is also known as Phong interpolation.) The vertex normals are computed just as in the intensity interpolation method. Referring to Figure 11.2 again, and representing the normal vector at a point P as NP , we have

NQ =  u*NB+(1-u)NA,  where  u=  AQ/AB 
NR =  w*NB+(1-w)NC,  where  w=  CR/CB 
NP =  v*NR+(1-v)NQ,  where  v=  QP/QR 

The normal vector can also be determined incrementally.

The result is much more realistic than the intensity interpolation method, especially when specular reflection is considered - the highlights are rendered more faithfully, and Mach banding is greatly reduced (but spheres and cylinders are still prone to Mach band effects).

The drawback of this method is that it requires more computations to process each interpolated point in a polygon; at each point the complete shading calculation must be performed using the interpolated normal at that point.

11.1.4 Dot-product interpolation This method is a compromise between intensity interpolation and normal vector interpolation. Here, the quantities interpolated from one end of the scan line to the other are the dot products N· L and (R· V)n . This method is sometimes called 'cheap Phong' interpolation.

You will probably have to read the whole article to get enough context to understand the equations.