fractals.clj

(ns physics-demos.fractals
  (:require
   [thi.ng.geom.core :as g]
   [thi.ng.geom.core.vector :as v :refer [vec2 vec3]]
   [thi.ng.geom.core.matrix :as mat :refer [M32 M44]]
   [thi.ng.geom.circle :as c]
   [thi.ng.geom.spatialtree :as accel]
   [thi.ng.geom.svg.core :as svg]
   [thi.ng.geom.physics.core :as phys]
   [thi.ng.geom.webgl.animator :refer [animate]]
   [thi.ng.domus.core :as dom]
   [shodan.inspection :as shodan]))


(def pi (.-PI js/Math))
(def cos #(.cos js/Math %))
(def sin #(.sin js/Math %))
(def sqrt #(.sqrt js/Math %))
(def sign #(if (>= % 0) 1 -1))
(def abs #(if (>= % 0) % (- %)))

(defn add-functional-meta [f]
  (fn [& args]
    (with-meta (apply f args)
      {:is f
       :args args
       })))

(def my-partial (add-functional-meta partial))


(defn merged-juxt[fs] (comp (partial reduce into []) (apply juxt fs)))

(defn rotate [a [x y]]
  [(- (* (cos a) x) (* (sin a) y))
   (+ (* (sin a) x) (* (cos a) y))])

(defn add [[x0 y0] [x1 y1]]
  [(+ x0 x1) (+ y0 y1)])

(defn minus [[x0 y0] [x1 y1]]
  [(- x0 x1) (- y0 y1)])

(defn rotate-around [center a p]
  (add center (rotate a (minus p center))))

(defn cross-product [[x0 y0] [x1 y1]]
  (+ (* x0 x1) (* y0 y1)))

(def norm-sq #(cross-product % %))

(defn scalar-multiply [k [x y]]
  [(* k x) (* k y)])



;; TODO macro !!!
;; (def comp (add-functional-meta comp))
;; (def juxt (add-functional-meta juxt))
;; (def apply (add-functional-meta apply))

;; (def add (add-functional-meta add))
;; (def rotate (add-functional-meta rotate))
;; (def scalar-multiply (add-functional-meta scalar-multiply))

;; (def mapv (add-functional-meta mapv))
;; (def merged-juxt (add-functional-meta merged-juxt))



(defn eltwise_op [op [x0 y0] [x1 y1]]
  [(op x0 x1) (op y0 y1)])


(defn barycenter [ps]
  (scalar-multiply (/ 1 (count ps)) (reduce add [0. 0.] ps)))

(defn regular-polygon [n]
  (take n (iterate (partial rotate (/ (* 2 pi) n)) [1. 0.])))

(defn fractal-step [[step-f step-elts] current-elts]
  (into step-elts (step-f current-elts)))

(defn fractal [[init-elts step-params] details]
  (nth (iterate (partial fractal-step step-params) init-elts) details))

(defn close-polygon [ps] (conj (into [] ps) (first ps)))

(defn sierpinski-params [n]
  (let[step-elt (close-polygon (regular-polygon n))
      scale-f (partial scalar-multiply (/ 1 (dec n)))]
  (condp = n
    3 [[]
       [(merged-juxt (for [i [0 1 2]](partial mapv (partial mapv (comp (partial add (rotate (+ pi (* i 2 (/ pi 3))) [1. 0.])) scale-f)))))
        [step-elt]]]
    4 [[]
       [(merged-juxt (let [d [-1 0 1]](for [dx d
                                            dy d
                                            :when (not= 0 dx dy)]
                                        (partial mapv (partial mapv (comp (partial add (scalar-multiply (sqrt 2.) [dx dy])) scale-f))))))
        [(map (partial rotate (/ pi 4)) step-elt)]]])))

(defn tree-params [angles]
  (let[branch [0 1]
       ratio (/ (+ 1 (sqrt 5.)) 2.)]
    [[]
     [(merged-juxt (for [a angles](partial mapv (partial mapv (comp (partial add branch)
                                                                    (partial scalar-multiply (/ 1 ratio))
                                                                    (partial rotate a))))))
      [[[0. 0] branch]]]]))

(def koch-params [[[[-0.5 0][0.5 0]]]
                  [(merged-juxt (for [[v a] [[[(/ -1 3) 0] 0]
                                             [[(/ 1 3) 0] 0]
                                             [(rotate (/ pi 3) [(/ 1 6) 0]) (/ pi -3)]
                                             [(rotate (/ pi -3) [(/ -1 6) 0]) (/ pi 3)]]]
                                  (partial mapv (partial mapv (comp (partial add v) (partial rotate a) (partial scalar-multiply (/ 1 3)))))))
                   []]])

(def positive-infinity (.-POSITIVE_INFINITY js/Number))
(def negative-infinity (.-NEGATIVE_INFINITY js/Number))

(defn bounding-box [ps]
  (reduce (fn [[[x-min y-min][x-max y-max]] [x y]]
            [ [(min x-min x) (min y-min y)] [(max x-max x) (max y-max y)]])
          [[positive-infinity positive-infinity] [negative-infinity negative-infinity]]
          ps))

(defn merge-bb [[[x-min y-min][x-max y-max]] [[X-min Y-min][X-max Y-max]]]
  [[(min x-min X-min) (min y-min Y-min)] [(max x-max X-max) (max y-max Y-max)]])

(defn adjust-bounding-box [preserve-ratio target-box scene]
  (let [scene-bb (reduce merge-bb (map bounding-box scene))
        center-f (fn[box] (let[center (barycenter box)][center (apply minus (mapv #(minus % center) box))]))
        [target-c target-centered] (center-f target-box)
        [scene-c scene-centered] (center-f scene-bb)
        scaling (partial eltwise_op *
                         (let[tmp (eltwise_op / target-centered scene-centered)]
                           (if preserve-ratio
                             (scalar-multiply (apply min (map abs tmp)) (map sign tmp))
                             tmp)))]
    (fn[v] (-> v (minus scene-c) scaling (add target-c)))))

(defn display-fractal [id params details]
 (let [root (dom/by-id id)
        display-params {:width 300 :height 300}
        color "#000"
        draw #(svg/line-strip % )
        display #(dom/create-dom! (let [[w h] ((juxt :width :height) display-params)
                                        f (adjust-bounding-box true [[0 h][w 0]] %)]
                                    (svg/svg display-params (svg/group {:stroke color} (map (comp draw (partial mapv f)) %))))
                                  root)
        scene (fractal params details)]
   (display scene)))

(defn -main
  []
  (do
    (let[ f (my-partial add [0 1])
         ]
      (shodan/inspect (meta f)))
  (doall (for [details (range 1 5)
               params [(sierpinski-params 3)
                       (sierpinski-params 4)
                       koch-params
                       (tree-params [(/ pi 6) (/ pi -3)])
                       ]]
           (display-fractal (str "fractals-" details) params details)))
  )
  )

(-main)