diff --git a/Generator/colour.py b/Generator/colour.py new file mode 100644 index 00000000..3843acfa --- /dev/null +++ b/Generator/colour.py @@ -0,0 +1,330 @@ +import numpy as np + + +biome_deep_water = { + "maxA": 0, + "minA": None, + "A": -6000, + "vA": 3000, + "T": 40, + "vT": 30, + "H": None, + "vH": None, + "r": 16, + "g": 25, + "b": 36, + "ref": 0.6 +} +biome_median_water = { + "maxA": 0, + "minA": None, + "A": -1000, + "vA": 2000, + "T": 40, + "vT": 30, + "H": None, + "vH": None, + "r": 16, + "g": 25, + "b": 36, + "ref": 0.6 +} +biome_shallow_water = { + "maxA": 0, + "minA": None, + "A": -100, + "vA": 300, + "T": 40, + "vT": 30, + "H": None, + "vH": None, + "r": 50, + "g": 79, + "b": 114, + "ref": 0.6 +} + + +biome_ice = { + "maxA": None, + "minA": 0, + "A": None, + "vA": None, + "T": -10, + "vT": 5, + "H": None, + "vH": None, + "r": 255, + "g": 255, + "b": 255, + "ref": 0.6 +} +biome_water_ice = { + "maxA": 0, + "minA": None, + "A": -500, + "vA": 500, + "T": -15, + "vT": 5, + "H": None, + "vH": None, + "r": 200, + "g": 200, + "b": 255, + "ref": 0.6 +} + +biome_cold_ice = { + "maxA": None, + "minA": 0, + "A": None, + "vA": None, + "T": -40, + "vT": 20, + "H": None, + "vH": None, + "r": 255, + "g": 255, + "b": 255, + "ref": 0.6 +} +biome_cold_water_ice = { + "maxA": 0, + "minA": None, + "A": None, + "vA": None, + "T": -45, + "vT": 20, + "H": None, + "vH": None, + "r": 200, + "g": 200, + "b": 255, + "ref": 0.6 +} + +biome_tundra = { + "maxA": None, + "minA": 0, + "A": 0, + "vA": 2000, + "T": -5, + "vT": 5, + "H": 0.5, + "vH": 2, + "r": 93, + "g": 80, + "b": 35, + "ref": 0 +} + +biome_boreal_forest = { + "maxA": None, + "minA": 0, + "A": 500, + "vA": 1500, + "T": 5, + "vT": 5, + "H": 2, + "vH": 2, + "r": 12, + "g": 63, + "b": 7, + "ref": 0 +} + +biome_forest = { + "maxA": None, + "minA": 0, + "A": 500, + "vA": 1500, + "T": 15, + "vT": 5, + "H": 2.5, + "vH": 3, + "r": 67, + "g": 88, + "b": 29, + "ref": 0 +} + +biome_cold_desert = { + "maxA": None, + "minA": 0, + "A": 500, + "vA": 1500, + "T": 10, + "vT": 7, + "H": 0, + "vH": 1, + "r": 245, + "g": 200, + "b": 151, + "ref": 0 +} +biome_warm_desert = { + "maxA": None, + "minA": 0, + "A": 500, + "vA": 1500, + "T": 25, + "vT": 7, + "H": 0, + "vH": 1, + "r": 245, + "g": 200, + "b": 151, + "ref": 0 +} + +biome_savanah = { + "maxA": None, + "minA": 0, + "A": 0, + "vA": 2000, + "T": 25, + "vT": 5, + "H": 1, + "vH": 1.5, + "r": 66, + "g": 105, + "b": 13, + "ref": 0 +} +biome_humid_savanah = { + "maxA": None, + "minA": 0, + "A": 0, + "vA": 2000, + "T": 25, + "vT": 5, + "H": 2, + "vH": 1.5, + "r": 53, + "g": 112, + "b": 46, + "ref": 0 +} +biome_jungle = { + "maxA": None, + "minA": 0, + "A": 500, + "vA": 1000, + "T": 25, + "vT": 5, + "H": 3, + "vH": 1.5, + "r": 5, + "g": 50, + "b": 5, + "ref": 0 +} + +biome_alpine_tundra = { + "maxA": None, + "minA": 0, + "A": 3000, + "vA": 1500, + "T": 0, + "vT": 5, + "H": 0.5, + "vH": 3, + "r": 131, + "g": 122, + "b": 75, + "ref": 0 +} +biome_high_moutain = { + "maxA": None, + "minA": 0, + "A": 6000, + "vA": 2000, + "T": None, + "vT": None, + "H": None, + "vH": None, + "r": 114, + "g": 97, + "b": 71, + "ref": 0 +} + +biome_list = [] + +biome_list.append(biome_deep_water) +biome_list.append(biome_median_water) +biome_list.append(biome_shallow_water) +biome_list.append(biome_ice) +biome_list.append(biome_cold_ice) +biome_list.append(biome_water_ice) +biome_list.append(biome_cold_water_ice) +biome_list.append(biome_tundra) +biome_list.append(biome_boreal_forest) +biome_list.append(biome_forest) +biome_list.append(biome_cold_desert) +biome_list.append(biome_warm_desert) +biome_list.append(biome_savanah) +biome_list.append(biome_humid_savanah) +biome_list.append(biome_jungle) +biome_list.append(biome_alpine_tundra) +biome_list.append(biome_high_moutain) + +T_eq = 30 +T_pole = -40 + + +def humidity(theta): + theta_deg = 90 * theta / (np.pi / 2) + + c1 = (5 / (1 + (theta_deg / 10)**2)) + c2 = 2.5 / (1 + ((theta_deg - 50) / 10)**2) + c3 = 2.5 / (1 + ((theta_deg + 50) / 10)**2) + + return 4 * (c1 + c2 + c3 - 0.5) / 4.692 + + +def colorize(x, y, z, c, vT, vH): + A = np.zeros_like(c) + A[c > 0.5] = ((c[c > 0.5] - 0.5) / 0.5) * 8000 # land + A[c < 0.5] = ((c[c < 0.5] - 0.5) / 0.5) * 12000 # ocean + + theta = np.arcsin(y / np.sqrt(x**2 + y**2 + z**2)) + + T = T_pole + np.cos(theta) * (T_eq - T_pole) + (vT - 0.5) * 20 + T[c > 0.5] -= 7 * A[c > 0.5] / 1000 + + H = humidity(theta) * (1 + (vH - 0.5)) + + rgbref = np.zeros((len(c), 4)) + pond_tot = np.zeros_like(c) + + for biome in biome_list: + ok = np.ones_like(pond_tot, dtype=bool)#np.bool) + + if biome['minA'] is not None: + ok = ok & (A >= biome['minA']) + if biome['maxA'] is not None: + ok = ok & (A <= biome['maxA']) + + pond = np.ones_like(pond_tot) + + if biome['A'] is not None: + pond[ok] *= np.exp(-((A[ok] - biome['A']) / biome['vA'])**2) + if biome['H'] is not None: + pond[ok] *= np.exp(-((H[ok] - biome['H']) / biome['vH'])**2) + if biome['T'] is not None: + pond[ok] *= np.exp(-((T[ok] - biome['T']) / biome['vT'])**2) + + rgbref[ok, 0] += pond[ok] * biome['r'] + rgbref[ok, 1] += pond[ok] * biome['g'] + rgbref[ok, 2] += pond[ok] * biome['b'] + rgbref[ok, 3] += pond[ok] * biome['ref'] + + pond_tot[ok] += pond[ok] + + rgbref[:, 0] /= (255 * pond_tot) + rgbref[:, 1] /= (255 * pond_tot) + rgbref[:, 2] /= (255 * pond_tot) + rgbref[:, 3] /= (pond_tot) + + return rgbref \ No newline at end of file diff --git a/Generator/examples/anim_example.gif b/Generator/examples/anim_example.gif new file mode 100644 index 00000000..9f534360 Binary files /dev/null and b/Generator/examples/anim_example.gif differ diff --git a/Generator/examples/out0.png b/Generator/examples/out0.png new file mode 100644 index 00000000..53f0f10a Binary files /dev/null and b/Generator/examples/out0.png differ diff --git a/Generator/examples/out1.png b/Generator/examples/out1.png new file mode 100644 index 00000000..442cbb82 Binary files /dev/null and b/Generator/examples/out1.png differ diff --git a/Generator/examples/out2.png b/Generator/examples/out2.png new file mode 100644 index 00000000..810db674 Binary files /dev/null and b/Generator/examples/out2.png differ diff --git a/Generator/geometry.py b/Generator/geometry.py new file mode 100644 index 00000000..1b01d5c8 --- /dev/null +++ b/Generator/geometry.py @@ -0,0 +1,60 @@ +import numpy as np + + +def dot(a, b): + x1, y1, z1 = a + x2, y2, z2 = b + return x1 * x2 + y1 * y2 + z1 * z2 + + +def cos_vec(a, b): + return dot(a, b) / np.sqrt(dot(a, a) * dot(b, b)) + + +def cross(a, b): + x1, y1, z1 = a + x2, y2, z2 = b + return (y1 * z2 - z1 * y2, z1 * x2 - x1 * z2, x1 * y2 - y1 * x2) + + +def reflection(d, n): + xd, yd, zd = d + xn, yn, zn = n + + u = dot(d, n) / dot(n, n) + + return (xd - 2 * u * xn, yd - 2 * u * yn, zd - 2 * u * zn) + + +def normal(x, y, z, h, dx, dy, dz, A=0.01): + ''' + x = R cos phi cos theta + y = R sin phi cos theta + z = R sin theta + ''' + + R = np.sqrt(x * x + y * y + z * z) + costheta = np.sqrt(x * x + y * y) / R + sintheta = z / R + cosphi = x / (R * costheta) + sinphi = y / (R * costheta) + ''' + r = R (1 + A*h) + ''' + r = R * (1 + A * h) + + drdphi = R * A * (dx * -y + dy * x) + drdtheta = R * A * (dx * (-z * cosphi) + dy * + (-z * cosphi) + dz * R * costheta) + + uphi = (r * -sinphi * costheta + drdphi * cosphi * costheta, + r * cosphi * costheta + drdphi * sinphi * costheta, 0) + + utheta = (r * cosphi * -sintheta + drdtheta * cosphi * costheta, + r * sinphi * -sintheta + drdtheta * sinphi * costheta, + r * costheta + drdtheta * sintheta) + + v = cross(uphi, utheta) + nv = dot(v, v) + v = v[0] / nv, v[1] / nv, v[2] / nv + return v \ No newline at end of file diff --git a/Generator/main.py b/Generator/main.py new file mode 100644 index 00000000..0b329cc9 --- /dev/null +++ b/Generator/main.py @@ -0,0 +1,198 @@ +import numpy as np +import matplotlib.pyplot as plt + +from tqdm import tqdm +from noise import double_perlin3D, background + + +from geometry import normal, cos_vec, reflection +from colour import colorize + +res = 2048 +radius = 0.4 + +multi_sampling = 1 + +eps = 1 / (100 * res) + +atmos_size = 0.02 +atmos_colour = (0.6, 0.851, 0.912) + +x, y = np.meshgrid(np.linspace(-0.5, 0.5, res), np.linspace(-0.5, 0.5, res)) + +dx, dy = np.meshgrid( + np.linspace(0, 1 / res, multi_sampling), + np.linspace(0, 1 / res, multi_sampling)) + + +def gen_rgbbase_and_normal(seed=69, rot=0): + im_rgbref_base = np.zeros((res, res, multi_sampling, multi_sampling, 4)) + normal_base = np.zeros((res, res, multi_sampling, multi_sampling, 3)) + cloud_map = np.zeros((res, res, multi_sampling, multi_sampling)) + + for i in tqdm(range(res)): + for d1 in range(multi_sampling): + for d2 in range(multi_sampling): + + z2 = radius**2 - (x[i, :] + dx[d1, d2])**2 - ( + y[i, :] + dy[d1, d2])**2 + ok = z2 >= 0 + + if np.any(ok): + x1, y1, z1 = (x[i, ok] + dx[d1, d2]), ( + y[i, ok] + dy[d1, d2]), np.sqrt(z2[ok]) + + ox1, oz1 = x1, z1 + + x1, z1, = x1 * np.cos(rot) - z1 * np.sin(rot), x1 * np.sin( + rot) + z1 * np.cos(rot) + + h = double_perlin3D(x1, y1, z1, seed=seed) + + vT = double_perlin3D( + x1, y1, z1, a=20, persistance=0.8, seed=seed + 132) + vH = double_perlin3D( + x1, y1, z1, a=20, persistance=0.8, seed=seed + 463) + + c = double_perlin3D( + x1, y1, z1, a=5, d=5, persistance=0.6, seed=seed + 212) + + c = .6 * c + .4 * double_perlin3D( + x1, + y1, + z1, + a=20, + d=10, + persistance=0.8, + seed=seed + 213) + + cloud_map[i, ok, d1, d2] = c + + rgbref = colorize(x1, y1, z1, h, vT, vH) + + dhdx = np.zeros_like(h) + dhdy = np.zeros_like(h) + dhdz = np.zeros_like(h) + + dhdx[h > 0.5] = (double_perlin3D( + x1[h > 0.5] + eps, y1[h > 0.5], z1[h > 0.5], seed=seed) + - h[h > 0.5]) / eps + dhdy[h > 0.5] = (double_perlin3D( + x1[h > 0.5], y1[h > 0.5] + eps, z1[h > 0.5], seed=seed) + - h[h > 0.5]) / eps + dhdz[h > 0.5] = (double_perlin3D( + x1[h > 0.5], y1[h > 0.5], z1[h > 0.5] + eps, seed=seed) + - h[h > 0.5]) / eps + + norm = normal(ox1, y1, oz1, h, dhdx, dhdy, dhdz, A=0.05) + + i_var = np.ones_like(h) + + i_var[h > 0.5] = (0.5 + double_perlin3D( + x1[h > 0.5], + y1[h > 0.5], + z1[h > 0.5], + a=20, + d=5, + persistance=0.8, + seed=seed + 214)) + + im_rgbref_base[i, ok, d1, d2, 0] = rgbref[:, 0] * i_var + im_rgbref_base[i, ok, d1, d2, 1] = rgbref[:, 1] * i_var + im_rgbref_base[i, ok, d1, d2, 2] = rgbref[:, 2] * i_var + im_rgbref_base[i, ok, d1, d2, 3] = rgbref[:, 3] + + normal_base[i, ok, d1, d2, 0] = norm[0] + normal_base[i, ok, d1, d2, 1] = norm[1] + normal_base[i, ok, d1, d2, 2] = norm[2] + + im_rgbref_base[i, ~ok, d1, d2, 0] = background( + x[i, ~ok], y[i, ~ok])[0, :] + im_rgbref_base[i, ~ok, d1, d2, 1] = background( + x[i, ~ok], y[i, ~ok])[1, :] + im_rgbref_base[i, ~ok, d1, d2, 2] = background( + x[i, ~ok], y[i, ~ok])[2, :] + + return im_rgbref_base, normal_base, cloud_map + + +def lum_perc(cos): + u = np.log(1 + 5 * np.abs(cos)) / np.log(1 + 5) + return u + + +def apply_illum(im_rgbref_base, normal_base, cloud_map, ilum=(1, 0.5, 1)): + im = np.zeros((res, res, 3)) + for i in range(res): + for d1 in range(multi_sampling): + for d2 in range(multi_sampling): + z2 = radius**2 - (x[i, :] + dx[d1, d2])**2 - ( + y[i, :] + dy[d1, d2])**2 + ok = z2 >= 0 + + if np.any(ok): + x1, y1, z1 = (x[i, ok] + dx[d1, d2]), ( + y[i, ok] + dy[d1, d2]), np.sqrt(z2[ok]) + + cos1 = cos_vec(ilum, (x1, y1, z1)) + norm = normal_base[i, ok, d1, d2, 0], normal_base[ + i, ok, d1, d2, 1], normal_base[i, ok, d1, d2, 2] + cos = np.maximum(cos_vec(ilum, norm), cos1 / 2) + + refvec = reflection((0, 0, -1), norm) + refillum = np.where( + cos1 >= 0, + np.maximum(cos_vec(refvec, ilum), 0)**100, 0) + + u = np.where(cos > 0, lum_perc(cos), 0) / multi_sampling**2 + q = refillum * im_rgbref_base[i, ok, d1, d2, + 3] / multi_sampling**2 + + r_sol = im_rgbref_base[i, ok, d1, d2, 0] * u + q + g_sol = im_rgbref_base[i, ok, d1, d2, 1] * u + q + b_sol = im_rgbref_base[i, ok, d1, d2, 2] * u + q + + u = np.where(cos1 > 0, lum_perc(cos1), + 0) / multi_sampling**2 + f = np.exp( + -14 * (np.clip(cloud_map[i, ok, d1, d2], 0.5, 1) - .5)) + + im[i, ok, 0] += r_sol * f + (1 - f) * u + im[i, ok, 1] += g_sol * f + (1 - f) * u + im[i, ok, 2] += b_sol * f + (1 - f) * u + + xx1 = x[i, ~ok] + dx[d1, d2] + yy1 = y[i, ~ok] + dy[d1, d2] + atmos_dens = np.exp(-(np.sqrt(xx1**2 + yy1**2) - radius) / + (atmos_size * radius)) + + atmos_illum = atmos_dens * (np.maximum( + cos_vec(ilum, (xx1, yy1, 0)), 0)) + + im[i, ~ok, 0] += ( + im_rgbref_base[i, ~ok, d1, d2, 0] + + atmos_illum * atmos_colour[0]) / multi_sampling**2 + im[i, ~ok, 1] += ( + im_rgbref_base[i, ~ok, d1, d2, 1] + + atmos_illum * atmos_colour[1]) / multi_sampling**2 + im[i, ~ok, 2] += ( + im_rgbref_base[i, ~ok, d1, d2, 2] + + atmos_illum * atmos_colour[2]) / multi_sampling**2 + + np.clip(im, 0, 1, out=im) + return im + + +def gen_image(seed): + a, b, c = gen_rgbbase_and_normal(seed=seed) + im = apply_illum(a, b, c) + plt.axis("off") + plt.imshow(im) + plt.imsave("examples/out"+str(seed)+".png", im) + plt.show() + # send img to ipfs + + +if __name__ == '__main__': + for seed in range(0, 10): + gen_image(seed) \ No newline at end of file diff --git a/Generator/main2D.py b/Generator/main2D.py new file mode 100644 index 00000000..388d5cbb --- /dev/null +++ b/Generator/main2D.py @@ -0,0 +1,68 @@ +import numpy as np +import matplotlib.pyplot as plt + +from colour import colorize +from noise import double_perlin3D + +res = 2048 +radius = 0.4 + +seed = 42 + +phi = np.linspace(0, 2 * np.pi, res) +theta = np.linspace(-np.pi / 2, np.pi / 2, res) + +phi, theta = np.meshgrid(phi, theta) + +x1 = radius * np.cos(theta) * np.cos(phi) +z1 = radius * np.cos(theta) * np.sin(phi) +y1 = radius * np.sin(theta) + +x1 = x1.flatten() +y1 = y1.flatten() +z1 = z1.flatten() + +h = double_perlin3D(x1, y1, z1, seed=seed) +vT = double_perlin3D(x1, y1, z1, a=20, persistance=0.8, seed=seed + 132) +vH = double_perlin3D(x1, y1, z1, a=20, persistance=0.8, seed=seed + 463) + +rgbref = colorize(x1, y1, z1, h, vT, vH) + +i_var = np.ones_like(h) + +i_var[h > 0.5] = (0.5 + double_perlin3D( + x1[h > 0.5], + y1[h > 0.5], + z1[h > 0.5], + a=20, + d=5, + persistance=0.8, + seed=seed + 214)) + +rgbref[:, 0] = rgbref[:, 0] * i_var +rgbref[:, 1] = rgbref[:, 1] * i_var +rgbref[:, 2] = rgbref[:, 2] * i_var +rgbref[:, 3] = rgbref[:, 3] + +heightmap = np.zeros((res, res, 3)) + +heightmap[:, :, 0] = np.reshape(np.where(h > 0.5, (h - 0.5) * 2, 0), + (res, res)) +heightmap[:, :, 1] = heightmap[:, :, 0] +heightmap[:, :, 2] = heightmap[:, :, 0] + +rgb = np.zeros((res, res, 3)) +rgb[:, :, 0] = np.reshape(rgbref[:, 0], (res, res)) +rgb[:, :, 1] = np.reshape(rgbref[:, 1], (res, res)) +rgb[:, :, 2] = np.reshape(rgbref[:, 2], (res, res)) + + +spec = np.zeros((res, res, 3)) +spec[:, :, 0] = np.reshape(rgbref[:, 3], (res, res)) +spec[:, :, 1] = spec[:, :, 0] +spec[:, :, 2] = spec[:, :, 0] + + +plt.imsave("heightmap.png", heightmap) +plt.imsave("rgb.png", rgb.clip(0, 1)) +plt.imsave("spec.png", spec) \ No newline at end of file diff --git a/Generator/main_gif.py b/Generator/main_gif.py new file mode 100644 index 00000000..48f2ea41 --- /dev/null +++ b/Generator/main_gif.py @@ -0,0 +1,35 @@ +import numpy as np +import matplotlib.pyplot as plt + +import os + +from tqdm import tqdm +from main import gen_rgbbase_and_normal, apply_illum +import shutil + +fig = plt.figure() + +N_seconds = 20 +FPS = 20 + +N_frame = FPS * N_seconds + +phi = np.linspace(0, 2 * np.pi, N_frame) + +ims = [] + +os.mkdir('temp') + +for i in tqdm(range(N_frame)): + rgbref_base, normal_base, cmap = gen_rgbbase_and_normal(seed=436, + rot=phi[i]) + im = apply_illum(rgbref_base, normal_base, cmap) + plt.imsave("temp/out" + str(i) + ".png", im) + +os.system( + "ffmpeg -r " + str(FPS) + + " -f image2 -s 1024x1024 -i temp/out%d.png -vcodec libx264 -crf 25 " + "-pix_fmt yuv420p test.mp4" +) + +shutil.rmtree('temp') \ No newline at end of file diff --git a/Generator/noise.py b/Generator/noise.py new file mode 100644 index 00000000..5a4faf9e --- /dev/null +++ b/Generator/noise.py @@ -0,0 +1,174 @@ +import numpy as np + +from geometry import dot + + +def fract(x): + return x - np.floor(x) + + +def hash(x, y=0, z=0): + return hash13(fract(x * .1031), fract(y * .1031), fract(z * .1031)) + + +def hash13(p3x, p3y, p3z): + s = p3x * (p3y + 19.19) + p3y * (p3z + 19.19) + p3z * (p3x + 19.19) + p3x = p3x + s + p3y = p3y + s + p3z = p3z + s + return fract((p3x + p3y) * p3z) + + +def hash83(p3x, p3y, p3z, p3X, p3Y, p3Z): + p3xx = (p3x + 19.19) + p3yy = (p3y + 19.19) + p3zz = (p3z + 19.19) + p3XX = (p3X + 19.19) + p3YY = (p3Y + 19.19) + p3ZZ = (p3Z + 19.19) + + uxy = p3x * p3yy + uyz = p3y * p3zz + uzx = p3z * p3xx + + uxY = p3x * p3YY + uyZ = p3y * p3ZZ + uzX = p3z * p3XX + + uXy = p3X * p3yy + uYz = p3Y * p3zz + uZx = p3Z * p3xx + + uXY = p3X * p3YY + uYZ = p3Y * p3ZZ + uZX = p3Z * p3XX + + sxyz = uxy + uyz + uzx + sXyz = uXy + uyz + uzX + sxYz = uxY + uYz + uzx + sXYz = uXY + uYz + uzX + sxyZ = uxy + uyZ + uZx + sXyZ = uXy + uyZ + uZX + sxYZ = uxY + uYZ + uZx + sXYZ = uXY + uYZ + uZX + + nxy = p3x + p3y + nXy = p3X + p3y + nxY = p3x + p3Y + nXY = p3X + p3Y + + vxyz = fract((nxy + sxyz + sxyz) * (p3z + sxyz)) + vXyz = fract((nXy + sXyz + sXyz) * (p3z + sXyz)) + vxYz = fract((nxY + sxYz + sxYz) * (p3z + sxYz)) + vXYz = fract((nXY + sXYz + sXYz) * (p3z + sXYz)) + vxyZ = fract((nxy + sxyZ + sxyZ) * (p3Z + sxyZ)) + vXyZ = fract((nXy + sXyZ + sXyZ) * (p3Z + sXyZ)) + vxYZ = fract((nxY + sxYZ + sxYZ) * (p3Z + sxYZ)) + vXYZ = fract((nXY + sXYZ + sXYZ) * (p3Z + sXYZ)) + + return vxyz, vXyz, vxYz, vXYz, vxyZ, vXyZ, vxYZ, vXYZ + + +def polynom_interpol(x): + x2 = x * x + x4 = x2 * x2 + return 6 * x4 * x - 15 * x4 + 10 * x2 * x + + +def noise3D(x, y, z, seed=42): + u = hash(seed) + x = x + u * 1431.26742 + y = y + u * 3262.42321 + z = z + u * 7929.27883 + + ix, iy, iz = np.floor(x), np.floor(y), np.floor(z) + + fx, fy, fz = x - ix, y - iy, z - iz + + wx = polynom_interpol(fx) + wy = polynom_interpol(fy) + wz = polynom_interpol(fz) + + p3x = fract(ix * .1031) + p3y = fract(iy * .1031) + p3z = fract(iz * .1031) + p3X = fract((ix + 1) * .1031) + p3Y = fract((iy + 1) * .1031) + p3Z = fract((iz + 1) * .1031) + + vxyz, vXyz, vxYz, vXYz, vxyZ, vXyZ, vxYZ, vXYZ = hash83( + p3x, p3y, p3z, p3X, p3Y, p3Z) + + a1 = vxyz + (vXyz - vxyz) * wx + a2 = vxYz + (vXYz - vxYz) * wx + a3 = vxyZ + (vXyZ - vxyZ) * wx + a4 = vxYZ + (vXYZ - vxYZ) * wx + + b1 = a1 + (a2 - a1) * wy + b2 = a3 + (a4 - a3) * wy + + c = b1 + (b2 - b1) * wz + + return c + + +def perlin3D(x, y, z, octave=7, persistance=0.5, seed=42): + p = 1 + totp = 0 + tot = 0 + + for _ in range(octave): + c = noise3D(x, y, z, seed) + + tot += p * c + totp += p + + p *= persistance + + x, y, z = x * 2, y * 2, z * 2 + + return tot / totp + + +def double_perlin3D(x, y, z, a=4, d=2, octave=8, persistance=0.6, seed=42): + dx = perlin3D(x, y, z, 5, 0.5, seed * hash(seed, 0, 10)) + dy = perlin3D(x, y, z, 5, 0.5, seed * hash(seed, 10, 0)) + dz = perlin3D(x, y, z, 5, 0.5, seed * hash(seed, 10, 10)) + + dx, dy, dz = d * (dx - 0.5), d * (dy - 0.5), d * (dz - 0.5) + + P = perlin3D(a * x + dx, a * y + dy, a * z + dz, octave, persistance, seed) + + return P + + +def field(x, y, z, t): + strength = 7. + .03 * np.log(1.e-6 + fract(np.sin(t) * 4373.11)) + accum = np.zeros_like(x) + prev = np.zeros_like(x) + tw = np.zeros_like(x) + p = x, y, z + for i in range(32): + mag = dot(p, p) + p = np.abs(p[0]) / mag, np.abs(p[1]) / mag, np.abs(p[2]) / mag + p = p[0] - .51, p[1] - .4, p[2] - 1.3 + w = np.exp(-i / 7.) + accum += w * np.exp(-strength * (np.abs(mag - prev)**2.3)) + tw += w + prev = mag + + return np.maximum(0., 5. * accum / tw - .7) + + +def background(x, y, t=437.2): + zeros = np.zeros_like(x) + p = np.array([x / 4 + 2, y / 4 - 1.3, -1 + zeros]) + p += 0.15 * (np.array([zeros + 1.0, zeros + 0.5, zeros + np.sin(t / 128.) + ])) + + t1 = field(p[0, :], p[1, :], p[2, :], t) + v = (1. - np.exp((np.abs(x) - 1.) * 6.)) * (1. - np.exp( + (np.abs(y) - 1.) * 6.)) + + c1 = (.4 + .6 * v) * np.array([1.8 * t1 * t1 * t1, 1.4 * t1 * t1, t1]) + return c1 \ No newline at end of file