Юля фрактальное рисование в C++


Мне бы код этой программы, написанной для рисования Юлия фракталы. Я специально искала отзывы на:

  • Стиль: я не писал много на C++, так что я заинтересован в зная лучше и более идиоматических способов форматирования и структурирования кода.
  • Производительность: я делаю все, что намеренное замедление мою программу, или что я мог бы оптимизировать выход?

julia.cpp

#include <iostream>
#include <complex>
#include <vector>
#include <string>
#include <algorithm>

#include <cmath>
#include <cstdlib>
#include <unistd.h>
#include <omp.h>

#include <png++/png.hpp>

// defining a few constants
#define MAX_ITERATIONS 2000
#define XPIXELS 1920
#define YPIXELS 1080

typedef std::complex<double> complex;

// code for this taken from
// https://stackoverflow.com/questions/865668/how-to-parse-command-line-arguments-in-c
class InputParser {
public:
    InputParser(int& argc, char **argv) {
        for (int i = 1; i < argc; ++i) {
            this->tokens.push_back(std::string(argv[i]));
        }
    }

    const std::string& getCmdOption(const std::string& option) const{
        std::vector<std::string>::const_iterator iter;
        iter = std::find(this->tokens.begin(), this->tokens.end(), option);
        if (iter != this->tokens.end() && ++iter != this->tokens.end()) {
            return *iter;
        }
        static const std::string empty_string("");
        return empty_string;
    }

    std::vector<std::string> getCmdOptionN(const std::string& option, int n) const{
        std::vector<std::string>::const_iterator iter;
        std::vector<std::string> results;
        std::vector<std::string> empty_vec;
        iter = std::find(this->tokens.begin(), this->tokens.end(), option);
        if (iter == this->tokens.end()) {
            return empty_vec;
        }
        for (int i = 1; i < n+1; i++) {
            if (iter+i != this->tokens.end()) {
                results.push_back(*(iter+i));
            } else {
                return empty_vec;
            }
        }
        return results;
    }

    bool cmdOptionExists(const std::string& option) const{
        return std::find(this->tokens.begin(), this->tokens.end(), option) != this->tokens.end();
    }

private:
    std::vector<std::string> tokens;
};

// struct representing vectors in R^3
// along with some basic operations on them
struct Vec3 {
    double x;
    double y;
    double z;

    Vec3() {}
    Vec3(double xx, double yy, double zz) : x(xx), y(yy), z(zz) {}

    Vec3 operator+(const Vec3& v) {
        return Vec3(x+v.x, y+v.y, z+v.z);
    }

    Vec3 operator-(const Vec3& v) {
        return Vec3(x-v.x, y-v.y, z-v.z);
    }

    Vec3 operator*(double& a) {
        return Vec3(a*x, a*y, a*z);
    }

    Vec3 operator/(double& a) {
        return Vec3(x/a, y/a, z/a);
    }

    bool operator==(const Vec3& v) {
        return x == v.x && y == v.y && z == v.z;
    }

    void print() {
        std:: cout << "(" << x << ", " << y << ", " << z << ")\n";
    }

    double dist(const Vec3& u, const Vec3& v) {
        return sqrt((u.x - v.x)*(u.x - v.x) +
                (u.y - v.y)*(u.y - v.y) +
                (u.z - v.z)*(u.z - v.z));
    }
};

Vec3 operator*(double& a, const Vec3& v)
{
    return Vec3(a*v.x, a*v.y, a*v.z);
};

// counts the number of iterations it takes for a complex function `f(z) = z^power + c0` evaluated iteratively
// at an initial point `init` to grow greater than 2 in magnitude
// normalized to achieve smoother coloring, look at this webpage for details:
// http://linas.org/art-gallery/escape/escape.html
double normalized_iterations(complex init, complex c0, int power)
{
    complex z = init;
    int iters = 0;
    while ((abs(z) <= 2) && (iters < MAX_ITERATIONS)) {
        z = std::pow(z,power) + c0;
        iters += 1;
    }
    double mu = iters;
    if ( iters < MAX_ITERATIONS ) {
        mu = iters + 1 - log(log(abs(z))) / log(power);
    }
    return mu;
}

// computes v + t(u - v)
// t should be a value between 0 and 1
Vec3 linear_interpolation(Vec3& v, Vec3& u, double t)
{
    return v + t*(u - v);
}

// creates a linear gradient of SIZE colours, using RGB values from PTS
// interspersed evenly
std::vector<Vec3> linear_interpolated_gradient(std::vector<Vec3> pts, int size)
{
    std::vector<Vec3> pal;
    int to_travel = size;
    int lines_left = pts.size();
    int pts_to_color;
    for (int i = 0; i < pts.size()-1; i++) {
        if (to_travel % lines_left != 0) {
            pts_to_color = (to_travel / lines_left)+1;
        } else {
            pts_to_color = to_travel / lines_left;
        }
        to_travel = to_travel - pts_to_color;
        lines_left--;
        double scaling = 1.0 / pts_to_color;
        Vec3 delta_vec = scaling*(pts[i+1] - pts[i]);
        Vec3 next_color = pts[i];
        for (int j = 0; j < pts_to_color; j++) {
            pal.push_back(next_color);
            next_color = next_color + delta_vec;
        }
    }
    return pal;
}


int main(int argc, char *argv[])
{
    const std::string& usage = "Usage: -f <filename> [-p <power>] -c <real_part> <imag_part> [-origin <x> <y>] [-z <zoom>] [-verbose]\nPower defaults to 2, origin defaults to (0,0)\n";

    // Parsing command line arguments
    InputParser input(argc, argv);
    const std::string& filename = input.getCmdOption("-f");
    if (filename.empty()) {
        std::cout << usage;
        return 0;
    }

    const std::string& power_string = input.getCmdOption("-p");
    int power = 2;
    if (!power_string.empty()) {
        power = stoi(power_string);
    }

    const std::vector<std::string>& complex_strings = input.getCmdOptionN("-c", 2);
    if (complex_strings.empty()) {
        std::cout << usage;
        return 0;
    }
    const double real_part = stod(complex_strings[0]);
    const double imag_part = stod(complex_strings[1]);

    double origin_x = 0.0, origin_y = 0.0;
    const std::vector<std::string>& origin_strings = input.getCmdOptionN("-origin", 2);
    if (!origin_strings.empty()) {
        origin_x = stod(origin_strings[0]);
    origin_y = stod(origin_strings[1]);
    }

    double zoom = 1.0;
    const std::string& zoom_string = input.getCmdOption("-z");
    if (!zoom_string.empty()) {
        zoom = stod(zoom_string);
    }

    bool verbose = input.cmdOptionExists("-verbose");

    // Setting up parameters
    const complex complex_constant(real_part, imag_part);

    // computing C -> pixel mapping
    double im_start = origin_y + 1.08/zoom;
    double re_start = origin_x - 1.92/zoom;
    double delta_y = 2*1.08/zoom / YPIXELS, delta_x = 2*1.92/zoom / XPIXELS;
    double im, re;

    if (verbose) {
        std::cout << "im_start = " << im_start << "\nre_start = " << re_start << std::endl;
        std::cout << "delta_y = " << delta_y << "\ndelta_x = " << delta_x << std::endl;
        std::cout << "zoom = " << zoom << std::endl;
        std::cout << "Running on " << omp_get_max_threads() << " threads" << std::endl;
    }

    // another thing that would be nice to add is allow the user to input a file
    // consisting of RGB triples to set up the color palette with
    std::vector<Vec3> colors;
    colors.push_back(Vec3(0, 0, 0));
    colors.push_back(Vec3(213, 67, 31));
    colors.push_back(Vec3(251, 255, 121));  
    colors.push_back(Vec3(62, 223, 89));
    colors.push_back(Vec3(43, 30, 218));
    colors.push_back(Vec3(0, 255, 247));

    std::vector<Vec3> palette = linear_interpolated_gradient(colors, 100);
    png::image<png::rgb_pixel> image(XPIXELS, YPIXELS);

    #pragma omp parallel for private(re) private(im)
    for (int y = 0; y < YPIXELS; y++) {
        if (verbose) {
            std::cout << "Computing row " << y+1 << '/' << YPIXELS << "...\n";
        }
        im = im_start - y*delta_y;
        for (int x = 0; x < XPIXELS; x++) {
            re = re_start + x*delta_x;
            complex init(re,im);
            double mu = normalized_iterations(init, complex_constant, power);
            // scale mu to be in the range of 1-100
            mu *= 100.0/MAX_ITERATIONS;
            double tmp;
            Vec3 color1 = palette[(int)floor(mu)];
            Vec3 color2 = palette[(int)ceil(mu)];
            Vec3 color = linear_interpolation(color1, color2, modf(mu, &tmp));
            image[y][x] = png::rgb_pixel(color.x, color.y, color.z);
        }
    }
    image.write(filename);
    return 0;
}


2613
18
задан 22 февраля 2018 в 10:02 Источник Поделиться
Комментарии
2 ответа

Неполные имена

Идентификатор пространства имен отсутствует много названий - например, std::sqrt, std::log, std::abs, std::stoi, std::stod. Это не портативный, чтобы полагаться на неполные имена определяются.

Анализатор входного сигнала

Многое из этого излишне многословен. Нет необходимости писать this->tokens все время, когда tokens это прекрасно видно. Например, я мог бы написать его как конструктор

InputParser(int argc, char **argv) {
for (int i = 1; i < argc; ++i) {
tokens.emplace_back(argv[i]);
}
}

Я изменилась argc быть int а не int&и построены строки в вектор.

getCmdOption может извлечь выгоду из использования auto объявить итератор (и я не уверена, что длинное имя помогает - это очевидно получая от аргументов, потому что его член наши доводы объект):

const std::string& optionValue(const std::string& option) const
{
static const std::string not_found{};
auto it = std::find(tokens.begin(), tokens.end(), option);
return it != tokens.end() && ++it != tokens.end() ? *it : not_found;
}


Это дает:

class InputParser
{
std::vector<std::string> tokens = {};

public:
InputParser(int argc, char **argv) {
for (int i = 1; i < argc; ++i) tokens.emplace_back(argv[i]);
}

const std::string& optionValue(const std::string& option) const
{
static const std::string not_found{};
auto it = std::find(tokens.begin(), tokens.end(), option);
return it != tokens.end() && ++it != tokens.end() ? *it : not_found;
}

std::vector<std::string> optionValues(const std::string& option, int n) const
{
static const std::vector<std::string> not_found{};
auto it = std::find(tokens.begin(), tokens.end(), option);
if (std::distance(it, tokens.end()) <= n) return not_found;

std::vector<std::string> results;
results.reserve(n);
while (n--) results.push_back(*++it);
return results;
}

bool contains(const std::string& option) const
{
return std::find(tokens.begin(), tokens.end(), option) != tokens.end();
}
};

Вектор

В Vec3D класс пропускает некоторые из очевидных операторы могут вам понадобиться в будущем (+=, /=и т. д.), Но вы можете добавить их, когда вы нуждаетесь в них. По умолчанию почленных == будет то же самое, что вы написали, Так вы на самом деле не нужно == или !=. Но они не очень полезны в любом случае, как сравнение значений с плавающей точкой на равенство, как правило, не так, как надо.

Операторы, что у нас может все быть объявлены const (и принять double аргументов по значению, а не по ссылке).

В print() метод нетрадиционный - мы обычно пишем что так:

friend std::ostream& operator<<(std::ostream& os, const Vec3& v)
{
return os << "(" << v.x << ", " << v.y << ", " << v.z << ")\n";
}

В dist() метод не использует this - мы могли бы сделать его статическим, или удалить один из аргументов. Я пишу без аргументов, а на два-аргумент один:

double dist() const {
return std::hypot(std::hypot(x, y), z);
}
double dist(const Vec3& other) const {
return (*this - other).dist();
}

(std::hypot() лучше вели себя и оптимизированной, чем почерк Пифагора формула)


Это дает:

struct Vec3 {
double x;
double y;
double z;

Vec3() : Vec3(0, 0, 0) {}
Vec3(double xx, double yy, double zz) : x(xx), y(yy), z(zz) {}

Vec3 operator+(const Vec3& v) const {
return Vec3(x+v.x, y+v.y, z+v.z);
}

Vec3 operator-(const Vec3& v) const {
return Vec3(x-v.x, y-v.y, z-v.z);
}

Vec3 operator*(double a) const {
return Vec3(a*x, a*y, a*z);
}

Vec3 operator/(double a) const {
return Vec3(x/a, y/a, z/a);
}

friend std::ostream& operator<<(std::ostream& os, const Vec3& v)
{
return os << "(" << v.x << ", " << v.y << ", " << v.z << ")\n";
}

double dist() const {
return std::hypot(std::hypot(x, y), z);
}
double dist(const Vec3& other) const {
return (*this - other).dist();
}
};

Vec3 operator*(double a, const Vec3& v)
{
return Vec3(a*v.x, a*v.y, a*v.z);
}

Написав все это, я позже обнаружил, что мы были только использовать этот класс для интерполяции значений пикселов (см. Ниже Ответ), и он получил вырезали из моей окончательной версии.

Подсчет итераций

Это будет немного более эффективным, чтобы вычислить квадрат величины z чем ее абсолютное значение: std::norm(z) <= 4. Я напишу подсчитывая петли как for цикл вроде этого:

double normalized_iterations(complex z, complex c0, int power)
{
int iters;
for (iters = 0; std::norm(z) <= 4; ++iters) {
if (iters == MAX_ITERATIONS) return iters;
z = std::pow(z,power) + c0;
}
return iters + 1 - std::log(std::log(std::abs(z))) / std::log(power);
}

Или даже

double normalized_iterations(complex z, const complex c0, const int power)
{
int iters;
for (iters = 0; std::norm(z) <= 4; z = std::pow(z,power) + c0) {
if (++iters == MAX_ITERATIONS) return iters;
}
return iters + 1 - std::log(std::log(std::abs(z))) / std::log(power);
}

Вы могли бы рассмотреть возможность написания собственных pow(complex, int) и тестирование, чтобы увидеть, является ли это более эффективным, чем std::pow(complex, double). Поиск "бинарное" по плану алгоритма (и помните, что отрицательные силы, как раз обратными величинами соответствующих позитивных).

Цветовой интерполяции

Я старался следовать логике linear_interpolated_gradient(). Тем не менее, хорошей новостью является то, что он нам не нужен вообще. Мы использовали его, чтобы создать 100-цветную палитру из 6-и цветовая палитра, но тогда мы использовали, что 100-цветовая палитра для интерполяции выходной цвет. Мы можем интерполировать прямо с 6-цветной палитрой и сохранить все это.

Далее, перепрофилирование Vec3 для представления цвета вводит в заблуждение. Давайте работать с png::rgb_pixel напрямую:

static const std::vector<png::rgb_pixel> colors{
{ 0, 0, 0},
{213, 67, 31},
{251, 255, 121},
{ 62, 223, 89},
{ 43, 30, 218},
{ 0, 255, 247}
};

png::rgb_pixel linear_interpolation(const png::rgb_pixel& v,
const png::rgb_pixel& u, double a)
{
auto const b = 1 - a;
return png::rgb_pixel(b*v.red + a*u.red,
b*v.green + a*u.green,
b*v.blue + a*u.blue);
}

Мы будем затем нужно настроить масштабирование для mu в main() чтобы поместиться в colors.size() а не palette.size() (что должно быть написано как 99, а не 100 - хитрая на одну ошибку).

Параметры обработки

Первая часть main() анализирует аргументы командной строки. Мы делаем проверку, как это:

if (filename.empty()) {
std::cout << usage;
return 0;
}

Мы должны печатать использование std::cerr и возвращают не ноль, чтобы указать отказ. Кроме того, мы могли бы сделать -f необязательно, и напишите std::cout если имя файла не указано - что позволяет нам трубы на дальнейшую переработку (например, pngcrush). Это, вероятно, хорошая идея, чтобы убедиться, что мы можем сохранить файл, прежде чем тратить время на подсчеты - лучше всего это сделать, открыв файл в этот момент, а позже через png::image::write_stream() для записи данных:

std::ofstream outfile;
const std::string& filename = input.optionValue("-f");
std::ostream& out = filename.empty()
? std::cout
: (outfile.open(filename, std::ios::out|std::ios::trunc|std::ios::binary), outfile);

if (!out) {
perror(filename.c_str());
return 1;
}

// ....

image.write_stream(out);
return 0;

Мы должны оставить открытием файла до проверив другие доводы, так что мы не перезаписывать существующий файл, если нам придется покинуть рано с сообщением об ошибке.

Расчета домен

Есть несколько констант 1.92 и 1.08 чье происхождение не очевидно. В конце концов я рассудил, что это XPIXELS / 1000.0 и YPIXELS / 1000.0 соответственно. Мы можем изменить, чтобы сделать это более ясным:

auto const x_scale = .001 * XPIXELS / zoom;
auto const y_scale = .001 * YPIXELS / zoom;

const double im_start = origin_y + y_scale;
const double re_start = origin_x - x_scale;
const double delta = .002 / zoom;

(нам не нужна отдельная delta_x и delta_y - они оценивали одно и то же значение).

Выходные данные журнала

В verbose выходной поток должен идти std::clogне std::cout. Это необходимо для конвейера для работы.

Уменьшить объем позиции пиксела

Мы можем двигаться im и re в for петли, и избежать необходимости объявлять их private с OpenMP.

Модель OpenMP планирование

Я подумал, что строки имеют различные вычисления раз друг другу, что #pragma omp parallel for schedule(dynamic) может улучшить исполнение Времен. Но мне кажется, что накладные расходы превышают выгоду, в моем тестировании.


Моя версия

#include <algorithm>
#include <array>
#include <complex>
#include <iostream>
#include <string>
#include <vector>

#include <cmath>
#include <cstdlib>

#include <omp.h>

#include <png++/png.hpp>

// defining a few constants - these should be user-specifiable
static constexpr auto max_iterations = 2000;
static constexpr auto width = 1920;
static constexpr auto height = 1080;

using complex = std::complex<double>;

class InputParser
{
std::vector<std::string> tokens = {};

public:
InputParser(int argc, char **argv) {
for (int i = 1; i < argc; ++i) tokens.emplace_back(argv[i]);
}

const std::string& optionValue(const std::string& option) const
{
static const std::string not_found{};
auto it = std::find(tokens.begin(), tokens.end(), option);
return it != tokens.end() && ++it != tokens.end() ? *it : not_found;
}

std::vector<std::string> optionValues(const std::string& option, int n) const
{
static const std::vector<std::string> not_found{};
auto it = std::find(tokens.begin(), tokens.end(), option);
if (std::distance(it, tokens.end()) <= n) return not_found;

std::vector<std::string> results;
results.reserve(n);
while (n--) results.push_back(*++it);
return results;
}

bool contains(const std::string& option) const
{
return std::find(tokens.begin(), tokens.end(), option) != tokens.end();
}
};

// counts the number of iterations it takes for a complex function `f(z) = z^power + c0` evaluated iteratively
// at an initial point `init` to grow greater than 2 in magnitude
// normalized to achieve smoother coloring, look at this webpage for details:
// http://linas.org/art-gallery/escape/escape.html
double normalized_iterations(complex z, const complex c0, const int power)
{
int iters;
for (iters = 0; std::norm(z) <= 4; z = std::pow(z, power) + c0) {
if (++iters == max_iterations) return 1;
}
return (iters + 1 - std::log(std::log(std::abs(z))) / std::log(power)) / max_iterations;
}

// computes v + t(u - v)
// t should be a value between 0 and 1
png::rgb_pixel linear_interpolation(const png::rgb_pixel& v, const png::rgb_pixel& u, double a)
{
auto const b = 1 - a;
return png::rgb_pixel(b*v.red + a*u.red,
b*v.green + a*u.green,
b*v.blue + a*u.blue);
}

int main(int argc, char *argv[])
{
static const auto usage =
"Usage: -f <filename> [-p <power>] -c <real_part> <imag_part>"
" [-origin <x> <y>] [-z <zoom>] [-verbose]\n"
"Power defaults to 2, origin defaults to (0,0)\n";

// Parsing command line arguments
InputParser input(argc, argv);

if (input.contains("-h") || input.contains("--help")) {
std::cout << usage;
return 0;
}

const bool verbose = input.contains("-verbose");
const auto filename = input.optionValue("-f");
const auto power_string = input.optionValue("-p");
const auto complex_strings = input.optionValues("-c", 2);
const auto origin_strings = input.optionValues("-origin", 2);

int power = 2;
if (!power_string.empty()) {
power = std::stoi(power_string);
}
if (power < 2) {
// a waste of time
std::cerr << "Power must be at least 2" << std::endl;
return 1;
}

if (complex_strings.empty()) {
std::cerr << usage;
return 1;
}
const complex complex_constant{std::stod(complex_strings[0]),
std::stod(complex_strings[1])};

double origin_x = 0.0, origin_y = 0.0;
if (!origin_strings.empty()) {
origin_x = std::stod(origin_strings[0]);
origin_y = std::stod(origin_strings[1]);
}

double zoom = 1.0;
const auto zoom_string = input.optionValue("-z");
if (!zoom_string.empty()) {
zoom = std::stod(zoom_string);
}

std::ofstream outfile;
std::ostream& out = filename.empty()
? std::cout
: (outfile.open(filename, std::ios::out|std::ios::trunc|std::ios::binary), outfile);

if (!out) {
perror(filename.c_str());
return 1;
}

// Julia set parameters
auto const x_scale = .001 * width / zoom;
auto const y_scale = .001 * height / zoom;

const double im_start = origin_y + y_scale;
const double re_start = origin_x - x_scale;
const double delta = .002 / zoom;

if (verbose) {
std::clog << "im_start = " << im_start << "\nre_start = " << re_start << std::endl;
std::clog << "delta = " << delta << std::endl;
std::clog << "zoom = " << zoom << std::endl;
std::clog << "Running on " << omp_get_max_threads() << " threads" << std::endl;
}

// Could we allow user input of these values?
static const std::vector<png::rgb_pixel> colors{
{ 0, 0, 0},
{213, 67, 31},
{251, 255, 121},
{ 62, 223, 89},
{ 43, 30, 218},
{ 0, 255, 247}
};
static const auto max_color = colors.size() - 1;

png::image<png::rgb_pixel> image(width, height);

#pragma omp parallel for
for (int y = 0; y < height; y++) {
if (verbose)
#pragma omp critical
{
std::clog << "Computing row " << y+1 << '/' << height << "...\n";
}
double im = im_start - y*delta;
for (int x = 0; x < width; x++) {
double re = re_start + x*delta;
double mu = normalized_iterations({re,im}, complex_constant, power);
// scale mu to be in the range of colors
mu *= max_color;
auto i_mu = static_cast<std::size_t>(mu);
auto color1 = colors[i_mu];
auto color2 = colors[std::min(i_mu+1, max_color)];
image[y][x] = linear_interpolation(color1, color2, mu-i_mu);
}
}

image.write_stream(out);
return 0;
}

26
ответ дан 22 февраля 2018 в 11:02 Источник Поделиться

Есть уже хороший ответ, поэтому я просто собираюсь поднять пару мелких точек.

У меня была такая же мысль, как Тоби Спейт о спасении квадратный корень abs(z)но я бы применить его в логарифм и избежать переоценки известным значением. Это может быть слишком далеко микро-оптимизации, но вы можете судить для себя, считаете ли вы попали в читабельности слишком много, чтобы быть стоит.

double normalized_iterations(complex z, const complex c0, const int power)
{
int iters = 0;
double norm = std::norm(z);
for ( ; norm <= 4; norm = std::norm(z)) {
if (++iters == MAX_ITERATIONS) return iters;
z = std::pow(z, power) + c0;
}
return iters + 1 - std::log(0.5 * std::log(norm)) / std::log(power);
}


И это немного вызов кадров. Линейная интерполяция в RGB создает мутные цвета:

RGB interpolation from green to red through muddy orange-brown

Используя цветовое пространство, которое предназначено, чтобы быть ближе к, как мы воспринимаем цвета, вы можете избежать потери вибрации. Е. Г. с ВСЛ, интерполируя между теми же конечными точками:

HSL interpolation from green to red through bright yellow

ВСЛ является, вероятно, хорошим компромиссом между скоростью и линейного восприятия, хотя если вы хотите, чтобы сосредоточиться на интерполяции качество за скорость можно посмотреть в XYZ и ЛАБ*.

14
ответ дан 22 февраля 2018 в 04:02 Источник Поделиться