{"id":1882,"date":"2023-11-29T10:20:07","date_gmt":"2023-11-29T15:20:07","guid":{"rendered":"https:\/\/blogs.mathworks.com\/matlab\/?p=1882"},"modified":"2023-11-29T10:20:07","modified_gmt":"2023-11-29T15:20:07","slug":"creating-natural-textures-with-power-law-noise-clouds-terrains-and-more","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/matlab\/2023\/11\/29\/creating-natural-textures-with-power-law-noise-clouds-terrains-and-more\/","title":{"rendered":"Creating natural textures with power-law noise: clouds, terrains, and more"},"content":{"rendered":"
\r\n
\"\"<\/a><\/div>\r\n
Today's guest writer is Adam Danz<\/a> whose name you might recognize from the MATLAB Central community. Adam is a developer in the MATLAB Graphics and Charting team and the host of the Graphics and App Building blog<\/span><\/a>. In this article he\u2019ll share insights to creating structured noise and applying it to MATLAB graphics.<\/div>\r\n
Noise is ubiquitous in nature. MATLAB has several tools to remove noise from data and images but in this article I'll review a method of generating<\/span> noise that is currently being used by participants in the Flipbook Mini Hack contest<\/a> to generate natural looking clouds, terrains, and textures including these examples below.<\/div>\r\n
\"CombinedAnimation.gif\"<\/div>\r\n
From left to right, <\/span>Rolling fog<\/span><\/a>, <\/span>Blowin'in the wind<\/span><\/a>, <\/span>Morning ascent<\/span><\/a>, <\/span>landScape<\/span><\/a>, from the 2023 MATLAB Flipbook Mini Hack contest.<\/span><\/div>\r\n
This type of structured noise has several classifications including 1\/f<\/span> noise, coherent noise, fractal noise, colored noise, pink noise, fBm, and Brownian noise, each with their own definitions. For this article, I will refer to this type of noise as power-law noise<\/span>, a term that describes the mathematical relationship between a signal's intensity and spatial or temporal frequency. Power-law noise can generate macro-scale variation such as hills and valleys while featuring micro-scale variation that adds texture. This type of noise is often found in natural phenomena from earthquakes to financial markets and it is quite easy to generate it in MATLAB.<\/div>\r\n

Five Simple steps to generate 2D power-law noise in MATLAB<\/h2>\r\n
Creating power-law noise is surprisingly simple and requires only three parameters and a few lines of code in MATLAB. I'll walk you through the five steps depicted in the summary image below.<\/div>\r\n
\"\"<\/div>\r\n
Summary of the five steps to creating 2D power-law noise<\/span><\/div>\r\n
1. Random noise<\/span>. Create a 2D matrix of random values. Some definitions of power-law noise require starting with Gaussian white noise which can be generated using randn()<\/span> but uniformly distributed random numbers using rand()<\/span> also produce fine results. The first two parameters are the size of the matrix and the random number generator seed (and generator type). Power-law noise is highly influenced by the initial values, so each seed and matrix size produces different results. The size parameter k<\/span> must be even.<\/div>\r\n
\r\n
\r\n
rng(0,'twister'<\/span>) % Parameter 1. random number generator seed<\/span><\/div>\r\n<\/div>\r\n
\r\n
k = 100; % Parameter 2. size of the square matrix, must be an even scalar<\/span><\/div>\r\n<\/div>\r\n
\r\n
m = randn(k); % k-by-k matrix<\/span><\/div>\r\n<\/div>\r\n<\/div>\r\n
2. Fourier transform<\/span>. Transform the 2D array from spatial to frequency domains using the 2D fast Fourier transform function fft2()<\/span>. This returns a matrix of frequency coefficients in a pyramid arrangement where coefficients for low frequencies are on the perimeter of the matrix and coefficients for high frequencies are in the middle of the matrix. Step four requires an inversion of this arrangement where coefficients for low frequencies are in the middle and coefficients for high frequencies are on the perimeter of the matrix. Use fftshift()<\/span> to invert the arrangement.<\/div>\r\n
\r\n
\r\n
mf = fftshift(fft2(m));<\/div>\r\n<\/div>\r\n<\/div>\r\n
How do I interpret the 2D frequency coefficients? <\/span>This paragraph is for curious readers and for my personal challenge to concisely explain it. 2D Fourier transforms decompose an image into a set of planar sinusoidal waves at various frequencies and orientations. The frequency coefficients are complex numbers that provide a set of instructions for how to combine or superimpose the planar waves to recreate the image. Consider the value at mf(u,v)<\/span>. The real and imaginary components of mf(u,v)<\/span> correspond to the amplitude and phase of the 2D sinusoidal plane that is defined by the frequencies assigned to row u<\/span> and column v<\/span>. Since MATLAB's image functions do not support complex values, the magnitudes of the frequency coefficients are shown in the second and fourth panels of the summary figure above. The magnitude of a complex numbers is computed using abs()<\/span>.<\/div>\r\n
3. Create the filter<\/span>. An important characteristic of power-law noise is the dominance of lower frequencies while still allowing higher frequencies to contribute to the variation. This creates macro patterns in the image such as 2D hills or cloud bellows while maintaining micro variations that add texture. In terms of frequency, sections of an image with fine details have high spatial frequency while smoother sections have low spatial frequency. The goal of the filter is to accentuate lower spatial frequencies in the image to create the macro patterns. Recall from step two that the shifted Fourier transform matrix mf<\/span> is organized such that coefficients for lower frequency waves are in the center. The 2D filter, therefore, should have an inverse relationship to frequency ($ 1\/f$ where $f = \\sqrt{f_u^2 + f_v^2}$ for 2D signals)<\/span> such that greater weight is placed on lower frequencies in the center and less weight on higher frequencies near the perimeter. The amount of emphasis placed on lower frequencies is controlled by the third parameter in this process, the exponent in $1\/f^\u03b1$\u00a0which is typically within the range of 0<\u03b1<=3<\/span>. Larger values of \u03b1<\/span> result in greater influence from low spatial frequencies which result in smoother noise transitions with less detailed texture as shown in the image below. The value of \u03b1<\/span> is used to define several classifications of noise such as white noise (\u03b1=0<\/span>) pink noise (\u03b1=1<\/span>) and Brownian noise (\u03b1=2)<\/span>.<\/div>\r\n
\"\"<\/div>\r\n
Exponent values are shown at the top of each image. Larger exponents result in smoother, lower frequency noise<\/span><\/div>\r\n
Every implementation of power-law noise I've seen creates the filter differently but follows the same pattern of weight distribution. This filter is quite robust at returning real number values in step five, though it requires that k<\/span> is even.<\/div>\r\n
\r\n
\r\n
a = 2; % Parameter 3. 1\/f\u1d45 exponent<\/span><\/div>\r\n<\/div>\r\n
\r\n
d = ((1:k)-(k\/2)-1).^2;<\/div>\r\n<\/div>\r\n
\r\n
dd = sqrt(d(:) + d(:)');<\/div>\r\n<\/div>\r\n
\r\n
filt = dd .^ -a; % frequency filter<\/span><\/div>\r\n<\/div>\r\n
\r\n
filt(isinf(filt))=1; % replace +\/-inf at DC or zero-frequency component<\/span><\/div>\r\n<\/div>\r\n<\/div>\r\n
4. Filtered frequencies<\/span>. Multiply the shifted Fourier transform matrix by the filter. Looking at the fourth panel of the 5-step summary image, it is difficult to imagine how this matrix contributes to the final image. We're used to thinking about matrices in the spatial domain where the value at coordinate (u,v)<\/span> determines the intensity of the pixel at (u,v)<\/span> of the image. Instead, this panel shows that lower spatial frequencies toward the middle of the matrix have higher coefficients which means that the lower spatial frequencies will have a greater effect on the final results.<\/div>\r\n
\r\n
\r\n
ff = mf .* filt;<\/div>\r\n<\/div>\r\n<\/div>\r\n
5. Inverse Fourier transform<\/span>. We started with a matrix of random values. A Fourier transform reconstructed the random matrix in the frequency domain and the filter accentuated lower frequencies. The final step is to convert the filtered frequencies back to the spatial domain. Shift the frequencies back to their original pyramidal arrangement with coefficients for lower frequencies on the perimeter and coefficients for higher frequencies at the center using ifftshift()<\/span>. Then apply a 2D inverse Fourier transform using iff2()<\/span>.<\/div>\r\n
\r\n
\r\n
b = ifft2(ifftshift(ff));<\/div>\r\n<\/div>\r\n<\/div>\r\n
The result (b<\/span>) is a matrix of real numbers forming structured noise the same size as the original image. View the results shown in panel five of the 5-step summary image using imagesc(b)<\/span> and apply a color bar to investigate the image's intensity values.<\/div>\r\n
It is common to scale the noise values to a range of [-1,1] which makes it easier to apply additional customized scaling. To scale the results, use rescale()<\/span>.<\/div>\r\n
\r\n
\r\n
bScaled = rescale(b,-1,1);<\/div>\r\n<\/div>\r\n<\/div>\r\n
This pattern of applying a filter to a Fourier transformed image and then inversing the process back to an image is common in image processing and provides a nice model for learning about image filters.<\/div>\r\n

Visualizing the power-law noise matrix.<\/h2>\r\n
The array of power-law noise values is commonly applied to color data, alpha value (transparency), or as height values in a 3D surface. The figure below shows the same noise matrix applied to four different visualizations.<\/div>\r\n
\"\"<\/div>\r\n
Applications of power-law noise in MATLAB graphics<\/span><\/div>\r\n
To produces these panels, create the power-law noise matrix (bScaled<\/span>) following the five steps above but replace the size parameter with k=500<\/span>. Then continue with the code block below. The second and third panels use a custom colormap terrainmap<\/span> which is attached to this article (download<\/a>). Note that the Mapping Toolbox has an improved terrain colormap, demcmap<\/span>.<\/div>\r\n
\r\n
\r\n
% Panel 1<\/span><\/div>\r\n<\/div>\r\n
\r\n
figure()<\/div>\r\n<\/div>\r\n
\r\n
imagesc(bScaled)<\/div>\r\n<\/div>\r\n
\r\n
set(gca,'YDir'<\/span>,'normal'<\/span>)<\/div>\r\n<\/div>\r\n
\r\n
colormap(gca,sky(256))<\/div>\r\n<\/div>\r\n
\r\n
colorbar<\/div>\r\n<\/div>\r\n
\r\n
%<\/span><\/div>\r\n<\/div>\r\n
\r\n
% Panel 2<\/span><\/div>\r\n<\/div>\r\n
\r\n
figure()<\/div>\r\n<\/div>\r\n
\r\n
imagesc(bScaled)<\/div>\r\n<\/div>\r\n
\r\n
set(gca,'YDir'<\/span>,'normal'<\/span>)<\/div>\r\n<\/div>\r\n
\r\n
colormap(gca,terrainmap()) % custom colormap<\/span><\/div>\r\n<\/div>\r\n
\r\n
colorbar<\/div>\r\n<\/div>\r\n
\r\n
%<\/span><\/div>\r\n<\/div>\r\n
\r\n
% Panel 3<\/span><\/div>\r\n<\/div>\r\n
\r\n
figure()<\/div>\r\n<\/div>\r\n
\r\n
x = linspace(-4, 4, width(bScaled));<\/div>\r\n<\/div>\r\n
\r\n
y = linspace(-4, 4, height(bScaled));<\/div>\r\n<\/div>\r\n
\r\n
surf(x,y,bScaled,EdgeColor='none'<\/span>)<\/div>\r\n<\/div>\r\n
\r\n
colormap(gca,terrainmap()) % custom colormap<\/span><\/div>\r\n<\/div>\r\n
\r\n
axis equal<\/span><\/div>\r\n<\/div>\r\n
\r\n
%<\/span><\/div>\r\n<\/div>\r\n
\r\n
% Panel 4<\/span><\/div>\r\n<\/div>\r\n
\r\n
figure()<\/div>\r\n<\/div>\r\n
\r\n
s = surf(zeros(k,k),FaceColor=[.28 .24 .54],EdgeColor='none'<\/span>,...<\/span><\/div>\r\n<\/div>\r\n
\r\n
\u00a0\u00a0\u00a0\u00a0AlphaData=bScaled,FaceAlpha='flat'<\/span>);<\/div>\r\n<\/div>\r\n
\r\n
view(2)<\/div>\r\n<\/div>\r\n<\/div>\r\n

Use light to reveal texture<\/h3>\r\n
The texture provided by power-law noise can go unnoticed when there isn't sufficient color variation in the graphics object. The use of lighting in MATLAB reveals the texture by adding detailed shading to objects. The image below shows the same mountain as in the surface plot above. Water level was added as a planar surface with color defined by a second noise matrix. Light position was selected manually using the cameratoolbar<\/a> tool and its values were extracted and placed in the code below for reproducibility. Camera properties were set to zoom in to the mountain.<\/div>\r\n
\"\"<\/div>\r\n
The use of lighting accentuates texture<\/span><\/div>\r\n
\r\n
rng(0,'twister'<\/span>)<\/div>\r\n
k = 500;<\/div>\r\n
mountainHeights = makenoise(k, 2); % mountain noise<\/span><\/div>\r\n
<\/div>\r\n
figure()<\/div>\r\n
x = linspace(-4, 4, k);<\/div>\r\n
y = linspace(-4, 4, k);<\/div>\r\n
mountainHeights = mountainHeights + 0.3; % raise the mountain<\/span><\/div>\r\n
surf(x,y,mountainHeights,EdgeColor='none'<\/span>)<\/div>\r\n
colormap(flipud(copper))<\/div>\r\n
<\/div>\r\n
hold on<\/span><\/div>\r\n
seaShading = makenoise(k, 1.5); % water color noise<\/span><\/div>\r\n
nColors = 500; % number of colors to use for the water<\/span><\/div>\r\n
seaColors = parula(nColors*3\/2);<\/div>\r\n
seaColors(nColors+1:end,:) = []; % Capture the blue-green end of parula<\/span><\/div>\r\n
% convert seaShading (-1:1) to colormap indices (1:500)<\/span><\/div>\r\n
seaShadingIdx = round((seaShading+1)*(nColors\/2-.5)+1);<\/div>\r\n
% Create mxnx3 matrix of CData<\/span><\/div>\r\n
cdata = reshape(seaColors(seaShadingIdx,:),nColors,nColors,3);<\/div>\r\n
S=surf(x*3, y*3, 0*seaShading,CData=cdata,EdgeColor='none'<\/span>);<\/div>\r\n
<\/div>\r\n
axis equal off<\/span><\/div>\r\n
<\/div>\r\n
clim([-0.5, 1.5])<\/div>\r\n
light(Position = [-.8 .1 .75]) % add lighting<\/span><\/div>\r\n
material dull<\/span><\/div>\r\n
camva(10) % camera viewing angle<\/span><\/div>\r\n
campos([-14 7 4]) % camera position<\/span><\/div>\r\n
camtarget([.4 .8 .4]) % camera target<\/span><\/div>\r\n
<\/div>\r\n
function <\/span>noise = makenoise(k, a)<\/div>\r\n
% Create the noise matrix<\/span><\/div>\r\n
% k is the matrix size (scalar)<\/span><\/div>\r\n
% a is the positive exponent (scalar)<\/span><\/div>\r\n
% noise are noise values with range [-1,1]<\/span><\/div>\r\n
m = randn(k);<\/div>\r\n
mf = fftshift(fft2(m));<\/div>\r\n
d = ((1:k)-(k\/2)-1).^2;<\/div>\r\n
dd = sqrt(d(:) + d(:)');<\/div>\r\n
filt = dd .^ -a;<\/div>\r\n
filt(isinf(filt))=1;<\/div>\r\n
ff = mf .* filt;<\/div>\r\n
b = ifft2(ifftshift(ff));<\/div>\r\n
noise = rescale(b,-1,1);<\/div>\r\n
end<\/span><\/div>\r\n<\/div>\r\n

Seamless tiling<\/h3>\r\n
Seamless tiling is another neat property of this style of noise. The edges of the matrix wrap seamlessly in a tiled arrangement as shown in the image below using the parula<\/span> colormap.<\/div>\r\n
\r\n
\r\n
figure()<\/div>\r\n<\/div>\r\n
\r\n
tiledlayout(3, 3, TileSpacing='none'<\/span>, Padding='tight'<\/span>)<\/div>\r\n<\/div>\r\n
\r\n
for <\/span>i = 1:9<\/div>\r\n<\/div>\r\n
\r\n
\u00a0\u00a0\u00a0\u00a0nexttile()<\/div>\r\n<\/div>\r\n
\r\n
\u00a0\u00a0\u00a0\u00a0imagesc(bScaled)<\/div>\r\n<\/div>\r\n
\r\n
\u00a0\u00a0\u00a0\u00a0set(gca,'XTickLabel'<\/span>,[],'ytickLabel'<\/span>,[])<\/div>\r\n<\/div>\r\n
\r\n
end<\/span><\/div>\r\n<\/div>\r\n<\/div>\r\n
\"\"<\/div>\r\n
3x3 tile of the same noise matrix shows circular wrapping around its edges<\/span><\/div>\r\n

Working in reverse: classifying noise by estimating the exponent<\/h2>\r\n
A 1D vector of power-law noise has the same macro-level hills and micro-level texture as the 2D case. This section below plots 1000 values of pink noise defined by an exponent of \u03b1=1<\/span>.<\/div>\r\n
\r\n
\r\n
rng default<\/span> % Parameter 1. random number generator seed<\/span><\/div>\r\n<\/div>\r\n
\r\n
k = 1000; % Parameter 2. size of the vector, must be even.<\/span><\/div>\r\n<\/div>\r\n
\r\n
m = randn(k,1);<\/div>\r\n<\/div>\r\n
\r\n
<\/div>\r\n<\/div>\r\n
\r\n
mf = fftshift(fft(m));<\/div>\r\n<\/div>\r\n
\r\n
<\/div>\r\n<\/div>\r\n
\r\n
a = 1; % Parameter 3. 1\/f\u1d45 exponent<\/span><\/div>\r\n<\/div>\r\n
\r\n
d = sqrt(((1:k).'-(k\/2)-1).^2);<\/div>\r\n<\/div>\r\n
\r\n
filt = d .^ (-a\/2); % frequency filter for 1D<\/span><\/div>\r\n<\/div>\r\n
\r\n
filt(isinf(filt))=1;<\/div>\r\n<\/div>\r\n
\r\n
<\/div>\r\n<\/div>\r\n
\r\n
ff = mf .* filt; % apply the filter<\/span><\/div>\r\n<\/div>\r\n
\r\n
v = ifft(ifftshift(ff)); % power-law noise vector<\/span><\/div>\r\n<\/div>\r\n
\r\n
vScaled = rescale(v,-1,1); % scaled power-law noise [-1,1]<\/span><\/div>\r\n<\/div>\r\n
\r\n
<\/div>\r\n<\/div>\r\n
\r\n
figure()<\/div>\r\n<\/div>\r\n
\r\n
plot(vScaled)<\/div>\r\n<\/div>\r\n
\r\n
grid on<\/span><\/div>\r\n<\/div>\r\n
\r\n
title('1D power-law noise'<\/span>)<\/div>\r\n<\/div>\r\n
\r\n
subtitle('\u03b1 = 1'<\/span>)<\/div>\r\n<\/div>\r\n<\/div>\r\n
\"\"<\/div>\r\n
Suppose you start with the noise vector and you want to classify the color of the noise by its exponent value. The exponent can be approximated as the slope of the $1\/f^\u03b1$\u00a0power spectrum in log-log coordinates. Fortunately, MATLAB's Signal Processing toolbox makes this task easy with the pspectrum<\/span> function.<\/div>\r\n
\r\n
\r\n
% Compute the power and frequency of the power-law noise vector<\/span><\/div>\r\n<\/div>\r\n
\r\n
% Requires signal processing toolbox<\/span><\/div>\r\n<\/div>\r\n
\r\n
[power, freq] = pspectrum(vScaled);<\/div>\r\n<\/div>\r\n
\r\n
%<\/span><\/div>\r\n<\/div>\r\n
\r\n
% Apply the log transform<\/span><\/div>\r\n<\/div>\r\n
\r\n
freqLog = log(freq);<\/div>\r\n<\/div>\r\n
\r\n
powerLog = log(power);<\/div>\r\n<\/div>\r\n
\r\n
%<\/span><\/div>\r\n<\/div>\r\n
\r\n
% Remove inf values<\/span><\/div>\r\n<\/div>\r\n
\r\n
infrm = isinf(freqLog) | isinf(powerLog);<\/div>\r\n<\/div>\r\n
\r\n
freqLog(infrm) = [];<\/div>\r\n<\/div>\r\n
\r\n
powerLog(infrm) = [];<\/div>\r\n<\/div>\r\n
\r\n
%<\/span><\/div>\r\n<\/div>\r\n
\r\n
% Fit the relationship between log-power and log-frequency<\/span><\/div>\r\n<\/div>\r\n
\r\n
% Requires the Curve Fitting toolbox<\/span><\/div>\r\n<\/div>\r\n
\r\n
F = fit(freqLog,powerLog,'poly1'<\/span>);<\/div>\r\n<\/div>\r\n
\r\n
%<\/span><\/div>\r\n<\/div>\r\n
\r\n
% Plot the power density<\/span><\/div>\r\n<\/div>\r\n
\r\n
figure()<\/div>\r\n<\/div>\r\n
\r\n
plot(freqLog,powerLog)<\/div>\r\n<\/div>\r\n
\r\n
%<\/span><\/div>\r\n<\/div>\r\n
\r\n
% Add the linear fit<\/span><\/div>\r\n<\/div>\r\n
\r\n
hold on<\/span><\/div>\r\n<\/div>\r\n
\r\n
h = plot(F);<\/div>\r\n<\/div>\r\n
\r\n
h.AffectAutoLimits = 'off'<\/span>;<\/div>\r\n<\/div>\r\n
\r\n
grid on<\/span><\/div>\r\n<\/div>\r\n
\r\n
axis tight<\/span><\/div>\r\n<\/div>\r\n
\r\n
%<\/span><\/div>\r\n<\/div>\r\n
\r\n
% Print the slope in the title<\/span><\/div>\r\n<\/div>\r\n
\r\n
coeffs = coeffvalues(F);<\/div>\r\n<\/div>\r\n
\r\n
xlabel('log frequency'<\/span>)<\/div>\r\n<\/div>\r\n
\r\n
ylabel('log power spectrum'<\/span>)<\/div>\r\n<\/div>\r\n
\r\n
title(sprintf('Slope = %.2f'<\/span>, coeffs(1)))<\/div>\r\n<\/div>\r\n<\/div>\r\n
\"\"<\/div>\r\n
The slope of the linear fit is -1.05 which would correctly classify this signal as pink noise (\u03b1=1<\/span>).<\/div>\r\n

Summary points<\/h2>\r\n