From 97a10386531589581cc283f91d3883612b0f719a Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Sun, 29 Jun 2008 12:04:00 +0000 Subject: improve greatly mandelbrot demo: - much better coloring - determine max number of iterations and choice between float and double at runtime based on zoom level - do draft renderings with increasing resolution before final rendering --- demos/mandelbrot/CMakeLists.txt | 4 -- demos/mandelbrot/README | 12 +++-- demos/mandelbrot/mandelbrot.cpp | 115 +++++++++++++++++++++++++++------------- demos/mandelbrot/mandelbrot.h | 26 +++------ 4 files changed, 93 insertions(+), 64 deletions(-) diff --git a/demos/mandelbrot/CMakeLists.txt b/demos/mandelbrot/CMakeLists.txt index 1441334c4..b7d948e88 100644 --- a/demos/mandelbrot/CMakeLists.txt +++ b/demos/mandelbrot/CMakeLists.txt @@ -7,10 +7,6 @@ IF (CMAKE_COMPILER_IS_GNUCXX) ADD_DEFINITIONS ( "-DNDEBUG" ) ENDIF (CMAKE_COMPILER_IS_GNUCXX) -IF (DOUBLE_PRECISION) - ADD_DEFINITIONS ( "-DREAL=double" ) -ENDIF (DOUBLE_PRECISION) - INCLUDE_DIRECTORIES( ${QT_INCLUDE_DIR} ) SET(mandelbrot_SRCS diff --git a/demos/mandelbrot/README b/demos/mandelbrot/README index 0b615dea1..a451d6551 100644 --- a/demos/mandelbrot/README +++ b/demos/mandelbrot/README @@ -1,6 +1,10 @@ *** Mandelbrot demo *** -Demonstrates: -* coefficient-wise operations on a matrix to perform general array computations which - may not have anything to do with linear algebra -* impact of vectorization on performance +Controls: +* Left mouse button to center view at a point. +* Drag vertically with left mouse button to zoom in and out. + +Be sure to enable SSE2 or AltiVec to improve performance. + +The number of iterations, and the choice between single and double precision, are +determined at runtime depending on the zoom level. diff --git a/demos/mandelbrot/mandelbrot.cpp b/demos/mandelbrot/mandelbrot.cpp index 255a90516..59b68256c 100644 --- a/demos/mandelbrot/mandelbrot.cpp +++ b/demos/mandelbrot/mandelbrot.cpp @@ -15,76 +15,115 @@ void MandelbrotWidget::resizeEvent(QResizeEvent *) } } -void MandelbrotWidget::paintEvent(QPaintEvent *) +template int MandelbrotWidget::render(int max_iter, int img_width, int img_height) { - QTime time; - time.start(); - int alignedWidth = (width()/packetSize)*packetSize; - real yradius = xradius * height() / width(); - vector2 start(center.x() - xradius, center.y() - yradius); - vector2 step(2*xradius/width(), 2*yradius/height()); + enum { packetSize = Eigen::ei_packet_traits::size }; // number of reals in a Packet + typedef Eigen::Matrix Packet; // wrap a Packet as a vector + + int alignedWidth = (img_width/packetSize)*packetSize; + float yradius = xradius * img_height / img_width; + Eigen::Vector2f start(center.x() - xradius, center.y() - yradius); + Eigen::Vector2f step(2*xradius/img_width, 2*yradius/img_height); int pix = 0, total_iter = 0; - static float max_speed = 0; - for(int y = 0; y < height(); y++) + for(int y = 0; y < img_height; y++) { // for each pixel, we're going to do the iteration z := z^2 + c where z and c are complex numbers, // starting with z = c = complex coord of the pixel. pzi and pzr denote the real and imaginary parts of z. // pci and pcr denote the real and imaginary parts of c. - packet pzi_start, pci_start; + Packet pzi_start, pci_start; for(int i = 0; i < packetSize; i++) pzi_start[i] = pci_start[i] = start.y() + y * step.y(); for(int x = 0; x < alignedWidth; x += packetSize, pix += packetSize) { - packet pcr, pci = pci_start, pzr, pzi = pzi_start, pzr_buf; + Packet pcr, pci = pci_start, pzr, pzi = pzi_start, pzr_buf; for(int i = 0; i < packetSize; i++) pzr[i] = pcr[i] = start.x() + (x+i) * step.x(); // do the iterations. Every 4 iterations we check for divergence, in which case we can stop iterating. - int j; - for(j = 0; j < iter/4 && (pzr.cwiseAbs2() + pzi.cwiseAbs2()).eval().minCoeff() < 4; j++) + int j = 0; + typedef Eigen::Matrix Packeti; + Packeti pix_iter = Packeti::zero(), pix_dont_diverge; + do { - total_iter += 4 * packetSize; for(int i = 0; i < 4; i++) { pzr_buf = pzr; pzr = pzr.cwiseAbs2() - pzi.cwiseAbs2() + pcr; pzi = 2 * pzr_buf.cwiseProduct(pzi) + pci; } + pix_dont_diverge = (pzr.cwiseAbs2() + pzi.cwiseAbs2()) + .cwiseLessThan(Packet::constant(4)) + .template cast(); + pix_iter += 4 * pix_dont_diverge; + j++; + total_iter += 4 * packetSize; } + while(j < max_iter/4 && pix_dont_diverge.any()); // compute arbitrary pixel colors - packet pblue, pgreen; - if(j == iter/4) - { - packet pampl = (pzr.cwiseAbs2() + pzi.cwiseAbs2()); - pblue = real(510) * (packet::constant(0.1) + pampl).cwiseInverse().cwiseMin(packet::ones()); - pgreen = real(2550) * (packet::constant(10) + pampl).cwiseInverse().cwiseMin(packet::constant(0.1)); - } - else pblue = pgreen = packet::zero(); - for(int i = 0; i < packetSize; i++) { - buffer[4*(pix+i)] = (unsigned char)(pblue[i]); - buffer[4*(pix+i)+1] = (unsigned char)(pgreen[i]); + + buffer[4*(pix+i)] = float(pix_iter[i])*255/max_iter; + buffer[4*(pix+i)+1] = 0; buffer[4*(pix+i)+2] = 0; } } // if the width is not a multiple of packetSize, fill the remainder in black - for(int x = alignedWidth; x < width(); x++, pix++) + for(int x = alignedWidth; x < img_width; x++, pix++) buffer[4*pix] = buffer[4*pix+1] = buffer[4*pix+2] = 0; } + return total_iter; +} + +void MandelbrotWidget::paintEvent(QPaintEvent *) +{ + float resolution = xradius*2/width(); + int max_iter = 64; + if(resolution < 1e-4f) max_iter += 32 * ( 4 - std::log10(resolution)); + max_iter = (max_iter/4)*4; + int img_width = width()/draft; + int img_height = height()/draft; + static float max_speed = 0; + int total_iter; + bool single_precision = resolution > 1e-6f; + + QTime time; + time.start(); + if(single_precision) + total_iter = render(max_iter, img_width, img_height); + else + total_iter = render(max_iter, img_width, img_height); int elapsed = time.elapsed(); - float speed = elapsed ? float(total_iter)*1000/elapsed : 0; - max_speed = std::max(max_speed, speed); - std::cout << elapsed << " ms elapsed, " - << total_iter << " iters, " - << speed << " iters/s (max " << max_speed << ")" << std::endl; - QImage image(buffer, width(), height(), QImage::Format_RGB32); + if(draft == 1) + { + float speed = elapsed ? float(total_iter)*1000/elapsed : 0; + max_speed = std::max(max_speed, speed); + std::cout << elapsed << " ms elapsed, " + << total_iter << " iters, " + << speed << " iters/s (max " << max_speed << ")" << std::endl; + int packetSize = single_precision ? int(Eigen::ei_packet_traits::size) + : int(Eigen::ei_packet_traits::size); + setWindowTitle(QString("resolution ")+QString::number(xradius*2/width(), 'e', 2) + +(single_precision ? QString(", single ") : QString(", double ")) + +QString("precision, ") + +(packetSize==1 ? QString("no vectorization") + : QString("vectorized (%1 per packet)").arg(packetSize))); + } + + QImage image(buffer, img_width, img_height, QImage::Format_RGB32); QPainter painter(this); - painter.drawImage(QPointF(0,0), image); + painter.drawImage(QPoint(0, 0), image.scaled(width(), height())); + + if(draft>1) + { + draft /= 2; + setWindowTitle(QString("recomputing at 1/%1 resolution...").arg(draft)); + update(); + } } void MandelbrotWidget::mousePressEvent(QMouseEvent *event) @@ -92,9 +131,10 @@ void MandelbrotWidget::mousePressEvent(QMouseEvent *event) if( event->buttons() & Qt::LeftButton ) { lastpos = event->pos(); - real yradius = xradius * height() / width(); - center = vector2(center.x() + (event->pos().x() - width()/2) * xradius * 2 / width(), - center.y() + (event->pos().y() - height()/2) * yradius * 2 / height()); + float yradius = xradius * height() / width(); + center = Eigen::Vector2f(center.x() + (event->pos().x() - width()/2) * xradius * 2 / width(), + center.y() + (event->pos().y() - height()/2) * yradius * 2 / height()); + draft = 16; update(); } } @@ -105,10 +145,11 @@ void MandelbrotWidget::mouseMoveEvent(QMouseEvent *event) lastpos = event->pos(); if( event->buttons() & Qt::LeftButton ) { - real t = 1 + 3 * real(delta.y()) / height(); + float t = 1 + 5 * float(delta.y()) / height(); if(t < 0.5) t = 0.5; if(t > 2) t = 2; xradius *= t; + draft = 16; update(); } } diff --git a/demos/mandelbrot/mandelbrot.h b/demos/mandelbrot/mandelbrot.h index de09a1ffb..ea2ee6084 100644 --- a/demos/mandelbrot/mandelbrot.h +++ b/demos/mandelbrot/mandelbrot.h @@ -5,43 +5,31 @@ #include #include -#ifdef REAL -typedef REAL real; -#else -typedef float real; -#endif - -enum { packetSize = Eigen::ei_packet_traits::size }; // number of reals in a packet -typedef Eigen::Matrix packet; // wrap a packet as a vector -typedef Eigen::Matrix vector2; // really just a complex number, but we're here to demo Eigen ! - -const int iter = 32; // the maximum number of iterations done per pixel. Must be a multiple of 4. - class MandelbrotWidget : public QWidget { Q_OBJECT - vector2 center; - real xradius; + Eigen::Vector2f center; + float xradius; int size; unsigned char *buffer; QPoint lastpos; + int draft; protected: void resizeEvent(QResizeEvent *); void paintEvent(QPaintEvent *); void mousePressEvent(QMouseEvent *event); void mouseMoveEvent(QMouseEvent *event); + template int render(int max_iter, int resx, int resy); public: - MandelbrotWidget() : QWidget(), center(real(0),real(0)), xradius(2), - size(0), buffer(0) + MandelbrotWidget() : QWidget(), center(0,0), xradius(2), + size(0), buffer(0), draft(16) { setAutoFillBackground(false); - setWindowTitle(QString("Mandelbrot/Eigen, sizeof(real)=")+QString::number(sizeof(real)) - +", sizeof(packet)="+QString::number(sizeof(packet))); } ~MandelbrotWidget() { if(buffer) delete[]buffer; } }; -#endif // MANDELBROT_H \ No newline at end of file +#endif // MANDELBROT_H -- cgit v1.2.3