柏林噪声实践(一) 海波

　　这篇文章用于记录柏林噪声的一些实践，在开始前，先看下维斯百科里对柏林噪声的一些说明.

　　用随机法产生的噪声图像和显然自然界物体的随机噪声有很大差别，不够真实。1985年Ken Perlin指出^[1]，一个理想的噪声应该具有以下性质：

1. 对旋转具有统计不变性；
2. 能量在频谱上集中于一个窄带，即：图像是连续的，高频分量受限；
3. 对变换具有统计不变性。

　　先来看下一张图：

　　

　　这二张图都是模仿海波(只是看二者波形，颜色啥的先不要看。。。),同样的模型密度，一张是随机点生成高度的，一张是经过柏林噪声生成的高度。通过二张图的对比，让我们对柏林噪声先有个初步的认识。

　　查找了一些资料，发现对柏林噪声的定义各有不同，有些定义是柏林噪声是由一系列的Coherent noise构成，而在维基里，定义的是前面的Coherent noise是柏林噪声，而前面的柏林噪声应定义为分形噪声，就是说，分形噪声是由柏林噪声来组成的。在这里，因为大部分文章都采用第一张说法，这里的说明也采用第一种定义。

　　那面上面的第二张图就是一个Coherent noise(后面直接简称noise)生成的高度点，网上关于noise生成方法很多，有着色器版的，有C语言版的，有改进版的，还有经典版的，好吧，刚开始看完全找不到北，最后下下心来，把Nvidia的一个关于柏林噪声的例子下来，因为想方便追踪一下执行过程，把其中Cg代码改为F#代码。Nvida例子 http://http.download.nvidia.com/developer/SDK/Individual_Samples/samples.html 其中Vertex Noise程序  vnoise地址 。

　　noise分别有一维，多维，最能说明情况的应该是二维，明白二维nose的生成过程，一维和多维的就很清楚了。先来看下noise的原理，网上讲解的非常多，下面我认为能清楚的一些说明摘抄下来：

　　首先是维基百科里说的。　　

]

Perlin在上述文章中提出了一种产生符合要求的一维噪声函数的简单方法，这是后续工作的基础：

1. 在一维坐标轴上，选择等距的无穷个点，将值空间划分为等长的线段（为方便起见，选择坐标值为整数的点），为每个点随机指定一个值和一个梯度（在一维的情况，梯度就是斜率）；
2. 对于坐标值为整数的点，将该点对应的值作为噪声图像上该点的值；对于坐标值不为整数的点，将相邻两点的值和根据梯度进行插值运算，获得该点的灰度。

　　插值使用的函数是一个在0处为1，在1处为0，在0.5处为0.5的连续单调减函数。例如，设  为左右两整数点的颜色， 为该点距离左边点的距离，使用  作为插值函数，则该点的值为 。

　　 是线性插值，得到的结果人工痕迹严重，且在整数点上不连续。Perlin建议使用  作为插值函数 ^[2]。对于二维的情况，可以将上述算法进行推广，即：

1. 为所有坐标为  且  都是整数的点指定一个值，同时指定一个梯度，这些点将空间分成方格；
2. 对于坐标轴为整数的点，即上述方格的顶点，将为它指定的值作为该点的值；对于某个方格内部的点 ，用所在方格四个顶点的值和梯度进行插值。

　　例如，对于点 ，令  它所在方格的四个顶点分别为： 、、、。令 ，这四个顶点对点  的贡献可以用它们的梯度（）和  点与这四个顶点的方向（、、、）进行点积获得。但是在二维的情况下，插值更为复杂。首先需要对  和  两点在  方向插值，得到点  的值；之后对  和  两点在x方向插值，得到点  的值；最后对  和  在 y 方向插值，得到  的值。

　　在三维的情况下，需要进行七次插值。可以证明，插值次数随着维数的增长指数增长。

　　经典Perlin噪声基本满足Perlin提出的噪声条件。但是由于  导数中含有线性分量，在计算相邻点差时会体现出人工效果，不够自然。经典Perlin噪声在进行下文的分形和运算后效果不够自然^[2]。

　　这个全是文字说明，下面再来段图解，可以和上面合在一起看。出自 http://webstaff.itn.liu.se/~stegu/simplexnoise/simplexnoise.pdf

　　建议大家看完这个PDF，我最开始看noise的代码，老是觉的不可理喻，至到看到这个图，就理解其中的几个比较主要的步骤。下面看下针对上面的那个例子改写的F# noise版本。

  1 type PerlinNoise() =
  2     let B = 256
  3     let BM,DBN,N = B-1,2*B+2,B*B
  4     //let B,BM,N,NM,NP,DBN=256,255,4096,4095,12,(2*256+2) 
  5     //permutation 排列表 实现一种伪随机的效果
  6     let p = Array.create DBN 0
  7     //gradient 梯度表 索引在排列表中
  8     let pg = Array.create DBN vec3.Zero
  9     let curve t =
 10         // t * t * (3. - 2. * t) ) 
 11         t * t * t * (t * (t * 6. - 15.) + 10.)         
 12     let lerp t a b = a + t*(b-a)
 13     //得到点vec在对应维度上(x,y,z)上的二边的索引，以及二边的贡献点，二阶平滑插值
 14     let pos vec =
 15         //去掉负值 1/B*floor(P)+ONEHALF; 
 16         let t = vec + float N + 0.5 // Math.Floor(float vec)/(float B) + 1.0/(float (2*B))
 17         //let t = Math.Floor(float vec)/(float B) + 1.0/(float (2*B))
 18         //点d的左边,前面,下面索引  (指的是排列表中的位置索引)
 19         let x0 = int (Math.Floor(t)) &&& BM
 20         //点d的右边,右面,上面索引
 21         let x1 = x0 + 1 &&& BM
 22         //点d在当前维度上,用于对应梯度中左边(x轴),前面(z轴),下面(y轴)各点贡献值
 23         let u0 = t - Math.Floor(t)
 24          //其中u0+(-u1)=1.0,-(1.0-u0) 右边,后边,上面各点贡献值 (
 25         let u1 = u0 - 1.0   
 26         //二阶平滑位置
 27         let cf = curve u0
 28         x0,x1,u0,u1,cf   
 29     let init() =
 30         let rand = Random(2)
 31         for i = 0 to BM do
 32             p.[i] <- i
 33             pg.[i] <- vec3(rand.NextDouble()*2.0 - 1.0,rand.NextDouble()*2.0 - 1.0,rand.NextDouble()*2.0 - 1.0)
 34             //pg.[i]  <- vec2( float (rand.Next()), float (rand.Next()))
 35             pg.[B+i] <- pg.[i]
 36         for i = 0 to BM do
 37             let j = (rand.Next() >>> 5) % BM
 38             let t = p.[i]
 39             p.[i] <- p.[j]
 40             p.[j] <- t
 41             p.[i + B] <- p.[i]           
 42     do
 43         init()        
 44     member this.Noise1 (x:float) =
 45         let x0,x1,u0,u1,cx = pos x
 46         //得到二点的梯度值.
 47         let ug0 = u0 * pg.[p.[x0]].X
 48         //-(1.0-f0) 负表示
 49         let ug1 = u1 * pg.[p.[x1]].X
 50         //u v 平滑位置 noise1位置
 51         lerp cx ug0 ug1    
 52     member this.Noise2 (x,y) =
 53         let x0,x1,u0,u1,cx = pos x
 54         let y0,y1,v0,v1,cy = pos y
 55         //2维中权重也需要混合 下标DBN，p.[x0] + y0 最大为2B
 56         let gx0y0 = pg.[p.[p.[x0] + y0]] 
 57         let gx1y0 = pg.[p.[p.[x1] + y0]] 
 58         let gx0y1 = pg.[p.[p.[x0] + y1]] 
 59         let gx1y1 = pg.[p.[p.[x1] + y1]] 
 60         let mixuv (v3:vec3) u v = v3.X * u + v3.Y * v
 61         //yo轴上x混合值
 62         let u0v0 = mixuv gx0y0 u0 v0 //gx0y0.X * u0 + gx0y0.Y * v0
 63         let u1v0 = mixuv gx1y0 u1 v0//gx1y0.X * u1 + gx1y0.Y * v0
 64         let ccx0 = lerp cx u0v0 u1v0
 65         //y1轴上x混合值
 66         let u0v1 = mixuv gx0y1 u0 v1//gx0y1.X * u0 + gx0y1.Y * v1
 67         let u1v1 = mixuv gx1y1 u1 v1//gx1y1.X * u1 + gx1y1.Y * v1
 68         let ccx1 = lerp cx u0v1 u1v1
 69         //混合二y轴
 70         let result = lerp cy ccx0 ccx1
 71         // -1 --- 1映射成 0 --- 1
 72         result //* 0.5 + 0.5
 73     member this.Noise3 (x,y,z) =
 74         let x0,x1,u0,u1,cx = pos x
 75         let y0,y1,v0,v1,cy = pos y
 76         let z0,z1,r0,r1,cz = pos z
 77         //权重与当前点(一维左右二个,二维上下左右四个,三维八个)上各分量的点积
 78         let mixuvr (v3:vec3) u v r = v3.X * u + v3.Y * v + v3.Z * r
 79         //3维中权重也需要混合 zo面 同二维
 80         let gx0y0z0 = pg.[p.[p.[x0] + y0] + z0]
 81         let gx1y0z0 = pg.[p.[p.[x1] + y0] + z0]
 82         let gx0y1z0 = pg.[p.[p.[x0] + y1] + z0]
 83         let gx1y1z0 = pg.[p.[p.[x1] + y1] + z0]   
 84 
 85         let u0v0r0 = mixuvr gx0y0z0 u0 v0 r0
 86         let u1v0r0 = mixuvr gx1y0z0 u1 v0 r0
 87         let ccx0z0 = lerp cx u0v0r0 u1v0r0
 88 
 89         let u0v1r0 = mixuvr gx0y1z0 u0 v1 r0
 90         let u1v1r0 = mixuvr gx1y1z0 u1 v1 r0
 91         let ccx1z0 = lerp cx u0v1r0 u1v1r0    
 92         
 93         let ccyz0 = lerp cy ccx0z0 ccx1z0
 94         //z1面 同二维
 95         let gx0y0r1 = pg.[p.[p.[x0] + y0] + z1]
 96         let gx1y0r1 = pg.[p.[p.[x1] + y0] + z1]
 97         let gx0y1r1 = pg.[p.[p.[x0] + y1] + z1]
 98         let gx1y1r1 = pg.[p.[p.[x1] + y1] + z1]       
 99         
100         let u0v0r1 = mixuvr gx0y0r1 u0 v0 r1
101         let u1v0r1 = mixuvr gx1y0r1 u1 v0 r1
102         let ccx0z1 = lerp cx u0v0r1 u1v0r1
103 
104         let u0v1r1 = mixuvr gx0y1r1 u0 v1 r1
105         let u1v1r1 = mixuvr gx1y1r1 u1 v1 r1
106         let ccx1z1 = lerp cx u0v1r1 u1v1r1 
107           
108         let ccyz1 = lerp cy ccx0z1 ccx1z1
109         //混合z轴
110         lerp cz ccyz0 ccyz1     
111     //octave 倍频,amplitude振幅,调整高低,frequency频率,越小间隔越平滑
112     member this.PerlinNoise2(x,y,octave,amplitude,frequency) =
113         let mutable total = 0.0
114         for i in [|0 .. octave - 1|] do
115             let freq  = Math.Pow(frequency,float i)
116             let amp = Math.Pow(amplitude,float i)
117             total <- total + this.Noise2(x*freq,y*freq)*amp
118         total
119     //生成一张图片RGBA源,用于测试
120     member this.Test() =        
121         init()
122         let arrayrgb = ArrayList<byte>()
123         let mutable maxvalue = 0.0
124         for i in [| 0 .. 255|] do
125             for j in [| 0 .. 255|] do
126                // let r = this.Noise2 (float i + float i/256.,float j + float j/256.) 
127                 //let r = this.Noise2 (float i,float j ) 
128                 let r = this.Noise2 (float i/256.,float j/256. ) 
129                 let color = byte (r*128.) + 128uy
130                 arrayrgb.AddRange([| color; color; color;color |])         
131         arrayrgb.ToArray()
132        
133 //如果我们想要坐标(i,j,k)的gradient g(i,j,k),而P里预存储了256个gradient, 那么:
134 //    g(i, j, k) = P[ ( i + P[ (j + P[k]) mod 256 ] ) mod 256 ]

PerlinNoise