Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP: implementing marching cubes and isosurface plotting #58

Draft
wants to merge 11 commits into
base: main
Choose a base branch
from
103 changes: 103 additions & 0 deletions src/emmy/mathbox/marching_cubes/marching_cubes.cljs
Original file line number Diff line number Diff line change
@@ -0,0 +1,103 @@
(ns emmy.mathbox.marching-cubes.marching-cubes
(:require [emmy.mathbox.marching-cubes.triangle-table :as tt]))

(defn CubeIndex
"Returns index into 256 value Triangle Lookup Table"
[vertex_vector rhs]
;; TODO double check the index implementation.
;; compares vertex position to rhs of surface function and
;; returns the 8 bit number associated with that edge
(let [bit_exp (map-indexed

(fn [index item]
(* (if (> item rhs) 1 0)
(Math/pow 2 index)))
vertex_vector)]
(reduce + bit_exp)))

(defn interp
[x y a b c]
(let [t (/ (- c a) (- b a))]
(+ (* x (- 1 t)) (* y t))))

(defn IsosurfaceEdgePoints
[x y z xStep yStep zStep v1 v2 v4 v8 v16 v32 v64 v128 lhs rhs index]
;;https://github.com/ChristopherChudzicki/math3d-react/blob/e352d15d275fc2e060441ecf11574e2fdec3211b/client/src/utils/marchingCubes/index.js#L54
(let [xInc (+ x xStep)
yInc (+ y yStep)
zInc (+ z zStep)]
(case index
0 [(interp x xInc v1 v2 rhs) y z]
1 [xInc y (interp z, zInc v2 v4 rhs)]
2 [(interp x xInc v4 v8 rhs) y zInc]
3 [x y (interp z zInc v1 v8 rhs)]
4 [(interp x xInc v16 v32 rhs) yInc z]
5 [xInc yInc (interp z zInc v32 v64 rhs)]
6 [(interp x xInc v128 v64 rhs) yInc zInc]
7 [x yInc (interp z zInc v16 v128 rhs)]
8 [x (interp y yInc v2 v16 rhs) z]
9 [xInc (interp y (+ y yInc) v2 v32 rhs) z]
10 [xInc (interp y yInc v4 v64 rhs)]
11 [x (interp y yInc v8 v128 rhs) zInc]
(println (str "Index must be between 0 and 11. Given: " index)) ;TODO how do I throw errors in clojure?
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Will update to throw the error as per the outdated comment.

)))

(defn GenRange [minVal maxVal numPoints]
(map (fn [x] (+ (* x (/ (- maxVal minVal) numPoints)) minVal))
;; (take (int (/ (- maxVal minVal) resolution)) (range)) ;;XXX this was wrong
(take numPoints (range))))

(defn GenerateGrid [xMin xMax yMin yMax zMin zMax resolution]
(println "generating grid")
(let [points ()
xVals (GenRange xMin xMax resolution)
yVals (GenRange yMin yMax resolution)
zVals (GenRange zMin zMax resolution)]
(for [x xVals y yVals z zVals]
[x y z])))

(defn MarchingCubesCore [lhs rhs x y z xStep yStep zStep]
(let [v1 (lhs x y z)
v2 (lhs (+ x xStep) y z)
v4 (lhs (+ x xStep) y (+ z zStep))
v8 (lhs x y (+ z zStep))
v16 (lhs x (+ y yStep) z)
v32 (lhs (+ x xStep) (+ y yStep) z)
v64 (lhs (+ x xStep) (+ y yStep) (+ z zStep))
v128 (lhs x (+ y yStep) (+ z zStep))]

(let [IsosurfaceEdgePointsPartial (partial IsosurfaceEdgePoints x y z xStep yStep zStep v1 v2 v4 v8 v16 v32 v64 v128 lhs rhs)
indexVal (CubeIndex [v1 v2 v4 v8 v16 v32 v64 v128] rhs)]

(for [triangle (tt/TriangleTable indexVal)]
(map IsosurfaceEdgePointsPartial triangle)))))

(defn MarchingCubes [lhs rhs
& {:keys [xMin xMax yMin yMax zMin zMax resolution]
:or {xMin -1 xMax 1 yMin -1 yMax 1 zMin -1 zMax 1 resolution 100}}]
(let [
xStep (/ 1 resolution)
yStep (/ 1 resolution)
zStep (/ 1 resolution)
grid (GenerateGrid xMin xMax yMin yMax zMin zMax resolution)
]

(filter not-empty ;; XXX somehow some of these only return 'points' with 2 entries....they seem to mostly evaluate to the RHS though!
(for [point grid]
(let [x (point 0)
y (point 1)
z (point 2)]
(MarchingCubesCore lhs rhs x y z xStep yStep zStep)))))
)

(comment
vnredmon marked this conversation as resolved.
Show resolved Hide resolved
;; Some repl functions and data

(defn TestLhs [x y z]
(reduce +
(map (fn [a] (Math/pow a 2)) [x y z])))

(do
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This block gives me a vector with a bunch of points who's magnitude is pretty close to 1.

I can add some tests as desired (probably would be good for me to do to learn clojure's testing facilities, but I can also do that independent of implementing this).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes please, we definitely want tests.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Will add. 👍

(def MarchingCubes TestLhs 1 )
nil)
)
265 changes: 265 additions & 0 deletions src/emmy/mathbox/marching_cubes/triangle_table.cljs
Original file line number Diff line number Diff line change
@@ -0,0 +1,265 @@
(ns emmy.mathbox.marching-cubes.triangle-table)


(def TriangleTable
"Lookup table that defines the isosurface facets for determining where the isosurface point sits relative to the isosurface facet"

[
[]
[[0 8 3]]
[[0 1 9]]
[[1 8 3] [9 8 1]]
[[1 2 10]]
[[0 8 3] [1 2 10]]
[[9 2 10] [0 2 9]]
[[2 8 3] [2 10 8] [10 9 8]]
[[3 11 2]]
[[0 11 2] [8 11 0]]
[[1 9 0] [2 3 11]]
[[1 11 2] [1 9 11] [9 8 11]]
[[3 10 1] [11 10 3]]
[[0 10 1] [0 8 10] [8 11 10]]
[[3 9 0] [3 11 9] [11 10 9]]
[[9 8 10] [10 8 11]]
[[4 7 8]]
[[4 3 0] [7 3 4]]
[[0 1 9] [8 4 7]]
[[4 1 9] [4 7 1] [7 3 1]]
[[1 2 10] [8 4 7]]
[[3 4 7] [3 0 4] [1 2 10]]
[[9 2 10] [9 0 2] [8 4 7]]
[[2 10 9] [2 9 7] [2 7 3] [7 9 4]]
[[8 4 7] [3 11 2]]
[[11 4 7] [11 2 4] [2 0 4]]
[[9 0 1] [8 4 7] [2 3 11]]
[[4 7 11] [9 4 11] [9 11 2] [9 2 1]]
[[3 10 1] [3 11 10] [7 8 4]]
[[1 11 10] [1 4 11] [1 0 4] [7 11 4]]
[[4 7 8] [9 0 11] [9 11 10] [11 0 3]]
[[4 7 11] [4 11 9] [9 11 10]]
[[9 5 4]]
[[9 5 4] [0 8 3]]
[[0 5 4] [1 5 0]]
[[8 5 4] [8 3 5] [3 1 5]]
[[1 2 10] [9 5 4]]
[[3 0 8] [1 2 10] [4 9 5]]
[[5 2 10] [5 4 2] [4 0 2]]
[[2 10 5] [3 2 5] [3 5 4] [3 4 8]]
[[9 5 4] [2 3 11]]
[[0 11 2] [0 8 11] [4 9 5]]
[[0 5 4] [0 1 5] [2 3 11]]
[[2 1 5] [2 5 8] [2 8 11] [4 8 5]]
[[10 3 11] [10 1 3] [9 5 4]]
[[4 9 5] [0 8 1] [8 10 1] [8 11 10]]
[[5 4 0] [5 0 11] [5 11 10] [11 0 3]]
[[5 4 8] [5 8 10] [10 8 11]]
[[9 7 8] [5 7 9]]
[[9 3 0] [9 5 3] [5 7 3]]
[[0 7 8] [0 1 7] [1 5 7]]
[[1 5 3] [3 5 7]]
[[9 7 8] [9 5 7] [10 1 2]]
[[10 1 2] [9 5 0] [5 3 0] [5 7 3]]
[[8 0 2] [8 2 5] [8 5 7] [10 5 2]]
[[2 10 5] [2 5 3] [3 5 7]]
[[7 9 5] [7 8 9] [3 11 2]]
[[9 5 7] [9 7 2] [9 2 0] [2 7 11]]
[[2 3 11] [0 1 8] [1 7 8] [1 5 7]]
[[11 2 1] [11 1 7] [7 1 5]]
[[9 5 8] [8 5 7] [10 1 3] [10 3 11]]
[[5 7 0] [5 0 9] [7 11 0] [1 0 10] [11 10 0]]
[[11 10 0] [11 0 3] [10 5 0] [8 0 7] [5 7 0]]
[[11 10 5] [7 11 5]]
[[10 6 5]]
[[0 8 3] [5 10 6]]
[[9 0 1] [5 10 6]]
[[1 8 3] [1 9 8] [5 10 6]]
[[1 6 5] [2 6 1]]
[[1 6 5] [1 2 6] [3 0 8]]
[[9 6 5] [9 0 6] [0 2 6]]
[[5 9 8] [5 8 2] [5 2 6] [3 2 8]]
[[2 3 11] [10 6 5]]
[[11 0 8] [11 2 0] [10 6 5]]
[[0 1 9] [2 3 11] [5 10 6]]
[[5 10 6] [1 9 2] [9 11 2] [9 8 11]]
[[6 3 11] [6 5 3] [5 1 3]]
[[0 8 11] [0 11 5] [0 5 1] [5 11 6]]
[[3 11 6] [0 3 6] [0 6 5] [0 5 9]]
[[6 5 9] [6 9 11] [11 9 8]]
[[5 10 6] [4 7 8]]
[[4 3 0] [4 7 3] [6 5 10]]
[[1 9 0] [5 10 6] [8 4 7]]
[[10 6 5] [1 9 7] [1 7 3] [7 9 4]]
[[6 1 2] [6 5 1] [4 7 8]]
[[1 2 5] [5 2 6] [3 0 4] [3 4 7]]
[[8 4 7] [9 0 5] [0 6 5] [0 2 6]]
[[7 3 9] [7 9 4] [3 2 9] [5 9 6] [2 6 9]]
[[3 11 2] [7 8 4] [10 6 5]]
[[5 10 6] [4 7 2] [4 2 0] [2 7 11]]
[[0 1 9] [4 7 8] [2 3 11] [5 10 6]]
[[9 2 1] [9 11 2] [9 4 11] [7 11 4] [5 10 6]]
[[8 4 7] [3 11 5] [3 5 1] [5 11 6]]
[[5 1 11] [5 11 6] [1 0 11] [7 11 4] [0 4 11]]
[[0 5 9] [0 6 5] [0 3 6] [11 6 3] [8 4 7]]
[[6 5 9] [6 9 11] [4 7 9] [7 11 9]]
[[10 4 9] [6 4 10]]
[[4 10 6] [4 9 10] [0 8 3]]
[[10 0 1] [10 6 0] [6 4 0]]
[[8 3 1] [8 1 6] [8 6 4] [6 1 10]]
[[1 4 9] [1 2 4] [2 6 4]]
[[3 0 8] [1 2 9] [2 4 9] [2 6 4]]
[[0 2 4] [4 2 6]]
[[8 3 2] [8 2 4] [4 2 6]]
[[10 4 9] [10 6 4] [11 2 3]]
[[0 8 2] [2 8 11] [4 9 10] [4 10 6]]
[[3 11 2] [0 1 6] [0 6 4] [6 1 10]]
[[6 4 1] [6 1 10] [4 8 1] [2 1 11] [8 11 1]]
[[9 6 4] [9 3 6] [9 1 3] [11 6 3]]
[[8 11 1] [8 1 0] [11 6 1] [9 1 4] [6 4 1]]
[[3 11 6] [3 6 0] [0 6 4]]
[[6 4 8] [11 6 8]]
[[7 10 6] [7 8 10] [8 9 10]]
[[0 7 3] [0 10 7] [0 9 10] [6 7 10]]
[[10 6 7] [1 10 7] [1 7 8] [1 8 0]]
[[10 6 7] [10 7 1] [1 7 3]]
[[1 2 6] [1 6 8] [1 8 9] [8 6 7]]
[[2 6 9] [2 9 1] [6 7 9] [0 9 3] [7 3 9]]
[[7 8 0] [7 0 6] [6 0 2]]
[[7 3 2] [6 7 2]]
[[2 3 11] [10 6 8] [10 8 9] [8 6 7]]
[[2 0 7] [2 7 11] [0 9 7] [6 7 10] [9 10 7]]
[[1 8 0] [1 7 8] [1 10 7] [6 7 10] [2 3 11]]
[[11 2 1] [11 1 7] [10 6 1] [6 7 1]]
[[8 9 6] [8 6 7] [9 1 6] [11 6 3] [1 3 6]]
[[0 9 1] [11 6 7]]
[[7 8 0] [7 0 6] [3 11 0] [11 6 0]]
[[7 11 6]]
[[7 6 11]]
[[3 0 8] [11 7 6]]
[[0 1 9] [11 7 6]]
[[8 1 9] [8 3 1] [11 7 6]]
[[10 1 2] [6 11 7]]
[[1 2 10] [3 0 8] [6 11 7]]
[[2 9 0] [2 10 9] [6 11 7]]
[[6 11 7] [2 10 3] [10 8 3] [10 9 8]]
[[7 2 3] [6 2 7]]
[[7 0 8] [7 6 0] [6 2 0]]
[[2 7 6] [2 3 7] [0 1 9]]
[[1 6 2] [1 8 6] [1 9 8] [8 7 6]]
[[10 7 6] [10 1 7] [1 3 7]]
[[10 7 6] [1 7 10] [1 8 7] [1 0 8]]
[[0 3 7] [0 7 10] [0 10 9] [6 10 7]]
[[7 6 10] [7 10 8] [8 10 9]]
[[6 8 4] [11 8 6]]
[[3 6 11] [3 0 6] [0 4 6]]
[[8 6 11] [8 4 6] [9 0 1]]
[[9 4 6] [9 6 3] [9 3 1] [11 3 6]]
[[6 8 4] [6 11 8] [2 10 1]]
[[1 2 10] [3 0 11] [0 6 11] [0 4 6]]
[[4 11 8] [4 6 11] [0 2 9] [2 10 9]]
[[10 9 3] [10 3 2] [9 4 3] [11 3 6] [4 6 3]]
[[8 2 3] [8 4 2] [4 6 2]]
[[0 4 2] [4 6 2]]
[[1 9 0] [2 3 4] [2 4 6] [4 3 8]]
[[1 9 4] [1 4 2] [2 4 6]]
[[8 1 3] [8 6 1] [8 4 6] [6 10 1]]
[[10 1 0] [10 0 6] [6 0 4]]
[[4 6 3] [4 3 8] [6 10 3] [0 3 9] [10 9 3]]
[[10 9 4] [6 10 4]]
[[4 9 5] [7 6 11]]
[[0 8 3] [4 9 5] [11 7 6]]
[[5 0 1] [5 4 0] [7 6 11]]
[[11 7 6] [8 3 4] [3 5 4] [3 1 5]]
[[9 5 4] [10 1 2] [7 6 11]]
[[6 11 7] [1 2 10] [0 8 3] [4 9 5]]
[[7 6 11] [5 4 10] [4 2 10] [4 0 2]]
[[3 4 8] [3 5 4] [3 2 5] [10 5 2] [11 7 6]]
[[7 2 3] [7 6 2] [5 4 9]]
[[9 5 4] [0 8 6] [0 6 2] [6 8 7]]
[[3 6 2] [3 7 6] [1 5 0] [5 4 0]]
[[6 2 8] [6 8 7] [2 1 8] [4 8 5] [1 5 8]]
[[9 5 4] [10 1 6] [1 7 6] [1 3 7]]
[[1 6 10] [1 7 6] [1 0 7] [8 7 0] [9 5 4]]
[[4 0 10] [4 10 5] [0 3 10] [6 10 7] [3 7 10]]
[[7 6 10] [7 10 8] [5 4 10] [4 8 10]]
[[6 9 5] [6 11 9] [11 8 9]]
[[3 6 11] [0 6 3] [0 5 6] [0 9 5]]
[[0 11 8] [0 5 11] [0 1 5] [5 6 11]]
[[6 11 3] [6 3 5] [5 3 1]]
[[1 2 10] [9 5 11] [9 11 8] [11 5 6]]
[[0 11 3] [0 6 11] [0 9 6] [5 6 9] [1 2 10]]
[[11 8 5] [11 5 6] [8 0 5] [10 5 2] [0 2 5]]
[[6 11 3] [6 3 5] [2 10 3] [10 5 3]]
[[5 8 9] [5 2 8] [5 6 2] [3 8 2]]
[[9 5 6] [9 6 0] [0 6 2]]
[[1 5 8] [1 8 0] [5 6 8] [3 8 2] [6 2 8]]
[[1 5 6] [2 1 6]]
[[1 3 6] [1 6 10] [3 8 6] [5 6 9] [8 9 6]]
[[10 1 0] [10 0 6] [9 5 0] [5 6 0]]
[[0 3 8] [5 6 10]]
[[10 5 6]]
[[11 5 10] [7 5 11]]
[[11 5 10] [11 7 5] [8 3 0]]
[[5 11 7] [5 10 11] [1 9 0]]
[[10 7 5] [10 11 7] [9 8 1] [8 3 1]]
[[11 1 2] [11 7 1] [7 5 1]]
[[0 8 3] [1 2 7] [1 7 5] [7 2 11]]
[[9 7 5] [9 2 7] [9 0 2] [2 11 7]]
[[7 5 2] [7 2 11] [5 9 2] [3 2 8] [9 8 2]]
[[2 5 10] [2 3 5] [3 7 5]]
[[8 2 0] [8 5 2] [8 7 5] [10 2 5]]
[[9 0 1] [5 10 3] [5 3 7] [3 10 2]]
[[9 8 2] [9 2 1] [8 7 2] [10 2 5] [7 5 2]]
[[1 3 5] [3 7 5]]
[[0 8 7] [0 7 1] [1 7 5]]
[[9 0 3] [9 3 5] [5 3 7]]
[[9 8 7] [5 9 7]]
[[5 8 4] [5 10 8] [10 11 8]]
[[5 0 4] [5 11 0] [5 10 11] [11 3 0]]
[[0 1 9] [8 4 10] [8 10 11] [10 4 5]]
[[10 11 4] [10 4 5] [11 3 4] [9 4 1] [3 1 4]]
[[2 5 1] [2 8 5] [2 11 8] [4 5 8]]
[[0 4 11] [0 11 3] [4 5 11] [2 11 1] [5 1 11]]
[[0 2 5] [0 5 9] [2 11 5] [4 5 8] [11 8 5]]
[[9 4 5] [2 11 3]]
[[2 5 10] [3 5 2] [3 4 5] [3 8 4]]
[[5 10 2] [5 2 4] [4 2 0]]
[[3 10 2] [3 5 10] [3 8 5] [4 5 8] [0 1 9]]
[[5 10 2] [5 2 4] [1 9 2] [9 4 2]]
[[8 4 5] [8 5 3] [3 5 1]]
[[0 4 5] [1 0 5]]
[[8 4 5] [8 5 3] [9 0 5] [0 3 5]]
[[9 4 5]]
[[4 11 7] [4 9 11] [9 10 11]]
[[0 8 3] [4 9 7] [9 11 7] [9 10 11]]
[[1 10 11] [1 11 4] [1 4 0] [7 4 11]]
[[3 1 4] [3 4 8] [1 10 4] [7 4 11] [10 11 4]]
[[4 11 7] [9 11 4] [9 2 11] [9 1 2]]
[[9 7 4] [9 11 7] [9 1 11] [2 11 1] [0 8 3]]
[[11 7 4] [11 4 2] [2 4 0]]
[[11 7 4] [11 4 2] [8 3 4] [3 2 4]]
[[2 9 10] [2 7 9] [2 3 7] [7 4 9]]
[[9 10 7] [9 7 4] [10 2 7] [8 7 0] [2 0 7]]
[[3 7 10] [3 10 2] [7 4 10] [1 10 0] [4 0 10]]
[[1 10 2] [8 7 4]]
[[4 9 1] [4 1 7] [7 1 3]]
[[4 9 1] [4 1 7] [0 8 1] [8 7 1]]
[[4 0 3] [7 4 3]]
[[4 8 7]]
[[9 10 8] [10 11 8]]
[[3 0 9] [3 9 11] [11 9 10]]
[[0 1 10] [0 10 8] [8 10 11]]
[[3 1 10] [11 3 10]]
[[1 2 11] [1 11 9] [9 11 8]]
[[3 0 9] [3 9 11] [1 2 9] [2 11 9]]
[[0 2 11] [8 0 11]]
[[3 2 11]]
[[2 3 8] [2 8 10] [10 8 9]]
[[9 10 2] [0 9 2]]
[[2 3 8] [2 8 10] [0 1 8] [1 10 8]]
[[1 10 2]]
[[1 3 8] [9 1 8]]
[[0 9 1]]
[[0 3 8]]
[]
]
)
7 changes: 7 additions & 0 deletions src/emmy/mathbox/plot.clj
Original file line number Diff line number Diff line change
Expand Up @@ -727,3 +727,10 @@
(let [[f-bind opts] (mc/compile-3d opts :f 3)]
(-> (c/wrap [f-bind] ['emmy.mathbox.plot/VectorField opts])
(ev/fragment scene))))
(defn implicit-surface
[opts]
(let [[f-bind opts] (mc/compile-3d opts :f 3)]
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@sritchie -- I'm getting caught up in this function. I know that mc/compile-3d is not what I want. But I'm murky on what's going on here -- these look like functions that take the clojure functions, turn them into clojurescript functions and return a fragment that has all the necessary details for a mathbox plot.

I think I need to implement a compile-implicit function that will take the Lhs and Rhs of the implicit function and do a bunch of the same stuff as compile-3d. Could you walk me through compile-3d and probably the compile-state-fn therein? I don't have a great high-level idea.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

so what is going on here is that we have a clojure function on the clj side that we need to turn into something like (js/Function. "input" "body") for the cljs side.

I can definitely help with this; I think the easiest thing to do now is to just IGNORE this part and then pass in a quoted function, like '(fn [in out] ....) with whatever signature you need for your code to work.

Then once everything is humming we can add this part in, which will remove the need for quotes and make everything faster.

wdyt?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Willing to try! I don't really know the tradeoffs, it seems like the development path will be speedier, and I'll probably learn something in the process.

(-> (c/wrap [f-bind] ['emmy.mathbox.plot/ImplicitSurface opts])
(ev/fragment scene)
)
))
Loading