97
97
but from now on we don't need to care.
100
int interleaving_count = mandelbrot->renderThreadCount();
102
100
// how many iterations we do on each sample before declaring that it doesn't diverge
103
const int max_iter = mandelbrot->maxIter();
101
max_iter = mandelbrot->maxIter();
105
103
// precompute some constants that will be used in the coloring computations
106
const float log_max_iter = std::log(float(max_iter));
104
log_max_iter = std::log(float(max_iter));
108
105
if(mandelbrot->min_iter_divergence() != 0 && mandelbrot->min_iter_divergence() != max_iter)
110
tmin = std::log(Real(mandelbrot->min_iter_divergence())) / log_max_iter;
107
tmin = std::log(float(mandelbrot->min_iter_divergence())) / log_max_iter;
113
const float gamma = mandelbrot->gamma();
114
const float tshift_to_gamma = 0.3f;
115
const float tshift = std::pow(tshift_to_gamma, 1.f/gamma);
116
const float tshift_plus_1_to_gamma = std::pow(tshift+1.f, gamma);
111
square_bailout_radius = 20; // value found in papers on continuous escape time formulas, changing it can degrade the smoothness.
112
log_log_square_bailout_radius = std::log(std::log(square_bailout_radius));
118
const float log_of_2 = std::log(2.f);
119
const float log_of_2log2 = std::log(2.f*log_of_2);
114
log_of_2 = std::log(2.f);
115
log_of_2log2 = std::log(2.f*log_of_2);
121
Color3 rgb1, rgb2, rgb3, hsv1, hsv2, hsv3;
122
117
mandelbrot->color1().getRgbF(&rgb1[0], &rgb1[1], &rgb1[2]);
123
118
mandelbrot->color1().getHsvF(&hsv1[0], &hsv1[1], &hsv1[2]);
124
119
mandelbrot->color2().getRgbF(&rgb2[0], &rgb2[1], &rgb2[2]);
125
120
mandelbrot->color2().getHsvF(&hsv2[0], &hsv2[1], &hsv2[2]);
126
121
mandelbrot->color3().getRgbF(&rgb3[0], &rgb3[1], &rgb3[2]);
127
122
mandelbrot->color3().getHsvF(&hsv3[0], &hsv3[1], &hsv3[2]);
129
/****** Beginning of rendering *********/
131
// iterate over scanlines. Notice how each thread works on a different set of scanlines -- since each thread uses
132
// a different value for interleaving_number
133
for(int y = interleaving_number; y < image->height(); y += interleaving_count)
135
// pointer to the first pixel in the scanline to render
136
unsigned char *pixel = const_cast<unsigned char*>(static_cast<const QImage *>(image)->scanLine(y));
138
// for each pixel, we're going to do the iteration z := z^2 + c where z and c are complex numbers,
139
// starting with z = c = complex coord of the pixel. pzr and pzi denote the real and imaginary parts of z.
140
// pcr and pci denote the real and imaginary parts of c.
141
Packet pzi_start, pci_start;
142
for(int i = 0; i < packet_size; i++) pzi_start[i] = pci_start[i] = imaginary_start + y * resolution;
144
// iterate over pixels in the current scanline, by steps of 'packet_size'.
145
for(int x = 0; x < image->width(); x += packet_size)
147
if(mandelbrot->abortRenderingAsSoonAsPossible()) return;
149
// initial values before we start iterating
150
Packet pcr, pci = pci_start, pzr, pzi = pzi_start, pzr_buf;
151
for(int i = 0; i < packet_size; i++) pzr[i] = pcr[i] = real_start + (x+i) * resolution;
153
/* first step: do iterations by batches of 'iter_before_test'. Testing only every 'iter_before_test'
154
* iterations allows to go faster.
158
Packeti pixel_iter = Packeti::Zero(), // number of iteration per pixel in the packet
159
pixel_diverge = Packeti::Zero(); // whether or not each pixel has already diverged
160
int count_not_yet_diverged = packet_size;
162
Packet pzr_previous, pzi_previous,
163
pzr_before_diverge = Packet::Zero(), pzi_before_diverge = Packet::Zero();
169
for(int i = 0; i < iter_before_test/4; i++) // we peel the inner loop by 4
171
LowlevelPacket lpzr, lpzi;
172
for(int repeat = 0; repeat < 4; repeat++)
125
template<typename Real>
126
void mandelbrot_render_tile_impl<Real>::computePacket(int x, int y, Color3 *pixels)
128
// for each pixel, we're going to do the iteration z := z^2 + c where z and c are complex numbers,
129
// starting with z = c = complex coord of the pixel. pzr and pzi denote the real and imaginary parts of z.
130
// pcr and pci denote the real and imaginary parts of c.
131
Packet pcr, pci, pzr, pzi, pzr_buf;
133
for(int i = 0; i < packet_size; i++) {
134
pzi[i] = pci[i] = tile.affix().y() + y * resolution;
135
pzr[i] = pcr[i] = tile.affix().x() + (x+i) * resolution;
138
/* first step: do iterations by batches of 'iter_before_test'. Testing only every 'iter_before_test'
139
* iterations allows to go faster.
143
Packeti pixel_iter = Packeti::Zero(), // number of iteration per pixel in the packet
144
pixel_diverge = Packeti::Zero(); // whether or not each pixel has already diverged
145
int count_not_yet_diverged = packet_size;
147
Packet pzr_previous, pzi_previous,
148
pzr_before_diverge = Packet::Zero(), pzi_before_diverge = Packet::Zero();
154
/* perform iter_before_test iterations */
155
for(int i = 0; i < iter_before_test/4; i++) // we peel the inner loop by 4
157
LowlevelPacket lpzr, lpzi;
158
for(int repeat = 0; repeat < 4; repeat++)
160
lpzr = Eigen::ei_pload(pzr.data());
161
lpzi = Eigen::ei_pload(pzi.data());
162
Eigen::ei_pstore(pzr.data(),
165
Eigen::ei_pmul(lpzr,lpzr),
166
Eigen::ei_pmul(lpzi,lpzi)
168
Eigen::ei_pload(pcr.data())
171
Eigen::ei_pstore(pzi.data(),
174
Eigen::ei_padd(lpzr,lpzr),
177
Eigen::ei_pload(pci.data())
183
/* test for divergence */
184
pzabs2 = pzr.cwise().square();
185
pzabs2 += pzi.cwise().square();
186
for(int i = 0; i < packet_size; i++) {
187
if(!(pixel_diverge[i])) {
188
if(pzabs2[i] > square_bailout_radius) {
189
pixel_diverge[i] = 1;
190
pzr_before_diverge[i] = pzr_previous[i];
191
pzi_before_diverge[i] = pzi_previous[i];
192
count_not_yet_diverged--;
194
else pixel_iter[i] += iter_before_test;
198
j += iter_before_test;
200
while(j < max_iter && count_not_yet_diverged);
202
/* Second step: we know the iteration count before divergence for each pixel but only up to precision
203
* 'iter_before_test'. We now want to get the exact iteration count before divergence,
204
* so we iterate again starting from where we were before divergence, and now we test at every iteration.
207
pzr = pzr_before_diverge;
208
pzi = pzi_before_diverge;
209
pixel_diverge = Packeti::Zero();
210
count_not_yet_diverged = packet_size;
211
typedef Eigen::Matrix<float,packet_size,1> Packet_to_float;
212
Packet_to_float square_escape_modulus = Packet_to_float::Zero();
216
pzr = pzr.cwise().square();
217
pzr -= pzi.cwise().square();
219
pzi = (2*pzr_buf).cwise()*pzi;
221
pzabs2 = pzr.cwise().square();
222
pzabs2 += pzi.cwise().square();
223
for(int i = 0; i < packet_size; i++) {
224
if(!(pixel_diverge[i])) {
225
if(pzabs2[i] > square_bailout_radius) {
226
pixel_diverge[i] = 1;
227
square_escape_modulus[i] = (float)pzabs2[i];
228
count_not_yet_diverged--;
230
else pixel_iter[i]++;
235
while(j < iter_before_test && count_not_yet_diverged);
237
if(count_not_yet_diverged < packet_size) found_exterior_point = true;
239
/* Now we know exactly the number of iterations before divergence, and the escape modulus. */
241
/* Third step: compute pixel colors. */
243
for(int i = 0; i < packet_size; i++)
245
Real log_log_escape_modulus = Real(0);
246
if(square_escape_modulus[i] > Real(1))
248
Real log_escape_modulus = std::log(square_escape_modulus[i]);
249
if(log_escape_modulus > Real(1))
250
log_log_escape_modulus = std::log(log_escape_modulus);
252
Real normalized_iter_count = pixel_iter[i] + (log_log_square_bailout_radius - log_log_escape_modulus) / log_of_2;
253
Real log_normalized_iter_count = (normalized_iter_count > Real(1)) ? std::log(normalized_iter_count) : Real(0);
254
Real t = log_normalized_iter_count / log_max_iter;
256
// Now, remember that we did a little statistical analysis on some samples
257
// to determine roughly what would be the smallest count of iterations before divergence.
258
// At the beginning of the present function, we used that to compute 'tmin'.
259
// Now we use it to make the gradient actually start at t=tmin. Lower values of t just give black.
260
// In other words, we are ensuring that no matter what the viewpoint, the darkest part of the image
261
// will always be black (at least if our statistical analysis went well).
262
// In other words, we are trying to get optimal contrast. This is important as otherwise there could be
263
// almost no contrast at all and an interesting viewpoint could give an almost blank image.
264
t = (t-tmin)/(1.f-tmin);
265
t = CLAMP(t, Real(0), Real(1));
267
float threshold1 = 0.09f;
268
float threshold2 = 0.3f;
270
pixels[i] = (t/threshold1)*rgb3;
272
else if(t < threshold2) {
273
pixels[i] = mix(rgb2, hsv2, rgb3, hsv3, (t - threshold1) / (threshold2 - threshold1));
276
pixels[i] = mix(rgb1, hsv1, rgb2, hsv2, (t - threshold2) / (1.f - threshold2));
282
/** This is the main rendering function. It renders only the current tile of the Mandelbrot wallpaper.
284
* The image parameter is the image to render the tile to. The size of this image determines how many pixels we render.
285
* So to get supersampled rendering, just pass a larger image here, and then scale it down to the real size using smooth
288
* This function only writes to 'image', it does not at all write to mandelbrot->image() directly. The mandelbrot parameter
289
* is only used to determine the parameters of the rendering such as viewpoint, colors, etc.
291
template<typename Real> void mandelbrot_render_tile(
292
Mandelbrot *mandelbrot,
293
const MandelbrotTile& tile
296
// if we're compiling the path with SSE2 explicitly enabled, since it's only used on x86 (and not even on x86-64 since it has SSE2 by
297
// default), let's only compile the code if it's going to be used!
299
#if defined(HAVE_PATH_WITH_SSE2_EXPLICTLY_ENABLED) || !defined(THIS_PATH_WITH_SSE2_EXPLICTLY_ENABLED)
301
enum { packet_size = Eigen::ei_packet_traits<Real>::size };
302
Color3 dummy_buffer[packet_size];
304
mandelbrot_render_tile_impl<Real> renderer(mandelbrot, tile);
306
int supersampling = renderer.supersampling;
307
int supersampled_packet_size = supersampling * packet_size;
308
int tile_x = tile.destination().x();
309
int tile_y = tile.destination().y();
310
int tile_width = tile.destination().width();
311
int tile_height = tile.destination().height();
312
int supersampled_tile_width = supersampling * tile_width;
313
int supersampled_tile_height = supersampling * tile_height;
315
// first render a part of the border to check if the tile is probably entirely inside the interior of the Mandelbrot set
316
for(int y = 1; y < supersampled_tile_height-1; y+=4) // render every 4th pixel on the border
318
renderer.computePacket(0,y,dummy_buffer);
319
renderer.computePacket(supersampled_tile_width - packet_size, y, dummy_buffer);
321
if(mandelbrot->abortRenderingAsSoonAsPossible()) return;
323
for(int x = 0; x < supersampled_tile_width; x+= 4*packet_size)
325
renderer.computePacket(x,0,dummy_buffer);
326
renderer.computePacket(x,supersampled_tile_height-1,dummy_buffer);
328
if(mandelbrot->abortRenderingAsSoonAsPossible()) return;
330
// render the bottom-right packet: due to our rendering only every 4-th packet on the border, we could be
331
// missing this corner of the tile.
332
renderer.computePacket(supersampled_tile_width-packet_size,supersampled_tile_height-1,dummy_buffer);
334
// now, if the tile looks like it's entirely inside the interior, just assume that's the case,
335
// so fill it (not using a QPainter so as to avoid having to lock a mutex in case the image is being accessed
336
// by another QPainter in the GUI thread) and return.
337
if(!(renderer.found_exterior_point))
339
for(int y = 0; y < tile_height; y++)
341
for(int x = 0; x < tile_width; x++)
344
= const_cast<unsigned char*>(
345
static_cast<const QImage *>(mandelbrot->image())->scanLine(tile_y+y)
347
pixel[0] = mandelbrot->color1().blue();
348
pixel[1] = mandelbrot->color1().green();
349
pixel[2] = mandelbrot->color1().red();
356
// ok now do the actual rendering. not much point trying to reuse the part of the border we've already rendered,
357
// it's few pixels and it would take some nontrivial code.
359
qreal one_over_supersampling_squared = qreal(1) / (supersampling*supersampling);
361
for(int y = 0; y < tile_height; y++)
363
for(int x = 0; x < tile_width; x += packet_size)
365
Color3 supersampled_buffer[MAX_SUPERSAMPLING][packet_size * MAX_SUPERSAMPLING];
367
for(int y2 = 0; y2 < supersampling; y2++)
369
for(int x2 = 0; x2 < supersampled_packet_size; x2 += packet_size)
371
renderer.computePacket(supersampling*x+x2,supersampling*y+y2,&supersampled_buffer[y2][x2]);
373
if(mandelbrot->abortRenderingAsSoonAsPossible()) return;
377
int pixels_to_write = std::min(int(packet_size), tile_width-x);
378
for(int i = 0; i < pixels_to_write; i++)
380
Color3 color = Color3::Zero();
381
for(int y2 = 0; y2 < supersampling; y2++)
383
for(int x2 = 0; x2 < supersampling; x2++)
174
lpzr = Eigen::ei_pload(pzr.data());
175
lpzi = Eigen::ei_pload(pzi.data());
176
Eigen::ei_pstore(pzr.data(),
179
Eigen::ei_pmul(lpzr,lpzr),
180
Eigen::ei_pmul(lpzi,lpzi)
182
Eigen::ei_pload(pcr.data())
185
Eigen::ei_pstore(pzi.data(),
188
Eigen::ei_padd(lpzr,lpzr),
191
Eigen::ei_pload(pci.data())
196
pzabs2 = pzr.cwise().square();
197
pzabs2 += pzi.cwise().square();
198
for(int i = 0; i < packet_size; i++) {
199
if(!(pixel_diverge[i])) {
201
pixel_diverge[i] = 1;
202
pzr_before_diverge[i] = pzr_previous[i];
203
pzi_before_diverge[i] = pzi_previous[i];
204
count_not_yet_diverged--;
206
else pixel_iter[i] += iter_before_test;
211
while(j < max_iter/iter_before_test && count_not_yet_diverged);
213
/* Second step: we know the iteration count before divergence for each pixel but only up to precision
214
* 'iter_before_test'. We now want to get the exact iteration count before divergence,
215
* so we iterate again starting from where we were before divergence, and now we test at every iteration.
219
pzr = pzr_before_diverge;
220
pzi = pzi_before_diverge;
221
pixel_diverge = Packeti::Zero();
222
count_not_yet_diverged = packet_size;
223
typedef Eigen::Matrix<float,packet_size,1> Packet_to_float;
224
Packet_to_float escape_modulus = Packet_to_float::Zero();
228
pzr = pzr.cwise().square();
229
pzr -= pzi.cwise().square();
231
pzi = (2*pzr_buf).cwise()*pzi;
233
pzabs2 = pzr.cwise().square();
234
pzabs2 += pzi.cwise().square();
235
for(int i = 0; i < packet_size; i++) {
236
if(!(pixel_diverge[i])) {
238
pixel_diverge[i] = 1;
239
escape_modulus[i] = (float)pzabs2[i];
240
count_not_yet_diverged--;
242
else pixel_iter[i]++;
247
while(j < iter_before_test && count_not_yet_diverged);
249
/* Now we know exactly the number of iterations before divergence, and the escape modulus. */
252
/* Third step: compute pixel colors. */
254
int pixels_to_write = std::min(int(packet_size), image->width()-x);
256
Packet_to_float log_log_escape_modulus = escape_modulus.template cast<float>();
257
log_log_escape_modulus = log_log_escape_modulus.cwise().log().cwise().log();
259
Packet_to_float shift, normalized_iter_count;
261
for(int i = 0; i < packet_size; i++)
263
shift[i] = (-log_log_escape_modulus[i]+log_of_2log2)/log_of_2;
264
shift[i]=CLAMP(shift[i],-1.f,0.f);
265
normalized_iter_count[i] = pixel_iter[i] + shift[i];
266
if(normalized_iter_count[i]<=1.f) normalized_iter_count[i]=1.f;
269
Packet_to_float log_normalized_iter_count = normalized_iter_count.cwise().log();
272
t = log_normalized_iter_count / log_max_iter;
274
for(int i = 0; i < pixels_to_write; i++)
276
// Now, remember that in MandelbrotRenderThread, we did a little statistical analysis on some samples
277
// to determine roughly what would be the smallest count of iterations before divergence.
278
// At the beginning of the present function, we used that to compute 'tmin'.
279
// Now we use it to make the gradient actually start at t=tmin. Lower values of t just give black.
280
// In other words, we are ensuring that no matter what the viewpoint, the darkest part of the image
281
// will always be black (at least if our statistical analysis went well).
282
// In other words, we are trying to get optimal contrast. This is important as otherwise there could be
283
// almost no contrast at all and an interesting viewpoint could give an almost blank image.
284
if(t[i] <= tmin) t[i] = 0.f;
285
else t[i] = (t[i]-tmin)/(1.f-tmin);
286
t[i] = CLAMP(t[i], 0.f, 1.f);
289
Packet_to_float t_plus_tshift = t.cwise()+tshift;
290
Packet_to_float t_plus_tshift_to_gamma = (t_plus_tshift.cwise().log() * gamma).cwise().exp();
293
// Another homemade formula. It seems that some amount of gamma correction is beneficial.
294
// However it must be avoided for t near zero as x^gamma is non differentiable at x=0.
295
// So we just shift t a little bit to avoid 0.
296
t = (t_plus_tshift_to_gamma.cwise()-tshift_to_gamma) / (tshift_plus_1_to_gamma - tshift_to_gamma);
298
for(int i = 0; i < pixels_to_write; i++)
304
else if(t[i] < 0.75f) {
305
rgb = mix(rgb2, hsv2, rgb3, hsv3, 4*t[i] - 2.f);
308
rgb = mix(rgb1, hsv1, rgb2, hsv2, 4*t[i] - 3.f);
311
pixel[0] = (unsigned char)(255*rgb[2]);
312
pixel[1] = (unsigned char)(255*rgb[1]);
313
pixel[2] = (unsigned char)(255*rgb[0]);
385
color += supersampled_buffer[y2][x2 + i*supersampling];
388
color *= one_over_supersampling_squared;
390
= const_cast<unsigned char*>(
391
static_cast<const QImage *>(mandelbrot->image())->scanLine(tile_y+y)
393
pixel[0] = qreal_to_uchar_color_channel(color[2]);
394
pixel[1] = qreal_to_uchar_color_channel(color[1]);
395
pixel[2] = qreal_to_uchar_color_channel(color[0]);
401
Q_UNUSED(mandelbrot);
321
406
template void mandelbrot_render_tile<float>(
322
407
Mandelbrot *mandelbrot,
323
int interleaving_number,
408
const MandelbrotTile& tile
327
411
template void mandelbrot_render_tile<double>(
328
412
Mandelbrot *mandelbrot,
329
int interleaving_number,
413
const MandelbrotTile& tile
added
removed
wallpapers/mandelbrot/render_impl.cpp
