[Boost-users] RFC: library for fractal algorithms

Subject: [Boost-users] RFC: library for fractal algorithms
From: Richard (legalize+jeeves_at_[hidden])
Date: 2014-12-16 22:21:17

[Please do not mail me a copy of your followup]

I'm interested in talking with other folks about the best way to
generalize various fractal algorithms in the style of the standard
library: iterators, containers and algorithms.  If there are any
libraries that already exist that apply to this sort of code, I'd love
to hear about them.

Most fractals consist of iterating a sequence for a point in the
complex plane.  The iterating sequence is parameterized for each point
in the complex plane and some characteristic of the iterating sequence
is used to color that point.  The result is an image.

For instance, the classic Mandelbrot fractal is:

    z <- z^2 + c

    z0 := 0

The complex number 'c' changes for each point in the complex plane.
This sequence is iterated until norm(z) > 4, or some maximum iteration
count is reached.  Pixels in the image correspond to each value of 'c'
and are typically colored according to the number of iterations needed
until norm(z) exceeded 4.

There are a couple things that make for interesting variation in the
algorithms used in rendering fractal images:

- choice of formula to iterate in the complex plane
- parameters for deciding the "bailout" criteria:
  - iteration counts
  - predicate function on z, i.e.
    [](complex<float> z){ return norm(z) > 4; }

Iterating until the bailout test passesis the simplest mechanism, but
there are gains to be made by applying more sophisticated mechanisms.
A significant optimization is to look for cycles in the sequences.
With very large iteration counts, detecting cycles in the iteration
can save significant amounts of time.  The cycle detection algorithm
is basically from Knuth and isn't specific to fractals but ends up
being mixed together with the Mandelbrot set computation in most
programs.  This seems to be another opportunity to use generic
composition in order to apply this orthogonal feature to any iterated
formula in the complex plane.

I did a little experiment casting the Mandelbrot formula as an
iterator and got the code at the end of this message.  Ideally I'd
like to explore approaches that result in code that is ignorant of the
underlying representation of a complex number in order to have
algorithms that can be applied to arbitrary precision numeric types.
My little experiment works, but even just looking at this little
experiment makes me think that the mandelbrot iterator is doing too
much: it is performing the z <- z^2 + c business, but it is also doing
the bailout test and these seem to be separate responsibilities.

#include <complex>
#include <iostream>
#include <iterator>
#include <vector>

class mandelbrot : public std::iterator<
    std::input_iterator_tag,
    std::complex<float>,
    unsigned>
{
public:
    mandelbrot()
        : c_{},
        iter_{},
        max_iter_{},
        bailout_{},
        z_{},
        bailed_out_{true}
    {
    }

    mandelbrot(std::complex<float> const& c,
            int max_iter = 256,
            float bailout = 2.0f)
        : c_{c},
        max_iter_{max_iter},
        bailout_{bailout*bailout},
        iter_{0},
        z_{0.f, 0.f},
        bailed_out_{false}
    {
    }

    mandelbrot& operator++()
    {
        z_ *= z_;
        z_ += c_;
        bailed_out_ = ++iter_ > max_iter_ ||
            std::norm(z_) >= bailout_;

        return *this;
    }

    bool operator==(mandelbrot const& rhs) const
    {
        if (bailed_out_ && rhs.bailed_out_)
        {
            return true;
        }

        if (bailed_out_ != rhs.bailed_out_)
        {
            return false;
        }

        return z_ == rhs.z_;
    }

    std::complex<float> const& operator*() const
    {
        return z_;
    }

private:
    std::complex<float> c_;
    int const max_iter_;
    float const bailout_;
    int iter_;
    std::complex<float> z_;
    bool bailed_out_;
};

bool operator!=(mandelbrot const &lhs, mandelbrot const &rhs)
{
    return !(lhs == rhs);
}

int main()
{
    float const r_min = -3.0f;
    float const r_max = 1.0f;
    float const i_min = -2.0f;
    float const i_max = 2.0f;

    std::vector<float> imag;
    imag.resize(24);
    for (int y = 0; y < 24; ++y)
    {
        imag[y] = i_min + (i_max - i_min)*(23 - y)/24.0f;
    }
    std::vector<float> real;
    real.resize(80);
    for (int x = 0; x < 80; ++x)
    {
        real[x] = r_min + (r_max - r_min)*x/80.0f;
    }

    std::string printables;
    for (int i = 32; i < 128; ++i)
    {
        printables.push_back(char(i));
    }
    for (int y = 0; y < 24; ++y)
    {
        float const im = imag[y];
        for (int x = 0; x < 80; ++x)
        {
            float const re = real[x];
            mandelbrot end;
            mandelbrot it{std::complex<float>{re, im}};
            unsigned count = std::distance(it, end);
            std::cout << printables[count % printables.size()];
        }
        std::cout << '\n';
    }
}
-- 
"The Direct3D Graphics Pipeline" free book <http://tinyurl.com/d3d-pipeline>
     The Computer Graphics Museum <http://computergraphicsmuseum.org>
         The Terminals Wiki <http://terminals.classiccmp.org>
  Legalize Adulthood! (my blog) <http://legalizeadulthood.wordpress.com>