summaryrefslogtreecommitdiffhomepage
path: root/eigen/demos/mandelbrot
diff options
context:
space:
mode:
Diffstat (limited to 'eigen/demos/mandelbrot')
-rw-r--r--eigen/demos/mandelbrot/CMakeLists.txt21
-rw-r--r--eigen/demos/mandelbrot/README10
-rw-r--r--eigen/demos/mandelbrot/mandelbrot.cpp213
-rw-r--r--eigen/demos/mandelbrot/mandelbrot.h71
4 files changed, 315 insertions, 0 deletions
diff --git a/eigen/demos/mandelbrot/CMakeLists.txt b/eigen/demos/mandelbrot/CMakeLists.txt
new file mode 100644
index 0000000..5c500e0
--- /dev/null
+++ b/eigen/demos/mandelbrot/CMakeLists.txt
@@ -0,0 +1,21 @@
+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)
+
+include_directories( ${QT_INCLUDE_DIR} )
+
+set(mandelbrot_SRCS
+ mandelbrot.cpp
+)
+
+qt4_automoc(${mandelbrot_SRCS})
+
+add_executable(mandelbrot ${mandelbrot_SRCS})
+add_dependencies(demos mandelbrot)
+
+target_link_libraries(mandelbrot ${QT_QTCORE_LIBRARY} ${QT_QTGUI_LIBRARY})
diff --git a/eigen/demos/mandelbrot/README b/eigen/demos/mandelbrot/README
new file mode 100644
index 0000000..a451d65
--- /dev/null
+++ b/eigen/demos/mandelbrot/README
@@ -0,0 +1,10 @@
+*** Mandelbrot demo ***
+
+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/eigen/demos/mandelbrot/mandelbrot.cpp b/eigen/demos/mandelbrot/mandelbrot.cpp
new file mode 100644
index 0000000..5d575d5
--- /dev/null
+++ b/eigen/demos/mandelbrot/mandelbrot.cpp
@@ -0,0 +1,213 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "mandelbrot.h"
+#include <iostream>
+#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];
+ }
+}
+
+template<typename T> struct iters_before_test { enum { ret = 8 }; };
+template<> struct iters_before_test<double> { enum { ret = 16 }; };
+
+template<typename Real> void MandelbrotThread::render(int img_width, int img_height)
+{
+ enum { packetSize = Eigen::internal::packet_traits<Real>::size }; // number of reals in a Packet
+ typedef Eigen::Array<Real, packetSize, 1> Packet; // wrap a Packet as a vector
+
+ enum { iters_before_test = iters_before_test<Real>::ret };
+ max_iter = (max_iter / iters_before_test) * iters_before_test;
+ const int alignedWidth = (img_width/packetSize)*packetSize;
+ unsigned char *const buffer = widget->buffer;
+ const double xradius = widget->xradius;
+ const double yradius = xradius * img_height / img_width;
+ const int threadcount = widget->threadcount;
+ typedef Eigen::Array<Real, 2, 1> Vector2;
+ Vector2 start(widget->center.x() - widget->xradius, widget->center.y() - yradius);
+ Vector2 step(2*widget->xradius/img_width, 2*yradius/img_height);
+ total_iter = 0;
+
+ for(int y = id; y < img_height; y += threadcount)
+ {
+ int pix = y * img_width;
+
+ // 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 iters_before_test iterations we check for divergence,
+ // in which case we can stop iterating.
+ int j = 0;
+ typedef Eigen::Matrix<int, packetSize, 1> Packeti;
+ Packeti pix_iter = Packeti::Zero(), // number of iteration per pixel in the packet
+ pix_dont_diverge; // whether or not each pixel has already diverged
+ do
+ {
+ for(int i = 0; i < iters_before_test/4; i++) // peel the inner loop by 4
+ {
+# define ITERATE \
+ pzr_buf = pzr; \
+ pzr = pzr.square(); \
+ pzr -= pzi.square(); \
+ pzr += pcr; \
+ pzi = (2*pzr_buf)*pzi; \
+ pzi += pci;
+ ITERATE ITERATE ITERATE ITERATE
+ }
+ pix_dont_diverge = ((pzr.square() + pzi.square())
+ .eval() // temporary fix as what follows is not yet vectorized by Eigen
+ <= Packet::Constant(4))
+ // the 4 here is not a magic value, it's a math fact that if
+ // the square modulus is >4 then divergence is inevitable.
+ .template cast<int>();
+ pix_iter += iters_before_test * pix_dont_diverge;
+ j++;
+ total_iter += iters_before_test * packetSize;
+ }
+ while(j < max_iter/iters_before_test && pix_dont_diverge.any()); // any() is not yet vectorized by Eigen
+
+ // compute pixel colors
+ for(int i = 0; i < packetSize; i++)
+ {
+ buffer[4*(pix+i)] = 255*pix_iter[i]/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 < img_width; x++, pix++)
+ buffer[4*pix] = buffer[4*pix+1] = buffer[4*pix+2] = 0;
+ }
+ return;
+}
+
+void MandelbrotThread::run()
+{
+ setTerminationEnabled(true);
+ double resolution = widget->xradius*2/widget->width();
+ max_iter = 128;
+ if(resolution < 1e-4f) max_iter += 128 * ( - 4 - std::log10(resolution));
+ int img_width = widget->width()/widget->draft;
+ int img_height = widget->height()/widget->draft;
+ single_precision = resolution > 1e-7f;
+
+ if(single_precision)
+ render<float>(img_width, img_height);
+ else
+ render<double>(img_width, img_height);
+}
+
+void MandelbrotWidget::paintEvent(QPaintEvent *)
+{
+ static float max_speed = 0;
+ long long total_iter = 0;
+
+ QTime time;
+ time.start();
+ for(int th = 0; th < threadcount; th++)
+ threads[th]->start(QThread::LowPriority);
+ for(int th = 0; th < threadcount; th++)
+ {
+ threads[th]->wait();
+ total_iter += threads[th]->total_iter;
+ }
+ 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 << threadcount << " threads, "
+ << elapsed << " ms, "
+ << speed << " iters/s (max " << max_speed << ")" << std::endl;
+ int packetSize = threads[0]->single_precision
+ ? int(Eigen::internal::packet_traits<float>::size)
+ : int(Eigen::internal::packet_traits<double>::size);
+ setWindowTitle(QString("resolution ")+QString::number(xradius*2/width(), 'e', 2)
+ +QString(", %1 iterations per pixel, ").arg(threads[0]->max_iter)
+ +(threads[0]->single_precision ? QString("single ") : QString("double "))
+ +QString("precision, ")
+ +(packetSize==1 ? QString("no vectorization")
+ : QString("vectorized (%1 per packet)").arg(packetSize)));
+ }
+
+ QImage image(buffer, width()/draft, height()/draft, QImage::Format_RGB32);
+ QPainter painter(this);
+ 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)
+{
+ if( event->buttons() & Qt::LeftButton )
+ {
+ lastpos = event->pos();
+ double yradius = xradius * height() / width();
+ center = Eigen::Vector2d(center.x() + (event->pos().x() - width()/2) * xradius * 2 / width(),
+ center.y() + (event->pos().y() - height()/2) * yradius * 2 / height());
+ draft = 16;
+ for(int th = 0; th < threadcount; th++)
+ threads[th]->terminate();
+ update();
+ }
+}
+
+void MandelbrotWidget::mouseMoveEvent(QMouseEvent *event)
+{
+ QPoint delta = event->pos() - lastpos;
+ lastpos = event->pos();
+ if( event->buttons() & Qt::LeftButton )
+ {
+ double t = 1 + 5 * double(delta.y()) / height();
+ if(t < 0.5) t = 0.5;
+ if(t > 2) t = 2;
+ xradius *= t;
+ draft = 16;
+ for(int th = 0; th < threadcount; th++)
+ threads[th]->terminate();
+ update();
+ }
+}
+
+int main(int argc, char *argv[])
+{
+ QApplication app(argc, argv);
+ MandelbrotWidget w;
+ w.show();
+ return app.exec();
+}
+
+#include "mandelbrot.moc"
diff --git a/eigen/demos/mandelbrot/mandelbrot.h b/eigen/demos/mandelbrot/mandelbrot.h
new file mode 100644
index 0000000..a687fd0
--- /dev/null
+++ b/eigen/demos/mandelbrot/mandelbrot.h
@@ -0,0 +1,71 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef MANDELBROT_H
+#define MANDELBROT_H
+
+#include <Eigen/Core>
+#include <QtGui/QApplication>
+#include <QtGui/QWidget>
+#include <QtCore/QThread>
+
+class MandelbrotWidget;
+
+class MandelbrotThread : public QThread
+{
+ friend class MandelbrotWidget;
+ MandelbrotWidget *widget;
+ long long total_iter;
+ int id, max_iter;
+ bool single_precision;
+
+ public:
+ MandelbrotThread(MandelbrotWidget *w, int i) : widget(w), id(i) {}
+ void run();
+ template<typename Real> void render(int img_width, int img_height);
+};
+
+class MandelbrotWidget : public QWidget
+{
+ Q_OBJECT
+
+ friend class MandelbrotThread;
+ Eigen::Vector2d center;
+ double xradius;
+ int size;
+ unsigned char *buffer;
+ QPoint lastpos;
+ int draft;
+ MandelbrotThread **threads;
+ int threadcount;
+
+ protected:
+ void resizeEvent(QResizeEvent *);
+ void paintEvent(QPaintEvent *);
+ void mousePressEvent(QMouseEvent *event);
+ void mouseMoveEvent(QMouseEvent *event);
+
+ public:
+ MandelbrotWidget() : QWidget(), center(0,0), xradius(2),
+ size(0), buffer(0), draft(16)
+ {
+ setAutoFillBackground(false);
+ threadcount = QThread::idealThreadCount();
+ threads = new MandelbrotThread*[threadcount];
+ for(int th = 0; th < threadcount; th++) threads[th] = new MandelbrotThread(this, th);
+ }
+ ~MandelbrotWidget()
+ {
+ if(buffer) delete[]buffer;
+ for(int th = 0; th < threadcount; th++) delete threads[th];
+ delete[] threads;
+ }
+};
+
+#endif // MANDELBROT_H
generated by cgit v1.2.3 (git 2.39.1) at 2025-04-07 15:02:20 +0000