mandelbrot.cpp

Go to the documentation of this file.

/******************************************************************************
** Copyright © 2002 by J.M.McGuiness, coder@hussar.me.uk
**
** This library is free software; you can redistribute it and/or
** modify it under the terms of the GNU Lesser General Public
** License as published by the Free Software Foundation; either
** version 2.1 of the License, or (at your option) any later version.
**
** This library is distributed in the hope that it will be useful,
** but WITHOUT ANY WARRANTY; without even the implied warranty of
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
** Lesser General Public License for more details.
**
** You should have received a copy of the GNU Lesser General Public
** License along with this library; if not, write to the Free Software
** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
 
#include "stdafx.h"
 
#define BOOST_TEST_MODULE libjmmcg_tests
#include <boost/test/included/unit_test.hpp>
 
#include "core/thread_pool_workers.hpp"
 
#include <cassert>
#include <stdexcept>
 
#include <boost/mpl/assert.hpp>
#include <boost/static_assert.hpp>
#include <boost/type_traits.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/tr1/array.hpp>
 
#include <algorithm>
#include <complex>
#include <iostream>
#include <iterator>
#include <limits>
#include <queue>
#include <vector>
 
using namespace libjmmcg;
 
struct point_type {
   typedef size_t element_type;
 
   element_type x, y;
 
   constexpr point_type() noexcept(true)
   : x(), y() {
   }
   constexpr point_type(element_type xi, element_type yi) noexcept(true)
   : x(xi), y(yi) {
   }
   bool operator==(point_type const &pt) const noexcept(true) {
      return x==pt.x && y==pt.y;
   }
};
 
template<class V, unsigned int W, unsigned int H, V Init>
struct bitmap {
   enum {
      width=W,
      height=H,
      initial_value=Init
   };
   struct element_type {
      V value;
      bool computed;
 
      constexpr element_type() noexcept(true)
      : value(Init), computed(false) {
      }
 
      friend std::ostream &__fastcall operator<<(std::ostream &os, element_type const &e) noexcept(false) {
         os<<std::hex<<e.value<<std::dec;
         return os;
      }
   };
   typedef typename std::tr1::array<element_type, height> column_type;
   /// The bitmaps is in the usual (x,y) cartesian co-ordinate format.
   typedef typename std::tr1::array<column_type, width> container_type;
   typedef typename column_type::const_iterator const_iterator;
   typedef typename column_type::iterator iterator;
   struct neighbourhood_item {
      point_type pt;
      iterator element;
 
      constexpr neighbourhood_item() noexcept(true)
      : pt(), element() {
      }
      explicit constexpr neighbourhood_item(iterator const &i) noexcept(true)
      : pt(), element(i) {
      }
      constexpr neighbourhood_item(point_type const &p, iterator const &i) noexcept(true)
      : pt(p), element(i) {
      }
   };
   typedef typename std::tr1::array<const neighbourhood_item, 8> const_neighbourhood_type;
   typedef typename std::tr1::array<neighbourhood_item, 8> neighbourhood_type;
 
   struct lower_iterations : std::binary_function<neighbourhood_item, neighbourhood_item, bool> {
      typedef std::binary_function<neighbourhood_item, neighbourhood_item, bool> operation_type;
      typedef typename operation_type::first_argument_type first_argument_type;
      typedef typename operation_type::second_argument_type second_argument_type;
      typedef typename operation_type::result_type result_type;
 
      result_type __fastcall operator()(first_argument_type const &f, second_argument_type const &s) const noexcept(true) {
         return f.element->value>s.element->value;
      }
   };
 
   container_type screen;
 
   const_iterator __fastcall find(point_type const &pt) const noexcept(true) {
      return ((screen.begin()+pt.x)->begin()+pt.y);
   }
   iterator __fastcall find(point_type const &pt) noexcept(true) {
      return ((screen.begin()+pt.x)->begin()+pt.y);
   }
 
   const_iterator __fastcall bottom_left() const noexcept(true) {
      return screen.begin()->begin();
   }
 
   const_iterator __fastcall top_right() const noexcept(true) {
      return ((screen.end()-1)->end()-1);
   }
 
   bool __fastcall inside(point_type const &pt) const noexcept(true) {
      return pt.x<width && pt.y<height;
   }
 
   const_neighbourhood_type neighbours(point_type const &p) const noexcept(true) {
      const const_neighbourhood_type ret={
         (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.
         (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.
         (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.
         (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.
         (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.
         (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.
         (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.
         (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.
      };
      return ret;
   }
   neighbourhood_type neighbours(point_type const &p) noexcept(true) {
      const neighbourhood_type ret={{
         (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.
         (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.
         (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.
         (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.
         (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.
         (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.
         (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.
         (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.
      }};
      return ret;
   }
 
   friend std::ostream &__fastcall operator<<(std::ostream &os, bitmap const &b) noexcept(false) {
      os<<"Width="<<b.width
         <<", height="<<b.height
         <<", initial value="<<static_cast<V>(b.initial_value)
         <<std::endl;
      for (typename container_type::size_type y=0; y<b.height; ++y) {
         for (typename container_type::size_type x=0; x<b.width; ++x) {
            os<<b.screen[x][y];
         }
         os<<"\n";
      }
      return os;
   }
};
 
template<class T>
struct complex_plane_t {
   typedef std::complex<T> element_type;
 
   const element_type bottom_left, top_right;
 
   __stdcall complex_plane_t(element_type const &bl, element_type const &tr) noexcept(true)
   : bottom_left(bl), top_right(tr) {
      assert(width()>0);
      assert(height()>0);
   }
 
   typename element_type::value_type __fastcall width() const noexcept(true) {
      return top_right.real()-bottom_left.real();
   }
 
   typename element_type::value_type __fastcall height() const noexcept(true) {
      return top_right.imag()-bottom_left.imag();
   }
 
   friend std::ostream &__fastcall operator<<(std::ostream &os, complex_plane_t const &cp) noexcept(false) {
      os<<cp.bottom_left<<cp.top_right;
      return os;
   }
};
 
template<class Op, class WkQ>
class greedy_render_t {
public:
   typedef Op operation_type;
   typedef typename operation_type::screen_type screen_type;
   typedef typename screen_type::lower_iterations::first_argument_type result_type;
   typedef WkQ work_queue_type;
 
private:
   screen_type &screen;
   work_queue_type &work_queue;
   operation_type const fn;
   result_type pt;
 
public:
   explicit __stdcall greedy_render_t(screen_type &scr, work_queue_type &wk, Op const &op, result_type const &p) noexcept(true)
   : screen(scr), work_queue(wk), fn(op), pt(p) {
      assert(pt.element);
   }
 
   void __fastcall start() noexcept(true) {
      process(pt);
   }
 
   void __fastcall process(result_type &pt_to_compute) noexcept(noexcept(fn.operator()(pt_to_compute.pt))) {
      pt_to_compute=pt;
      assert(pt_to_compute.element);
      if (!pt_to_compute.element->computed) {
         pt_to_compute.element->value=fn.operator()(pt_to_compute.pt);
         pt_to_compute.element->computed=true;
         add_neighbours(pt_to_compute);
      }
   }
 
   bool __fastcall operator<(greedy_render_t const &gr) const noexcept(true) {
      return typename screen_type::lower_iterations().operator()(pt, gr.pt);
   }
 
   friend std::ostream &__fastcall operator<<(std::ostream &os, greedy_render_t const &cp) noexcept(false) {
      os<<cp.fn;
      return os;
   }
 
private:
   /**
      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.
   */
   void __fastcall add_neighbours(result_type const &cur_pt) noexcept(true) {
      const typename screen_type::neighbourhood_type neighbours(screen.neighbours(cur_pt.pt));
      std::for_each(
         neighbours.begin(),
         neighbours.end(),
         [&screen, &work_queue, &fn, &cur_pt](typename screen_type::neighbourhood_type::const_iterator i) {
            assert(i->element);
            if (i->element!=screen.screen.begin()->end() && !i->element->computed && i->element->value==screen.initial_value) {
               auto const &computed_pt=work_queue<<typename work_queue_type::joinable()<<greedy_render_t<Op, WkQ>(screen, work_queue, fn, *i);
               *computed_pt;
            }
         }
      }
   }
};
 
template<class CP, class Scr>
struct Mandelbrot : public std::unary_function<point_type, unsigned long> {
   typedef CP complex_plane_type;
   typedef Scr screen_type;
   typedef std::unary_function<point_type, unsigned long>::argument_type argument_type;
   typedef std::unary_function<point_type, unsigned long>::result_type result_type;
   typedef result_type iterations_type;
   typedef typename complex_plane_type::element_type::value_type bailout_type;
 
   const iterations_type max_iters;
   const bailout_type bailout_sqrd;
 
   __stdcall Mandelbrot(complex_plane_type &cp, screen_type &sc, iterations_type mi, bailout_type b=2) noexcept(true)
   : max_iters(mi), bailout_sqrd(b*b), plane(cp), d_x(plane.width()/sc.width), d_y(plane.height()/sc.height) {
   }
 
   iterations_type __fastcall operator()(argument_type const &p) const noexcept(true) {
      iterations_type iter=0;
      const typename complex_plane_type::element_type c=pt_to_plane(p);
      typename complex_plane_type::element_type z(0);
      while (++iter<max_iters && std::norm(z)<bailout_sqrd) {
         z=sqr(z)+c;
      };
      return iter;
   }
 
   friend std::ostream &__fastcall operator<<(std::ostream &os, Mandelbrot const &m) noexcept(false) {
      os<<"Max. iterations="<<m.max_iters
         <<", bailout="<<std::sqrt(m.bailout_sqrd)
         <<", step=("<<m.d_x
         <<", "<<m.d_y<<")";
      return os;
   }
 
private:
   complex_plane_type &plane;
   const typename complex_plane_type::element_type::value_type d_x;
   const typename complex_plane_type::element_type::value_type d_y;
 
   typename complex_plane_type::element_type __fastcall pt_to_plane(argument_type const &p) const noexcept(true) {
      return typename complex_plane_type::element_type(
         plane.bottom_left.real()+d_x*p.x,
         plane.bottom_left.imag()+d_y*p.y
      );
   };
 
   static typename complex_plane_type::element_type __fastcall sqr(typename complex_plane_type::element_type const &c) noexcept(true) {
      return typename complex_plane_type::element_type(c.real()*c.real()-c.imag()*c.imag(), 2*c.real()*c.imag());
   }
};
 
BOOST_AUTO_TEST_SUITE(fractals)
 
BOOST_AUTO_TEST_CASE(mandelbrot)
{
   typedef bitmap<unsigned long, 170, 170, -1> screen_t;
   typedef complex_plane_t<long double> plane_t;
   typedef Mandelbrot<plane_t, screen_t> Mandelbrot_t;
   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;
 
   typedef ppd::thread_pool<
      ppd::pool_traits::work_distribution_mode_t::worker_threads_get_work<pool_traits::work_distribution_mode_t::queue_model_t::pool_owns_queue>,
      ppd::pool_traits::size_mode_t::fixed_size,
      ppd::pool_aspects<
         ppd::generic_traits::return_data::joinable,
         ppd::platform_api,
         ppd::heavyweight_threading,
         ppd::pool_traits::prioritised_queue,
         std::less,
         1,
         ppd::basic_statistics/*,
         ppd::control_flow_graph*/
      >
   > pool_type;
   typedef greedy_render_t<Mandelbrot_t, pool_type> renderer_t;
 
   pool_type pool(2);
   plane_t complex_plane(plane_t::element_type(-2.0, -1.5), plane_t::element_type(0.9, 1.5));
   std::cout<<complex_plane<<std::endl;
   screen_t screen;
   const point_type start(0, 0);
   const screen_t::lower_iterations::first_argument_type pt(start, screen.find(start));
   renderer_t renderer(screen, pool, Mandelbrot_t(complex_plane, screen, 15), pt);
   std::cout<<renderer<<std::endl;
   renderer.start();
   std::cout<<screen<<std::endl;
}
 
BOOST_AUTO_TEST_SUITE_END()