以下内容翻译自Mathworks官方软件工程师Adam Danz的博客文章《Creating natural textures with power-law noise: clouds, terrains, and more》，值得一提的是他同时也是 MATLAB Graphics and Charting团队的开发人员。
我发现这篇文章是因为这篇文章中举例子阶段恰好用到了我在比赛中的一张图，因此被@。此篇文章很好的讲解了如何利用幂律噪声制作一些很真实自然的云、地形和纹理。感觉很值得分享出来，于是将其翻译了一下。
正文
您是否曾对 Flipbook Mini Hack 比赛中创意无限的参赛者所打造的令人叹为观止的自然风貌动画赞叹不已？例如@Tim，@Jenny Bosten和@Zhaoxu Liu/slandarer，在比赛中的非凡作品：
@Tim
https://www.mathworks.com/matlabcentral/communitycontests/contests/6/entries/13177
@Jenny Bosten
https://www.mathworks.com/matlabcentral/communitycontests/contests/6/entries/13695
@Tim
https://www.mathworks.com/matlabcentral/communitycontests/contests/6/entries/13820
@Zhaoxu Liu/slandarer
https://www.mathworks.com/matlabcentral/communitycontests/contests/6/entries/12757
从前到后，分别是来自2023年 MATLAB Flipbook Mini Hack 比赛的四个作品 Rolling fog, Blowin’in the wind, Morning ascent, landScape。
这种类型的结构化噪声有几种分类，包括 1/f 噪声，相干噪声，分形噪声，有色噪声，粉红噪声，分形布朗运动等，每种都有其自己的定义。在本文中，我将把这种类型的噪声称为幂律噪声，这个术语描述了信号强度与空间或时间频率之间的数学关系。幂律噪声可以生成宏观尺度的变化，如丘陵和山谷，同时具有微观尺度的变化，增加了纹理。这种噪声经常出现在从地震到金融市场等自然现象中，而且在 MATLAB 中生成它相当容易。
五个简单步骤在MATLAB生成二维幂律噪声
创建幂律噪声出奇地简单，只需要在 MATLAB 中使用三个参数和几行代码。我将为您讲解下面总结图中所描述的五个步骤。
5步骤总结图
第1步：随机噪声
创建一个二维矩阵，其中填充随机值。一些幂律噪声的定义要求从高斯白噪声开始，可以使用randn()生成，但使用rand()生成均匀分布的随机数也能产生良好的结果。前两个参数是矩阵的大小和随机数生成器的种子（以及生成器类型）。幂律噪声受初始值的影响很大，因此每个种子和矩阵大小都会产生不同的结果。大小参数 k 必须是偶数。
rng(0,'twister') % Parameter 1. random number generator seed
k = 100; % Parameter 2. size of the square matrix, must be an even scalar
m = randn(k); % k-by-k matrix
第2步：傅里叶变换
使用二维快速傅里叶变换函数fft2()将二维数组从空间域转换为频率域。这将返回一个金字塔排列的频率系数矩阵，其中低频系数位于矩阵的周边，而高频系数位于矩阵的中间。第四步需要对这种排列进行反转，使得低频系数位于中间，而高频系数位于矩阵的周边。使用fftshift()来反转排列。
mf = fftshift(fft2(m));
如何解释二维频率系数？二维傅里叶变换将图像分解为一组在不同频率和方向上的平面正弦波。频率系数是复数，提供了一组指令，说明如何组合或叠加这些平面波以重新创建图像。考虑 mf(u,v) 处的值。mf(u,v) 的实部和虚部对应于由分配给行 u 和列 v 的频率定义的二维正弦平面的幅度和相位。由于 MATLAB 的图像函数不支持复数值，因此在上面的总结图的第二和第四个面板中显示了频率系数的幅度。复数的幅度使用abs()计算。
第3步：创建滤波器
幂律噪声的一个重要特征是低频率的主导性，同时仍允许高频率对变化做出贡献。这在图像中产生宏观模式，比如二维山丘或云层，同时保持添加纹理的微观变化。就频率而言，具有细节的图像部分具有高空间频率，而较平滑的部分具有低空间频率。滤波器的目标是强调图像中的较低空间频率，以创建宏观模式。回想一下第二步，移位后的傅里叶变换矩阵 mf 的组织方式，使得较低频率波的系数位于中心。因此，2D滤波器应该与频率呈反比关系 ，其中 ， 使得在中心对较低频率施加更大的权重，在周边对较高频率施加较小的权重。在这个过程中，对较低频率的强调程度由第三个参数控制，即 中的指数，通常在 的范围内。较大的 值导致低空间频率的更大影响，从而产生平滑的噪声过渡，纹理细节较少，如下图所示。 的值用于定义几种噪声分类，如白噪声 、粉红噪声 和布朗噪声 。
指数值显示在每个图像的上方。指数越大，噪声越平滑、频率越低
每个我见过的幂律噪声实现都以不同方式创建滤波器，但遵循相同的权重分布模式。这个滤波器在第五步中返回实数值时非常稳健，尽管需要 k 为偶数。
a = 2; % Parameter 3. 1/fᵅ exponent
d = ((1:k)-(k/2)-1).^2;
dd = sqrt(d(:) + d(:)');
filt = dd .^ -a; % frequency filter
filt(isinf(filt))=1; % replace +/-inf at DC or zero-frequency component
第4步：滤波后的频率
将移位后的傅里叶变换矩阵与滤波器相乘。观察5步骤总结图的第四个面板，很难想象这个矩阵如何对最终图像有所贡献。我们习惯于在空间域中思考矩阵，其中坐标 (u,v) 处的值决定了图像中坐标 (u,v) 的像素强度。然而，该面板显示较低的空间频率朝着矩阵中心具有较高的系数，这意味着较低的空间频率将对最终结果产生更大的影响。
ff = mf .* filt;
第5步：逆傅里叶变换
我们从一个随机值矩阵开始。傅里叶变换将随机矩阵重建为频率域，并且滤波器强调了较低的频率。最后一步是将经过滤波的频率转换回空间域。使用ifftshift()将频率移回到其原始的金字塔排列，其中较低频率的系数位于周边，较高频率的系数位于中心。然后使用ifft2()进行二维逆傅里叶变换。
b = ifft2(ifftshift(ff));
结果(b)是一个由实数构成的矩阵，形成了与原始图像相同大小的结构化噪声。使用imagesc(b)查看5步骤总结图中的结果，并应用一个色条来调查图像的强度值。通常将噪声值缩放到[-1,1]范围是很常见的，这样更容易应用额外的缩放。为了缩放结果，可以使用rescale()函数。
bScaled = rescale(b,-1,1);
这种将滤波器应用于傅里叶变换图像，然后将过程逆转回图像的模式在图像处理中很常见，且为学习图像滤波提供了一个很好的模型。
幂律噪声矩阵的可视化
幂律噪声数值数组通常应用于颜色数据、alpha值(透明度)，或作为3D表面的高度数值。下图展示了相同噪声矩阵应用于四种不同可视化效果的情况。
幂律噪声在MATLAB图形中的应用
为了生成这些图像，按照上述的五个步骤创建幂律噪声矩阵(bScaled)，但将大小参数替换为k=500。然后继续执行下面的代码块。第二和第三个面板使用自定义的配色 "terrainmap"，该配色可在该仓库被找到(https://github.com/mathworks/matlab-blog/tree/main/2023/powerLawNoise)。同时，Mapping Toolbox 提供了改进的地形配色 "demcmap"。
% Panel 1
figure()
imagesc(bScaled)
set(gca,'YDir','normal')
colormap(gca,sky(256))
colorbar
%
% Panel 2
figure()
imagesc(bScaled)
set(gca,'YDir','normal')
colormap(gca,terrainmap()) % custom colormap
colorbar
%
% Panel 3
figure()
x = linspace(-4, 4, width(bScaled));
y = linspace(-4, 4, height(bScaled));
surf(x,y,bScaled,EdgeColor='none')
colormap(gca,terrainmap()) % custom colormap
axis equal
%
% Panel 4
figure()
s = surf(zeros(k,k),FaceColor=[.28 .24 .54],EdgeColor='none',...
AlphaData=bScaled,FaceAlpha='flat');
view(2)
使用灯光显示纹理
由幂律噪声提供的纹理在图形对象中颜色变化不足时可能会被忽略。在 MATLAB 中使用光照可以通过为对象添加详细的阴影来显示纹理。下图显示了与上面的曲面图中相同的山脉。水位被添加为一个平面表面，其颜色由第二个噪声矩阵定义。使用 cameratoolbar 工具手动选择了光源位置，并将其数值提取并放置在下面的代码中以便复现。摄像机属性被设置为放大到山脉。
使用灯光突出纹理
rng(0,'twister')
k = 500;
mountainHeights = makenoise(k, 2); % mountain noise
figure()
x = linspace(-4, 4, k);
y = linspace(-4, 4, k);
mountainHeights = mountainHeights + 0.3; % raise the mountain
surf(x,y,mountainHeights,EdgeColor='none')
colormap(flipud(copper))
hold on
seaShading = makenoise(k, 1.5); % water color noise
nColors = 500; % number of colors to use for the water
seaColors = parula(nColors*3/2);
seaColors(nColors+1:end,:) = []; % Capture the blue-green end of parula
% convert seaShading (-1:1) to colormap indices (1:500)
seaShadingIdx = round((seaShading+1)*(nColors/2-.5)+1);
% Create mxnx3 matrix of CData
cdata = reshape(seaColors(seaShadingIdx,:),nColors,nColors,3);
S=surf(x*3, y*3, 0*seaShading,CData=cdata,EdgeColor='none');
axis equal off
clim([-0.5, 1.5])
light(Position = [-.8 .1 .75]) % add lighting
material dull
camva(10) % camera viewing angle
campos([-14 7 4]) % camera position
camtarget([.4 .8 .4]) % camera target
function noise = makenoise(k, a)
% Create the noise matrix
% k is the matrix size (scalar)
% a is the positive exponent (scalar)
% noise are noise values with range [-1,1]
m = randn(k);
mf = fftshift(fft2(m));
d = ((1:k)-(k/2)-1).^2;
dd = sqrt(d(:) + d(:)');
filt = dd .^ -a;
filt(isinf(filt))=1;
ff = mf .* filt;
b = ifft2(ifftshift(ff));
noise = rescale(b,-1,1);
end
无缝拼接
这种类型的噪声还具有无缝平铺的特性。如下图所示，矩阵的边缘在使用 parula 颜色地图的瓦片排列中无缝地连接在一起。
figure()
tiledlayout(3, 3, TileSpacing='none', Padding='tight')
for i = 1:9
nexttile()
imagesc(bScaled)
set(gca,'XTickLabel',[],'ytickLabel',[])
end
逆向工作：通过估计指数对噪声进行分类
一个一维的幂律噪声向量具有与二维情况相同的宏观层面的起伏和微观层面的纹理。下面的部分绘制了由指数 定义的 1000 个粉红噪声数值。
rng default % Parameter 1. random number generator seed
k = 1000; % Parameter 2. size of the vector, must be even.
m = randn(k,1);
mf = fftshift(fft(m));
a = 1; % Parameter 3. 1/fᵅ exponent
d = sqrt(((1:k).'-(k/2)-1).^2);
filt = d .^ (-a/2); % frequency filter for 1D
filt(isinf(filt))=1;
ff = mf .* filt; % apply the filter
v = ifft(ifftshift(ff)); % power-law noise vector
vScaled = rescale(v,-1,1); % scaled power-law noise [-1,1]
figure()
plot(vScaled)
grid on
title('1D power-law noise')
subtitle('α = 1')
假设您从噪声向量开始，并希望通过其指数值对噪声的颜色进行分类。指数可以近似为对数-对数坐标中 功率谱的斜率。幸运的是，MATLAB 的信号处理工具箱通过 pspectrum 函数使得这一任务变得简单。
% Compute the power and frequency of the power-law noise vector
% Requires signal processing toolbox
[power, freq] = pspectrum(vScaled);
%
% Apply the log transform
freqLog = log(freq);
powerLog = log(power);
%
% Remove inf values
infrm = isinf(freqLog) | isinf(powerLog);
freqLog(infrm) = [];
powerLog(infrm) = [];
%
% Fit the relationship between log-power and log-frequency
% Requires the Curve Fitting toolbox
F = fit(freqLog,powerLog,'poly1');
%
% Plot the power density
figure()
plot(freqLog,powerLog)
%
% Add the linear fit
hold on
h = plot(F);
h.AffectAutoLimits = 'off';
grid on
axis tight
%
% Print the slope in the title
coeffs = coeffvalues(F);
xlabel('log frequency')
ylabel('log power spectrum')
title(sprintf('Slope = %.2f', coeffs(1)))
线性拟合的斜率为 ，这将正确地将该信号分类为粉红噪声 。
总结
幂律噪声（如粉红噪声和布朗噪声）可用于创建具有自然噪声变化的二维和三维对象。 在 MATLAB 中创建幂律噪声矩阵非常简单，只需五个步骤。 幂律噪声创建的基本原理涉及通过加权频率系数的逆频率来过滤随机值的傅立叶变换的二维矩阵，然后应用逆傅立叶变换。 二维幂律噪声矩阵可以作为颜色数据、alpha 数据（透明度）、高度数据或强度值应用于图像和曲面。 在 MATLAB 中使用光照可以凸显曲面的纹理。 几位 MATLAB 用户已经利用这一技术在 Flipbook Mini Hack 比赛中创建了令人印象深刻的图像。
参考及附加资源
colored_noise(https://www.mathworks.com/matlabcentral/fileexchange/121108-colored-noise) is a wonderful tool on the File Exchange by Abhranil Das that creates power-law noise in any number of dimensions (also on https://github.com/abhranildas/coloured_noise). Abhranil Das' 2022 thesis on Camouflage Detection & Signal Discrimination (http://dx.doi.org/10.13140/RG.2.2.10585.80487 Ch. 5) contains a terrific interpretation of 2D Fourier transform coefficients. Paul Bourke's 1998 articles(https://paulbourke.net/fractals/noise/) on generating noise and applications in graphics is a source of inspiration and creativity.
注
本文翻译自：Creating natural textures with power-law noise: clouds, terrains, and more(https://blogs.mathworks.com/matlab/2023/11/29/creating-natural-textures-with-power-law-noise-clouds-terrains-and-more/) 本文所提到的竞赛：MATLAB Flipbook Mini Hack(https://www.mathworks.com/matlabcentral/contests/2023-matlab-mini-hack.html) 原作者的博客：Graphics and App Building blog(https://blogs.mathworks.com/graphics-and-apps/)
