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 (at) 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" 214
