From 541c8cb42e7b6862ec05f13a743e5f57f629fb5e Mon Sep 17 00:00:00 2001 From: Crizomb <62544756+Crizomb@users.noreply.github.com> Date: Fri, 19 Aug 2022 19:28:22 +0200 Subject: [PATCH] create --- mandelbulb/__pycache__/camera.cpython-310.pyc | Bin 0 -> 1705 bytes mandelbulb/__pycache__/tools.cpython-310.pyc | Bin 0 -> 2129 bytes mandelbulb/camera.py | 60 + mandelbulb/controls.py | 1 + mandelbulb/main.py | 109 ++ mandelbulb/resources/programs/fragment.glsl | 175 +++ mandelbulb/resources/programs/hg_sdf.glsl | 1004 +++++++++++++++++ mandelbulb/resources/programs/vertex.glsl | 8 + mandelbulb/resources/textures/fur.jpg | 0 mandelbulb/resources/textures/noise.png | 0 mandelbulb/resources/textures/wall.jpg | 0 mandelbulb/test.py | 5 + mandelbulb/tools.py | 48 + 13 files changed, 1410 insertions(+) create mode 100644 mandelbulb/__pycache__/camera.cpython-310.pyc create mode 100644 mandelbulb/__pycache__/tools.cpython-310.pyc create mode 100644 mandelbulb/camera.py create mode 100644 mandelbulb/controls.py create mode 100644 mandelbulb/main.py create mode 100644 mandelbulb/resources/programs/fragment.glsl create mode 100644 mandelbulb/resources/programs/hg_sdf.glsl create mode 100644 mandelbulb/resources/programs/vertex.glsl create mode 100644 mandelbulb/resources/textures/fur.jpg create mode 100644 mandelbulb/resources/textures/noise.png create mode 100644 mandelbulb/resources/textures/wall.jpg create mode 100644 mandelbulb/test.py create mode 100644 mandelbulb/tools.py diff --git a/mandelbulb/__pycache__/camera.cpython-310.pyc b/mandelbulb/__pycache__/camera.cpython-310.pyc new file mode 100644 index 0000000000000000000000000000000000000000..46a5db81c89191abbe6bad0dd23ec1d1f9aca2e1 GIT binary patch literal 1705 zcmZux&u=3&6t+E{$s|pG&|SK%io^vDvxkkAUN|7MPy}ctqG~CZa2O3^w@H}E1beb< zqUnY7Z{WN=!as!ME2sVoT;O|7+HET=z0daZv;E%po?UHjb{N`UzyJHwpOUe^X>zv( zn7n}4X8^?%FIdS3oYPno;*|QJ2nO;SrlbllnF^O;5USvNmPFUs!xHY=Dyy@+$gIVL zKketFzuO{AUO?=(0L=zmv4K!r2V<#(#%Ve{i~RX=5UEJDF4>@^V%0_;>$d8sjZ3!V zgO2K|P3Sh%7Vh76(Q#I4lR06YyI?vTxv(B;xPnM;TZj=wc?_{H0SzA`&Yy*Cy=1@f zhBqwbN&tllfPqKZ&^T;JA2sZ4>xx}S^l>AvXh&Nycp&=;zox^Jz;JXAAnCZ!4zsB? zLhg1 z>f+2FZir&jBA|r-!OSCc|NPh+R;BLcg)a4^?w!o@VP?t`Q;l_A+uqs8*t&Q0RS#Y& zZ7kfw++_2#%uGHUP0mv_dVNpP!{wm-*>vvWG#yPwb(+$(b{DY4BN5$7EOz)}8y=JW zjcqPqchkse5}%t-F{qK9CBG0gt9gwIQIKWOh^2rws5$ZxEX9Sa1LQ|SA6$e@h)e{1 z&;2!JA(5^1ZYpb?e68~%C(EopN+Ltfdfo;S@qF#^mB$|jWRrUOR`8f? zcK}Zyme`H?9`Eq3kbJSZ4sVs3oUXrV6d2{*2RFt5!Lcz{{9x6K1BOr#CY9hAv|T@2vt@kwplJyp^n5jZ`ZeLg_7U=#YK0ha0bnY)?!Qdj>vPRic=_GMlut zL9j*8B{%@A*>PgN7(^5OYdT*z@n$`Kvv2{3e8|o=_I(Jkbjlvz<&PozLI#VSH9xmi zd@uj>bVxi!#B+oED}KoK(I2u4(Fo*L^vxmG{H00YTC=pQW_UQ=S!%Scidj9XCXSa6 zymDP{c3OQgIWP26oM7o+*rBvNTDx`E;$HG7k++_GfC=6{y8@<2)?54vWA6Wf!2bih zpGZR__{Zmq#k_}U!dI0>Z}D9*xo9@UkEwgmdY3w2>$Z|sYyG54?`wj`lcuhdB4jly zXxvTHUuIdcUJ;Y23p0Z!>%z05%4+jI-9_c<9dccZ?_-pgx*k@_d_?+>3EsKnN5063 z9;#YP6%n(29`D4l*g&aELu*OoHg3OqZ4FfyzTB!RtP3ZzaymB&RNnUgE^zT@t1vz% OdTk8o;xWZySNsQk;d4&_ literal 0 HcmV?d00001 diff --git a/mandelbulb/__pycache__/tools.cpython-310.pyc b/mandelbulb/__pycache__/tools.cpython-310.pyc new file mode 100644 index 0000000000000000000000000000000000000000..847e9798fb17e851abc07cf6ce6eafc6160a11c7 GIT binary patch literal 2129 zcmZWqO=}}Z5bf$7eap6N%a%VjOUM#fg~Zy}fxu!+LVQUM#_T$<3@m0mJ(A_jXtX`; z#WLedWJoSK1o8uRkNzk54Rg)OkXvrKBvmEJmbC-(`gKorSHJ42o_KDqEb#gJkN-}# z93lRq&h(SV&Nuibefg1&dmXRS7AOXxP_?b7P;@BUCAvm zhp_CH(IyvGqN-a35y07;JFoN9)`hze)Y_m6x_AxlBJ+|iBVW=Lv{=?v#yP@vUeD`= zYsanYnqIsX?uuU0%h;~!x?aJyp;vVS+cmwW*RfsK8+sGl4Gk|v^U*D;5eH2tQ+}Le zD(UtzWj%!8qTiI6(~C`%DIM4*WI5l&NrI%LJ?r_Qmn63mr~8T)89&4^e&RcPlO+Tz z&L9;jv}Gza{3c(D8+kyhX*=i#H!_pGewOo$@dizm!9@lGQ^k2V!qD`8-PUOw1uZ`e zqM&D6hXeoAGtr@m&jQ~jt#;QWwzae0N=`i;nB>?7*vG|=BhUD!-CpNdci%o3XxTUp zldb+BbKVC2_B?v|hjc7~sbk~Q_9qYVJ2~`qA}5t?m8y{-Q^5~r4x#Co8u3Mh3@%v+ z%NW5j4U^LyGDyg29ZZ9by$nt=@IMG<4!=$X^Bo71yuXJt3ML#hl_3pgf!8_0?izL1 z2{s5Qj?){6Ji9|$%nHFr2tFW89^zCC?DaX46w`g3P%PcXZP}6Uv6jpx8squj9fhIeHHt+wFat%Q^TI%qPfwnmc)j&ul4$A?86{SRh6lxj&(b7|Ob62#mfun2M{q z-w@|HmN)cAsWOuwvj2jQ;0Zc5pCUZi4wXXzUyjf#BovaMQ zB^@GOS<-%q?{?g82fw4}?*iNgN$nvWj5E zF9I)O(DJ&N;?{>KkfZt_VE?+wE)x1MESZ(u6; zd^41H!){{zIO^|xMvKwPb#eJ5IuGq}aAC({wsynV_rhck$;7@1hLRNvhHJlxeB|;O zF1=G7?G!BA3SZY#Xd?)8nFmFUOeh26{G zb2MT%RWM)Cgv~71JCFK -res2.x) ? res1 : vec3(-res2.x, res2.yz); +} + +vec3 fOpUnionId(in vec3 res1, in vec3 res2){ + return (res1.x < res2.x) ? res1 : res2; +} + + +vec3 map(in vec3 p){ + //pMod3(p, vec3(10.0)); + + float planeId = 1; + float planeColor = packColor(vec3(50, 50, 50)); + float planeDist = fPlane(p, vec3(0, 1, 0), 1.0); + vec3 plane = vec3(planeDist, planeColor, planeId); + + float sphereId = 2; + float sphereColor = packColor(vec3(200, 0, 0)); + vec3 spherePos = vec3(2.0, 0.0, 0.0); + float sphereDist = fSphere(p-spherePos, 1); + vec3 sphere = vec3(sphereDist, sphereColor, sphereId); + + float playerId = 3; + float playerColor = packColor(vec3(0, 200, 0));; + vec3 playerPos = ro+vec3(0, -1, DIST_CAM); + float playerDist = fCapsule(p-playerPos, 0.2, 0.2); + vec3 player = vec3(playerDist, playerColor, playerId); + + float mandelBulbId = 3; + float mandelBulbColor = packColor(vec3(50, 50, 50));; + float mandelBulbDist = fmandelBulbe(p);; + vec3 mandelBulb = vec3(mandelBulbDist, mandelBulbColor, mandelBulbId);; + + //vec3 res = fOpUnionId(plane, sphere); + vec3 res = mandelBulb; + //res = fOpUnionId(res, mandelBulb); + + + return res; +} + +vec4 rayMarch(vec3 ro, vec3 rd) { + vec3 hit = vec3(0.0); + vec3 object = vec3(0.0); + vec3 p = vec3(0.0); + float compt = 0; + for (float i = 0; i < MAX_STEPS; i++) { + p = ro + object.x * rd; + hit = map(p); + object.x += hit.x; + object.yz = hit.yz; + if (hit.z == 2 || hit.z == 3) compt += 1; + + if (abs(hit.x) < EPSILON || object.x > MAX_DIST){ + if (hit.z != 2) return vec4(object, compt); + return vec4(object, 0); + } + } + return vec4(object, compt); +} + +vec3 getNormal(vec3 p) { + vec2 e = vec2(EPSILON, 0.0); + vec3 n = vec3(map(p).x) - vec3(map(p - e.xyy).x, map(p - e.yxy).x, map(p - e.yyx).x); + return normalize(n); +} + +vec3 getLight(vec3 p, vec3 rd, vec3 color) { + vec3 lightPos = vec3(10.0, 55.0, -20.0); + vec3 L = normalize(lightPos - p); + vec3 N = getNormal(p); + vec3 V = -rd; + vec3 R = reflect(-L, N); + + vec3 specColor = vec3(0.5); + vec3 specular = specColor * pow(clamp(dot(R, V), 0.0, 1.0), 10.0); + vec3 diffuse = color * clamp(dot(L, N), 0.0, 1.0); + vec3 ambient = color * 0.05; + vec3 fresnel = 0.25 * color * pow(1.0 + dot(rd, N), 3.0); + + // shadows + float d = rayMarch(p + N * 0.02, normalize(lightPos)).x; + if (d < length(lightPos - p)) return ambient + fresnel; + + return diffuse + ambient + specular + fresnel; +} + +vec3 getMaterial(vec3 p, float color_pack, float id) { + vec3 m = vec3(0.0); + vec3 color = unpackColor(color_pack); + switch (int(id)) { + case 1: + m = vec3(0.2 + 0.4 * mod(floor(p.x) + floor(p.z), 2.0)); break; + } + return color+m; + return color; +} + + + + +void render(inout vec3 col, in vec3 ro, in mat3 cam, in vec2 uv) { + vec3 r1 = vec3(0, 0, 1); + vec3 lookAt = ro + r1; + vec3 rd = cam*normalize(vec3(uv, FOV)); + vec4 result = rayMarch(ro, rd); + vec3 object = result.xyz; + vec3 glowing_color = vec3(0.5, 0, 0); + float nb_itter = result.w; + + + + vec3 background = vec3(0, 0, 0); + if (object.x < MAX_DIST) { + vec3 p = ro + object.x * rd; + vec3 material = getMaterial(p, object.y, object.z); + col += getLight(p, rd, material); + // fog + col = mix(col, background, 1.0 - exp(-0.00002 * object.x * object.x)); + + } else { + // col += background - max(0.9 * rd.y, 0.0); + col = col; + } + float glow_fact = smoothstep(0, MAX_STEPS, nb_itter); + col += glowing_color*glow_fact; +} + +void main() { + vec2 uv = (2.0 * gl_FragCoord.xy - u_resolution.xy) / u_resolution.y; + + vec3 col; + render(col, ro, cam, uv); + + // gamma correction + col = pow(col, vec3(0.4545)); + fragColor = vec4(col, 1.0); +} diff --git a/mandelbulb/resources/programs/hg_sdf.glsl b/mandelbulb/resources/programs/hg_sdf.glsl new file mode 100644 index 0000000..312715a --- /dev/null +++ b/mandelbulb/resources/programs/hg_sdf.glsl @@ -0,0 +1,1004 @@ +//////////////////////////////////////////////////////////////// +// +// HG_SDF +// +// GLSL LIBRARY FOR BUILDING SIGNED DISTANCE BOUNDS +// +// version 2021-07-28 +// +// Check https://mercury.sexy/hg_sdf for updates +// and usage examples. Send feedback to spheretracing@mercury.sexy. +// +// Brought to you by MERCURY https://mercury.sexy/ +// +// +// +// Released dual-licensed under +// Creative Commons Attribution-NonCommercial (CC BY-NC) +// or +// MIT License +// at your choice. +// +// SPDX-License-Identifier: MIT OR CC-BY-NC-4.0 +// +// ///// +// +// CC-BY-NC-4.0 +// https://creativecommons.org/licenses/by-nc/4.0/legalcode +// https://creativecommons.org/licenses/by-nc/4.0/ +// +// ///// +// +// MIT License +// +// Copyright (c) 2011-2021 Mercury Demogroup +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in all +// copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. +// +// ///// +// +//////////////////////////////////////////////////////////////// +// +// How to use this: +// +// 1. Build some system to #include glsl files in each other. +// Include this one at the very start. Or just paste everywhere. +// 2. Build a sphere tracer. See those papers: +// * "Sphere Tracing" https://link.springer.com/article/10.1007%2Fs003710050084 +// * "Enhanced Sphere Tracing" http://diglib.eg.org/handle/10.2312/stag.20141233.001-008 +// * "Improved Ray Casting of Procedural Distance Bounds" https://www.bibsonomy.org/bibtex/258e85442234c3ace18ba4d89de94e57d +// The Raymnarching Toolbox Thread on pouet can be helpful as well +// http://www.pouet.net/topic.php?which=7931&page=1 +// and contains links to many more resources. +// 3. Use the tools in this library to build your distance bound f(). +// 4. ??? +// 5. Win a compo. +// +// (6. Buy us a beer or a good vodka or something, if you like.) +// +//////////////////////////////////////////////////////////////// +// +// Table of Contents: +// +// * Helper functions and macros +// * Collection of some primitive objects +// * Domain Manipulation operators +// * Object combination operators +// +//////////////////////////////////////////////////////////////// +// +// Why use this? +// +// The point of this lib is that everything is structured according +// to patterns that we ended up using when building geometry. +// It makes it more easy to write code that is reusable and that somebody +// else can actually understand. Especially code on Shadertoy (which seems +// to be what everybody else is looking at for "inspiration") tends to be +// really ugly. So we were forced to do something about the situation and +// release this lib ;) +// +// Everything in here can probably be done in some better way. +// Please experiment. We'd love some feedback, especially if you +// use it in a scene production. +// +// The main patterns for building geometry this way are: +// * Stay Lipschitz continuous. That means: don't have any distance +// gradient larger than 1. Try to be as close to 1 as possible - +// Distances are euclidean distances, don't fudge around. +// Underestimating distances will happen. That's why calling +// it a "distance bound" is more correct. Don't ever multiply +// distances by some value to "fix" a Lipschitz continuity +// violation. The invariant is: each fSomething() function returns +// a correct distance bound. +// * Use very few primitives and combine them as building blocks +// using combine opertors that preserve the invariant. +// * Multiply objects by repeating the domain (space). +// If you are using a loop inside your distance function, you are +// probably doing it wrong (or you are building boring fractals). +// * At right-angle intersections between objects, build a new local +// coordinate system from the two distances to combine them in +// interesting ways. +// * As usual, there are always times when it is best to not follow +// specific patterns. +// +//////////////////////////////////////////////////////////////// +// +// FAQ +// +// Q: Why is there no sphere tracing code in this lib? +// A: Because our system is way too complex and always changing. +// This is the constant part. Also we'd like everyone to +// explore for themselves. +// +// Q: This does not work when I paste it into Shadertoy!!!! +// A: Yes. It is GLSL, not GLSL ES. We like real OpenGL +// because it has way more features and is more likely +// to work compared to browser-based WebGL. We recommend +// you consider using OpenGL for your productions. Most +// of this can be ported easily though. +// +// Q: How do I material? +// A: We recommend something like this: +// Write a material ID, the distance and the local coordinate +// p into some global variables whenever an object's distance is +// smaller than the stored distance. Then, at the end, evaluate +// the material to get color, roughness, etc., and do the shading. +// +// Q: I found an error. Or I made some function that would fit in +// in this lib. Or I have some suggestion. +// A: Awesome! Drop us a mail at spheretracing@mercury.sexy. +// +// Q: Why is this not on github? +// A: Because we were too lazy. If we get bugged about it enough, +// we'll do it. +// +// Q: Your license sucks for me. +// A: Oh. What should we change it to? +// +// Q: I have trouble understanding what is going on with my distances. +// A: Some visualization of the distance field helps. Try drawing a +// plane that you can sweep through your scene with some color +// representation of the distance field at each point and/or iso +// lines at regular intervals. Visualizing the length of the +// gradient (or better: how much it deviates from being equal to 1) +// is immensely helpful for understanding which parts of the +// distance field are broken. +// +//////////////////////////////////////////////////////////////// + + + + + + +//////////////////////////////////////////////////////////////// +// +// HELPER FUNCTIONS/MACROS +// +//////////////////////////////////////////////////////////////// + +#define PI 3.14159265 +#define TAU (2*PI) +#define PHI (sqrt(5)*0.5 + 0.5) + +// Clamp to [0,1] - this operation is free under certain circumstances. +// For further information see +// http://www.humus.name/Articles/Persson_LowLevelThinking.pdf and +// http://www.humus.name/Articles/Persson_LowlevelShaderOptimization.pdf +#define saturate(x) clamp(x, 0, 1) + + +// Sign function that doesn't return 0 +float sgn(float x) { + return (x<0)?-1.0:1.0; +} + +vec2 sgn(vec2 v) { + return vec2((v.x<0)?-1:1, (v.y<0)?-1:1); +} + +float square (float x) { + return x*x; +} + +vec2 square (vec2 x) { + return x*x; +} + +vec3 square (vec3 x) { + return x*x; +} + +float lengthSqr(vec3 x) { + return dot(x, x); +} + + +// Maximum/minumum elements of a vector +float vmax(vec2 v) { + return max(v.x, v.y); +} + +float vmax(vec3 v) { + return max(max(v.x, v.y), v.z); +} + +float vmax(vec4 v) { + return max(max(v.x, v.y), max(v.z, v.w)); +} + +float vmin(vec2 v) { + return min(v.x, v.y); +} + +float vmin(vec3 v) { + return min(min(v.x, v.y), v.z); +} + +float vmin(vec4 v) { + return min(min(v.x, v.y), min(v.z, v.w)); +} + +float sigmoid(in float x, in float c1, in float c2){ + return 1/(1+exp(-c1*(x-c2))); +} + + + + +//////////////////////////////////////////////////////////////// +// +// PRIMITIVE DISTANCE FUNCTIONS +// +//////////////////////////////////////////////////////////////// +// +// Conventions: +// +// Everything that is a distance function is called fSomething. +// The first argument is always a point in 2 or 3-space called

. +// Unless otherwise noted, (if the object has an intrinsic "up" +// side or direction) the y axis is "up" and the object is +// centered at the origin. +// +//////////////////////////////////////////////////////////////// + +float fSphere(vec3 p, float r) { + return length(p) - r; +} + +// Plane with normal n (n is normalized) at some distance from the origin +float fPlane(vec3 p, vec3 n, float distanceFromOrigin) { + return dot(p, n) + distanceFromOrigin; +} + +// Cheap Box: distance to corners is overestimated +float fBoxCheap(vec3 p, vec3 b) { //cheap box + return vmax(abs(p) - b); +} + +// Box: correct distance to corners +float fBox(vec3 p, vec3 b) { + vec3 d = abs(p) - b; + return length(max(d, vec3(0))) + vmax(min(d, vec3(0))); +} + +// Same as above, but in two dimensions (an endless box) +float fBox2Cheap(vec2 p, vec2 b) { + return vmax(abs(p)-b); +} + +float fBox2(vec2 p, vec2 b) { + vec2 d = abs(p) - b; + return length(max(d, vec2(0))) + vmax(min(d, vec2(0))); +} + + +// Endless "corner" +float fCorner (vec2 p) { + return length(max(p, vec2(0))) + vmax(min(p, vec2(0))); +} + +// Blobby ball object. You've probably seen it somewhere. This is not a correct distance bound, beware. +float fBlob(vec3 p) { + p = abs(p); + if (p.x < max(p.y, p.z)) p = p.yzx; + if (p.x < max(p.y, p.z)) p = p.yzx; + float b = max(max(max( + dot(p, normalize(vec3(1, 1, 1))), + dot(p.xz, normalize(vec2(PHI+1, 1)))), + dot(p.yx, normalize(vec2(1, PHI)))), + dot(p.xz, normalize(vec2(1, PHI)))); + float l = length(p); + return l - 1.5 - 0.2 * (1.5 / 2)* cos(min(sqrt(1.01 - b / l)*(PI / 0.25), PI)); +} + +// Cylinder standing upright on the xz plane +float fCylinder(vec3 p, float r, float height) { + float d = length(p.xz) - r; + d = max(d, abs(p.y) - height); + return d; +} + +// Capsule: A Cylinder with round caps on both sides +float fCapsule(vec3 p, float r, float c) { + return mix(length(p.xz) - r, length(vec3(p.x, abs(p.y) - c, p.z)) - r, step(c, abs(p.y))); +} + +// Distance to line segment between and , used for fCapsule() version 2below +float fLineSegment(vec3 p, vec3 a, vec3 b) { + vec3 ab = b - a; + float t = saturate(dot(p - a, ab) / dot(ab, ab)); + return length((ab*t + a) - p); +} + +// Capsule version 2: between two end points and with radius r +float fCapsule(vec3 p, vec3 a, vec3 b, float r) { + return fLineSegment(p, a, b) - r; +} + +// Torus in the XZ-plane +float fTorus(vec3 p, float smallRadius, float largeRadius) { + return length(vec2(length(p.xz) - largeRadius, p.y)) - smallRadius; +} + +// A circle line. Can also be used to make a torus by subtracting the smaller radius of the torus. +float fCircle(vec3 p, float r) { + float l = length(p.xz) - r; + return length(vec2(p.y, l)); +} + +// A circular disc with no thickness (i.e. a cylinder with no height). +// Subtract some value to make a flat disc with rounded edge. +float fDisc(vec3 p, float r) { + float l = length(p.xz) - r; + return l < 0 ? abs(p.y) : length(vec2(p.y, l)); +} + +// Hexagonal prism, circumcircle variant +float fHexagonCircumcircle(vec3 p, vec2 h) { + vec3 q = abs(p); + return max(q.y - h.y, max(q.x*sqrt(3)*0.5 + q.z*0.5, q.z) - h.x); + //this is mathematically equivalent to this line, but less efficient: + //return max(q.y - h.y, max(dot(vec2(cos(PI/3), sin(PI/3)), q.zx), q.z) - h.x); +} + +// Hexagonal prism, incircle variant +float fHexagonIncircle(vec3 p, vec2 h) { + return fHexagonCircumcircle(p, vec2(h.x*sqrt(3)*0.5, h.y)); +} + +// Cone with correct distances to tip and base circle. Y is up, 0 is in the middle of the base. +float fCone(vec3 p, float radius, float height) { + vec2 q = vec2(length(p.xz), p.y); + vec2 tip = q - vec2(0, height); + vec2 mantleDir = normalize(vec2(height, radius)); + float mantle = dot(tip, mantleDir); + float d = max(mantle, -q.y); + float projected = dot(tip, vec2(mantleDir.y, -mantleDir.x)); + + // distance to tip + if ((q.y > height) && (projected < 0)) { + d = max(d, length(tip)); + } + + // distance to base ring + if ((q.x > radius) && (projected > length(vec2(height, radius)))) { + d = max(d, length(q - vec2(radius, 0))); + } + return d; +} + +// +// "Generalized Distance Functions" by Akleman and Chen. +// see the Paper at https://www.viz.tamu.edu/faculty/ergun/research/implicitmodeling/papers/sm99.pdf +// +// This set of constants is used to construct a large variety of geometric primitives. +// Indices are shifted by 1 compared to the paper because we start counting at Zero. +// Some of those are slow whenever a driver decides to not unroll the loop, +// which seems to happen for fIcosahedron und fTruncatedIcosahedron on nvidia 350.12 at least. +// Specialized implementations can well be faster in all cases. +// + +const vec3 GDFVectors[19] = vec3[]( + normalize(vec3(1, 0, 0)), + normalize(vec3(0, 1, 0)), + normalize(vec3(0, 0, 1)), + + normalize(vec3(1, 1, 1 )), + normalize(vec3(-1, 1, 1)), + normalize(vec3(1, -1, 1)), + normalize(vec3(1, 1, -1)), + + normalize(vec3(0, 1, PHI+1)), + normalize(vec3(0, -1, PHI+1)), + normalize(vec3(PHI+1, 0, 1)), + normalize(vec3(-PHI-1, 0, 1)), + normalize(vec3(1, PHI+1, 0)), + normalize(vec3(-1, PHI+1, 0)), + + normalize(vec3(0, PHI, 1)), + normalize(vec3(0, -PHI, 1)), + normalize(vec3(1, 0, PHI)), + normalize(vec3(-1, 0, PHI)), + normalize(vec3(PHI, 1, 0)), + normalize(vec3(-PHI, 1, 0)) +); + +// Version with variable exponent. +// This is slow and does not produce correct distances, but allows for bulging of objects. +float fGDF(vec3 p, float r, float e, int begin, int end) { + float d = 0; + for (int i = begin; i <= end; ++i) + d += pow(abs(dot(p, GDFVectors[i])), e); + return pow(d, 1/e) - r; +} + +// Version with without exponent, creates objects with sharp edges and flat faces +float fGDF(vec3 p, float r, int begin, int end) { + float d = 0; + for (int i = begin; i <= end; ++i) + d = max(d, abs(dot(p, GDFVectors[i]))); + return d - r; +} + +// Primitives follow: + +float fOctahedron(vec3 p, float r, float e) { + return fGDF(p, r, e, 3, 6); +} + +float fDodecahedron(vec3 p, float r, float e) { + return fGDF(p, r, e, 13, 18); +} + +float fIcosahedron(vec3 p, float r, float e) { + return fGDF(p, r, e, 3, 12); +} + +float fTruncatedOctahedron(vec3 p, float r, float e) { + return fGDF(p, r, e, 0, 6); +} + +float fTruncatedIcosahedron(vec3 p, float r, float e) { + return fGDF(p, r, e, 3, 18); +} + +float fOctahedron(vec3 p, float r) { + return fGDF(p, r, 3, 6); +} + +float fDodecahedron(vec3 p, float r) { + return fGDF(p, r, 13, 18); +} + +float fIcosahedron(vec3 p, float r) { + return fGDF(p, r, 3, 12); +} + +float fTruncatedOctahedron(vec3 p, float r) { + return fGDF(p, r, 0, 6); +} + +float fTruncatedIcosahedron(vec3 p, float r) { + return fGDF(p, r, 3, 18); +} + + +//////////////////////////////////////////////////////////////// +// +// DOMAIN MANIPULATION OPERATORS +// +//////////////////////////////////////////////////////////////// +// +// Conventions: +// +// Everything that modifies the domain is named pSomething. +// +// Many operate only on a subset of the three dimensions. For those, +// you must choose the dimensions that you want manipulated +// by supplying e.g. or +// +// is always the first argument and modified in place. +// +// Many of the operators partition space into cells. An identifier +// or cell index is returned, if possible. This return value is +// intended to be optionally used e.g. as a random seed to change +// parameters of the distance functions inside the cells. +// +// Unless stated otherwise, for cell index 0,

is unchanged and cells +// are centered on the origin so objects don't have to be moved to fit. +// +// +//////////////////////////////////////////////////////////////// + + + +// Rotate around a coordinate axis (i.e. in a plane perpendicular to that axis) by angle . +// Read like this: R(p.xz, a) rotates "x towards z". +// This is fast if is a compile-time constant and slower (but still practical) if not. +void pR(inout vec2 p, float a) { + p = cos(a)*p + sin(a)*vec2(p.y, -p.x); +} + +// Shortcut for 45-degrees rotation +void pR45(inout vec2 p) { + p = (p + vec2(p.y, -p.x))*sqrt(0.5); +} + +// Repeat space along one axis. Use like this to repeat along the x axis: +// - using the return value is optional. +float pMod1(inout float p, float size) { + float halfsize = size*0.5; + float c = floor((p + halfsize)/size); + p = mod(p + halfsize, size) - halfsize; + return c; +} + +// Same, but mirror every second cell so they match at the boundaries +float pModMirror1(inout float p, float size) { + float halfsize = size*0.5; + float c = floor((p + halfsize)/size); + p = mod(p + halfsize,size) - halfsize; + p *= mod(c, 2.0)*2 - 1; + return c; +} + +// Repeat the domain only in positive direction. Everything in the negative half-space is unchanged. +float pModSingle1(inout float p, float size) { + float halfsize = size*0.5; + float c = floor((p + halfsize)/size); + if (p >= 0) + p = mod(p + halfsize, size) - halfsize; + return c; +} + +// Repeat only a few times: from indices to (similar to above, but more flexible) +float pModInterval1(inout float p, float size, float start, float stop) { + float halfsize = size*0.5; + float c = floor((p + halfsize)/size); + p = mod(p+halfsize, size) - halfsize; + if (c > stop) { //yes, this might not be the best thing numerically. + p += size*(c - stop); + c = stop; + } + if (c = (repetitions/2)) c = abs(c); + return c; +} + +// Repeat in two dimensions +vec2 pMod2(inout vec2 p, vec2 size) { + vec2 c = floor((p + size*0.5)/size); + p = mod(p + size*0.5,size) - size*0.5; + return c; +} + +// Same, but mirror every second cell so all boundaries match +vec2 pModMirror2(inout vec2 p, vec2 size) { + vec2 halfsize = size*0.5; + vec2 c = floor((p + halfsize)/size); + p = mod(p + halfsize, size) - halfsize; + p *= mod(c,vec2(2))*2 - vec2(1); + return c; +} + +// Same, but mirror every second cell at the diagonal as well +vec2 pModGrid2(inout vec2 p, vec2 size) { + vec2 c = floor((p + size*0.5)/size); + p = mod(p + size*0.5, size) - size*0.5; + p *= mod(c,vec2(2))*2 - vec2(1); + p -= size/2; + if (p.x > p.y) p.xy = p.yx; + return floor(c/2); +} + +// Repeat in three dimensions +vec3 pMod3(inout vec3 p, vec3 size) { + vec3 c = floor((p + size*0.5)/size); + p = mod(p + size*0.5, size) - size*0.5; + return c; +} + +// Mirror at an axis-aligned plane which is at a specified distance from the origin. +float pMirror (inout float p, float dist) { + float s = sgn(p); + p = abs(p)-dist; + return s; +} + +// Mirror in both dimensions and at the diagonal, yielding one eighth of the space. +// translate by dist before mirroring. +vec2 pMirrorOctant (inout vec2 p, vec2 dist) { + vec2 s = sgn(p); + pMirror(p.x, dist.x); + pMirror(p.y, dist.y); + if (p.y > p.x) + p.xy = p.yx; + return s; +} + +// Reflect space at a plane +float pReflect(inout vec3 p, vec3 planeNormal, float offset) { + float t = dot(p, planeNormal)+offset; + if (t < 0) { + p = p - (2*t)*planeNormal; + } + return sgn(t); +} + + +//////////////////////////////////////////////////////////////// +// +// OBJECT COMBINATION OPERATORS +// +//////////////////////////////////////////////////////////////// +// +// We usually need the following boolean operators to combine two objects: +// Union: OR(a,b) +// Intersection: AND(a,b) +// Difference: AND(a,!b) +// (a and b being the distances to the objects). +// +// The trivial implementations are min(a,b) for union, max(a,b) for intersection +// and max(a,-b) for difference. To combine objects in more interesting ways to +// produce rounded edges, chamfers, stairs, etc. instead of plain sharp edges we +// can use combination operators. It is common to use some kind of "smooth minimum" +// instead of min(), but we don't like that because it does not preserve Lipschitz +// continuity in many cases. +// +// Naming convention: since they return a distance, they are called fOpSomething. +// The different flavours usually implement all the boolean operators above +// and are called fOpUnionRound, fOpIntersectionRound, etc. +// +// The basic idea: Assume the object surfaces intersect at a right angle. The two +// distances and constitute a new local two-dimensional coordinate system +// with the actual intersection as the origin. In this coordinate system, we can +// evaluate any 2D distance function we want in order to shape the edge. +// +// The operators below are just those that we found useful or interesting and should +// be seen as examples. There are infinitely more possible operators. +// +// They are designed to actually produce correct distances or distance bounds, unlike +// popular "smooth minimum" operators, on the condition that the gradients of the two +// SDFs are at right angles. When they are off by more than 30 degrees or so, the +// Lipschitz condition will no longer hold (i.e. you might get artifacts). The worst +// case is parallel surfaces that are close to each other. +// +// Most have a float argument to specify the radius of the feature they represent. +// This should be much smaller than the object size. +// +// Some of them have checks like "if ((-a < r) && (-b < r))" that restrict +// their influence (and computation cost) to a certain area. You might +// want to lift that restriction or enforce it. We have left it as comments +// in some cases. +// +// usage example: +// +// float fTwoBoxes(vec3 p) { +// float box0 = fBox(p, vec3(1)); +// float box1 = fBox(p-vec3(1), vec3(1)); +// return fOpUnionChamfer(box0, box1, 0.2); +// } +// +//////////////////////////////////////////////////////////////// + + +// The "Chamfer" flavour makes a 45-degree chamfered edge (the diagonal of a square of size ): +float fOpUnionChamfer(float a, float b, float r) { + return min(min(a, b), (a - r + b)*sqrt(0.5)); +} + +// Intersection has to deal with what is normally the inside of the resulting object +// when using union, which we normally don't care about too much. Thus, intersection +// implementations sometimes differ from union implementations. +float fOpIntersectionChamfer(float a, float b, float r) { + return max(max(a, b), (a + r + b)*sqrt(0.5)); +} + +// Difference can be built from Intersection or Union: +float fOpDifferenceChamfer (float a, float b, float r) { + return fOpIntersectionChamfer(a, -b, r); +} + +// The "Round" variant uses a quarter-circle to join the two objects smoothly: +float fOpUnionRound(float a, float b, float r) { + vec2 u = max(vec2(r - a,r - b), vec2(0)); + return max(r, min (a, b)) - length(u); +} + +float fOpIntersectionRound(float a, float b, float r) { + vec2 u = max(vec2(r + a,r + b), vec2(0)); + return min(-r, max (a, b)) + length(u); +} + +float fOpDifferenceRound (float a, float b, float r) { + return fOpIntersectionRound(a, -b, r); +} + + +// The "Columns" flavour makes n-1 circular columns at a 45 degree angle: +float fOpUnionColumns(float a, float b, float r, float n) { + if ((a < r) && (b < r)) { + vec2 p = vec2(a, b); + float columnradius = r*sqrt(2)/((n-1)*2+sqrt(2)); + pR45(p); + p.x -= sqrt(2)/2*r; + p.x += columnradius*sqrt(2); + if (mod(n,2) == 1) { + p.y += columnradius; + } + // At this point, we have turned 45 degrees and moved at a point on the + // diagonal that we want to place the columns on. + // Now, repeat the domain along this direction and place a circle. + pMod1(p.y, columnradius*2); + float result = length(p) - columnradius; + result = min(result, p.x); + result = min(result, a); + return min(result, b); + } else { + return min(a, b); + } +} + +float fOpDifferenceColumns(float a, float b, float r, float n) { + a = -a; + float m = min(a, b); + //avoid the expensive computation where not needed (produces discontinuity though) + if ((a < r) && (b < r)) { + vec2 p = vec2(a, b); + float columnradius = r*sqrt(2)/n/2.0; + columnradius = r*sqrt(2)/((n-1)*2+sqrt(2)); + + pR45(p); + p.y += columnradius; + p.x -= sqrt(2)/2*r; + p.x += -columnradius*sqrt(2)/2; + + if (mod(n,2) == 1) { + p.y += columnradius; + } + pMod1(p.y,columnradius*2); + + float result = -length(p) + columnradius; + result = max(result, p.x); + result = min(result, a); + return -min(result, b); + } else { + return -m; + } +} + +float fOpIntersectionColumns(float a, float b, float r, float n) { + return fOpDifferenceColumns(a,-b,r, n); +} + +// The "Stairs" flavour produces n-1 steps of a staircase: +// much less stupid version by paniq +float fOpUnionStairs(float a, float b, float r, float n) { + float s = r/n; + float u = b-r; + return min(min(a,b), 0.5 * (u + a + abs ((mod (u - a + s, 2 * s)) - s))); +} + +// We can just call Union since stairs are symmetric. +float fOpIntersectionStairs(float a, float b, float r, float n) { + return -fOpUnionStairs(-a, -b, r, n); +} + +float fOpDifferenceStairs(float a, float b, float r, float n) { + return -fOpUnionStairs(-a, b, r, n); +} + + +// Similar to fOpUnionRound, but more lipschitz-y at acute angles +// (and less so at 90 degrees). Useful when fudging around too much +// by MediaMolecule, from Alex Evans' siggraph slides +float fOpUnionSoft(float a, float b, float r) { + float e = max(r - abs(a - b), 0); + return min(a, b) - e*e*0.25/r; +} + + +// produces a cylindical pipe that runs along the intersection. +// No objects remain, only the pipe. This is not a boolean operator. +float fOpPipe(float a, float b, float r) { + return length(vec2(a, b)) - r; +} + +// first object gets a v-shaped engraving where it intersect the second +float fOpEngrave(float a, float b, float r) { + return max(a, (a + r - abs(b))*sqrt(0.5)); +} + +// first object gets a capenter-style groove cut out +float fOpGroove(float a, float b, float ra, float rb) { + return max(a, min(a + ra, rb - abs(b))); +} + +// first object gets a capenter-style tongue attached +float fOpTongue(float a, float b, float ra, float rb) { + return min(a, max(a - ra, abs(b) - rb)); +} + +// +const int Iterations = 100; +const int Power = 6; +const float Bailout = 8; + +float fmandelBulbe(vec3 pos) { + vec3 z = pos*0.3; + float dr = 1.0; + float r = 0.0; + for (int i = 0; i < Iterations ; i++) { + r = length(z); + if (r>Bailout) break; + + // convert to polar coordinates + float theta = acos(z.z/r); + float phi = atan(z.y,z.x); + dr = pow( r, Power-1.0)*Power*dr + 1.0; + + // scale and rotate the point + float zr = pow( r,Power); + theta = theta*Power; + phi = phi*Power; + + // convert back to cartesian coordinates + z = zr*vec3(sin(theta)*cos(phi), sin(phi)*sin(theta), cos(theta)); + z+=pos; + } + return 0.5*log(r)*r/dr; +} + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +// propre \ No newline at end of file diff --git a/mandelbulb/resources/programs/vertex.glsl b/mandelbulb/resources/programs/vertex.glsl new file mode 100644 index 0000000..5789a3c --- /dev/null +++ b/mandelbulb/resources/programs/vertex.glsl @@ -0,0 +1,8 @@ +#version 440 core + + +in vec3 in_position; + +void main() { + gl_Position = vec4(in_position, 1); +} diff --git a/mandelbulb/resources/textures/fur.jpg b/mandelbulb/resources/textures/fur.jpg new file mode 100644 index 0000000..e69de29 diff --git a/mandelbulb/resources/textures/noise.png b/mandelbulb/resources/textures/noise.png new file mode 100644 index 0000000..e69de29 diff --git a/mandelbulb/resources/textures/wall.jpg b/mandelbulb/resources/textures/wall.jpg new file mode 100644 index 0000000..e69de29 diff --git a/mandelbulb/test.py b/mandelbulb/test.py new file mode 100644 index 0000000..7ce8c90 --- /dev/null +++ b/mandelbulb/test.py @@ -0,0 +1,5 @@ +import numpy as np +import tools + +m = np.array(((1, 2, 3), (4, 5, 6), (7, 8, 9))) +print(tools.mat3ToTuple(m)) diff --git a/mandelbulb/tools.py b/mandelbulb/tools.py new file mode 100644 index 0000000..2334850 --- /dev/null +++ b/mandelbulb/tools.py @@ -0,0 +1,48 @@ +import numpy as np +from numpy import cos, sin, tan, exp + +from numpy.linalg import norm, det +from numpy import cross +import numba +from dataclasses import dataclass + +PI = np.pi +TAU = 2*np.pi + + + +def vec2(x, y, dtype=float): + return np.array((x, y), dtype=dtype) + +def vec3(x, y, z, dtype=float): + return np.array((x, y, z), dtype=dtype) + +def mat3(a : np.array, b : np.array, c : np.array): + return np.array((a, b, c)).T + +def mat3ToTuple(m): + return m[0, 0], m[0, 1], m[0, 2], m[1, 0], m[1, 1], m[1, 2], m[2, 0], m[2, 1], m[2, 2] + #return m[0, 0], m[1, 0], m[2, 0], m[0, 1], m[1, 1], m[2, 1], m[0, 2], m[1, 2], m[2, 2] + +def getRot2D(theta): + return np.array( [( cos(theta), -sin(theta) ), (sin(theta), cos(theta))] ) + +def pR(p, a): + return cos(a)*p + sin(a)*vec2(p[1], -p[0]) + + +def add_tuples(t1, t2): + assert len(t1) == len(t2) + temp = [0]*len(t1) + for i in range(len(t1)): + temp[i] += t1[i] + t2[i] + return tuple(temp) + +def mult_ext_tuples(a, t1): + return tuple([a*t for t in t1]) + +def normalize(vec): + return vec/norm(vec) + + +