How can I make a Worley Noise algorithm faster?

Viewed 482

I wrote an algorithm, that produces 2D Worley Noise in python. But it takes several seconds for a resolution with 400x400 Pixels. I tried to subdivide the space like here, but I am not able to implement it in Python. I also tried a different method to calculate the nearest point, but it did not work faster... How can I make it work faster?

class Point2D:
        def __init__(self, x, y):
            self.p = [x,y]
        
        def x(self):
            return(self.p[0])
        def y(self):
            return(self.p[1])
        def distance(point1, point2):
            import math
            return(math.sqrt((point1.x()-point2.x())**2 + (point1.y()-point2.y())**2))
        def getDistances(origin, li : list):
            distances = []
            for ll in li:
                distances.append(Point2D.distance(origin, Point2D(ll[0], ll[1])))
            return(distances)
        
class WorleyNoise:
      
    def __init__(self, height, width, density):
        
        self.height = height
        self.width = width
        self.density = density
    def auto(self, option):

        self.generatePoints()
        self.calculateNoise(option)
        self.showNoise()
    def generatePoints(self):
        import numpy as np
        self.points = []
        for _ in range(self.density):
            self.points.append([np.random.randint(0, self.width,1)[0], np.random.randint(0, self.height,1)[0]])
    def showPoints(self):
        import matplotlib.pyplot as plt
        plt.scatter([self.points[i][0] for i in range(len(self.points))], [self.points[l][1] for l in range(len(self.points))])
        plt.show()
    def calculateNoise(self, option):
        import time
        self.data = [[0] * self.width for _ in range(self.height)]
        
        for h in range(self.height):
            start = time.time()
            for w in range(self.width):
                self.distances = Point2D.getDistances(Point2D(w, h), self.points) 
                self.distances.sort()
                self.data[h][w] = self.distances[option]
            print(round(h/(self.height)*100), "%", f" {(time.time()-start) * (self.height - h)} Seconds")
        
    def showNoise(self):
        import matplotlib.pyplot as plt
        plt.imshow(self.data, cmap = "gray")
        plt.show()
        
w = WorleyNoise(height = 200, width = 200, density = 40)
w.auto(0)

Edit:

That‘s the fastest solution I have found until now:


def worleynoise(width, height, density):
    from numpy import random, mgrid, dstack, empty
    from scipy.spatial import cKDTree

    points = [[random.randint(0, height), random.randint(0, width)] for _ in range(density)]  # Generates Points(y, x)
   
    coord = dstack(mgrid[0:height, 0:width])  # Makes array with coordinates as values

    tree = cKDTree(points)  # Build Tree
    distances = tree.query(coord, workers=-1)[0]  # Calculate distances (workers=-1: Uses all CPU Cores)
    return distances


if __name__ == "__main__":
    from matplotlib.pyplot import imshow, show
    import time
    start = time.time()
    w = worleynoise(1000, 1000, 225)
    print(time.time() - start)
    imshow(w)
    show()

The difference to brute force becomes more apparent with a higher amount of points.

4 Answers

Looping over width and height is inefficient if you can use numpy to broadcast your operations, here is your modified code where a new broadcastCalculateNoise is used with a benchmarking of both old and new methods:

import time
import numpy as np
import matplotlib.pyplot as plt


class Point2D:
    def __init__(self, x, y):
        self.p = [x, y]

    def x(self):
        return (self.p[0])

    def y(self):
        return (self.p[1])

    def distance(point1, point2):
        import math
        return (math.sqrt((point1.x() - point2.x()) ** 2 + (point1.y() - point2.y()) ** 2))

    def getDistances(origin, li: list):
        distances = []
        for ll in li:
            distances.append(Point2D.distance(origin, Point2D(ll[0], ll[1])))
        return (distances)


class WorleyNoise:

    def __init__(self, height, width, density, use_broadcast_ops):

        self.height = height
        self.width = width
        self.density = density
        self.use_broadcast_ops = use_broadcast_ops

    def auto(self, option):

        self.generatePoints()
        start = time.time()
        if self.use_broadcast_ops:
            self.broadcastCalculateNoise(option)
        else:
            self.calculateNoise(option)
        end = time.time()
        print("total time : " + str(end - start) + " seconds")
        self.showNoise()

    def generatePoints(self):
        self.points = []
        for _ in range(self.density):
            self.points.append([np.random.randint(0, self.width, 1)[0], np.random.randint(0, self.height, 1)[0]])

    def showPoints(self):
        plt.scatter([self.points[i][0] for i in range(len(self.points))],
                    [self.points[l][1] for l in range(len(self.points))])
        plt.show()

    def calculateNoise(self, option):
        self.data = [[0] * self.width for _ in range(self.height)]
        for h in range(self.height):
            for w in range(self.width):
                self.distances = Point2D.getDistances(Point2D(w, h), self.points)
                self.distances.sort()
                self.data[h][w] = self.distances[option]

    def broadcastCalculateNoise(self, option):
        # casting points to numpy, it is of shape (nb_point, 2)
        points = np.array(self.points)
        # simple array of x and y coordinates for each coordinate
        xs = np.arange(self.width)
        ys = np.arange(self.height)
        # use the previously computed xs to get point.x - x for each x
        # notice the use of np.newaxis to control the broadcasting of the result to
        # an array of shape (nb_point, width)
        x_dist = np.power(points[:, 0, np.newaxis] - xs, 2)
        # same for ys, giving a (nb_point, height) shaped array
        y_dist = np.power(points[:, 1, np.newaxis] - ys, 2)
        # use the two last array to compute distance : sqrt((p.x - x) ** 2 + (.y - y ) ** 2))
        d = np.sqrt(x_dist[:, :, np.newaxis] + y_dist[:, np.newaxis, :])
        # d is of shape (nb_point, width, height), but we must sort it along the first axis
        distances = np.sort(d, axis=0)
        self.data = distances[option]

    def showNoise(self):
        import matplotlib.pyplot as plt
        plt.imshow(self.data, cmap="gray")
        plt.show()


print("Using old fashioned loops :")
w = WorleyNoise(height=200, width=200, density=40, use_broadcast_ops=False)
w.auto(0)
print("Harnessing numpy power :")
w = WorleyNoise(height=200, width=200, density=40, use_broadcast_ops=True)
w.auto(0)

results in :

Using old fashioned loops :
total time : 2.945978879928589 seconds
Harnessing numpy power :
total time : 0.03304648399353027 seconds

Process finished with exit code 0

That‘s the fastest solution I found until now:


def worleynoise(width, height, density):
    from numpy import random, mgrid, dstack, empty
    from scipy.spatial import cKDTree

    points = [[random.randint(0, height), random.randint(0, width)] for _ in range(density)]  # Generates Points(y, x)
   
    coord = dstack(mgrid[0:height, 0:width])  # Makes array with coordinates as values

    tree = cKDTree(points)  # Build Tree
    distances = tree.query(coord, workers=-1)[0]  # Calculate distances (workers=-1: Uses all CPU Cores)
    return distances


if __name__ == "__main__":
    from matplotlib.pyplot import imshow, show
    import time
    start = time.time()
    w = worleynoise(1000, 1000, 225)
    print(time.time() - start)
    imshow(w)
    show()

The difference to brute force becomes more apparent with a higher amount of points.

The accepted answer shows off "numpy power", something I'm unfamiliar with. However, this process can be sped up with the power of standard Python instead/as well. Here is what I've been able to come up with, by removing conditional imports, list initializations (comprehensions are used), numpy.random.randint (random.randint seems faster), and a semi-unnecessary class. This isn't beautiful, but it also isn't ugly and it saves a lot of time.

import random
import math
import matplotlib.pyplot as plt
        
class WorleyNoise:     
    def __init__(self, height, width, npoints):       
        self.height = height
        self.width = width
        self.npoints = npoints
    def generatePoints(self):
        self.points = [(random.randint(0, self.width), random.randint(0, self.height)) for x in range(self.npoints)]
    def showPoints(self):
        plt.scatter([self.points[x][0] for x in range(self.npoints)], [self.points[x][1] for x in range(self.npoints)])
        plt.show()
    def calculateNoise(self, option):
        self.data = [[sorted([math.sqrt((p[0]-x)**2 + (p[1]-y)**2) for p in self.points])[option] for x in range(self.width)] for y in range(self.height)]
    def showNoise(self):
        plt.imshow(self.data, cmap = "gray")
        plt.show()
        
w = WorleyNoise(200, 200, 40)
w.generatePoints()
w.calculateNoise(0)
w.showNoise()

By my calculations (CPython, Unix, x86_64), the original (minus the time checking, printing, and matplotlib display) takes 6.037s, and the modified version (minus the matplotlib display) takes 2.556 seconds.

Apart from the above answers, there is another great tool( Numba ) to make python code faster by compiling complete code before executing it. It makes python act as a compiled language rather than as an interpreter.

Related