/
utils.clj
154 lines (135 loc) · 5.45 KB
/
utils.clj
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
(ns fastmath.fields.utils
(:require [fastmath.random :as r]
[fastmath.core :as m]
[fastmath.vector :as v])
(:import [fastmath.vector Vec2 Vec3]))
#_(set! *warn-on-reflection* true)
(set! *unchecked-math* :warn-on-boxed)
(m/use-primitive-operators)
(def unitx (Vec2. 1.0 0.0))
(def zerov (Vec2. 0.0 0.0))
(deftype Pair [a b])
(deftype Triplet [a b c])
(deftype Tetrad [a b c d])
(defn sdrand
"Symetric random from [-mx -mn] and [mn mx]"
^double [^double mn ^double mx]
(let [rand (r/drand mn mx)]
(r/randval rand (* -1.0 rand))))
(defn sirand
"Symmetric irand"
(^long [^long mx] (sirand 1 mx))
(^long [^long mn ^long mx]
(let [rand (r/lrand mn mx)]
(r/randval rand (* -1 rand)))))
(defn spow
^double [^double x ^double p]
(if (neg? x)
(- (m/pow (- x) p))
(m/pow x p)))
;;
(def offsets (vec (for [x (range -1 2) y (range -1 2)] (Vec2. x y))))
(defn closest
^long [P ^long n U]
(loop [i (long 0)
d2min Double/MAX_VALUE
j (long 0)]
(if (< i n)
(let [d2 (v/dist-sq (P i) U)
low? (< d2 d2min)]
(recur (inc i) (if low? d2 d2min) (if low? i j)))
j)))
(defn vratio
^double [^Vec2 P ^Vec2 Q ^Vec2 U]
(let [PmQ (v/sub P Q)]
(if (v/is-zero? PmQ)
1.0
(-> (v/sub U Q)
(v/emult PmQ)
(v/sum)
(/ (v/magsq PmQ))
(* 2.0)))))
(defn voronoi
^double [P ^long n ^long q U]
(let [Pq (P q)]
(loop [i (long 0)
ratiomax Double/MIN_VALUE]
(if (< i n)
(if (== i q)
(recur (inc i) ratiomax)
(recur (inc i) (max ratiomax (vratio (P i) Pq U))))
ratiomax))))
;;
(defrecord JacobiData [^double sn ^double cn ^double dn])
(defn- jacobi-data [^Vec3 v] (JacobiData. (.x v) (.y v) (.z v)))
;; order: sn,cn,dn
;; emmc = 1-k^2 (m=k^2, m=1-emmc)
(defn jacobi-elliptic
(^JacobiData [^double uu ^double emmc] (jacobi-elliptic uu emmc 8 0.003))
(^JacobiData [^double uu ^double emmc ^long iters ^double accuracy]
(jacobi-data (if-not (zero? emmc)
(let [bo? (neg? emmc)
^Vec3 emcud (if bo?
(let [d (- 1.0 emmc)
emc (/ (- emmc) d)
d (m/sqrt d)]
(Vec3. emc (* d uu) d))
(Vec3. emmc uu 0.0))
em (double-array iters)
en (double-array iters)
^Vec3 emccl (loop [a 1.0
l (long 0)
emc (.x emcud)]
(if (< l iters)
(do
(aset ^doubles em l a)
(let [emc (m/sqrt emc)]
(aset ^doubles en l emc)
(let [c (* 0.5 (+ a emc))]
(if (< (m/abs (- a emc)) (* accuracy a))
(Vec3. emc c l)
(recur c (inc l) (* a emc))))))
(Vec3. emc a (dec iters))))
u (* (.y emccl) (.y emcud))
sn (m/sin u)
cn (m/cos u)
^Vec3 scd (if-not (zero? sn)
(loop [dn 1.0
i (long (.z emccl))
a (/ cn sn)
c (* (.y emccl) a)]
(if-not (neg? i)
(let [b (aget ^doubles em i)
a (* c a)
c (* dn c)]
(recur (/ (+ (aget ^doubles en i) a)
(+ a b)) (dec i) (/ c b) c))
(let [a (/ (m/sqrt (inc (* c c))))
sn (if (neg? sn) (- a) a)]
(Vec3. sn (* c sn) dn))))
(Vec3. sn cn 1.0))]
(if bo?
(Vec3. (/ (.x scd) (.z emcud)) (.z scd) (.y scd))
scd))
(let [cn (/ (m/cosh uu))]
(Vec3. (m/tanh uu) cn cn))))))
;; Local noise generating functions (for various noise implementations)
(defn make-noise-variation
"Calculate shift by angle taken from noise"
([^double amount ^double scale noise-fn] (make-noise-variation amount scale 1.0 noise-fn))
([^double amount ^double scale ^double m noise-fn]
(let [m (* m/TWO_PI m)]
(fn [v]
(let [^Vec2 vv (v/mult v scale)
angle (* m ^double (noise-fn (.x vv) (.y vv)))]
(v/add v (Vec2. (* amount (m/cos angle))
(* amount (m/sin angle)))))))))
(defn make-noise-variation2
"Calculate shift by vector taken from noise"
[^double amount scale noise-fn]
(fn [v]
(let [^Vec2 vv (v/mult v scale)
x1 (- ^double (noise-fn (.x vv) (.y vv)) 0.5)
y1 (- ^double (noise-fn (.y vv) m/E (.x vv)) 0.5)]
(v/add v (Vec2. (* amount x1)
(* amount y1))))))