mandelbrot.cpp
Go to the documentation of this file.
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@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 
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;
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;
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 
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
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);
209  w.show();
210  return app.exec();
211 }
212 
213 #include "mandelbrot.moc"
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
gtsam.examples.DogLegOptimizerExample.int
int
Definition: DogLegOptimizerExample.py:111
MandelbrotThread::max_iter
int max_iter
Definition: mandelbrot.h:25
relicense.update
def update(text)
Definition: relicense.py:46
MandelbrotWidget::mousePressEvent
void mousePressEvent(QMouseEvent *event)
Definition: mandelbrot.cpp:173
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
MandelbrotThread::total_iter
long long total_iter
Definition: mandelbrot.h:24
gtsam::Vector2
Eigen::Vector2d Vector2
Definition: Vector.h:42
MandelbrotWidget::mouseMoveEvent
void mouseMoveEvent(QMouseEvent *event)
Definition: mandelbrot.cpp:188
Packet
internal::packet_traits< Scalar >::type Packet
Definition: benchmark-blocking-sizes.cpp:54
benchmark.elapsed
elapsed
Definition: benchmark.py:86
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
Eigen::internal::packet_traits
Definition: GenericPacketMath.h:106
MandelbrotWidget::buffer
unsigned char * buffer
Definition: mandelbrot.h:42
MandelbrotWidget::paintEvent
void paintEvent(QPaintEvent *)
Definition: mandelbrot.cpp:127
Eigen::Array
General-purpose arrays with easy API for coefficient-wise operations.
Definition: Array.h:45
buffer
Definition: pytypes.h:2270
MandelbrotWidget::resizeEvent
void resizeEvent(QResizeEvent *)
Definition: mandelbrot.cpp:17
Vector2::y
float y
Definition: test_operator_overloading.cpp:94
MandelbrotThread::render
void render(int img_width, int img_height)
Definition: mandelbrot.cpp:31
MandelbrotWidget::draft
int draft
Definition: mandelbrot.h:44
MandelbrotWidget::threads
MandelbrotThread ** threads
Definition: mandelbrot.h:45
iters_before_test
Definition: mandelbrot.cpp:28
mandelbrot.h
gtsam.examples.PlanarManipulatorExample.delta
def delta(g0, g1)
Definition: PlanarManipulatorExample.py:45
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
MandelbrotThread::widget
MandelbrotWidget * widget
Definition: mandelbrot.h:23
MandelbrotWidget::size
int size
Definition: mandelbrot.h:41
MandelbrotWidget::threadcount
int threadcount
Definition: mandelbrot.h:46
gtsam::utils.visual_isam.step
def step(data, isam, result, truth, currPoseIndex, isamArgs=())
Definition: visual_isam.py:82
Vector2
Definition: test_operator_overloading.cpp:18
Eigen::Triplet< double >
arg
EIGEN_DEVICE_FUNC const EIGEN_STRONG_INLINE ArgReturnType arg() const
Definition: ArrayCwiseUnaryOps.h:66
main
int main(int argc, char *argv[])
Definition: mandelbrot.cpp:205
time
#define time
Definition: timeAdaptAutoDiff.cpp:31
iters_before_test::ret
@ ret
Definition: mandelbrot.cpp:28
MandelbrotThread::single_precision
bool single_precision
Definition: mandelbrot.h:26
y
Scalar * y
Definition: level1_cplx_impl.h:124
tree::f
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
Definition: testExpression.cpp:218
MandelbrotWidget::center
Eigen::Vector2d center
Definition: mandelbrot.h:39
ITERATE
#define ITERATE
MandelbrotThread::run
void run()
Definition: mandelbrot.cpp:111
MandelbrotWidget
Definition: mandelbrot.h:34
gtsam.examples.DogLegOptimizerExample.float
float
Definition: DogLegOptimizerExample.py:113
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
Vector2::x
float x
Definition: test_operator_overloading.cpp:94
align_3::t
Point2 t(10, 10)
max
#define max(a, b)
Definition: datatypes.h:20
MandelbrotWidget::xradius
double xradius
Definition: mandelbrot.h:40
MandelbrotWidget::lastpos
QPoint lastpos
Definition: mandelbrot.h:43
eval
internal::nested_eval< T, 1 >::type eval(const T &xpr)
Definition: sparse_permutations.cpp:38
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
log10
const EIGEN_DEVICE_FUNC Log10ReturnType log10() const
Definition: ArrayCwiseUnaryOps.h:158


gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:02:44