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

@@ -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

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.

demos/mandelbrot/mandelbrot.cpp

@@ -15,76 +15,115 @@ void MandelbrotWidget::resizeEvent(QResizeEvent *)
   }
 }
 
-void MandelbrotWidget::paintEvent(QPaintEvent *)
+template<typename Real> 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());
-  int pix = 0, total_iter = 0;
-  static float max_speed = 0;
+  enum { packetSize = Eigen::ei_packet_traits<Real>::size }; // number of reals in a Packet
+  typedef Eigen::Matrix<Real, packetSize, 1> Packet; // wrap a Packet as a vector
 
-  for(int y = 0; y < height(); y++)
+  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;
+
+  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<int, packetSize, 1> 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<int>();
+        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;
   }
-  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;
+  return total_iter;
+}
 
-  QImage image(buffer, width(), height(), QImage::Format_RGB32);
+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<float>(max_iter, img_width, img_height);
+  else
+    total_iter = render<double>(max_iter, img_width, img_height);
+  int elapsed = time.elapsed();
+
+  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<float>::size)
+                                      : int(Eigen::ei_packet_traits<double>::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();
   }
 }

demos/mandelbrot/mandelbrot.h

@@ -5,43 +5,31 @@
 #include <QtGui/QApplication>
 #include <QtGui/QWidget>
 
-#ifdef REAL
-typedef REAL real;
-#else
-typedef float real;
-#endif
-
-enum { packetSize = Eigen::ei_packet_traits<real>::size }; // number of reals in a packet
-typedef Eigen::Matrix<real, packetSize, 1> packet; // wrap a packet as a vector
-typedef Eigen::Matrix<real, 2, 1> 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<typename Real> 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
+#endif // MANDELBROT_H