From d8591efbaea4e2d251802ae17c2a11f8ba363695 Mon Sep 17 00:00:00 2001 From: gdamms Date: Thu, 26 Jan 2023 11:15:52 +0100 Subject: [PATCH] levelset 3d --- src/levelset.py | 69 +++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 64 insertions(+), 5 deletions(-) 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)