newton.hpp

Go to the documentation of this file.

#pragma once
#ifndef OGLPLUS_IMAGES_NEWTON_1107121519_HPP
#define OGLPLUS_IMAGES_NEWTON_1107121519_HPP

#include <oglplus/images/image.hpp>
#include <oglplus/math/vector.hpp>

#include <cassert>

namespace oglplus {
namespace images {


class NewtonFractal
 : public Image
{
private:
    // complex number division
    static Vec2f _cdiv(Vec2f a, Vec2f b)
    {
        float d = Dot(b, b);
        if(d == 0.0f) return a;
        else return Vec2f(
            (a.x()*b.x() + a.y()*b.y()) / d,
            (a.y()*b.x() - a.x()*b.y()) / d
        );
    }

    template <typename T>
    static T _mix(T a, T b, float coef)
    {
        return a*(1.0f - coef) + b*coef;
    }

    template <typename Function, typename Mixer, std::size_t N>
    void _make(
        GLsizei width,
        GLsizei height,
        Function,
        Mixer mixer,
        Vec2f lb,
        Vec2f rt,
        Vector<float, N> c1,
        Vector<float, N> c2
    )
    {
        auto p = this->_begin<GLfloat>();

        for(GLsizei i=0; i!=width; ++i)
        for(GLsizei j=0; j!=height; ++j)
        {
            Vec2f z(
                _mix(lb.x(), rt.x(), float(i)/float(width-1)),
                _mix(lb.y(), rt.y(), float(j)/float(height-1))
            );
            std::size_t n, max = 256;
            for(n = 0; n != max; ++n)
            {
                Vec2f zn = z - _cdiv(
                    Function::f(z),
                    Function::df(z)
                );
                if(Distance(zn, z) < 0.00001f) break;
                z = zn;
            }
            Vector<float, N> c = _mix(
                c1,
                c2,
                mixer(float(n) / float(max-1))
            );
            for(n=0; n!=N; ++n)
            {
                assert(p != this->_end<GLfloat>());
                *p = c.At(n);
                ++p;
            }
        }
        assert(p == this->_end<GLfloat>());
    }
public:
    struct X3Minus1
    {
        static Vec2f f(Vec2f n)
        {
            return Vec2f(
                n.x()*n.x()*n.x() - 3.f*n.x()*n.y()*n.y() - 1.f,
                -n.y()*n.y()*n.y() + 3.f*n.x()*n.x()*n.y()
            );
        }

        static Vec2f df(Vec2f n)
        {
            return 3.0f * Vec2f(
                n.x()*n.x() - n.y()*n.y(),
                2.0 * n.x() * n.y()
            );
        }
    };

    struct X4Minus1
    {
        static Vec2f f(Vec2f n)
        {
            return Vec2f(
                n.x()*n.x()*n.x()*n.x() +
                n.y()*n.y()*n.y()*n.y() -
                6.f*n.x()*n.x()*n.y()*n.y() - 1.f,
                4.f*n.x()*n.x()*n.x()*n.y() -
                4.f*n.x()*n.y()*n.y()*n.y()
            );
        }

        static Vec2f df(Vec2f n)
        {
            return 4.0f * Vec2f(
                n.x()*n.x()*n.x() - 3.f*n.x()*n.y()*n.y(),
                -n.y()*n.y()*n.y() + 3.f*n.x()*n.x()*n.y()
            );
        }
    };

    typedef X3Minus1 DefaultFunction;

    struct NoopMixer
    {
        template <typename T>
        T operator()(T value) const
        {
            return value;
        }
    };

    struct PowMixer
    {
        float _exponent;

        PowMixer(float exponent)
         : _exponent(exponent)
        { }

        template <typename T>
        T operator()(T value) const
        {
            return std::pow(value, T(_exponent));
        }
    };

    typedef NoopMixer DefaultMixer;

#if OGLPLUS_DOCUMENTATION_ONLY || !OGLPLUS_NO_FUNCTION_TEMPLATE_DEFAULT_ARGS

    template <
        typename Function = DefaultFunction,
        typename Mixer = DefaultMixer
    >
    NewtonFractal(
        GLsizei width,
        GLsizei height,
        Vec3f c1,
        Vec3f c2,
        Vec2f lb = Vec2f(-1.0f, -1.0f),
        Vec2f rt = Vec2f( 1.0f,  1.0f),
        Function func = Function(),
        Mixer mixer = Mixer()
    ): Image(width, height, 1, 3, (GLfloat*)0)
    {
        _make(width, height, func, mixer, lb, rt, c1, c2);
    }


    template <
        typename Function = DefaultFunction,
        typename Mixer = DefaultMixer
    >
    NewtonFractal(
        GLsizei width,
        GLsizei height,
        Function func = Function(),
        Mixer mixer = Mixer()
    ): Image(width, height, 1, 1, (GLfloat*)0)
    {
        _make(
            width, height,
            func,
            mixer,
            Vec2f(-1.0f, -1.0f), Vec2f(1.0f, 1.0f),
            Vec1f(0.0f), Vec1f(1.0f)
        );
    }
#else
    template <typename Function, typename Mixer>
    NewtonFractal(
        GLsizei width,
        GLsizei height,
        Vec3f c1,
        Vec3f c2,
        Vec2f lb,
        Vec2f rt,
        Function func,
        Mixer mixer
    ): Image(width, height, 1, 3, (GLfloat*)0)
    {
        _make(width, height, func, mixer, lb, rt, c1, c2);
    }

    template <typename Function, typename Mixer>
    NewtonFractal(
        GLsizei width,
        GLsizei height,
        Function func,
        Mixer mixer
    ): Image(width, height, 1, 1, (GLfloat*)0)
    {
        _make(
            width, height,
            func,
            mixer,
            Vec2f(-1.0f, -1.0f), Vec2f(1.0f, 1.0f),
            Vec1f(0.0f), Vec1f(1.0f)
        );
    }
#endif
};

} // images
} // oglplus

#endif // include guard