Circle map bifurcation generator

main.cpp

#include <cv.h>
#include <highgui.h>

#include <cmath>
#include <string>
#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <algorithm> 

const int POINTS = 20000;
const int ITERS = 1000;
const int WIDTH = 600;
const int HEIGHT = 1200;
const double OMEGA = 1.0 / 3;

const double X_LIM = 1;
const double Y_LIM = 4 * M_PI;

double calculate(double omega, double K, std::vector<double> &xs) {
    double theta = 0;
 
    int sz = xs.size();
    for (int i = 0; i < ITERS + sz; ++i) {
        theta = (theta + omega - K / (2 * M_PI) * std::sin(2 * M_PI * theta));
        theta = theta - floor(theta / X_LIM) * X_LIM;
        
        int j = i - ITERS;
        if (j >= 0) xs[j] = theta;
    }

    return theta;
}

cv::Vec3d color(double x) {
    if (4 * x < 1) return cv::Vec3d(0, 0, 4*x);
    if (4 * x < 2) return cv::Vec3d(0, 4*x - 1, 2 - 4*x);
    if (4 * x < 3) return cv::Vec3d(4*x-2, 1, 0);
    return cv::Vec3d(1, 4 - 4*x, 0);
}

int main(int argc, char **argv) {

    cv::Mat out = cv::Mat(HEIGHT, WIDTH, CV_8UC3);

    unsigned **count = new unsigned*[HEIGHT];
    unsigned *max = new unsigned[HEIGHT];
    for (int i = 0; i < HEIGHT; ++i) {
        count[i] = new unsigned[WIDTH];
        for (int j = 0; j < WIDTH; ++j) {
            count[i][j] = 0;
        }
        max[i] = 0;
    }

    std::cout << "Generating data..." << std::endl;
    
    for (int i = 0; i < HEIGHT; ++i) {
        double K = Y_LIM / HEIGHT * i;
        std::cout << "Processing " << (100 * i / HEIGHT) << "%...\r";

        std::vector<double> xs(POINTS);
        calculate(OMEGA, K, xs);

        for (int k = 0; k < xs.size(); ++k) {
            double x = xs[k];
            int j = floor(x / X_LIM * WIDTH);
            ++count[i][j];
            if (count[i][j] > max[i]) max[i] = count[i][j];
        }
    }
    std::cout << "Data generated." << std::endl;

    for (int i = 0; i < HEIGHT; ++i) {
        for (int j = 0; j < WIDTH; ++j) {
            if (max[i] > 0) {
                double v = double(count[i][j]) / max[i];
                cv::Vec3d c = color(v);
                out.at<cv::Vec3b>(HEIGHT - i - 1, j) = cv::Vec3b(255 * c[2], 255 * c[1], 255 * c[0]);

            } else {
                out.at<cv::Vec3b>(HEIGHT - i - 1, j) = cv::Vec3b(0, 0, 0);
            }
        }
    }

    // выводим изображение
    cv::imwrite("out.png", out);

    return 0;
}

Сделал возможность менять 1/4, 1/2, 3/4 в функции расчёта цвета
Добавил логарифмическую шкалу. Обычный логарифм оказался уныло-ровным, поэтому можно сделать несколько итераций, чем больше — тем логарифмичнее. Но приблизиться к оригиналу так и не удалось, надо играться с итерациями и светом, либо придумывать новую шкалу.

#include <opencv/cv.h>
#include <opencv/highgui.h>

#include <cmath>
#include <string>
#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <algorithm> 

#define M_PI 3.14159265358

const int POINTS = 20000;
const int ITERS = 1000;
const int WIDTH = 450;
const int HEIGHT = 900;
const double OMEGA = 1.0 / 3;

const double X_LIM = 1;
const double Y_LIM = 4 * M_PI;

const bool LOG_SCALE = true;
const int LOG_ITERATIONS = 10;

double calculate(double omega, double K, std::vector<double> &xs) {
    double theta = 0;

    int sz = xs.size();
    for (int i = 0; i < ITERS + sz; ++i) {
        theta = (theta + omega - K / (2 * M_PI) * std::sin(2 * M_PI * theta));
        theta = theta - floor(theta / X_LIM) * X_LIM;

        int j = i - ITERS;
        if (j >= 0) xs[j] = theta;
    }

    return theta;
}

cv::Vec3d color(double x) {
    double p1 = 1.0 / 4;
    double p2 = 2.0 / 4;
    double p3 = 3.0 / 4;
    if (x < p1) return cv::Vec3d(0, 0, x / p1);
    if (x < p2) return cv::Vec3d(0, (x - p1) / (p2 - p1), (p2 - x) / (p2 - p1));
    if (x < p3) return cv::Vec3d((x - p2) / (p3 - p2), 1, 0);
    return cv::Vec3d(1, (1 - x) / (1 - p3), 0);
}

double logScale(double x) {
    double result = x;
    for (int i = 0; i < LOG_ITERATIONS; i++) {
        result = log(1 + result) / log(2);
    }
    return result;
}

int main(int argc, char **argv) {

    cv::Mat out = cv::Mat(HEIGHT, WIDTH, CV_8UC3);

    unsigned **count = new unsigned*[HEIGHT];
    unsigned *max = new unsigned[HEIGHT];
    for (int i = 0; i < HEIGHT; ++i) {
        count[i] = new unsigned[WIDTH];
        for (int j = 0; j < WIDTH; ++j) {
            count[i][j] = 0;
        }
        max[i] = 0;
    }

    std::cout << "Generating data..." << std::endl;

    for (int i = 0; i < HEIGHT; ++i) {
        double K = Y_LIM / HEIGHT * i;
        std::cout << "Processing " << (100 * i / HEIGHT) << "%...\r";

        std::vector<double> xs(POINTS);
        calculate(OMEGA, K, xs);

        for (size_t k = 0; k < xs.size(); ++k) {
            double x = xs[k];
            int j = static_cast<int>(floor(x / X_LIM * WIDTH));
            ++count[i][j];
            if (count[i][j] > max[i]) max[i] = count[i][j];
        }
    }
    std::cout << "Data generated." << std::endl;

    for (int i = 0; i < HEIGHT; ++i) {
        for (int j = 0; j < WIDTH; ++j) {
            if (max[i] > 0) {
                double v = static_cast<double>(count[i][j]) / max[i];
                if (LOG_SCALE) {
                    v = logScale(v);
                }
                cv::Vec3d c = color(v);
                out.at<cv::Vec3b>(HEIGHT - i - 1, j) = cv::Vec3b(255 * c[2], 255 * c[1], 255 * c[0]);

            } else {
                out.at<cv::Vec3b>(HEIGHT - i - 1, j) = cv::Vec3b(0, 0, 0);
            }
        }
    }

    cvSaveImage("out.png", &(IplImage(out)));

    return 0;
}