root/examples/mandelbrot.cpp

DEFINITIONS

1. pt
2. start
3. process
4. add_neighbours
5. BOOST_AUTO_TEST_SUITE

   1 /******************************************************************************
   2 ** $Header: svn+ssh://jmmcg@svn.code.sf.net/p/libjmmcg/code/trunk/libjmmcg/examples/mandelbrot.cpp 2185 2017-10-13 10:14:17Z jmmcg $
   3 **
   4 ** Copyright (C) 2002 by J.M.McGuiness, coder@hussar.me.uk
   5 **
   6 ** This library is free software; you can redistribute it and/or
   7 ** modify it under the terms of the GNU Lesser General Public
   8 ** License as published by the Free Software Foundation; either
   9 ** version 2.1 of the License, or (at your option) any later version.
  10 **
  11 ** This library is distributed in the hope that it will be useful,
  12 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  14 ** Lesser General Public License for more details.
  15 **
  16 ** You should have received a copy of the GNU Lesser General Public
  17 ** License along with this library; if not, write to the Free Software
  18 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  19 */
  20 
  21 #include "stdafx.h"
  22 
  23 #define BOOST_TEST_MODULE libjmmcg_tests
  24 #include <boost/test/included/unit_test.hpp>
  25 
  26 #include "core/thread_pool_workers.hpp"
  27 
  28 #include <cassert>
  29 #include <stdexcept>
  30 
  31 #include <boost/bind.hpp>
  32 #include <boost/mpl/assert.hpp>
  33 #include <boost/static_assert.hpp>
  34 #include <boost/type_traits.hpp>
  35 #include <boost/scoped_ptr.hpp>
  36 #include <boost/shared_ptr.hpp>
  37 #include <boost/tr1/array.hpp>
  38 
  39 #include <algorithm>
  40 #include <complex>
  41 #include <iostream>
  42 #include <iterator>
  43 #include <limits>
  44 #include <queue>
  45 #include <vector>
  46 
  47 struct point_type {
  48         typedef size_t element_type;
  49 
  50         element_type x, y;
  51 
  52         constexpr point_type() noexcept(true)
  53         : x(), y() {
  54         }
  55         constexpr point_type(element_type xi, element_type yi) noexcept(true)
  56         : x(xi), y(yi) {
  57         }
  58         bool operator==(point_type const &pt) const noexcept(true) {
  59                 return x==pt.x && y==pt.y;
  60         }
  61 };
  62 
  63 template<class V, unsigned int W, unsigned int H, V Init>
  64 struct bitmap {
  65         enum {
  66                 width=W,
  67                 height=H,
  68                 initial_value=Init
  69         };
  70         struct element_type {
  71                 V value;
  72                 bool computed;
  73 
  74                 constexpr element_type() noexcept(true)
  75                 : value(Init), computed(false) {
  76                 }
  77 
  78                 friend std::ostream &__fastcall operator<<(std::ostream &os, element_type const &e) noexcept(false) {
  79                         os<<std::hex<<e.value<<std::dec;
  80                         return os;
  81                 }
  82         };
  83         typedef typename std::tr1::array<element_type, height> column_type;
  84         /// The bitmaps is in the usual (x,y) cartesian co-ordinate format.
  85         typedef typename std::tr1::array<column_type, width> container_type;
  86         typedef typename column_type::const_iterator const_iterator;
  87         typedef typename column_type::iterator iterator;
  88         struct neighbourhood_item {
  89                 point_type pt;
  90                 iterator element;
  91 
  92                 constexpr neighbourhood_item() noexcept(true)
  93                 : pt(), element() {
  94                 }
  95                 explicit constexpr neighbourhood_item(iterator const &i) noexcept(true)
  96                 : pt(), element(i) {
  97                 }
  98                 constexpr neighbourhood_item(point_type const &p, iterator const &i) noexcept(true)
  99                 : pt(p), element(i) {
 100                 }
 101         };
 102         typedef typename std::tr1::array<const neighbourhood_item, 8> const_neighbourhood_type;
 103         typedef typename std::tr1::array<neighbourhood_item, 8> neighbourhood_type;
 104 
 105         struct lower_iterations : std::binary_function<neighbourhood_item, neighbourhood_item, bool> {
 106                 typedef std::binary_function<neighbourhood_item, neighbourhood_item, bool> operation_type;
 107                 typedef typename operation_type::first_argument_type first_argument_type;
 108                 typedef typename operation_type::second_argument_type second_argument_type;
 109                 typedef typename operation_type::result_type result_type;
 110 
 111                 result_type __fastcall operator()(first_argument_type const &f, second_argument_type const &s) const noexcept(true) {
 112                         return f.element->value>s.element->value;
 113                 }
 114         };
 115 
 116         container_type screen;
 117 
 118         const_iterator __fastcall find(point_type const &pt) const noexcept(true) {
 119                 return ((screen.begin()+pt.x)->begin()+pt.y);
 120         }
 121         iterator __fastcall find(point_type const &pt) noexcept(true) {
 122                 return ((screen.begin()+pt.x)->begin()+pt.y);
 123         }
 124 
 125         const_iterator __fastcall bottom_left() const noexcept(true) {
 126                 return screen.begin()->begin();
 127         }
 128 
 129         const_iterator __fastcall top_right() const noexcept(true) {
 130                 return ((screen.end()-1)->end()-1);
 131         }
 132 
 133         bool __fastcall inside(point_type const &pt) const noexcept(true) {
 134                 return pt.x<width && pt.y<height;
 135         }
 136 
 137         const_neighbourhood_type neighbours(point_type const &p) const noexcept(true) {
 138                 const const_neighbourhood_type ret={
 139                         (inside(point_type(p.x, p.y+1)) ? typename const_neighbourhood_type::value_type(point_type(p.x, p.y+1), find(point_type(p.x, p.y+1))) : typename const_neighbourhood_type::value_type(screen.begin()->end())),  // North.
 140                         (inside(point_type(p.x+1, p.y+1)) ? typename const_neighbourhood_type::value_type(point_type(p.x+1, p.y+1), find(point_type(p.x+1, p.y+1))) : typename const_neighbourhood_type::value_type(screen.begin()->end())),    // NE.
 141                         (inside(point_type(p.x+1, p.y)) ? typename const_neighbourhood_type::value_type(point_type(p.x+1, p.y), find(point_type(p.x+1, p.y))) : typename const_neighbourhood_type::value_type(screen.begin()->end())),  // East.
 142                         (inside(point_type(p.x+1, p.y-1)) ? typename const_neighbourhood_type::value_type(point_type(p.x+1, p.y-1), find(point_type(p.x+1, p.y-1))) : typename const_neighbourhood_type::value_type(screen.begin()->end())),    // SE.
 143                         (inside(point_type(p.x, p.y-1)) ? typename const_neighbourhood_type::value_type(point_type(p.x, p.y-1), find(point_type(p.x, p.y-1))) : typename const_neighbourhood_type::value_type(screen.begin()->end())),  // South.
 144                         (inside(point_type(p.x-1, p.y-1)) ? typename const_neighbourhood_type::value_type(point_type(p.x-1, p.y-1), find(point_type(p.x-1, p.y-1))) : typename const_neighbourhood_type::value_type(screen.begin()->end())),    // SW.
 145                         (inside(point_type(p.x-1, p.y)) ? typename const_neighbourhood_type::value_type(point_type(p.x-1, p.y), find(point_type(p.x-1, p.y))) : typename const_neighbourhood_type::value_type(screen.begin()->end())),  // West.
 146                         (inside(point_type(p.x-1, p.y+1)) ? typename const_neighbourhood_type::value_type(point_type(p.x-1, p.y+1), find(point_type(p.x-1, p.y+1))) : typename const_neighbourhood_type::value_type(screen.begin()->end()))     // NW.
 147                 };
 148                 return ret;
 149         }
 150         neighbourhood_type neighbours(point_type const &p) noexcept(true) {
 151                 const neighbourhood_type ret={{
 152                         (inside(point_type(p.x, p.y+1)) ? typename neighbourhood_type::value_type(point_type(p.x, p.y+1), find(point_type(p.x, p.y+1))) : typename neighbourhood_type::value_type(screen.begin()->end())),      // North.
 153                         (inside(point_type(p.x+1, p.y+1)) ? typename neighbourhood_type::value_type(point_type(p.x+1, p.y+1), find(point_type(p.x+1, p.y+1))) : typename neighbourhood_type::value_type(screen.begin()->end())),        // NE.
 154                         (inside(point_type(p.x+1, p.y)) ? typename neighbourhood_type::value_type(point_type(p.x+1, p.y), find(point_type(p.x+1, p.y))) : typename neighbourhood_type::value_type(screen.begin()->end())),      // East.
 155                         (inside(point_type(p.x+1, p.y-1)) ? typename neighbourhood_type::value_type(point_type(p.x+1, p.y-1), find(point_type(p.x+1, p.y-1))) : typename neighbourhood_type::value_type(screen.begin()->end())),        // SE.
 156                         (inside(point_type(p.x, p.y-1)) ? typename neighbourhood_type::value_type(point_type(p.x, p.y-1), find(point_type(p.x, p.y-1))) : typename neighbourhood_type::value_type(screen.begin()->end())),      // South.
 157                         (inside(point_type(p.x-1, p.y-1)) ? typename neighbourhood_type::value_type(point_type(p.x-1, p.y-1), find(point_type(p.x-1, p.y-1))) : typename neighbourhood_type::value_type(screen.begin()->end())),        // SW.
 158                         (inside(point_type(p.x-1, p.y)) ? typename neighbourhood_type::value_type(point_type(p.x-1, p.y), find(point_type(p.x-1, p.y))) : typename neighbourhood_type::value_type(screen.begin()->end())),      // West.
 159                         (inside(point_type(p.x-1, p.y+1)) ? typename neighbourhood_type::value_type(point_type(p.x-1, p.y+1), find(point_type(p.x-1, p.y+1))) : typename neighbourhood_type::value_type(screen.begin()->end())) // NW.
 160                 }};
 161                 return ret;
 162         }
 163 
 164         friend std::ostream &__fastcall operator<<(std::ostream &os, bitmap const &b) noexcept(false) {
 165                 os<<"Width="<<b.width
 166                         <<", height="<<b.height
 167                         <<", initial value="<<static_cast<V>(b.initial_value)
 168                         <<std::endl;
 169                 for (typename container_type::size_type y=0; y<b.height; ++y) {
 170                         for (typename container_type::size_type x=0; x<b.width; ++x) {
 171                                 os<<b.screen[x][y];
 172                         }
 173                         os<<"\n";
 174                 }
 175                 return os;
 176         }
 177 };
 178 
 179 template<class T>
 180 struct complex_plane_t {
 181         typedef std::complex<T> element_type;
 182 
 183         const element_type bottom_left, top_right;
 184 
 185         __stdcall complex_plane_t(element_type const &bl, element_type const &tr) noexcept(true)
 186         : bottom_left(bl), top_right(tr) {
 187                 assert(width()>0);
 188                 assert(height()>0);
 189         }
 190 
 191         typename element_type::value_type __fastcall width() const noexcept(true) {
 192                 return top_right.real()-bottom_left.real();
 193         }
 194 
 195         typename element_type::value_type __fastcall height() const noexcept(true) {
 196                 return top_right.imag()-bottom_left.imag();
 197         }
 198 
 199         friend std::ostream &__fastcall operator<<(std::ostream &os, complex_plane_t const &cp) noexcept(false) {
 200                 os<<cp.bottom_left<<cp.top_right;
 201                 return os;
 202         }
 203 };
 204 
 205 template<class Op, class WkQ>
 206 class greedy_render_t {
 207 public:
 208         typedef Op operation_type;
 209         typedef typename operation_type::screen_type screen_type;
 210         typedef typename screen_type::lower_iterations::first_argument_type result_type;
 211         typedef WkQ work_queue_type;
 212 
 213 private:
 214         screen_type &screen;
 215         work_queue_type &work_queue;
 216         operation_type const fn;
 217         result_type pt;
 218 
 219 public:
 220         explicit __stdcall greedy_render_t(screen_type &scr, work_queue_type &wk, Op const &op, result_type const &p) noexcept(true)
 221         : screen(scr), work_queue(wk), fn(op), pt(p) {
 222                 assert(pt.element);
 223         }
 224 
 225         void __fastcall start() noexcept(true) {
 226                 process(pt);
 227         }
 228 
 229         void __fastcall process(result_type &pt_to_compute) noexcept(noexcept(fn.operator()(pt_to_compute.pt))) {
 230                 pt_to_compute=pt;
 231                 assert(pt_to_compute.element);
 232                 if (!pt_to_compute.element->computed) {
 233                         pt_to_compute.element->value=fn.operator()(pt_to_compute.pt);
 234                         pt_to_compute.element->computed=true;
 235                         add_neighbours(pt_to_compute);
 236                 }
 237         }
 238 
 239         bool __fastcall operator<(greedy_render_t const &gr) const noexcept(true) {
 240                 return typename screen_type::lower_iterations().operator()(pt, gr.pt);
 241         }
 242 
 243         friend std::ostream &__fastcall operator<<(std::ostream &os, greedy_render_t const &cp) noexcept(false) {
 244                 os<<cp.fn;
 245                 return os;
 246         }
 247 
 248 private:
 249         /**
 250                 TODO JMG: this doesn't work as it recurses, adding new items to the queue, yet waiting for each item to be computed. The big clue this is balmy is that it's using a for_each (which has side-effects), with a wait in the functor. The side effect is the setting of all iteration values into the bitmap, when finished, the render is completed. This is not reflected in this function.
 251         */
 252         void __fastcall add_neighbours(result_type const &cur_pt) noexcept(true) {
 253                 const typename screen_type::neighbourhood_type neighbours(screen.neighbours(cur_pt.pt));
 254                 std::for_each(
 255                         neighbours.begin(),
 256                         neighbours.end(),
 257                         [&screen, &work_queue, &fn, &cur_pt](typename screen_type::neighbourhood_type::const_iterator i) {
 258                                 assert(i->element);
 259                                 if (i->element!=screen.screen.begin()->end() && !i->element->computed && i->element->value==screen.initial_value) {
 260                                         auto const &computed_pt=work_queue<<typename work_queue_type::joinable()<<greedy_render_t<Op, WkQ>(screen, work_queue, fn, *i);
 261                                         *computed_pt;
 262                                 }
 263                         }
 264                 }
 265         }
 266 };
 267 
 268 template<class CP, class Scr>
 269 struct Mandelbrot : public std::unary_function<point_type, unsigned long> {
 270         typedef CP complex_plane_type;
 271         typedef Scr screen_type;
 272         typedef std::unary_function<point_type, unsigned long>::argument_type argument_type;
 273         typedef std::unary_function<point_type, unsigned long>::result_type result_type;
 274         typedef result_type iterations_type;
 275         typedef typename complex_plane_type::element_type::value_type bailout_type;
 276 
 277         const iterations_type max_iters;
 278         const bailout_type bailout_sqrd;
 279 
 280         __stdcall Mandelbrot(complex_plane_type &cp, screen_type &sc, iterations_type mi, bailout_type b=2) noexcept(true)
 281         : max_iters(mi), bailout_sqrd(b*b), plane(cp), d_x(plane.width()/sc.width), d_y(plane.height()/sc.height) {
 282         }
 283 
 284         iterations_type __fastcall operator()(argument_type const &p) const noexcept(true) {
 285                 iterations_type iter=0;
 286                 const typename complex_plane_type::element_type c=pt_to_plane(p);
 287                 typename complex_plane_type::element_type z(0);
 288                 while (++iter<max_iters && std::norm(z)<bailout_sqrd) {
 289                         z=sqr(z)+c;
 290                 };
 291                 return iter;
 292         }
 293 
 294         friend std::ostream &__fastcall operator<<(std::ostream &os, Mandelbrot const &m) noexcept(false) {
 295                 os<<"Max. iterations="<<m.max_iters
 296                         <<", bailout="<<std::sqrt(m.bailout_sqrd)
 297                         <<", step=("<<m.d_x
 298                         <<", "<<m.d_y<<")";
 299                 return os;
 300         }
 301 
 302 private:
 303         complex_plane_type &plane;
 304         const typename complex_plane_type::element_type::value_type d_x;
 305         const typename complex_plane_type::element_type::value_type d_y;
 306 
 307         typename complex_plane_type::element_type __fastcall pt_to_plane(argument_type const &p) const noexcept(true) {
 308                 return typename complex_plane_type::element_type(
 309                         plane.bottom_left.real()+d_x*p.x,
 310                         plane.bottom_left.imag()+d_y*p.y
 311                 );
 312         };
 313 
 314         static typename complex_plane_type::element_type __fastcall sqr(typename complex_plane_type::element_type const &c) noexcept(true) {
 315                 return typename complex_plane_type::element_type(c.real()*c.real()-c.imag()*c.imag(), 2*c.real()*c.imag());
 316         }
 317 };
 318 
 319 BOOST_AUTO_TEST_SUITE(fractals)
 320 
 321 BOOST_AUTO_TEST_CASE(mandelbrot)
 322 {
 323         typedef bitmap<unsigned long, 170, 170, -1> screen_t;
 324         typedef complex_plane_t<long double> plane_t;
 325         typedef Mandelbrot<plane_t, screen_t> Mandelbrot_t;
 326         typedef std::priority_queue<screen_t::lower_iterations::first_argument_type, std::vector<screen_t::lower_iterations::first_argument_type>, screen_t::lower_iterations> work_queue_type;
 327 
 328         typedef jmmcg::ppd::thread_pool<
 329                 jmmcg::ppd::pool_traits::work_distribution_mode_t::worker_threads_get_work<pool_traits::work_distribution_mode_t::queue_model_t::pool_owns_queue>,
 330                 jmmcg::ppd::pool_traits::size_mode_t::fixed_size,
 331                 jmmcg::ppd::pool_aspects<
 332                         jmmcg::ppd::generic_traits::return_data::joinable,
 333                         jmmcg::ppd::platform_api,
 334                         jmmcg::ppd::heavyweight_threading,
 335                         jmmcg::ppd::pool_traits::prioritised_queue,
 336                         std::less,
 337                         1,
 338                         jmmcg::ppd::basic_statistics/*,
 339                         jmmcg::ppd::control_flow_graph*/
 340                 >
 341         > pool_type;
 342         typedef greedy_render_t<Mandelbrot_t, pool_type> renderer_t;
 343 
 344         pool_type pool(2);
 345         plane_t complex_plane(plane_t::element_type(-2.0, -1.5), plane_t::element_type(0.9, 1.5));
 346         std::cout<<complex_plane<<std::endl;
 347         screen_t screen;
 348         const point_type start(0, 0);
 349         const screen_t::lower_iterations::first_argument_type pt(start, screen.find(start));
 350         renderer_t renderer(screen, pool, Mandelbrot_t(complex_plane, screen, 15), pt);
 351         std::cout<<renderer<<std::endl;
 352         renderer.start();
 353         std::cout<<screen<<std::endl;
 354 }
 355 
 356 BOOST_AUTO_TEST_SUITE_END()