julia renderer in C; render worker in Haskell; test the score mechanism
[maximus:mandulia.git] / julia.c
1 #include <errno.h>
2 #include <limits.h>
3 #include <math.h>
4 #include <stdio.h>
5 #include <stdlib.h>
6 #include <string.h>
7
8 #include "julia.h"
9
10 #ifndef PRECISION
11 #define PRECISION 1
12 #endif
13
14 #if PRECISION == 1
15 typedef float R;
16 #define cos cosf
17 #define sqrt sqrtf
18 #define log logf
19 #define log2 log2f
20 #define fmin fminf
21 #define fmax fmaxf
22 #define SCAN "%f"
23 #else
24 typedef double R;
25 #define SCAN "%lf"
26 #endif
27
28 static const int channels = 4;
29 static const R escapeRadius = 65536.0;
30 static const R escapeRadius2 = 65536.0 * 65536.0;
31 static const int escapeLimit = 64;
32
33 static const R sqrt5 = 2.23606797749979;
34 static const R sqrt6 = 2.449489742783178;
35 static const R sqrt7 = 2.6457513110645907;
36
37 struct point {
38   unsigned char *pixel;
39   R x;
40   R y;
41 };
42
43 struct context {
44   struct point *points;
45   int width;
46   int height;
47 };
48
49 struct context *julia_new(int width, int height) {
50   if (width > 0 && height > 0) {
51     struct context *ctx = calloc(1, sizeof(struct context));
52     if (ctx) {
53       ctx->points = calloc(width * height, sizeof(struct point));
54       if (ctx->points) {
55         ctx->width = width;
56         ctx->height = height;
57         return ctx;
58       }
59       free(ctx);
60     }
61   }
62   return 0;
63 }
64
65 void julia_delete(struct context *ctx) {
66   if (ctx) {
67     if (ctx->points) {
68       free(ctx->points);
69     }
70     free(ctx);
71   }
72 }
73
74 static inline int min(int x, int y) {
75   return x < y ? x : y;
76 }
77
78 static inline int max(int x, int y) {
79   return x > y ? x : y;
80 }
81
82 static inline int colour(R a) {
83   int channel = 128 * cos(a) + 128;
84   return min(max(channel, 0), 255);
85 }
86
87 void julia(struct context *ctx, unsigned char *image, double dcx, double dcy) {
88   if (! ctx || ! image) {
89     return;
90   }
91   struct point *points = ctx->points;
92   const int width = ctx->width;
93   const int height = ctx->height;
94   const int stride = width * channels;
95   const R cx = dcx;
96   const R cy = dcy;
97   const R hw = height * 3.0 / width;
98   const R wh = width  * 3.0 / height;
99   const R sx = width < height ? 3 : wh;
100   const R sy = width < height ? hw : 3;
101   const int ntotal = width * height;
102   int npoints;
103   { // initialize
104     npoints = ntotal;
105     struct point *p = points;
106     for (int j = 0; j < height; ++j) {
107       R y = j * sy / height - sy/2;
108       for (int i = 0; i < width; ++i) {
109         R x = i * sx / width - sx/2;
110         p->pixel = image + j * stride + i * channels;
111         p->x = x;
112         p->y = y;
113         p->pixel[0] = 255;
114         p->pixel[1] = 255;
115         p->pixel[2] = 255;
116         p->pixel[3] = 255;
117         ++p;
118       }
119     }
120   }
121   int ndone = 0;
122   int n = 0;
123   int escapes;
124   do {
125     escapes = 0;
126     { // iterate
127       struct point *p = points;
128       for (int k = 0; k < npoints; ++k) {
129         R x = p->x;
130         R y = p->y;
131         int ok = 1;
132         for (int e = 0; ok && (e < escapeLimit); ++e) {
133           R zx2 = x * x - y * y;
134           R zy2 = 2 * x * y;
135           R zx = zx2 + cx;
136           R zy = zy2 + cy;
137           R zd2 = zx * zx + zy * zy;
138           if (zd2 >= escapeRadius2) {
139             R d = sqrt(zd2);
140             // colourify pixel
141             unsigned char *pixel = p->pixel;
142             R v = 0.025 * ((n+e) - log2(log2(d)/log2(escapeRadius)));
143             R a = fmin(fmax(8 * (v - 0.1) * v, 0), 1);
144             pixel[0] = a * colour(sqrt5 * v    );
145             pixel[1] = a * colour(sqrt6 * v + 1);
146             pixel[2] = a * colour(sqrt7 * v + 2);
147             pixel[3] = a * 255;
148             p->pixel = 0;
149             ++escapes;
150             ok = 0;
151           } else {
152             x = zx;
153             y = zy;
154           }
155         }
156         p->x = x;
157         p->y = y;
158         ++p;
159       }
160       n += escapeLimit;
161     }
162     { // compact memory
163       struct point *src = points;
164       struct point *dst = points;
165       int ncopied = 0;
166       for (int k = 0; k < npoints; ++k) {
167         if (! (src->pixel)) {
168           ++src;
169           ++ndone;
170         } else {
171           unsigned char *pixel = src->pixel;
172           R x = src->x;
173           R y = src->y;
174           dst->pixel = pixel;
175           dst->x = x;
176           dst->y = y;
177           ++src;
178           ++dst;
179           ++ncopied;
180         }
181       }
182       npoints = ncopied;
183     }
184   } while(escapes);
185 }