add mandelbrot demo

5 files changed, 204 insertions, 0 deletions

diff --git a/demos/CMakeLists.txt b/demos/CMakeLists.txt
new file mode 100644
index 000000000..b1e206096
--- /dev/null
+++ b/demos/CMakeLists.txt
@@ -0,0 +1,3 @@
+IF(BUILD_DEMOS)
+ADD_SUBDIRECTORY(mandelbrot)
+ENDIF(BUILD_DEMOS)
diff --git a/demos/mandelbrot/CMakeLists.txt b/demos/mandelbrot/CMakeLists.txt
new file mode 100644
index 000000000..1441334c4
--- /dev/null
+++ b/demos/mandelbrot/CMakeLists.txt
@@ -0,0 +1,24 @@
+FIND_PACKAGE(Qt4 REQUIRED)
+
+SET(CMAKE_INCLUDE_CURRENT_DIR ON)
+
+IF (CMAKE_COMPILER_IS_GNUCXX)
+   SET ( CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2")
+   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
+    mandelbrot.cpp
+)
+
+QT4_AUTOMOC(${mandelbrot_SRCS})
+
+ADD_EXECUTABLE(mandelbrot ${mandelbrot_SRCS})
+
+TARGET_LINK_LIBRARIES(mandelbrot ${QT_QTCORE_LIBRARY} ${QT_QTGUI_LIBRARY})
diff --git a/demos/mandelbrot/README b/demos/mandelbrot/README
new file mode 100644
index 000000000..0b615dea1
--- /dev/null
+++ b/demos/mandelbrot/README
@@ -0,0 +1,6 @@
+*** 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
diff --git a/demos/mandelbrot/mandelbrot.cpp b/demos/mandelbrot/mandelbrot.cpp
new file mode 100644
index 000000000..255a90516
--- /dev/null
+++ b/demos/mandelbrot/mandelbrot.cpp
@@ -0,0 +1,124 @@
+#include "mandelbrot.h"
+#include<QtGui/QPainter>
+#include<QtGui/QImage>
+#include<QtGui/QMouseEvent>
+#include<QtCore/QTime>
+
+void MandelbrotWidget::resizeEvent(QResizeEvent *)
+{
+  if(size < width() * height())
+  {
+    std::cout << "reallocate buffer" << std::endl;
+    size = width() * height();
+    if(buffer) delete[]buffer;
+    buffer = new unsigned char[4*size];
+  }
+}
+
+void MandelbrotWidget::paintEvent(QPaintEvent *)
+{
+  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;
+
+  for(int y = 0; y < 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;
+    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;
+      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++)
+      {
+        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;
+        }
+      }
+
+      // 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)+2] = 0;
+      }
+    }
+
+    // if the width is not a multiple of packetSize, fill the remainder in black
+    for(int x = alignedWidth; x < 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;
+
+  QImage image(buffer, width(), height(), QImage::Format_RGB32);
+  QPainter painter(this);
+  painter.drawImage(QPointF(0,0), image);
+}
+
+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());
+    update();
+  }
+}
+
+void MandelbrotWidget::mouseMoveEvent(QMouseEvent *event)
+{
+  QPoint delta = event->pos() - lastpos;
+  lastpos = event->pos();
+  if( event->buttons() & Qt::LeftButton )
+  {
+    real t = 1 + 3 * real(delta.y()) / height();
+    if(t < 0.5) t = 0.5;
+    if(t > 2) t = 2;
+    xradius *= t;
+    update();
+  }
+}
+
+int main(int argc, char *argv[])
+{
+  QApplication app(argc, argv);
+  MandelbrotWidget w;
+  w.show();
+  return app.exec();
+}
+
+#include "mandelbrot.moc"
diff --git a/demos/mandelbrot/mandelbrot.h b/demos/mandelbrot/mandelbrot.h
new file mode 100644
index 000000000..de09a1ffb
--- /dev/null
+++ b/demos/mandelbrot/mandelbrot.h
@@ -0,0 +1,47 @@
+#ifndef MANDELBROT_H
+#define MANDELBROT_H
+
+#include <Eigen/Array>
+#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;
+    int size;
+    unsigned char *buffer;
+    QPoint lastpos;
+
+  protected:
+    void resizeEvent(QResizeEvent *);
+    void paintEvent(QPaintEvent *);
+    void mousePressEvent(QMouseEvent *event);
+    void mouseMoveEvent(QMouseEvent *event);
+
+  public:
+    MandelbrotWidget() : QWidget(), center(real(0),real(0)), xradius(2),
+                         size(0), buffer(0)
+    {
+      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