Julia Fractal Drawing in C++

Question

I'd like a code review for this program written to draw Julia fractals. I'm specifically looking for feedback on:

• Style: I haven't written much C++ before, so I'm interested in knowing better or more idiomatic ways of formatting or structuring my code.
• Performance: Am I doing anything that's needlessly slowing down my program, or that I could optimize out?

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;
}

Answer 1

Unqualified names

The namespace identifier is missing from a lot of names - e.g. std::sqrt, std::log, std::abs, std::stoi, std::stod. It's not portable to rely on the unqualified names being defined.

Input parser

A lot of this is unnecessarily verbose. There's no need to write this->tokens all the time when tokens is perfectly clear. For example, I'd write its constructor as

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

I've changed argc to be int rather than int&, and constructed the strings in-place in the vector.

getCmdOption could benefit from using auto to declare the iterator (and I'm not convinced the long name helps - it's obviously getting from command arguments because its a member of our arguments object):

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;
}

---

This gives:

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();
    }
};

Vector

The Vec3D class misses some of the obvious operators you might need in future (+=, /=, etc), but you can add them when you need them. The default memberwise == would be the same as what you've written, so you don't actually need == or !=. But they aren't very useful anyway, as comparing floating-point values for equality doesn't usually work as we want.

The operators that we do have can all be declared const (and accept double arguments by value rather than by reference).

The print() method is unconventional - we normally write that as:

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

The dist() method isn't using this - we could make it static, or remove one of the arguments. I'd write a no-argument form, and use that for the two-argument one:

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

(std::hypot() is better behaved and optimised than hand-writing Pythagoras' formula)

---

This gives:

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);
}

Having written all this, I later discovered that we were only using this class for interpolating pixel values (see later in answer), and it got excised from my final version.

Counting iterations

It's going to be a bit more efficient to compute the squared magnitude of z than its absolute value: std::norm(z) <= 4. I'd write the counting loop as a for loop like this:

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);
}

Or even

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);
}

You might consider writing your own pow(complex, int) and testing to see whether that's more efficient than std::pow(complex, double). Search for "binary exponentiation" for an outline of the algorithm (and remember that negative powers are just the reciprocals of the corresponding positive ones).

Colour interpolation

I struggled to follow the logic of linear_interpolated_gradient(). However, the good news is that we don't need it at all. We were using it to generate a 100-colour palette from an initial 6-colour palette, but then we were using that 100-colour palette to interpolate an output colour. We can interpolate directly from the 6-colour palette and save all that.

Further, repurposing Vec3 to represent colours is misleading. Let's work with png::rgb_pixel directly:

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);
}

We'll then need to adjust the scaling for mu in main() to fit in colors.size() rather than palette.size() (which should have been written as 99, not 100 - a tricky off-by-one error).

Options processing

The first part of main() examines the command-line arguments. We do some checking, like this:

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

We ought to print the usage to std::cerr and return non-zero to indicate failure. Alternatively, we could make -f optional, and write to std::cout if a filename is not specified - that allows us to pipe into further processing (e.g. pngcrush). It's probably a good idea to ensure we can write the file before spending time calculating - this is best done by opening the file at this point, and later using png::image::write_stream() to write the data:

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;

We should leave opening the file until after checking the other arguments, so that we don't overwrite an existing file if we have to exit early with an error message.

Calculating the domain

There's a couple of constants 1.92 and 1.08 whose origins aren't obvious. I eventually reasoned that these are XPIXELS / 1000.0 and YPIXELS / 1000.0 respectively. We can change to make that more explicit:

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;

(we don't need separate delta_x and delta_y - they evaluated to the same value).

Logging output

The verbose output stream needs to go to std::clog, not std::cout. That's necessary for pipelining to work.

Reduce scope of pixel position

We can move im and re into the for loop, and avoid needing to declare them private to OpenMP.

OpenMP scheduling

I thought that because rows have different computation times to each other, that #pragma omp parallel for schedule(dynamic) might improve execution times. But it seems that the overhead outweighs the benefit, in my testing.

---

My version

#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;
}

Answer 2

There's already a good answer, so I'm just going to raise a couple of small points.

I had the same thought as Toby Speight about saving the square root of abs(z), but I would apply it also in the logarithm and avoid re-evaluating a known value. This may be taking micro-optimisation too far, but you can judge for yourself whether you think the hit to readability is too much to be worth it.

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);
}

---

And this is a bit of a frame challenge. Linear interpolation in RGB creates muddy colours:

Using a colour space which is designed to be closer to the way we perceive colours you can avoid the loss of vibrancy. E.g. with HSL, interpolating between the same endpoints:

HSL is probably a good compromise between speed and linear perception, although if you want to focus on interpolation quality over speed you could look into XYZ and Lab*.