2D Perlin noise generaton needs Perfomance

Question

I'm trying to use a Perlin noise generator to make the tiles of a map, but I think the code is too slow, I mean, I took some hours to complete the generation of 1000 x 1000 chart.
(I didn't used timeit.timeit() to check the time, I simply executed the script in the afternoon and in the evening I got the result...)

Note: My code is based in this pseudocode that I through it was a Perlin noise generator, but pikalek inform me that the page was wrong and it was a value noise generator.

My code is divided in 3 classes and in "2 parts of execution": map generation and result's shows. Also, it use only vanilla libraries to work, but to show the results you can use matplotlib and NumPy or pillow.

Classes:

• D:
  Inside this class it's the cubic interpolation that I use to make the interpolation of the values. This class inherit in another two classes: D1 (1D Perlin) and D2 (2D Perlin). Here I will only use D2, so I won't post the D1 class, but if you want it, you just have to ask in comment and I will publish it.
• D2: Make 2D noise and show it.

And yes, I know that cubic interpolation is much slower and I shall use a cosine interpolation. Using time.time() and the line a = D2(10000, 10) I gatherer this values:

Cubic:  46.46249842643738 seconds.
Cosine: 11.931403160095215 seconds.

But, just look at the picture, it talks by itself.
Cubic is awesome.

Parts of execution: (In the code you won't find it classified in "parts")

• Map generation: It use several functions to work: Cubic_Interpolate, Noise, Smooth_Noise, Interpolate_Noise and Perlin.

  • The generation start executing the Perlin function, which manages everything.
    Perlin iterate over the y and x axis of the grid (or map) and then iterate again but over the o (octaves). For each octave it doubles the frequency and amplitude and get a value from Interpolate_Noise, that value is added to the previously values. When the octaves end, the final value is appended to a list called line, and when the x axis iteration end, the line list is appended to the result list of list (the chart). Also, there is a line result list, that is similar to result but instead of been a list of lists (grid), it's just a linear list.

  • The magics happens inside Interpolate_Noise. It interpolate of the surroundings of the cord (x, y) by this way. This is the slowest part of the code, but I don't think it can be changed.

    Note that the image is in 1D...

                                                                        _ 
    (x-1, y+2) (x  , y+2) (x+1, y+2) (x+2, y+2) ---> Interpolate i4 (v3) |
    (x-1, y+1) (x  , y+1) (x+1, y+1) (x+2, y+1) ---> Interpolate i3 (v2) | Interpolate Complete!
    (x-1, y  ) (x  , y  ) (x+1, y  ) (x+2, y  ) ---> Interpolate i2 (v1) |
    (x-1, y-1) (x  , y-1) (x+1, y-1) (x+2, y-1) ---> Interpolate i1 (v0)_|
    

  • Interpolated_Noise get all its values from Smooth_Noise, a function that takes the average values from the surrondings of (x,y). So, it make a average with:

    (x-1,y+1) (x  , y+1) (x+1,y+1)
    (x-1,y  ) (x  , y  ) (x+1,y  )
    (x-1,y-1) (x  , y-1) (x+1,y-1)
    

  • And finally, Smooth_Noise make an average with the results of Noise, which generate coherents number by a wreid way that I don't understand much.

• Show results: Once the map is made, it only rest how to show the results to a human. By that I make 5 ways, make a photo with pillow in black and white or in colors, or get a graph in 2D, 3D or both with matplotlib. To make the code smaller I will only left here: the 2D + 3D from matplotlib, and the photo from pillow (Black and white, and color).

Code

import math, random                         # Map Generation
import matplotlib.pyplot as plt             # 2D and 3D graph
from mpl_toolkits.mplot3d import Axes3D     # 3D graph
import numpy as np                          # 2D and 3D graph
from PIL import Image                       # Picture

class D():
    def Cubic_Interpolate(self, v0, v1, v2, v3, x):
        P = (v3 - v2) - (v0 - v1)
        Q = (v0 - v1) - P
        R = v2 - v0
        S = v1
        return P * x**3 + Q * x**2 + R * x + S

    #def Linear_Interpolate(self, a, b, x):
    #    ''' Strongly not recomend. '''
    #    return a * (1 - x) + b * x

    #def Cosine_Interpolate(self, a, b, x):
    #    ''' Faster but ugly. '''
    #    ft = x * math.pi # 3.1415927
    #    f = (1 - math.cos(ft)) * 0.5
    #    return a * (1 - f) + b * f

class D2(D):
    def __init__(self, lenght, octaves = 1):

        self.lenght_axes = round(lenght ** 0.5)
        self.lenght = self.lenght_axes ** 2                
        self.result, self.line_result = self.Perlin(self.lenght_axes, octaves)

    def Noise(self, x, y):
        n = x + y * 57
        n = (n<<13) ^ n
        return ( 1.0 - ( (n * (n * n * 15731 + 789221) + 1376312589) & 0x7fffffff) / 1073741824.0)

    def Smooth_Noise(self, x, y, smooth = 5 ):    # I plant to make a re-work here, so if possible I don't want a modification in this function in order to answer this question.
        corners = (self.Noise(x - 1, y - 1) + self.Noise(x + 1, y - 1) + self.Noise(x - 1, y + 1) + self.Noise(x + 1, y + 1) ) / 16
        sides   = (self.Noise(x - 1, y) + self.Noise(x + 1, y) + self.Noise(x, y - 1)  + self.Noise(x, y + 1) ) /  8
        center  =  self.Noise(x, y) / 4
        return corners + sides + center

    def Interpolate_Noise(self, x, y):

        round_x = math.floor(x)
        frac_x  = x - round_x

        round_y = math.floor(y)
        frac_y  = y - round_y

        v11 = self.Smooth_Noise(round_x - 1, round_y - 1)
        v12 = self.Smooth_Noise(round_x    , round_y - 1)
        v13 = self.Smooth_Noise(round_x + 1, round_y - 1)
        v14 = self.Smooth_Noise(round_x + 2, round_y - 1)
        i1 = self.Cubic_Interpolate(v11, v12, v13, v14, frac_x)

        v21 = self.Smooth_Noise(round_x - 1, round_y)
        v22 = self.Smooth_Noise(round_x    , round_y)
        v23 = self.Smooth_Noise(round_x + 1, round_y)
        v24 = self.Smooth_Noise(round_x + 2, round_y)
        i2 = self.Cubic_Interpolate(v21, v22, v23, v24, frac_x)

        v31 = self.Smooth_Noise(round_x - 1, round_y + 1)
        v32 = self.Smooth_Noise(round_x    , round_y + 1)
        v33 = self.Smooth_Noise(round_x + 1, round_y + 1)
        v34 = self.Smooth_Noise(round_x + 2, round_y + 1)
        i3 = self.Cubic_Interpolate(v31, v32, v33, v34, frac_x)

        v41 = self.Smooth_Noise(round_x - 1, round_y + 2)
        v42 = self.Smooth_Noise(round_x    , round_y + 2)
        v43 = self.Smooth_Noise(round_x + 1, round_y + 2)
        v44 = self.Smooth_Noise(round_x + 2, round_y + 2)
        i4 = self.Cubic_Interpolate(v41, v42, v43, v44, frac_x)

        return self.Cubic_Interpolate(i1, i2, i3, i4, frac_y)

    #def Interpolate_Noise(self, x, y):
    #    ''' In case you want linear or cosine interpolation. '''
    #
    #    Interpolation = self.Linear_Interpolate or self.Cosine_Interpolate
    #
    #    round_x = math.floor(x)
    #    frac_x  = x - round_x
    #
    #    round_y = math.floor(y)
    #    frac_y  = y - round_y
    #
    #    a1 = self.Smooth_Noise(round_x    , round_y)
    #    b1 = self.Smooth_Noise(round_x + 1, round_y)
    #    a2 = self.Smooth_Noise(round_x    , round_y + 1)
    #    b2 = self.Smooth_Noise(round_x + 1, round_y + 1)
    # 
    #    i1 = Interpolation(a1, b1, frac_x)
    #    i2 = Interpolation(a2, b2, frac_x)
    # 
    #    return self.Cubic_Interpolate(i1, i2, frac_x)

    def Perlin(self, lenght_axes, octaves, zoom = 0.01, amplitude_base = 0.5):
        result = []
        line_result = []

        for y in range(lenght_axes):
            line = []
            for x in range(lenght_axes):
                value = 0
                for o in range(octaves):
                    frequency = 2 ** o
                    amplitude = amplitude_base ** o
                    value += self.Interpolate_Noise(x * frequency * zoom, y * frequency * zoom) * amplitude
                line.append(value)
                line_result.append(value)
            result.append(line)
            print(f"{y:5} / {lenght_axes} ({y/lenght_axes*100:.2f}%): {round(y/lenght_axes*20) * '#'} {(20-round(y/lenght_axes*20)) * ' '}. Remaining {lenght_axes-y}.")
        return result, line_result             

    def graph_2d(self, color = 'viridis'):
        # Other colors: https://matplotlib.org/examples/color/colormaps_reference.html
        Z = np.array(self.result)     
        fig, ax1 = plt.subplots()
        pos = ax1.imshow(Z, cmap=color, interpolation='none')
        fig.colorbar(pos)
        plt.show()

    def graph_3d(self, color = 'viridis'):
        # Other colors: https://matplotlib.org/examples/color/colormaps_reference.html
        X = np.arange(self.lenght_axes)
        Y = np.arange(self.lenght_axes)
        X, Y = np.meshgrid(X, Y)        
        Z = np.array(self.result)
        fig = plt.figure()
        ax = Axes3D(fig)
        ax1 = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=color)
        fig.colorbar(ax1, shrink=0.5, aspect=5)
        plt.show()

    def graph(self, color = 'viridis'):
        # Other colors: https://matplotlib.org/examples/color/colormaps_reference.html
        fig = plt.figure()
        Z = np.array(self.result)

        vmin = np.amin(Z)
        vmax = np.amax(Z)

        ax = fig.add_subplot(1, 2, 1, projection='3d')
        X = np.arange(self.lenght_axes)
        Y = np.arange(self.lenght_axes)
        X, Y = np.meshgrid(X, Y)        
        d3 = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=color, linewidth=0, antialiased=False, vmin=vmin, vmax=vmax)
        fig.colorbar(d3)

        ax = fig.add_subplot(1, 2, 2)
        d2 = ax.imshow(Z, cmap=color, interpolation='none', vmin=vmin, vmax=vmax)
        fig.colorbar(d2)

        plt.show()

    def image(self, color = True):
        line = []
        vmax = max(self.line_result)
        if not color:
            vmin = abs(min(self.line_result))
        for v in self.line_result:
            r = g = b = 0
            if color:
                value = v / vmax * 255
                if value > 0:                    
                    b = value
                else:
                    r = abs(value)
            else:
                r = g = b = 127 + (v / vmax * 128)

            line.append((round(r), round(g), round(b)))  
        img = Image.new('RGB', (self.lenght_axes, self.lenght_axes))
        img.putdata(line)
        img.save('chart1.png')
        img.show()

    test = D2(total_titles, octaves)
    test.image(color_bolean)     # Photo
    test.graph()     # Graph

I tried to make the code faster, but my attemps failed.

• In my first attempt I made a thread for each iteration of x and y in Perlin, but the code result even slower. (I deleted the code)

• My second attempt was made a thread for each iteration of only y in Perlin, but the code was about only 0.09% faster but much more complex (I use time.time() to check it). It was:

  def Perlin(self, lenght_axes, octaves, zoom = 0.01, amplitude_base = 0.5):
  
      lines_queue = Queue(maxsize=0)
      _queue = Queue(maxsize=0)
      self.threads = []
      results = Queue(maxsize=0)
  
      def do_line(y, lenght_axes = lenght_axes, octaves = octaves, amplitude_base = amplitude_base, zoom = zoom):
          line = []
          for x in range(lenght_axes):
              value = 0
              for o in range(octaves):
                  frequency = 2 ** o
                  amplitude = amplitude_base ** o
                  value += self.Interpolate_Noise(x * frequency * zoom, y * frequency * zoom) * amplitude
              line.append(value)                
          lines_queue.put((y, line))
          _queue.put(1)
  
      def manage(lenght_axes = lenght_axes):
          result = [0] * lenght_axes
          works = lenght_axes
          while True:
              if not lines_queue.empty():
                  y, values = lines_queue.get()
                  result[y] = values
              if not _queue.empty():
                  works -= _queue.get()
              print(works)
              if works == 0:
                  break
          line_result = []
          for line in result:
              line_result.extend(line)
  
          results.put(result)
          results.put(line_result)
  
      managing = threading.Thread(target = manage)
      managing.start()
  
      for y in range(lenght_axes):
          self.threads.append(threading.Thread(target = do_line, args = (y, lenght_axes)))
          self.threads[-1].start()
  
      managing.join()
  
      return results.get(), results.get() 
  

CProfile

It's the first time I use it, So I am not sure how to read it, I found this interesting (I crop it with paint).

Profile

It's the first time I use it, So I am not sure how to read it, I found this interesting (I crop it with paint).

Answer 1

Since you are already mentioning numpy for displaying, you should use it for the calculation as well. It excels in fast numeric calculations (by being efficiently implemented in C).

All you need to change is an import at the top, two lines in Interpolate_Noise and your Perlin function, the other functions work with numpy.arrays without any modification:

import numpy as np

...

def Interpolate_Noise(x, y):

    round_x = np.array(np.floor(x), dtype=int)
    frac_x = x - round_x

    round_y = np.array(np.floor(y), dtype=int)
    frac_y = y - round_y
    ...


def Perlin(length_axes, octaves, zoom=0.01, amplitude_base=0.5):
    x = np.arange(length_axes)
    y = np.arange(length_axes)
    o = np.arange(octaves)
    xx, yy, oo = np.meshgrid(x, y, o)

    frequency = 2 ** oo * zoom
    amplitude = amplitude_base ** oo
    result = (Interpolate_Noise(xx * frequency, yy * frequency)
              * amplitude).sum(axis=-1)

    return result, result.flatten()

This gives a huge speed-boost, as can be seen from this graph, where I scanned length_axes from 2 to 20, for octaves = 50 (commenting out your progress printing code, to avoid bias due to printing).

The output is the same as for your code, within numerical precision (differences are on the order of 10^-17).

It relies on numpy.meshgrid to give matrices with all combinations of x, y, and o.

The overall algorithmic complexity is still the same, as can be seen in this second graph, where I continued increasing length_axes to 100 for the numpy function. It just takes longer until the execution time gets large:

---

Note that I removed all traces of this being a class. This is because you don't need it to be a class, as far as I can tell.

Here are some additional notes about your code:

1. D, D2 and so on are bad class names. Try something more descriptive like NoiseGenerator2D or so.
2. Python has an official style-guide, PEP8, which recommends using lower_case for functions and variables, PascalCase for classes and UPPER_CASE for constants.
3. You should avoid magic numbers. Put the numbers in Noise into constants, defined at the module level.