1// Copyright 2019 The Marl Authors.
2//
3// Licensed under the Apache License, Version 2.0 (the "License");
4// you may not use this file except in compliance with the License.
5// You may obtain a copy of the License at
6//
7// https://www.apache.org/licenses/LICENSE-2.0
8//
9// Unless required by applicable law or agreed to in writing, software
10// distributed under the License is distributed on an "AS IS" BASIS,
11// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12// See the License for the specific language governing permissions and
13// limitations under the License.
14
15// This is an example application that uses Marl to parallelize the calculation
16// of a Julia fractal.
17
18#include "marl/defer.h"
19#include "marl/scheduler.h"
20#include "marl/thread.h"
21#include "marl/waitgroup.h"
22
23#include <fstream>
24
25#include <math.h>
26#include <stdint.h>
27
28// A color formed from a red, green and blue component.
29template <typename T>
30struct Color {
31 T r, g, b;
32
33 inline Color<T>& operator+=(const Color<T>& rhs) {
34 r += rhs.r;
35 g += rhs.g;
36 b += rhs.b;
37 return *this;
38 }
39
40 inline Color<T>& operator/=(T rhs) {
41 r /= rhs;
42 g /= rhs;
43 b /= rhs;
44 return *this;
45 }
46};
47
48// colorize returns a 'rainbow-color' for the scalar v.
49inline Color<float> colorize(float v) {
50 constexpr float PI = 3.141592653589793f;
51 constexpr float PI_2_THIRDS = 2.0f * PI / 3.0f;
52 return Color<float>{
53 0.5f + 0.5f * cosf(v + 0 * PI_2_THIRDS),
54 0.5f + 0.5f * cosf(v + 1 * PI_2_THIRDS),
55 0.5f + 0.5f * cosf(v + 2 * PI_2_THIRDS),
56 };
57}
58
59// lerp returns the linear interpolation between min and max using the weight x.
60inline float lerp(float x, float min, float max) {
61 return min + x * (max - min);
62}
63
64// julia calculates the Julia-set fractal value for the given coordinate and
65// constant. See https://en.wikipedia.org/wiki/Julia_set for more information.
66Color<float> julia(float x, float y, float cx, float cy) {
67 for (int i = 0; i < 1000; i++) {
68 if (x * x + y * y > 4) {
69 return colorize(sqrtf(static_cast<float>(i)));
70 }
71
72 auto xtemp = x * x - y * y;
73 y = 2 * x * y + cy;
74 x = xtemp + cx;
75 }
76
77 return {};
78}
79
80// writeBMP writes the given image as a bitmap to the given file, returning
81// true on success and false on error.
82bool writeBMP(const Color<uint8_t>* texels,
83 int width,
84 int height,
85 const char* path) {
86 auto file = fopen(path, "wb");
87 if (!file) {
88 fprintf(stderr, "Could not open file '%s'\n", path);
89 return false;
90 }
91 defer(fclose(file));
92
93 bool ok = true;
94 auto put1 = [&](uint8_t val) { ok = ok && fwrite(&val, 1, 1, file) == 1; };
95 auto put2 = [&](uint16_t val) { put1(static_cast<uint8_t>(val));
96 put1(static_cast<uint8_t>(val >> 8)); };
97 auto put4 = [&](uint32_t val) { put2(static_cast<uint16_t>(val));
98 put2(static_cast<uint16_t>(val >> 16)); };
99
100 const uint32_t padding = -(3 * width) & 3U; // in bytes
101 const uint32_t stride = 3 * width + padding; // in bytes
102 const uint32_t offset = 54;
103
104 // Bitmap file header
105 put1('B'); // header field
106 put1('M');
107 put4(offset + stride * height); // size in bytes
108 put4(0); // reserved
109 put4(offset);
110
111 // BITMAPINFOHEADER
112 put4(40); // size of header in bytes
113 put4(width); // width in pixels
114 put4(height); // height in pixels
115 put2(1); // number of color planes
116 put2(24); // bits per pixel
117 put4(0); // compression scheme (none)
118 put4(0); // size
119 put4(72); // horizontal resolution
120 put4(72); // vertical resolution
121 put4(0); // color pallete size
122 put4(0); // 'important colors' count
123
124 for (int y = height - 1; y >= 0; y--) {
125 for (int x = 0; x < width; x++) {
126 auto& texel = texels[x + y * width];
127 put1(texel.b);
128 put1(texel.g);
129 put1(texel.r);
130 }
131 for (uint32_t i = 0; i < padding; i++) {
132 put1(0);
133 }
134 }
135
136 return ok;
137}
138
139// Constants used for rendering the fractal.
140constexpr uint32_t imageWidth = 2048;
141constexpr uint32_t imageHeight = 2048;
142constexpr int samplesPerPixelW = 3;
143constexpr int samplesPerPixelH = 3;
144constexpr float windowMinX = -0.5f;
145constexpr float windowMaxX = +0.5f;
146constexpr float windowMinY = -0.5f;
147constexpr float windowMaxY = +0.5f;
148constexpr float cx = -0.8f;
149constexpr float cy = 0.156f;
150
151int main() {
152 // Create a marl scheduler using the full number of logical cpus.
153 // Bind this scheduler to the main thread so we can call marl::schedule()
154 marl::Scheduler scheduler(marl::Scheduler::Config::allCores());
155 scheduler.bind();
156 defer(scheduler.unbind()); // unbind before destructing the scheduler.
157
158 // Allocate the image.
159 auto pixels = new Color<uint8_t>[imageWidth * imageHeight];
160 defer(delete[] pixels); // free memory before returning.
161
162 // Create a wait group that will be used to synchronize the tasks.
163 // The wait group is constructed with an initial count of imageHeight as
164 // there will be a total of imageHeight tasks.
165 marl::WaitGroup wg(imageHeight);
166
167 // For each line of the image...
168 for (uint32_t y = 0; y < imageHeight; y++) {
169 // Schedule a task to calculate the image for this line.
170 // These may run concurrently across hardware threads.
171 marl::schedule([=] {
172 // Before this task returns, decrement the wait group counter.
173 // This is used to indicate that the task is done.
174 defer(wg.done());
175
176 for (uint32_t x = 0; x < imageWidth; x++) {
177 // Calculate the fractal pixel color.
178 Color<float> color = {};
179 // Take a number of sub-pixel samples.
180 for (int sy = 0; sy < samplesPerPixelH; sy++) {
181 auto fy = float(y) + (sy / float(samplesPerPixelH));
182 auto dy = float(fy) / float(imageHeight);
183 for (int sx = 0; sx < samplesPerPixelW; sx++) {
184 auto fx = float(x) + (sx / float(samplesPerPixelW));
185 auto dx = float(fx) / float(imageWidth);
186 color += julia(lerp(dx, windowMinX, windowMaxX),
187 lerp(dy, windowMinY, windowMaxY), cx, cy);
188 }
189 }
190 // Average the color.
191 color /= samplesPerPixelW * samplesPerPixelH;
192 // Write the pixel out to the image buffer.
193 pixels[x + y * imageWidth] = {static_cast<uint8_t>(color.r * 255),
194 static_cast<uint8_t>(color.g * 255),
195 static_cast<uint8_t>(color.b * 255)};
196 }
197 });
198 }
199
200 // Wait until all image lines have been calculated.
201 wg.wait();
202
203 // Write the image to "fractal.bmp".
204 if (!writeBMP(pixels, imageWidth, imageHeight, "fractal.bmp")) {
205 return 1;
206 }
207
208 // All done.
209 return 0;
210}
211