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
17void 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
28template<typename T> struct iters_before_test { enum { ret = 8 }; };
29template<> struct iters_before_test<double> { enum { ret = 16 }; };
30
31template<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
111void 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
127void 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
173void 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
188void 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
205int 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