1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10 #include "mandelbrot.h"
11 #include <iostream>
12 #include<QtGui/QPainter>
13 #include<QtGui/QImage>
14 #include<QtGui/QMouseEvent>
15 #include<QtCore/QTime>
16
resizeEvent(QResizeEvent *)17 void MandelbrotWidget::resizeEvent(QResizeEvent *)
18 {
19 if(size < width() * height())
20 {
21 std::cout << "reallocate buffer" << std::endl;
22 size = width() * height();
23 if(buffer) delete[]buffer;
24 buffer = new unsigned char[4*size];
25 }
26 }
27
28 template<typename T> struct iters_before_test { enum { ret = 8 }; };
29 template<> struct iters_before_test<double> { enum { ret = 16 }; };
30
render(int img_width,int img_height)31 template<typename Real> void MandelbrotThread::render(int img_width, int img_height)
32 {
33 enum { packetSize = Eigen::internal::packet_traits<Real>::size }; // number of reals in a Packet
34 typedef Eigen::Array<Real, packetSize, 1> Packet; // wrap a Packet as a vector
35
36 enum { iters_before_test = iters_before_test<Real>::ret };
37 max_iter = (max_iter / iters_before_test) * iters_before_test;
38 const int alignedWidth = (img_width/packetSize)*packetSize;
39 unsigned char *const buffer = widget->buffer;
40 const double xradius = widget->xradius;
41 const double yradius = xradius * img_height / img_width;
42 const int threadcount = widget->threadcount;
43 typedef Eigen::Array<Real, 2, 1> Vector2;
44 Vector2 start(widget->center.x() - widget->xradius, widget->center.y() - yradius);
45 Vector2 step(2*widget->xradius/img_width, 2*yradius/img_height);
46 total_iter = 0;
47
48 for(int y = id; y < img_height; y += threadcount)
49 {
50 int pix = y * img_width;
51
52 // for each pixel, we're going to do the iteration z := z^2 + c where z and c are complex numbers,
53 // starting with z = c = complex coord of the pixel. pzi and pzr denote the real and imaginary parts of z.
54 // pci and pcr denote the real and imaginary parts of c.
55
56 Packet pzi_start, pci_start;
57 for(int i = 0; i < packetSize; i++) pzi_start[i] = pci_start[i] = start.y() + y * step.y();
58
59 for(int x = 0; x < alignedWidth; x += packetSize, pix += packetSize)
60 {
61 Packet pcr, pci = pci_start, pzr, pzi = pzi_start, pzr_buf;
62 for(int i = 0; i < packetSize; i++) pzr[i] = pcr[i] = start.x() + (x+i) * step.x();
63
64 // do the iterations. Every iters_before_test iterations we check for divergence,
65 // in which case we can stop iterating.
66 int j = 0;
67 typedef Eigen::Matrix<int, packetSize, 1> Packeti;
68 Packeti pix_iter = Packeti::Zero(), // number of iteration per pixel in the packet
69 pix_dont_diverge; // whether or not each pixel has already diverged
70 do
71 {
72 for(int i = 0; i < iters_before_test/4; i++) // peel the inner loop by 4
73 {
74 # define ITERATE \
75 pzr_buf = pzr; \
76 pzr = pzr.square(); \
77 pzr -= pzi.square(); \
78 pzr += pcr; \
79 pzi = (2*pzr_buf)*pzi; \
80 pzi += pci;
81 ITERATE ITERATE ITERATE ITERATE
82 }
83 pix_dont_diverge = ((pzr.square() + pzi.square())
84 .eval() // temporary fix as what follows is not yet vectorized by Eigen
85 <= Packet::Constant(4))
86 // the 4 here is not a magic value, it's a math fact that if
87 // the square modulus is >4 then divergence is inevitable.
88 .template cast<int>();
89 pix_iter += iters_before_test * pix_dont_diverge;
90 j++;
91 total_iter += iters_before_test * packetSize;
92 }
93 while(j < max_iter/iters_before_test && pix_dont_diverge.any()); // any() is not yet vectorized by Eigen
94
95 // compute pixel colors
96 for(int i = 0; i < packetSize; i++)
97 {
98 buffer[4*(pix+i)] = 255*pix_iter[i]/max_iter;
99 buffer[4*(pix+i)+1] = 0;
100 buffer[4*(pix+i)+2] = 0;
101 }
102 }
103
104 // if the width is not a multiple of packetSize, fill the remainder in black
105 for(int x = alignedWidth; x < img_width; x++, pix++)
106 buffer[4*pix] = buffer[4*pix+1] = buffer[4*pix+2] = 0;
107 }
108 return;
109 }
110
run()111 void MandelbrotThread::run()
112 {
113 setTerminationEnabled(true);
114 double resolution = widget->xradius*2/widget->width();
115 max_iter = 128;
116 if(resolution < 1e-4f) max_iter += 128 * ( - 4 - std::log10(resolution));
117 int img_width = widget->width()/widget->draft;
118 int img_height = widget->height()/widget->draft;
119 single_precision = resolution > 1e-7f;
120
121 if(single_precision)
122 render<float>(img_width, img_height);
123 else
124 render<double>(img_width, img_height);
125 }
126
paintEvent(QPaintEvent *)127 void MandelbrotWidget::paintEvent(QPaintEvent *)
128 {
129 static float max_speed = 0;
130 long long total_iter = 0;
131
132 QTime time;
133 time.start();
134 for(int th = 0; th < threadcount; th++)
135 threads[th]->start(QThread::LowPriority);
136 for(int th = 0; th < threadcount; th++)
137 {
138 threads[th]->wait();
139 total_iter += threads[th]->total_iter;
140 }
141 int elapsed = time.elapsed();
142
143 if(draft == 1)
144 {
145 float speed = elapsed ? float(total_iter)*1000/elapsed : 0;
146 max_speed = std::max(max_speed, speed);
147 std::cout << threadcount << " threads, "
148 << elapsed << " ms, "
149 << speed << " iters/s (max " << max_speed << ")" << std::endl;
150 int packetSize = threads[0]->single_precision
151 ? int(Eigen::internal::packet_traits<float>::size)
152 : int(Eigen::internal::packet_traits<double>::size);
153 setWindowTitle(QString("resolution ")+QString::number(xradius*2/width(), 'e', 2)
154 +QString(", %1 iterations per pixel, ").arg(threads[0]->max_iter)
155 +(threads[0]->single_precision ? QString("single ") : QString("double "))
156 +QString("precision, ")
157 +(packetSize==1 ? QString("no vectorization")
158 : QString("vectorized (%1 per packet)").arg(packetSize)));
159 }
160
161 QImage image(buffer, width()/draft, height()/draft, QImage::Format_RGB32);
162 QPainter painter(this);
163 painter.drawImage(QPoint(0, 0), image.scaled(width(), height()));
164
165 if(draft>1)
166 {
167 draft /= 2;
168 setWindowTitle(QString("recomputing at 1/%1 resolution...").arg(draft));
169 update();
170 }
171 }
172
mousePressEvent(QMouseEvent * event)173 void MandelbrotWidget::mousePressEvent(QMouseEvent *event)
174 {
175 if( event->buttons() & Qt::LeftButton )
176 {
177 lastpos = event->pos();
178 double yradius = xradius * height() / width();
179 center = Eigen::Vector2d(center.x() + (event->pos().x() - width()/2) * xradius * 2 / width(),
180 center.y() + (event->pos().y() - height()/2) * yradius * 2 / height());
181 draft = 16;
182 for(int th = 0; th < threadcount; th++)
183 threads[th]->terminate();
184 update();
185 }
186 }
187
mouseMoveEvent(QMouseEvent * event)188 void MandelbrotWidget::mouseMoveEvent(QMouseEvent *event)
189 {
190 QPoint delta = event->pos() - lastpos;
191 lastpos = event->pos();
192 if( event->buttons() & Qt::LeftButton )
193 {
194 double t = 1 + 5 * double(delta.y()) / height();
195 if(t < 0.5) t = 0.5;
196 if(t > 2) t = 2;
197 xradius *= t;
198 draft = 16;
199 for(int th = 0; th < threadcount; th++)
200 threads[th]->terminate();
201 update();
202 }
203 }
204
main(int argc,char * argv[])205 int main(int argc, char *argv[])
206 {
207 QApplication app(argc, argv);
208 MandelbrotWidget w;
209 w.show();
210 return app.exec();
211 }
212
213 #include "mandelbrot.moc"
214