diff --git a/src/levelset.py b/src/levelset.py
index 3345ef2..9d0dcc3 100644
--- a/src/levelset.py
+++ b/src/levelset.py
@@ -5,12 +5,71 @@ import imageio.v2 as imageio
 import perlin_noise
 
 
-# X, Y, Z = np.mgrid[:100, :100, :100]
-# V = np.sqrt((X - 50)**2 + (Y - 50)**2 + (Z - 50)**2)
+def generate_perlin_noise_3d(shape, res):
+    def f(t):
+        return 6*t**5 - 15*t**4 + 10*t**3
+
+    delta = (res[0] / shape[0], res[1] / shape[1], res[2] / shape[2])
+    d = (shape[0] // res[0], shape[1] // res[1], shape[2] // res[2])
+    grid = np.mgrid[0:res[0]:delta[0],0:res[1]:delta[1],0:res[2]:delta[2]]
+    grid = grid.transpose(1, 2, 3, 0) % 1
+    # Gradients
+    theta = 2*np.pi*np.random.rand(res[0]+1, res[1]+1, res[2]+1)
+    phi = 2*np.pi*np.random.rand(res[0]+1, res[1]+1, res[2]+1)
+    gradients = np.stack((np.sin(phi)*np.cos(theta), np.sin(phi)*np.sin(theta), np.cos(phi)), axis=3)
+    gradients[-1] = gradients[0]
+    g000 = gradients[0:-1,0:-1,0:-1].repeat(d[0], 0).repeat(d[1], 1).repeat(d[2], 2)
+    g100 = gradients[1:  ,0:-1,0:-1].repeat(d[0], 0).repeat(d[1], 1).repeat(d[2], 2)
+    g010 = gradients[0:-1,1:  ,0:-1].repeat(d[0], 0).repeat(d[1], 1).repeat(d[2], 2)
+    g110 = gradients[1:  ,1:  ,0:-1].repeat(d[0], 0).repeat(d[1], 1).repeat(d[2], 2)
+    g001 = gradients[0:-1,0:-1,1:  ].repeat(d[0], 0).repeat(d[1], 1).repeat(d[2], 2)
+    g101 = gradients[1:  ,0:-1,1:  ].repeat(d[0], 0).repeat(d[1], 1).repeat(d[2], 2)
+    g011 = gradients[0:-1,1:  ,1:  ].repeat(d[0], 0).repeat(d[1], 1).repeat(d[2], 2)
+    g111 = gradients[1:  ,1:  ,1:  ].repeat(d[0], 0).repeat(d[1], 1).repeat(d[2], 2)
+    # Ramps
+    n000 = np.sum(np.stack((grid[:,:,:,0]  , grid[:,:,:,1]  , grid[:,:,:,2]  ), axis=3) * g000, 3)
+    n100 = np.sum(np.stack((grid[:,:,:,0]-1, grid[:,:,:,1]  , grid[:,:,:,2]  ), axis=3) * g100, 3)
+    n010 = np.sum(np.stack((grid[:,:,:,0]  , grid[:,:,:,1]-1, grid[:,:,:,2]  ), axis=3) * g010, 3)
+    n110 = np.sum(np.stack((grid[:,:,:,0]-1, grid[:,:,:,1]-1, grid[:,:,:,2]  ), axis=3) * g110, 3)
+    n001 = np.sum(np.stack((grid[:,:,:,0]  , grid[:,:,:,1]  , grid[:,:,:,2]-1), axis=3) * g001, 3)
+    n101 = np.sum(np.stack((grid[:,:,:,0]-1, grid[:,:,:,1]  , grid[:,:,:,2]-1), axis=3) * g101, 3)
+    n011 = np.sum(np.stack((grid[:,:,:,0]  , grid[:,:,:,1]-1, grid[:,:,:,2]-1), axis=3) * g011, 3)
+    n111 = np.sum(np.stack((grid[:,:,:,0]-1, grid[:,:,:,1]-1, grid[:,:,:,2]-1), axis=3) * g111, 3)
+    # Interpolation
+    t = f(grid)
+    n00 = n000*(1-t[:,:,:,0]) + t[:,:,:,0]*n100
+    n10 = n010*(1-t[:,:,:,0]) + t[:,:,:,0]*n110
+    n01 = n001*(1-t[:,:,:,0]) + t[:,:,:,0]*n101
+    n11 = n011*(1-t[:,:,:,0]) + t[:,:,:,0]*n111
+    n0 = (1-t[:,:,:,1])*n00 + t[:,:,:,1]*n10
+    n1 = (1-t[:,:,:,1])*n01 + t[:,:,:,1]*n11
+    return ((1-t[:,:,:,2])*n0 + t[:,:,:,2]*n1)
+
+V = 10 * generate_perlin_noise_3d((100, 100, 100), (10, 10, 10))
+
+X, Y, Z = np.mgrid[:100, :100, :100]
+V += np.sqrt((X-50)**2 + (Y-50)**2 + (Z-50)**2)
+V = (V - V.min()) / (V.max() - V.min())
+
+frame_list = []
+
+fig = plt.figure()
+ax = fig.add_subplot(111, projection='3d')
+for i in np.linspace(0.05, 0.5, 100):
+    ax.clear()
+    vertices, triangles = mcubes.marching_cubes(V, i)
+    ax.plot_trisurf(vertices[:,0], vertices[:,1], vertices[:,2], triangles=triangles)
+    ax.set_xlim(0, 100)
+    ax.set_ylim(0, 100)
+    ax.set_zlim(0, 100)
+    ax.set_xticklabels([])
+    ax.set_yticklabels([])
+    ax.set_zticklabels([])
+    plt.savefig(f"/tmp/frame.png", bbox_inches='tight', pad_inches=0, dpi=300, transparent=True)
+    frame_list.append(imageio.imread(f"/tmp/frame.png"))
+
+imageio.mimsave('picture.gif', frame_list + frame_list[::-1], fps=60)
 
-# for i in range(5, 45):
-#     vertices, triangles = mcubes.marching_cubes(V, i)
-#     mcubes.export_obj(vertices, triangles, f"cube_{i}.obj")
 
 noise = perlin_noise.PerlinNoise(octaves=6, seed=1)