-
Notifications
You must be signed in to change notification settings - Fork 13
/
mandelbrot.fut
115 lines (97 loc) · 3.93 KB
/
mandelbrot.fut
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
-- Port of Accelerate's Mandelbrot example.
import "lib/github.com/diku-dk/complex/complex"
import "lib/github.com/athas/matte/colour"
module complex = mk_complex f32
type complex = complex.complex
def dot (c: complex): f32 =
let (r, i) = (complex.re c, complex.im c)
in (r * r + i * i)
def divergence (limit: i32) (radius: f32) (c0: complex): (complex,i32) =
loop (c, i) = (c0, 0) while i < limit && dot c < radius do
(c0 complex.+ c complex.* c,
i + 1)
def mandelbrot (screenX: i64) (screenY: i64) (limit: i32) (radius: f32)
((xmin, ymin, xmax, ymax): (f32,f32,f32,f32))
: [screenX][screenY](complex,i32) =
let sizex = (xmax - xmin)
let sizey = (ymax - ymin)
in map (\x ->
map (\y ->
let c0 = complex.mk (xmin + (f32.i64 x * sizex) / f32.i64 screenX)
(ymin + (f32.i64 y * sizey) / f32.i64 screenY)
in divergence limit radius c0)
(iota screenY))
(iota screenX)
-- Remaining code is about how to turn a divergence into a pretty
-- colour. That is, the most important part!
-- | Cubic interpolation.
def cubic (x0:f32,x1:f32) (y0:f32,y1:f32) (m0:f32,m1:f32) (x:f32) =
let h = x1 - x0
let t = (x - x0) / h
let h_00 = (1.0 + 2.0*t) * (1.0 - t) ** 2.0
let h_10 = t * (1.0 - t) ** 2.0
let h_01 = t ** 2.0 * (3.0 - 2.0 * t)
let h_11 = t ** 2.0 * (t - 1.0)
in y0 * h_00 + h * m0 * h_10 + y1 * h_01 + h * m1 * h_11
def interp (x0:f32, x1:f32)
(y0:argb.colour, y1:argb.colour)
((mr0,mg0,mb0):(f32,f32,f32),
(mr1,mg1,mb1):(f32,f32,f32))
(x: f32) =
let (r0,g0,b0,_) = argb.to_rgba y0
let (r1,g1,b1,_) = argb.to_rgba y1
in argb.from_rgba (cubic (x0,x1) (r0,r1) (mr0,mr1) x)
(cubic (x0,x1) (g0,g1) (mg0,mg1) x)
(cubic (x0,x1) (b0,b1) (mb0,mb1) x)
0.0
-- the ultraPalette from Accelerate.
def mk_palette (points: i32) (ix: i32): argb.colour =
let p = r32 ix / r32 points
let p0 = 0.0
let p1 = 0.16
let p2 = 0.42
let p3 = 0.6425
let p4 = 0.8575
let p5 = 1.0
let rgb8 (r: i32) (g: i32) (b: i32) =
argb.from_rgba (r32 r / 255.0) (r32 g / 255.0) (r32 b / 255.0) 0.0
let c0 = rgb8 0 7 100
let c1 = rgb8 32 107 203
let c2 = rgb8 237 255 255
let c3 = rgb8 255 170 0
let c4 = rgb8 0 2 0
let c5 = c0
let m0 = (0.7843138, 2.4509804, 2.52451)
let m1 = (1.93816, 2.341629, 1.6544118)
let m2 = (1.7046283, 0.0, 0.0)
let m3 = (0.0, -2.2812111, 0.0)
let m4 = (0.0, 0.0, 0.0)
let m5 = m0
in
if p <= p1 then interp (p0,p1) (c0,c1) (m0,m1) p else
if p <= p2 then interp (p1,p2) (c1,c2) (m1,m2) p else
if p <= p3 then interp (p2,p3) (c2,c3) (m2,m3) p else
if p <= p4 then interp (p3,p4) (c3,c4) (m3,m4) p else
interp (p4,p5) (c4,c5) (m4,m5) p
def log2 (x: f32) = f32.log x / f32.log 2.0
def escape_to_colour (limit: i32) (points: i32)
(z: complex, n: i32): argb.colour =
if limit == n then argb.black
else let smooth = log2 (log2 (complex.mag z))
let scale = 256.0
let shift = 1664.0
let ix = t32 (f32.sqrt (r32 n + 1.0 - smooth) * scale + shift)
in mk_palette points (ix %% points)
def mandelbrot_render (screenX: i64) (screenY: i64)
(xcentre: f32) (ycentre: f32) (width: f32)
(limit: i32) (radius: f32): [screenX][screenY]u32 =
let aspect_ratio = (f32.i64 screenX / f32.i64 screenY)
let (xmin,ymin) = ((xcentre - width/2),
(ycentre - (1.0/aspect_ratio)*width/2.0))
let (xmax,ymax) = ((xcentre + width/2),
(ycentre + (1.0/aspect_ratio)*width/2.0))
let escapes = mandelbrot screenX screenY limit radius (xmin, ymin, xmax, ymax)
let points = 2048
in map (\row -> map (escape_to_colour limit points) row) escapes
def main (n: i64) =
mandelbrot_render n n (-0.7f32) 0f32 3.067f32 100 16f32