jaja

libmandel/src/opencl/float.cl

@@ -42,20 +42,32 @@ __kernel void iterate_vec4(__global float* A, const int width, float xl, float y
    int4 count = (int4)(0);
 
    int n = 0;
-   while (n < max) {
-       float4 ab = a * b;
-       float4 cmpVal = fma(a, a, b * b);
-       int4 cmp = isless(cmpVal, (float4)(16.0f));
-       if (!any(cmp)) break;
-       a = fma(a, a, -fma(b, b, -ca));
-       b = fma(2, ab, cb);
-       if (smooth) {
+   if (smooth) {
+       while (n < max) {
+           float4 ab = a * b;
+           float4 cmpVal = fma(a, a, b * b);
+           int4 cmp = isless(cmpVal, (float4)(16.0f));
+           if (!any(cmp)) break;
+           a = fma(a, a, -fma(b, b, -ca));
+           b = fma(2, ab, cb);
            resa = as_float4(as_int4(a) & cmp | (as_int4(resa) & ~cmp));
            resb = as_float4(as_int4(b) & cmp | (as_int4(resb) & ~cmp));
+           count += cmp & (int4)(1);
+           n++;
        }
-       count += cmp & (int4)(1);
-       n++;
    }
+    else {
+       while (n < max) {
+           float4 ab = a * b;
+           float4 cmpVal = fma(a, a, b * b);
+           int4 cmp = isless(cmpVal, (float4)(16.0f));
+           if (!any(cmp)) break;
+           a = fma(a, a, -fma(b, b, -ca));
+           b = fma(2, ab, cb);
+           count += cmp & (int4)(1);
+           n++;
+       }
+    }
    for (int i = 0; i < 4 && i + x < width; i++) {
        if (smooth != 0)
            A[index + i] = ((float) count[i]) + 1 - log(log(fma(resa[i], resa[i], resb[i] * resb[i])) / 2) / log(2.0f);

mandelvid/CMakeLists.txt

@@ -0,0 +1,80 @@
+cmake_minimum_required(VERSION 3.9)
+
+project(mandelvid VERSION 1.0.0 DESCRIPTION "mandel video generator")
+
+set(ARCH "X86_64" CACHE STRING "Target Architecture")
+
+
+find_package(PkgConfig)
+pkg_check_modules(LIBAVCODEC REQUIRED libavcodec)
+pkg_check_modules(LIBAVUTIL REQUIRED libavutil)
+pkg_check_modules(LIBAVFORMAT REQUIRED libavformat)
+pkg_check_modules(LIBSWSCALE REQUIRED libswscale)
+
+
+add_subdirectory(../libmandel ./libmandel)
+
+find_package(OpenCL)
+find_package(OpenMP)
+#set(Boost_DEBUG 1)
+set(Boost_USE_STATIC_LIBS ON)
+find_package(Boost 1.53)
+
+
+set(CMAKE_CXX_STANDARD 17)
+
+
+FILE(GLOB mvidsources src/*.cpp)
+FILE(GLOB mvidheaders include/*.h)
+
+add_executable(mvg ${mvidsources})
+target_include_directories(mvg PUBLIC "include")
+target_include_directories(mvg PUBLIC ${LIBAVCODEC_INCLUDE_DIRS})
+target_include_directories(mvg PUBLIC ${LIBAVUTIL_INCLUDE_DIRS})
+target_include_directories(mvg PUBLIC ${LIBAVFORMAT_INCLUDE_DIRS})
+target_include_directories(mvg PUBLIC ${LIBSWSCALE_INCLUDE_DIRS})
+target_link_libraries(mvg PUBLIC mandel)
+target_link_libraries(mvg PUBLIC qd)
+target_link_libraries(mvg PUBLIC ${LIBAVCODEC_LIBRARIES})
+target_link_libraries(mvg PUBLIC ${LIBAVUTIL_LIBRARIES})
+target_link_libraries(mvg PUBLIC ${LIBAVFORMAT_LIBRARIES})
+target_link_libraries(mvg PUBLIC ${LIBSWSCALE_LIBRARIES})
+
+include(CheckIPOSupported)
+check_ipo_supported(RESULT LTO_SUPPORTED)
+
+if (LTO_SUPPORTED AND WITH_LTO)
+    message("Enabling link-time optimization.")
+    set_property(TARGET mvg PROPERTY INTERPROCEDURAL_OPTIMIZATION TRUE)
+endif()
+
+
+#if(OPENCL_FOUND)
+#    target_compile_definitions(mvg PUBLIC WITH_OPENCL)
+#    target_include_directories(mvg PUBLIC
+#        "include"
+#        ${OpenCL_INCLUDE_DIRS}
+#    )
+#    link_directories(${OpenCL_LIBRARY})
+#else(OPENCL_FOUND)
+#    include_directories("include")
+#endif(OPENCL_FOUND)
+#
+#if (APPLE AND OpenCL_FOUND)
+#    SET(CMAKE_CXX_FLAGS  "${CMAKE_CXX_FLAGS} -framework OpenCL")
+#endif()
+#
+#    target_compile_definitions(mandel PUBLIC WITH_BOOST)
+#if(Boost_FOUND)
+#    target_compile_definitions(mandel PUBLIC WITH_BOOST)
+#    target_include_directories(mandel PRIVATE ${Boost_INCLUDE_DIRS})
+#    target_link_libraries(mandel PRIVATE ${Boost_LIBRARIES})
+#endif(Boost_FOUND)
+#
+#if(OpenMP_CXX_FOUND)
+#    target_link_libraries(mandel PUBLIC OpenMP::OpenMP_CXX)
+#endif()
+#if(OpenCL_FOUND)
+#    target_link_libraries(mandel PUBLIC OpenCL::OpenCL)
+#endif()
+

mandelvid/include/Bitmap.h

@@ -0,0 +1,62 @@
+#pragma once
+#ifndef BITMAP_H_
+#define BITMAP_H_
+
+#include "Color.h"
+#include <memory>
+#include <string>
+#include <functional>
+
+template<typename Pixel>
+struct Bitmap
+{
+    long width, height;
+    std::unique_ptr<Pixel[]> pixels;
+public:
+    Bitmap(void) :
+        width{ 0 },
+        height{ 0 },
+        pixels{ 0 }
+    {
+    }
+
+
+    Bitmap(long width, long height) :
+        width{ width }, height{ height },
+        pixels{ std::make_unique<Pixel[]>(width * height) }
+    {
+    }
+
+    Bitmap(Bitmap&&) = default;
+    Bitmap& operator = (Bitmap&&) = default;
+
+    ~Bitmap() = default;
+
+    template<typename T>
+    Bitmap<T> map(std::function<T(Pixel)> f) const {
+        Bitmap<T> b{ width, height };
+        for (long i = 0; i < width * height; i++) {
+            b.pixels[i] = f(pixels[i]);
+        }
+        return b;
+    }
+
+    Pixel& get(long x, long y)
+    {
+        return pixels[x + y * width];
+    }
+
+    const Pixel& get(long x, long y) const
+    {
+        return pixels[x + y * width];
+    }
+
+    void print(void)
+    {
+        for (size_t i = 0; i < width * height; i++) {
+            printf("%03d ", int(pixels[i].r));
+        }
+    }
+};
+
+#endif // BITMAP_H_

mandelvid/include/Color.h

@@ -0,0 +1,48 @@
+#pragma once
+#ifndef COLOR_H_
+#define COLOR_H_
+
+#include <cinttypes>
+
+
+struct RGBColor;
+struct RGBColorf;
+
+/*!
+ * \brief Represents a color in the sRGB color space with 8-bit channels
+ */
+struct RGBColor
+{
+    uint8_t r, g, b;
+
+    inline RGBColor(void) :
+        r{ 0 }, g{ 0 }, b{ 0 }
+    {
+    }
+
+    inline RGBColor(uint8_t r, uint8_t g, uint8_t b) :
+        r{ r }, g{ g }, b{ b }
+    {
+    }
+
+    RGBColor(const RGBColorf& rgb);
+};
+
+
+/*!
+ * \brief Represents a color in a linear RGB color space with 32-bit floating
+ *        point numbers as channels.
+ */
+struct RGBColorf
+{
+    float r, g, b;
+
+    inline RGBColorf(float r, float g, float b) :
+        r{ r }, g{ g }, b{ b }
+    {
+    }
+
+    RGBColorf(const RGBColor& srgb);
+};
+
+#endif // COLOR_H_

mandelvid/include/CubicSpline.h

@@ -0,0 +1,19 @@
+#ifndef CUBICSPLINE_H
+#define CUBICSPLINE_H
+
+#include <vector>
+#include <utility>
+#include <tuple>
+
+class CubicSpline
+{
+    /// contains x, y and y' of each interpolation point
+    std::vector<std::tuple<float, float, float>> points;
+    bool useSlopes;
+public:
+    CubicSpline(const std::vector<std::pair<float, float>>& dataPoints, bool useSlopes);
+
+    float interpolateAt(float x);
+};
+
+#endif // CUBICSPLINE_H

mandelvid/include/Gradient.h

@@ -0,0 +1,38 @@
+#ifndef GRADIENT_H
+#define GRADIENT_H
+
+#include <vector>
+#include <string>
+#include "Color.h"
+#include <tuple>
+#include <cinttypes>
+
+
+class Gradient
+{
+    /// the colors of this gradient stored in linear RGB format
+    /// so they can be easily interpolated
+    std::vector<RGBColorf> colors;
+    float max;
+    bool repeat;
+public:
+    Gradient(void);
+    Gradient(std::vector<std::pair<RGBColor, float>> colors, bool repeat = false, int precalcSteps = -1);
+
+    static Gradient defaultGradient(void);
+
+    //static Gradient readXml(const std::string& xml);
+
+    /*!
+     * \brief get a color at a specific position in this gradient
+     * \param x the position
+     * \return the color in sRGB format
+     */
+    RGBColor get(float x) const;
+private:
+    static RGBColorf lerpColors(RGBColorf a, RGBColorf b, float val);
+    static RGBColor lerpColors(RGBColor a, RGBColor b, float val);
+    std::tuple<RGBColor, RGBColor, float> getNeighbors(float x) const;
+};
+
+#endif // GRADIENT_H

mandelvid/include/MandelVideoGenerator.h

@@ -0,0 +1,43 @@
+#ifndef MANDELVIDEOGENERATOR_H
+#define MANDELVIDEOGENERATOR_H
+
+#include "MandelUtil.h"
+#include "Gradient.h"
+#include "Bitmap.h"
+
+struct ExportVideoInfo
+{
+    mnd::MandelViewport start;
+    mnd::MandelViewport end;
+    Gradient gradient;
+
+    int width;
+    int height;
+    int maxIterations;
+    int fps;
+    double zoomSpeed;
+
+    std::string path;
+
+    /// bitrate in kbps
+    int bitrate;
+
+    std::string preset;
+};
+
+
+class MandelVideoGenerator
+{
+    const ExportVideoInfo evi;
+public:
+    MandelVideoGenerator(const ExportVideoInfo& evi);
+
+    void generate(void);
+
+private:
+    Bitmap<RGBColor> overlay(const Bitmap<RGBColor>& outer,
+                             const Bitmap<RGBColor>& inner,
+                             double scale);
+};
+
+#endif // MANDELVIDEOGENERATOR_H

mandelvid/include/VideoStream.h

@@ -0,0 +1,50 @@
+#pragma once
+#define FFMPEG_ENABLED
+#ifdef FFMPEG_ENABLED
+
+
+#ifndef VIDEO_STREAM_H_
+#define VIDEO_STREAM_H_
+
+#include <string>
+#include "Bitmap.h"
+
+extern "C" {
+#   include <libavformat/avformat.h>
+#   include <libavformat/avio.h>
+#   include <libavcodec/avcodec.h>
+#   include <libavformat/avformat.h>
+#   include <libavutil/imgutils.h>
+#   include <libavutil/opt.h>
+#   include <libswscale/swscale.h>
+}
+
+class VideoStream
+{
+    const AVCodec* codec;
+    AVCodecContext* codecContext;
+    AVFormatContext* formatContext;
+    AVCodecParameters* params;
+    //FILE* file;
+    AVFrame* picture;
+    AVPacket* pkt;
+    AVStream* stream;
+    SwsContext* swsContext;
+    static const uint8_t endcode[];
+
+    int width;
+    int height;
+
+    int64_t frameIndex = 0;
+public:
+    VideoStream(int width, int height, const std::string& filename, int bitrate, int fps, const char* preset);
+    ~VideoStream(void);
+
+    void encode(AVFrame* frame);
+
+    void addFrame(const Bitmap<RGBColor>& frame);
+};
+
+#endif // VIDEO_STREAM_H_
+
+#endif // FFMPEG_ENABLED

mandelvid/src/Bitmap.cpp

@@ -0,0 +1,2 @@
+#include "Bitmap.h"
+

mandelvid/src/Color.cpp

@@ -0,0 +1,22 @@
+#include "Color.h"
+#include <cmath>
+#include <algorithm>
+
+
+RGBColor::RGBColor(const RGBColorf& rgb)
+{
+    float cr = std::clamp(rgb.r, 0.0f, 1.0f);
+    float cg = std::clamp(rgb.g, 0.0f, 1.0f);
+    float cb = std::clamp(rgb.b, 0.0f, 1.0f);
+    r = uint8_t(cr * cr * 255.0f);
+    g = uint8_t(cg * cg * 255.0f);
+    b = uint8_t(cb * cb * 255.0f);
+}
+
+
+RGBColorf::RGBColorf(const RGBColor& srgb)
+{
+    r = ::sqrtf(srgb.r / 255.0f);
+    g = ::sqrtf(srgb.g / 255.0f);
+    b = ::sqrtf(srgb.b / 255.0f);
+}

mandelvid/src/CubicSpline.cpp

@@ -0,0 +1,66 @@
+#include "CubicSpline.h"
+
+CubicSpline::CubicSpline(const std::vector<std::pair<float, float> >& dataPoints, bool useSlopes) :
+    useSlopes{ useSlopes }
+{
+    if (dataPoints.size() < 2) {
+        return;
+    }
+
+    points.push_back({ dataPoints[0].first, dataPoints[0].second,
+                       (dataPoints[1].second - dataPoints[0].second) /
+                       (dataPoints[1].first - dataPoints[0].first) });
+    for (size_t i = 1; i < dataPoints.size() - 1; i++) {
+
+        auto& dp1 = dataPoints[i - 1];
+        auto& dp2 = dataPoints[i];
+        auto& dp3 = dataPoints[i + 1];
+
+        float w1 = dp2.first - dp1.first;
+        float w2 = dp3.first - dp2.first;
+        float h1 = dp2.second - dp1.second;
+        float h2 = dp3.second - dp2.second;
+
+        float s1 = h1 / w1;
+        float s2 = h2 / w2;
+
+        float avgSlope = (s1 + s2) / 2;
+        points.push_back({ dp2.first, dp2.second, avgSlope });
+    }
+    points.push_back({ dataPoints[dataPoints.size() - 1].first, dataPoints[dataPoints.size() - 1].second,
+                       (dataPoints[dataPoints.size() - 2].second - dataPoints[dataPoints.size() - 1].second) /
+                       (dataPoints[dataPoints.size() - 2].first - dataPoints[dataPoints.size() - 1].first) });
+}
+
+
+float CubicSpline::interpolateAt(float x)
+{
+    const static auto h00 = [] (float t) { return (1 + 2 * t) * (1 - t) * (1 - t); };
+    const static auto h01 = [] (float t) { return t * t * (3 - 2 * t); };
+    const static auto h10 = [] (float t) { return t * (1 - t) * (1 - t); };
+    const static auto h11 = [] (float t) { return t * t * (t - 1); };
+    for (auto it = points.begin(); it != points.end() && (it + 1) != points.end(); ++it) {
+        auto& left = *it;
+        auto& right = *(it + 1);
+        float xleft = std::get<0>(left);
+        float xright = std::get<0>(right);
+        if (xleft < x && xright >= x) {
+            float w = (xright - xleft);
+            float t = (x - xleft) / w;
+            float yleft = std::get<1>(left);
+            float yright = std::get<1>(right);
+            float sleft = std::get<2>(left);
+            float sright = std::get<2>(right);
+
+            float inter = h00(t) * yleft +
+                          h01(t) * yright;
+
+            if (useSlopes)
+                inter += h10(t) * w * sleft +
+                         h11(t) * w * sright;
+
+            return inter;
+        }
+    }
+    return std::get<1>(points[points.size() - 1]);
+}

mandelvid/src/Gradient.cpp

@@ -0,0 +1,152 @@
+#include "Gradient.h"
+
+#include "CubicSpline.h"
+
+#include <cmath>
+#include <algorithm>
+#include <functional>
+
+
+Gradient::Gradient(void) :
+    max{ 1.0 }
+{
+}
+
+
+Gradient::Gradient(std::vector<std::pair<RGBColor, float>> colors, bool repeat, int precalcSteps) :
+    repeat{ repeat }
+{
+    if(colors.empty() || colors.size() < 2)
+        return;
+    std::sort(colors.begin(), colors.end(),
+        [] (const auto& a, const auto& b) {
+            return a.second < b.second;
+        });
+
+    max = colors.at(colors.size() - 1).second;
+
+    std::vector<std::pair<RGBColorf, float>> linearColors;
+    std::transform(colors.begin(), colors.end(), std::back_inserter(linearColors),
+                   [] (auto c) { return c; });
+
+    std::vector<std::pair<float, float>> rs;
+    std::vector<std::pair<float, float>> gs;
+    std::vector<std::pair<float, float>> bs;
+
+    std::transform(linearColors.begin(), linearColors.end(), std::back_inserter(rs),
+                   [] (auto p) { return std::pair{ p.second, p.first.r }; });
+    std::transform(linearColors.begin(), linearColors.end(), std::back_inserter(gs),
+                   [] (auto p) { return std::pair{ p.second, p.first.g }; });
+    std::transform(linearColors.begin(), linearColors.end(), std::back_inserter(bs),
+                   [] (auto p) { return std::pair{ p.second, p.first.b }; });
+
+    CubicSpline rsp(rs, false);
+    CubicSpline gsp(gs, false);
+    CubicSpline bsp(bs, false);
+
+    if(precalcSteps <= 0) {
+        precalcSteps = int(max * 15) + 10;
+    }
+
+    for (int i = 0; i < precalcSteps; i++) {
+        float position = i * max / precalcSteps;
+        RGBColorf at = {
+            rsp.interpolateAt(position),
+            gsp.interpolateAt(position),
+            bsp.interpolateAt(position)
+        };
+        this->colors.push_back(at);
+    }
+}
+
+
+Gradient Gradient::defaultGradient(void)
+{
+    /*QFile res(":/gradients/default");
+    res.open(QIODevice::ReadOnly);
+    QString str = QString::fromUtf8(res.readAll());
+    return readXml(str);*/
+    return Gradient({
+        { RGBColor{ 0, 0, 0 }, 0.0f },
+        { RGBColor{ 0, 255, 255 }, 30.0f },
+        { RGBColor{ 50, 100, 170 }, 60.0f },
+        { RGBColor{ 180, 140, 20 }, 90.0f },
+        { RGBColor{ 255, 255, 0 }, 120.0f },
+        { RGBColor{ 143, 67, 0 }, 150.0f },
+        { RGBColor{ 255, 255, 255 }, 180.0f },
+        { RGBColor{ 20, 30, 180 }, 210.0f },
+        { RGBColor{ 20, 190, 30 }, 240.0f },
+        { RGBColor{ 120, 240, 120 }, 270.0f },
+        { RGBColor{ 40, 40, 40 }, 300.0f },
+    }, true);
+}
+
+
+RGBColor Gradient::get(float x) const
+{
+    if (colors.empty() || std::isnan(x) || std::isinf(x))
+        return RGBColor();
+    /*const auto [left, right, lerp] = getNeighbors(x);
+    RGBColor lerped = lerpColors(left, right, lerp);
+    return lerped;*/
+
+    if (x < 0)
+        return colors[0];
+    if (x > this->max) {
+        if (repeat)
+            x = ::fmodf(x, this->max);
+        else
+            x = this->max;
+    }
+    float pos = x * colors.size() / max;
+    if (pos < 0) {
+        pos = 0;
+    }
+    if (pos > colors.size() - 1) {
+        pos = colors.size() - 1;
+    }
+
+    int left = int(pos);
+    int right = int(pos + 1);
+    float lerp = pos - left;
+
+    if (lerp < 1e-5f) {
+        return colors[left];
+    }
+    else {
+        return lerpColors(colors[left], colors[right], lerp);
+    }
+}
+
+
+RGBColorf Gradient::lerpColors(RGBColorf a, RGBColorf b, float val)
+{
+    return RGBColorf {
+        b.r * val + a.r * (1 - val),
+        b.g * val + a.g * (1 - val),
+        b.b * val + a.b * (1 - val)
+    };
+}
+
+
+RGBColor Gradient::lerpColors(RGBColor a, RGBColor b, float val)
+{
+    return RGBColor{ lerpColors(RGBColorf{ a }, RGBColorf{ b }, val) };
+}
+
+
+/*std::tuple<RGBColor, RGBColor, float> Gradient::getNeighbors(float x) const
+{
+    for (auto it = colors.begin(); it != colors.end(); ++it) {
+        if (it->second > x) {
+            if (it == colors.begin()) {
+                return { it->first, it->first, 0 };
+            }
+            else {
+                float lerp = (x - (it - 1)->second) / (it->second - (it - 1)->second);
+                return { (it - 1)->first, it->first, lerp };
+            }
+        }
+    }
+    return { (colors.end() - 1)->first, (colors.end() - 1)->first, 0 };
+}*/

mandelvid/src/MandelVideoGenerator.cpp

@@ -0,0 +1,147 @@
+#include "MandelVideoGenerator.h"
+#include "VideoStream.h"
+#include "Mandel.h"
+#include <thread>
+#include <cmath>
+
+MandelVideoGenerator::MandelVideoGenerator(const ExportVideoInfo& evi) :
+    evi{ evi }
+{
+}
+
+
+void MandelVideoGenerator::generate(void)
+{
+    mnd::MandelContext ctxt = mnd::initializeContext();
+    mnd::Generator& gen = ctxt.getDefaultGenerator();
+    mnd::MandelInfo mi;
+    mi.bWidth = evi.width * 2;
+    mi.bHeight = evi.height * 2;
+    mi.maxIter = evi.maxIterations;
+
+    VideoStream vs(evi.width, evi.height, evi.path, evi.bitrate, evi.fps, evi.preset.c_str());
+
+    mnd::Real x = evi.end.x + evi.end.width / 2;
+    mnd::Real y = evi.end.y + evi.end.height / 2;
+    mnd::Real w = evi.start.width;
+    mnd::Real h = evi.start.height;
+
+    mnd::Real bigW = 10000000000000000.0;
+    double bigFac = 1.0;
+    Bitmap<RGBColor> big;
+    Bitmap<RGBColor> small;
+
+    while(w > evi.end.width || h > evi.end.height) {
+        mi.view = mnd::MandelViewport{ x - w/2, y - h/2, w, h };
+
+        if (bigW > 2 * w) {
+            Bitmap<float> raw{ evi.width * 2, evi.height * 2 };
+            gen.generate(mi, raw.pixels.get());
+            //auto before = std::chrono::high_resolution_clock::now();
+            big = raw.map<RGBColor>([&mi, this] (float i) {
+                return i >= mi.maxIter ? RGBColor{ 0, 0, 0 } : evi.gradient.get(i);
+            });
+            /*mi.view.zoomCenter(0.5);
+            gen.generate(mi, raw.pixels.get());
+            small = raw.map<RGBColor>([] (float x) { return
+                RGBColor{ uint8_t(::sin(x / 100) * 127 + 127), uint8_t(::sin(x / 213) * 127 + 127), uint8_t(::cos(x / 173) * 127 + 127) };
+            });*/
+            bigW = w;
+            bigFac = 1.0;
+        }
+
+        vs.addFrame(overlay(big, small, bigFac));
+
+        w *= ::pow(0.99, evi.zoomSpeed);
+        h *= ::pow(0.99, evi.zoomSpeed);
+        bigFac *= ::pow(0.99, evi.zoomSpeed);
+    }
+}
+
+
+inline RGBColor biliniear(const Bitmap<RGBColor>& img, double x, double y)
+{
+    int xfloor = int(::floor(x));
+    int yfloor = int(::floor(y));
+    int xceil = int(::ceil(x));
+    int yceil = int(::ceil(y));
+
+    double xLerp = x - xfloor;
+    double yLerp = y - yfloor;
+
+    RGBColor samples[2][2] = {
+        {
+            img.get(xfloor, yfloor),
+            img.get(xfloor, yceil),
+        },
+        {
+            img.get(xceil, yfloor),
+            img.get(xceil, yceil),
+        }
+    };
+
+    double r = 0, g = 0, b = 0;
+
+    auto mklin = [] (double x) {
+        return x;
+    };
+    auto unlin = [] (double x) {
+        return x;
+    };
+
+    r += (1 - xLerp) * (1 - yLerp) * mklin(samples[0][0].r);
+    r += (1 - xLerp) * yLerp * mklin(samples[0][1].r);
+    r += xLerp * (1 - yLerp) * mklin(samples[1][0].r);
+    r += xLerp * yLerp * mklin(samples[1][1].r);
+
+    g += (1 - xLerp) * (1 - yLerp) * mklin(samples[0][0].g);
+    g += (1 - xLerp) * yLerp * mklin(samples[0][1].g);
+    g += xLerp * (1 - yLerp) * mklin(samples[1][0].g);
+    g += xLerp * yLerp * mklin(samples[1][1].g);
+
+    b += (1 - xLerp) * (1 - yLerp) * mklin(samples[0][0].b);
+    b += (1 - xLerp) * yLerp * mklin(samples[0][1].b);
+    b += xLerp * (1 - yLerp) * mklin(samples[1][0].b);
+    b += xLerp * yLerp * mklin(samples[1][1].b);
+
+    return RGBColor{ uint8_t(unlin(r)), uint8_t(unlin(g)), uint8_t(unlin(b)) };
+}
+
+
+inline RGBColor nearest(const Bitmap<RGBColor>& img, double x, double y)
+{
+    int xfloor = int(::floor(x));
+    int yfloor = int(::floor(y));
+    return img.get(xfloor, yfloor);
+}
+
+
+Bitmap<RGBColor> MandelVideoGenerator::overlay(const Bitmap<RGBColor>& outer,
+                         const Bitmap<RGBColor>& inner, double scale)
+{
+    printf("%lf\n", scale);
+    Bitmap<RGBColor> ret{ outer.width / 2, outer.height / 2 };
+    double newW = outer.width * scale * 2;
+    double newH = outer.height * scale * 2;
+    double newX = outer.width * (1 - scale) / 2;
+    double newY = outer.height * (1 - scale) / 2;
+
+    auto before = std::chrono::high_resolution_clock::now();
+#pragma omp parallel for schedule(static, 1)
+    for (int i = 0; i < ret.height; i++) {
+        for (int j = 0; j < ret.width; j++) {
+            double newJ = newX + j * newW / outer.width;
+            double newI = newY + i * newH / outer.height;
+            RGBColor a = biliniear(outer, newJ, newI);
+            ret.get(j, i) = a;
+        }
+    }
+    auto after = std::chrono::high_resolution_clock::now();
+    printf("gradient applied in: %lld microseconds\n", std::chrono::duration_cast<std::chrono::microseconds>(after - before).count());
+    fflush(stdout);
+    /*for (int i = 0; i < ret.height * ret.width; i++) {
+        ret.pixels[i] = outer.pixels[i];
+    }*/
+
+    return ret;
+}

mandelvid/src/VideoStream.cpp

@@ -0,0 +1,232 @@
+#include "VideoStream.h"
+
+#include <iostream>
+
+
+
+#if LIBAVCODEC_VERSION_INT < AV_VERSION_INT(55,28,1)
+#define av_frame_alloc  avcodec_alloc_frame
+#define av_frame_free  avcodec_free_frame
+#endif
+
+const uint8_t VideoStream::endcode[] = { 0, 0, 1, 0xb7 };
+
+
+VideoStream::VideoStream(int width, int height, const std::string& filename, int bitrate, int fps, const char* preset) :
+    width{ width & (~1) }, height{ height & (~1) }
+{
+    // only needed with ffmpeg version < 4
+    //avcodec_register_all();
+
+    codec = avcodec_find_encoder(AV_CODEC_ID_H264);
+    if (!codec) {
+        fprintf(stderr, "invalid codec\n");
+        exit(1);
+    }
+
+    AVOutputFormat* oformat = av_guess_format(nullptr, filename.c_str(), nullptr);
+    if (!oformat)
+        oformat = av_guess_format("mp4", nullptr, nullptr);
+    if (oformat == nullptr)
+        throw "invalid format";
+
+    codecContext = avcodec_alloc_context3(codec);
+
+    pkt = av_packet_alloc();
+    if (!pkt)
+        exit(1);
+
+    codecContext->bit_rate = bitrate * 1000;
+    codecContext->width = width;
+    codecContext->height = height;
+    codecContext->time_base = AVRational{ 1, fps };
+    codecContext->framerate = AVRational{ fps, 1 };
+
+    codecContext->gop_size = 5; /* emit one intra frame every five frames */
+    codecContext->max_b_frames = 1;
+    codecContext->pix_fmt = AV_PIX_FMT_YUV420P;
+
+    formatContext = avformat_alloc_context();
+    formatContext->oformat = oformat;
+    formatContext->video_codec_id = oformat->video_codec;
+
+    stream = avformat_new_stream(formatContext, codec);
+    if (!stream)
+        throw "error";
+
+    params = avcodec_parameters_alloc();
+    avcodec_parameters_from_context(params, codecContext);
+    stream->codecpar = params;
+
+    /*AVCPBProperties *props;
+    props = (AVCPBProperties*) av_stream_new_side_data(
+        stream, AV_PKT_DATA_CPB_PROPERTIES, sizeof(*props));
+    props->buffer_size = 1024 * 1024;
+    props->max_bitrate = 0;
+    props->min_bitrate = 0;
+    props->avg_bitrate = 0;
+    props->vbv_delay = UINT64_MAX;*/
+
+    if (codec->id == AV_CODEC_ID_H264)
+        av_opt_set(codecContext->priv_data, "preset", preset, 0);
+
+    if (avcodec_open2(codecContext, codec, nullptr) < 0) {
+        fprintf(stderr, "could not open codec\n");
+        exit(1);
+    }
+    avio_open(&formatContext->pb, filename.c_str(), AVIO_FLAG_WRITE);
+
+    if (avformat_write_header(formatContext, NULL) < 0) {
+        throw "error";
+    }
+    /*file = fopen(filename.c_str(), "wb");
+    if (!file) {
+        fprintf(stderr, "could not open %s\n", filename.c_str());
+        exit(1);
+    }*/
+
+    picture = av_frame_alloc();
+    av_frame_make_writable(picture);
+    picture->format = codecContext->pix_fmt;
+    picture->width  = codecContext->width;
+    picture->height = codecContext->height;
+
+    int retval = av_frame_get_buffer(picture, 0);
+    if (retval < 0) {
+        fprintf(stderr, "could not alloc the frame data\n");
+        exit(1);
+    }
+    //av_image_alloc(picture->data, picture->linesize, width, height, codecContext->pix_fmt, 32);
+
+    swsContext = sws_getContext(width, height,
+        AV_PIX_FMT_RGB24, width, height,
+        AV_PIX_FMT_YUV420P, 0, 0, 0, 0);
+}
+
+
+void VideoStream::encode(AVFrame* frame)
+{
+    int ret;
+
+    /* send the frame to the encoder */
+    ret = avcodec_send_frame(codecContext, frame);
+    if (ret < 0) {
+        fprintf(stderr, "error sending a frame for encoding\n");
+        exit(1);
+    }
+
+    while (ret >= 0) {
+        ret = avcodec_receive_packet(codecContext, pkt);
+        //ret = avcodec_encode_video2(codecContext, pkt, picture, &gotPacket);
+        if (ret == AVERROR(EAGAIN) || ret == AVERROR_EOF)
+            return;
+        else if (ret < 0) {
+            fprintf(stderr, "error during encoding\n");
+            exit(1);
+        }
+
+        printf("encoded frame %3d\"PRId64\" (size=%5d)\n", pkt->pts, pkt->size);
+        //fwrite(pkt->data, 1, pkt->size, outfile);
+        //av_interleaved_write_frame(formatContext, pkt);
+
+        av_packet_rescale_ts(pkt, AVRational{1, 60}, stream->time_base);
+        pkt->stream_index = stream->index;
+
+        av_write_frame(formatContext, pkt);
+        av_packet_unref(pkt);
+    }
+}
+
+
+VideoStream::~VideoStream()
+{
+    /* flush the encoder */
+    encode(nullptr);
+    av_write_trailer(this->formatContext);
+
+    /* add sequence end code to have a real MPEG file */
+    //fwrite(endcode, 1, sizeof(endcode), file);
+    //fclose(file);
+
+    avcodec_close(codecContext);
+    avio_close(formatContext->pb);
+    av_frame_unref(picture);
+    //av_free(codecContext);
+    avcodec_parameters_free(&params);
+    avcodec_free_context(&codecContext);
+    av_frame_free(&picture);
+    av_packet_free(&pkt);
+
+/*
+    AVPacket pkt;
+    av_init_packet(&pkt);
+    pkt.data = nullptr;
+    pkt.size = 0;
+
+    for (;;) {
+        avcodec_send_frame(codecContext, NULL);
+        if (avcodec_receive_packet(codecContext, &pkt) == 0) {
+            av_interleaved_write_frame(codecContext, &pkt);
+            av_packet_unref(&pkt);
+        }
+        else {
+            break;
+        }
+    }
+
+    av_write_trailer();
+    if (!(oformat->flags & AVFMT_NOFILE)) {
+        int err = avio_close(ofctx->pb);
+        if (err < 0) {
+            Debug("Failed to close file", err);
+        }
+    }*/
+
+}
+
+
+void VideoStream::addFrame(const Bitmap<RGBColor>& frame)
+{
+    int retval = av_frame_make_writable(picture);
+    if (retval < 0)
+        exit(1);
+
+    /* prepare a dummy image */
+    /* Y */
+    /*for(int y = 0; y < height; y++) {
+        for(int x = 0; x < width; x++) {
+            picture->data[0][y * picture->linesize[0] + x] = frame.get(x, y).r / 2;
+        }
+    }*/
+
+    /* Cb and Cr */
+    /*for(int y=0;y<height / 2;y++) {
+        for(int x=0;x<width / 2;x++) {
+            picture->data[1][y * picture->linesize[1] + x] = frame.get(x * 2, y * 2).g / 2;
+            picture->data[2][y * picture->linesize[2] + x] = frame.get(x * 2, y * 2).b / 2;
+        }
+    }*/
+
+    /*auto gammaCorrect = [] (const RGBColor& rgb) {
+        const float gamma = 2.2f;
+        return RGBColor {
+            uint8_t(::powf(rgb.r / 255.0f, 1.0f / gamma) * 255),
+            uint8_t(::powf(rgb.g / 255.0f, 1.0f / gamma) * 255),
+            uint8_t(::powf(rgb.b / 255.0f, 1.0f / gamma) * 255),
+        };
+    };
+
+    Bitmap<RGBColor> gammaCorrected = frame.map<RGBColor>(gammaCorrect);*/
+
+    const uint8_t* pixelPointer[] = { reinterpret_cast<const uint8_t*>(frame.pixels.get()), 0 };
+    const int linesizeIn[] = { int(frame.width * sizeof(RGBColor)) };
+
+    sws_scale(swsContext, pixelPointer, linesizeIn, 0,
+        frame.height, picture->data, picture->linesize);
+
+    picture->pts = frameIndex++;
+
+    /* encode the image */
+    encode(picture);
+}
+

mandelvid/src/main.cpp

@@ -0,0 +1,40 @@
+#include "MandelVideoGenerator.h"
+#include "Gradient.h"
+#include "Mandel.h"
+#include "Fixed.h"
+
+
+int main() {
+    //mnd::MandelContext mndCtxt = mnd::initializeContext();
+
+
+    ExportVideoInfo evi;
+    
+    evi.start = mnd::MandelViewport::standardView();
+    evi.end = mnd::MandelViewport {
+        mnd::Real("-1.5016327722130767973008541252724123393337183519056236025189105693015282429244791506194548898968185999262221668435271537932672968559900159142085320685031"),
+        mnd::Real("9.1949171527697821768939276268368163504538591789778359909730511642378316080598664365235178721745031546786105261407973733873085119833457073054327967448264e-06"),
+        mnd::Real("1.6236294899543021550377844129369984149698872979955210084321757728274664401182171658849308001321609757279087031477100527629814577654596624031152718524352e-28"),
+        mnd::Real("1.2246019034401093377903721086780361028058704962292211685926779200766324399350798858672587301860274703389823933260119617558370004128301410779021141722617e-28")
+    };
+    //evi.end.zoomCenter(1.0e+27);
+    evi.gradient = Gradient::defaultGradient();
+
+    evi.width = 2560;
+    evi.height = 1440;
+    evi.maxIterations = 20000;
+    evi.fps = 60;
+    evi.zoomSpeed = 0.9;
+    evi.path = "video.avi";
+    evi.bitrate = 15000;
+    evi.preset = "slow";
+
+    evi.start.adjustAspectRatio(evi.width, evi.height);
+
+    MandelVideoGenerator mvg(evi);
+
+    mvg.generate();
+    return 0;
+}
+
+