ENSEEIHT/projet-probleme-inverse-3D
1
0
Fork 0
Code Issues Pull requests Packages Projects Releases Wiki Activity

Merge branch 'bordel_laurent' of git.inpt.fr:tocard-inc/enseeiht/projet-be into bordel_laurent

Browse source
This commit is contained in:
Laureηt 2023-01-26 20:05:39 +01:00
parent bbeeee6aee d8591efbae
commit 2ae473a89f
Signed by: Laurent
SSH key fingerprint: SHA256:kZEpW8cMJ54PDeCvOhzreNr4FSh6R13CMGH/POoO8DI
1 changed files with 64 additions and 5 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

69
src/levelset.py
View file

@ -5,12 +5,71 @@ import numpy as np
import perlin_noise
from rich.progress import track
# 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)