[136] | 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 |
|
---|
| 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 |
|
---|
| 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 |
|
---|
| 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 |
|
---|
| 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 |
|
---|
| 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 |
|
---|
| 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 |
|
---|
| 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"
|
---|